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

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.

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.

1 Like

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.

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.

1 Like

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?

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.

1 Like

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?

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.

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?