GridFunction to numpy array to new GridFunction

Main Forum Old-Forum - do not post here
1

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

2

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:

k_gfu.Interpolate ( f ( c_gfu) )

for interpolation, or

k_gfu.Set ( f ( c_gfu) )

for local L2-projection.

Joachim

3

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

4

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

Attachment: 3366d64d2e56164d4ff06c23e5099587

5

Hi Joachim,

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

Alex

6

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

7

Hi,

the fes can be a product space consisting of different component spaces.

the fesir can be a vectorial IRSpace, with all components having the same points.

When going from the polynomials to the points, you form one vector coefficient function as argument for the gfuir.Interpolate().
When going back to the polynomial space, you can extract components of the ir-gridfunction: gfuir[2:4]

Joachim

8

Hi Joachim,

Thank you for your response. Can you please elaborate on vectorial IRSpace? Do you mean:

fesir = IntegrationRuleSpace([fesir_u, fesir_p])?

9

If you want a vectorial IRSpace with vector dimension 17 you define it as

fes = IntegrationRuleSpace(...) ** 17
10

Dear Joachim,

Thank you for the response. This works. I do have an additional question.

Now that I’ve used IntegrationRuleSpace the data in gfuir.vec.data is equal to the actual u,p values as opposed to polynomial coefficients in the original gfu.vec.data.

However, I am not sure where the spatial data corresponding to these values is stored.

Normally with the original gfu I’d evaluate it like:

gfu(mesh(p1,p2)) or gfu.components0

For gfuir this won’t work.

Is there another way to obtain the corresponding spatial locations of the values?

11

you may create another GridFunction gfxyz in an IntegrationRuleSpace of dimension 3, and then you interpolate the coordinate coefficientfunctions as

gfxzy.Interpolate( CF( (x,y,z) ) )

then the vector entries are the coordinates at the points

12

Thank you, this worked.

13

Dear Joachim,

I just have a follow up question on this. How come the values in gfuir are sampled at the points inside of the coordinate coefficient function? Is there any way I can chose my own points, or, can they be the same points as the original gfu?

14

@joachim

Do .Interpolate and .Set make any differences?
In the following example, they are completely the same.

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=10)

J1 = GridFunction(fes)
J1.Interpolate(sin(x))
J1 = J1.vec.FV()

J2 = GridFunction(fes)
J2.Set(sin(x))
J2 = J2.vec.FV()
15

they are similar, but not the same. both do some kind of interpolation:

gfu.Set does element-wise L2 projection, plus averaging
gfu.Interpolate uses canonical degrees of freedom, if the element provides them.

The usual spaces provide both functions, some may provide only the one or the other. Performance can also differ, but this should rarely be an issue here.

And, there is also a historic difference:
gfu.Interpolate(func, mesh.Boundaries("bnd1"))
does not zero all other boundaries, the Set does so.

Joachim

1 Like