I am working with an advection-diffusion equation with source terms to model biological growth.

the source terms is of the type f_i = k_i * c_i

where i is the species, f is the whole source term, c is the concentration, and k is the coefficient.

That coefficient k is spatially varying. It is the result of solving a linear optimization problem for a given set of concentrations, so we do not know it ahead of time or have a closed-form expression for it.

We are solving for k explicitly at each point in time, using the previous time-stepâ€™s concentration to calculate itâ€™s values. Thus at each time-step it is spatially varying, but constant. The concentrations, c_i, however are the current timestepâ€™s concentration, and thus f is solved implictly.

The concentrations are stored on a H1 3rd order FES. I have chosen to store the kâ€™s on the same finite element space so that they are defined at the same points, in order to minimize any interpolation or anything similar. I am doing this because the source term will be re-calculated at every time-step and I want it to be as fast as possible.

I have the concentrations for the previous time-step as gridfunctions (c_gfu), and I can access their data in the form of a numpy array using either c_np = np.array(c_gfu.vec.data[:]) or c_np = c_gfu.FV().NumPy(). Both seem to do the same thing.

The issue Iâ€™m having is that I think that c_np actually has coefficients for the interpolating polynomials rather than the values at the points where itâ€™s solved. So on a very coarse mesh I set c_gfu to a initial condition of c=4, I find that only c_np[0:12] (or some small subset of all of the values in c_np) is equal to 4, the rest are equal to ~0 (really 1e-16). If I use a spatially varying initial condition than all of the values are non-zero.

My plan was to write a function â€śfâ€ť which performed the following:

k_np = f(c_np)

and then k_gfu.vec.data[:] = k_np.

But this does not work if c_np and k_np do not contain the values of the nodes, and instead contain coefficients.

Am I trying to work with these in a wrong fashion?

What can I do in order to achieve what Iâ€™m trying to do?

Alex