GridFunction to numpy array to new GridFunction

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

Hi Alex

your observations are right, the vector coefficients are not the point values of the function, i.e. we are not using a nodal basis (for p >= 2).

You may want to use:
for interpolation, or

for local L2-projection.

Joachim

Hi Joachim,

Thanks for the confirmation and the suggestion.

Unfortunately the underlying optimization problem will only accept individual floats, or numpy array of floats. I.e. f(c_np) where c has to be floats or bumpy arrays.

It is not a simple algebraic equation in terms of c_np so I don't believe that I can simply replace c_np with c_gfu.

For a bit more context:

From above: rxn_i = k_i * c_i

At every point in space, k_i(x) = F(c_1(x), c_2(x), ..., c_N(x))

As you said, if I were working with p=1, then I could read nodal values directly and calculate k_i at each node by converting c_gfu_i into c_np_i and then working directly with c_np_i.

However I am using higher order elements, and thus c_gfu_i does not contain point values.

I am attempting to create a spatially varying, even within each mesh element, field for the k_i's in order to best capture the spatial variation in reaction rate.



Perhaps I am not taking the best approach here.

Is there any way to go from c_gfu -> c_np (which contains values) -> k_gfu?

Alex

Hi Alex,

a relatively new option is to use IntegrationRuleSpace. It provides one value per integration point, but you can choose the points arbitrary within the reference element.

Going from H1 to IRSpace is very simple by interpolation. One possibility to go back is global L2 - projection (but could be optimized to have a local operation).

Hope this gives further ideas to solve your problem.

Joachim

Hi Joachim,

I realized I never replied. Thank you! That does the job perfectly!

Alex

Hi Joachim, thank you for providing this example code. I was wondering if this can be used for Gridfunctions with mixed finite element spaces.