Expressions for boundary condition as a function of neural-network outputs

[boundaryconditions]

Hi, I am trying to reproduce a study of active flow control of the flow over a 2D cylinder with 2 jets positioned at ±90 degrees. I am solving the compressible NS equations for now. As a start, I set the boundary conditions for the jet as, for example,

[soln-bcs-jet_top]
type = sub-in-frv
theta = atan2(y,x) ; theta, not used for now
omega = 10*pi/180 ; jet width in radians
theta0_top = 90*pi/180 ; center of jet
R = 0.5 ; radius of cylinder

rho = rhoc
u = 0.09 / (%(omega)s*%(R)s*%(R)s) * x
v = 0.09 / (%(omega)s*%(R)s*%(R)s) * y
p = rhoc/(M*M*gamma) ; T = 1

where the “0.09” value is supposed to be the mass flow rate. With steady blowing, I get some result like the following:

The next step is to change the “0.09” value to be a function that will change at different time steps, involving one or two values inferred from a neural network. For now, say instead of 0.09, we obtain mdot (mass flow rate) from a neural network as an output. Furthermore, at the same time, we also need to feed some input to the neural network in form of primitive variables at some chosen locations (essentially similar to what the sampler plugin does). Can someone point me to the files that I could begin to modify in the source code to achieve this? Any additional advice is welcome.

TL;DR : I am seeking advice in modifying the source code to make expressions used in the setting the boundary condition to be not only a function of space and time, but a function of some parameters that would be inferred from a neural network.

[fdw]

This is possible with a custom boundary condition. In particular, boundary conditions have a prepare method which is called on the host-side in advance of a time-step stage. This gives them an opportunity to compute any quantities they like and modify any arrays they like. These arrays can then be used by the boundary condition kernel.

Unfortunately, while the functionality is there none of our current boundary conditions make use of this feature and so I do not have any examples.

Regards, Freddie.

[boundaryconditions]

Line 254 of PyFR/pyfr/solvers/base/system.py has the following

    def _prepare_kernels(self, t, uinbank, foutbank):
        _, binders = self._get_kernels(uinbank, foutbank)

        for b in self._bc_inters:
            b.prepare(t)

        for b in binders:
            b(t=t)

Is this the one called before every-time step? Furthermore, would the next step regarding custom BCs would be to add some code in PyFR/pyfr/solvers/navstokes/inters.py ? For instance,

class  NavierStokesSubInflowFrvNeuralBCInters(NavierStokesBaseBCInters):
    type = 'sub-in-frv-neural'
    cflux_state = 'ghost'

    def prepare(self, t):
         # do things here

    def __init__(self, be, lhs, elemap, cfgsect, cfg):
        super().__init__(be, lhs, elemap, cfgsect, cfg)

        self.c |= self._exp_opts(
            ['rho', 'u', 'v', 'w'][:self.ndims + 1], lhs,
            default={'u': 0, 'v': 0, 'w': 0}
        )

Thank you for your time.

[fdw]

Pretty much. Although inside the constructor you’ll need to allocate some backend matrices (arrays) to hold whatever data you’re interested in passing to the kernel.

Then you’ll need to call _set_external to make those arrays visible to the kernel.

Next, inside the kernel definition itself you’ll need something like

<%pyfr:macro name='bc_rsolve_state' params='ul, nl, ur' externs='ploc, t, your_array'>

where your_array is a fill-in for whatever name you give the extern. Then, inside prepare you will need to call the set method of the matrix you’ve allocated to change its value on the backend.

Regards, Freddie.

[boundaryconditions]

So I put the following in PyFR/pyfr/solvers/navstokes/inters.py :

class  NavierStokesSubInflowFrvNeuralBCInters(NavierStokesBaseBCInters):
    type = 'sub-in-frv-neural'
    cflux_state = 'ghost'

    def __init__(self, be, lhs, elemap, cfgsect, cfg):
        self.backend = be
        super().__init__(be, lhs, elemap, cfgsect, cfg)

        self.c |= self._exp_opts(
            ['rho', 'u', 'v', 'w'][:self.ndims + 1], lhs,
            default={'u': 0, 'v': 0, 'w': 0}
        )

        self.nn_params = self.backend.matrix((1, 2),tags={'align'},dtype=self.backend.fpdtype)
        self._set_external('nn_params', 'fpdtype_t[1][2]',value=self.nn_params)

    def prepare(self, t):
        nn_params = self.nn_params.get()
        nn_params[0][0] = 0.09 # mdot
        #nn_params[0][1] = 0.00 # reserved for future purposes
        self.nn_params.set(nn_params)

and created PyFR/pyfr/solvers/navstokes/kernels/bcs/sub-in-frv-neural.mako with the following content :

<%namespace module='pyfr.backends.base.makoutil' name='pyfr'/>
<%include file='pyfr.solvers.navstokes.kernels.bcs.common'/>

<%pyfr:macro name='bc_rsolve_state' params='ul, nl, ur' externs='ploc, t, nn_params'>
    ur[0] = ${c['rho']};
% for i, v in enumerate('uvw'[:ndims]):
    ur[${i + 1}] = nn_params[0][0]*(${c['rho']}) * (${c[v]});
% endfor
    ur[${nvars - 1}] = ul[${nvars - 1}]
                     - 0.5*(1.0/ul[0])*${pyfr.dot('ul[{i}]', i=(1, ndims + 1))}
                     + 0.5*(1.0/ur[0])*${pyfr.dot('ur[{i}]', i=(1, ndims + 1))};
</%pyfr:macro>

<%pyfr:alias name='bc_ldg_state' func='bc_rsolve_state'/>
<%pyfr:alias name='bc_ldg_grad_state' func='bc_common_grad_zero'/>

Here I was hoping that by multiplying nn_params[0][0] to the conserved variables corresponding to the momentum equation, I am effectively multiplying velocity by this variable. So, in the .ini file, I changed the definition of velocities to

u = 1.0 / (%(omega)s*%(R)s*%(R)s) * x
v = 1.0 / (%(omega)s*%(R)s*%(R)s) * y

as compared to

u = 0.09 / (%(omega)s*%(R)s*%(R)s) * x
v = 0.09 / (%(omega)s*%(R)s*%(R)s) * y

from before. However, I am not getting the expected results. I think the nn_params[0][0] that I set is being accessed by the kernel because if I set this to a large value like 1000, the program crashes. Then, I think the value is being altered somewhere along the way or I put the variable at the wrong location in the mako file. Are you able to see any obvious errors in my code?

[fdw]

@boundaryconditions

fpdtype_t[1][2]

This means that at each flux point there are two distinct values. However, you want every flux point to see the same value. As such I suspect you want broadcast fpdtype_t[1][2]. (I think there are a couple of other topics on the forum regarding this which may be helpful references.)

Regards, Freddie.

[boundaryconditions]

Yes, in broadcast fpdtype_t[1][2] did the trick. For my purpose, the next step is to extract the primitive variables at certain locations to feed it to the machine learning library. Is it OK if I continue discussing things in this thread?

[fdw]

The point sampler plugin contains all of the required machinery for this. In future (unreleased) versions of PyFR this capability will be extracted out into its own class so anyone can use it, but for now you will need to hack the code a little.

Also, do not forget that prepare is being called for all MPI ranks at the same time, so when running in parallel you may need to broadcast/gather/scatter quantities. Note, however, that every ranks gets its own copy of the nn_params array.

Regards, Freddie.

[boundaryconditions]

Thank you, I will look into MPI once I get a serial version running. I managed to get something working here. However, I realize that the point sampling needs to be done only once per timestep (or once every 5000 time steps or so). What currently happens in my code is that the sampling is done once each for each “sub-in-frv-neural” boundaries. My intention is to sample the points once for a time step, feed the samples to a neural network, and then obtain some output which will be stored in the nn_params array. Then, each BC of type sub-in-frv-neural will share the values that have been stored in the nn_params array, but use the shared values in a slightly different way, say one uses +nn_params[0] and the other uses -nn_params[0] (I think this is easily done by flipping the signs of velocities in my .ini file). Therefore, I think I need to do the sampling outside of the boundary condition code (to be more precise, I need to probe values at certain points, and also calculate aerodynamic forces and moments like in the fluidforces plugin). Do you have any suggestions how I could go about doing this?

[fdw]

Here you want a singleton class in Python which multiple boundary condition instances can make use of.

Regards, Freddie.

[boundaryconditions]

So I decided to create a plugin as a common ground for the boundary conditions to get values from. I got some decent result after my neural network model was trained:

My latest boundary condition code looks like this:

class NavierStokesSubInflowFrvNeuralBCInters(NavierStokesBaseBCInters):
    type = 'sub-in-frv-neural'
    cflux_state = 'ghost'

    def __init__(self, intg, be, lhs, elemap, cfgsect, cfg):
        self.backend = be
        self.intg = intg
        super().__init__(be, lhs, elemap, cfgsect, cfg)
        
        # Basic initialization
        self.c |= self._exp_opts(
            ['rho', 'u', 'v', 'w'][:self.ndims + 1], lhs,
            default={'u': 0, 'v': 0, 'w': 0}
        )

        # some config parameters
        self.jcenter = self.backend.matrix((1,1))
        self.polarity = self.backend.matrix((1,1))
        self._set_external('jcenter', 'broadcast fpdtype_t[1][1]', 
                         value=self.jcenter)
        self._set_external('polarity', 'broadcast fpdtype_t[1][1]', 
                         value=self.polarity)

        # Neural network parameters
        self.nn_params = self.backend.matrix((1, 2), tags={'align'})
        self._set_external('nn_params', 'broadcast fpdtype_t[1][2]', 
                         value=self.nn_params)

        # Initial value
        nn_params = self.nn_params.get()
        nn_params[0][0] = 0.0
        self.nn_params.set(nn_params)

        # Fixed values
        jcenter = self.jcenter.get()
        jcenter[0][0] = self.cfg.getfloat(cfgsect, 'jcenter')
        self.jcenter.set(jcenter)
        polarity = self.polarity.get()
        polarity[0][0] = self.cfg.getfloat(cfgsect, 'polarity')
        self.polarity.set(polarity)
        #print(f"jcenter: {jcenter[0][0]}, polarity: {polarity[0][0]}")

        # Time constant for control signal smoothing
        self.alpha = 5e-2

    def prepare(self, t):
        # Direct access to control signal from solver environment
        target = self.intg.system.env.current_control
        
        # Update parameters with smoothing
        nn_params = self.nn_params.get()
        nn_params[0][0] += self.intg._dt/self.alpha * (target - nn_params[0][0])
        self.nn_params.set(nn_params)
        #print(f"Control signal: {target}, Updated nn_params: {nn_params[0][0]}")

Unfortunately, the two lines in the prepare method that have the .get() and .set() commands for the externed variable “nn_params” cause the entire program to slow down by a factor of atleast 4.4 times ! To my understanding this is because of the high number of calls (since prepare is called every time step) and the I/O overhead between the CPU and GPU (I have tested CUDA and HIP backends). I can think of having a CPU-side copy of nn_params to avoid using the .get() command but I don’t see any way to avoid a .set() call (which causes the most overhead as per my measurements). Is there any obvious way I could overcome this overhead?

[fdw]

You can certainly avoid the .get() calls. Secondly, the overhead is likely large because your simulation is proportionally quite small. In a 3D simulation where time steps are a lot more expensive it may well be insignificant. Thirdly, you’ll want to see if you really do need to update the parameters every step. For an explicit simulation it is almost certainly excessive.

Regards, Freddie,

[boundaryconditions]

For now, I have sped things up by calling the .set() function less frequently:

class NavierStokesSubInflowFrvNeuralBCInters(NavierStokesBaseBCInters):
    type = 'sub-in-frv-neural'
    cflux_state = 'ghost'

    def __init__(self, intg, be, lhs, elemap, cfgsect, cfg):
.
.
.
        self._set_external('nn_params', 'broadcast fpdtype_t[1][1]', 
                         value=self.nn_params)
.
.
.
    def prepare(self, t):
        # Direct access to control signal from solver environment
        target = self.intg.system.env.current_control

        new_param = self._current_param + self.intg._dt/self.alpha * (target - self._current_param)

        # Only update backend if change is significant
        if abs(new_param - self._current_param) > 1e-2:
            self.nn_params.set(np.array([[new_param]]))

        self._current_param = new_param

However I wish to change the way the parameters are in the near future. In particular, the parameter is changed by this line

new_param = self._current_param + self.intg._dt/self.alpha * (target - self._current_param)

and this has to be done every timestep. One significant way of speeding things up would be to extern target (needs to be updated less frequently) and alpha (just a constant, broadcast once), and that way I can only broadcast target once in a while. To do this, I need to do the calculation shown in the above equation within the Mako kernel. The only piece I am missing is the time-step size dt. Can this be easily obtained in a manner similar to how the current time t is available to the kernel?

[fdw]

Any use of \Delta t is almost certainly in error.

The prepare method is called for every stage of an RK scheme. These stages occur at times between t and t + \Delta t.

Regards, Freddie.

[boundaryconditions]

Would storing t inside a variable within the kernel and calculating \Delta t based off of it be a good workaround?

[fdw]

No. Any use of \Delta t is fundamentally, numerically, wrong. This is why we prefix it with an _ and do not make it available to kernels. Relaxation time scales need to be based on physics.

Regards, Freddie.

[boundaryconditions]

I see. I do have another option proposed here which doesn’t depend on \Delta t:
\displaystyle Q = \frac{Q_1-Q_0}{T_a}\left(t-t_0\right) + Q_0
This is supposed to change the value Q at different t (in my Mako kernel Q corresponds to control), like described here:

I implemented this using:

class NavierStokesSubInflowFrvNeuralBCInters(NavierStokesBaseBCInters):
    type = 'sub-in-frv-neural'
    cflux_state = 'ghost'

    def __init__(self, intg, be, lhs, elemap, cfgsect, cfg):
        self.backend = be
        self.intg = intg
        super().__init__(be, lhs, elemap, cfgsect, cfg)
        
        # Basic initialization
        self.c |= self._exp_opts(
            ['rho', 'u', 'v', 'w'][:self.ndims + 1], lhs,
            default={'u': 0, 'v': 0, 'w': 0}
        )

        # some config parameters
        self.t_act_interval = self.backend.matrix((1,1))
        self._set_external('t_act_interval', 'broadcast fpdtype_t[1][1]', 
                         value=self.t_act_interval)

        # Neural network + control parameters
        self.control_params = self.backend.matrix((1,3))
        self._set_external('control_params', 'broadcast fpdtype_t[1][3]', 
                         value=self.control_params)

        # Initial value
        self.control_params.set(np.array([[0.0, 0.0, 0.0]])) #(Q0,Q1,t0)

        # Cache current parameter value 
        self._current_target = 0.0

    def prepare(self, t):
        # Direct access to control signal from solver environment
        new_target = self.intg.system.env.current_control

        # Only update backend if there is a new new_target
        if new_target != self._current_target:
            self.t_act_interval.set(np.array([[self.intg.system.env.action_interval]]))
            self.control_params.set(np.array([[self._current_target, new_target, t]]))
            self._current_target = new_target

and the corresponding Mako kernel is:

<%namespace module='pyfr.backends.base.makoutil' name='pyfr'/>
<%include file='pyfr.solvers.navstokes.kernels.bcs.common'/>

<%pyfr:macro name='bc_rsolve_state' params='ul, nl, ur' externs='ploc, t, control_params, t_act_interval'>
    <% control = "(control_params[0][1]-control_params[0][0])/t_act_interval[0][0]*(t-control_params[0][2]) + control_params[0][0]"%>
    ur[0] = ${c['rho']};
% for i, v in enumerate('uvw'[:ndims]):
    ur[${i + 1}] = ${control}*(${c['rho']}) * (${c[v]});
% endfor
    ur[${nvars - 1}] = ul[${nvars - 1}]
                     - 0.5*(1.0/ul[0])*${pyfr.dot('ul[{i}]', i=(1, ndims + 1))}
                     + 0.5*(1.0/ur[0])*${pyfr.dot('ur[{i}]', i=(1, ndims + 1))};
</%pyfr:macro>

<%pyfr:alias name='bc_ldg_state' func='bc_rsolve_state'/>
<%pyfr:alias name='bc_ldg_grad_state' func='bc_common_grad_zero'/>

But does the Mako kernel for boundary condition run at all stages of Runge-Kutta schemes as well? In that case me updating the interface values might be wrong.

[fdw]

Yes, the BC kernels run at every single RK stage.

Regards, Freddie.