DOFs of higher order (>1) basis functions being ordered differently

Hello everyone,

I have a problem with DOFs corresponding to higher order basis functions not getting the same index in the gridfunction each time. In my application, computations can take quite some time, so I’d like to save my solutions, so that they can be visualised later.

The way I’m trying to achieve this is as follows: when I compute the solution, the mesh is saved to a .vol-file and the solution gridfunction is saved to a file containing all its coefficients in the right order. Later, when I’d like to load this solution again, I create a new finite element space (which is a product space in my application) using the mesh loaded from the .vol-file, correct dirichlet boundary settings and correct order. Details about how I save the mesh and solution, and how the finite element space is generated are given below.

The trouble is, the DOFs in the new finite element space corresponding to higher (>1) order basis functions seem to be reshuffled. This makes the coefficients corresponding to those DOFs saved in the gridfunction incorrect.

Is there a way to make sure the DOFs in the new finite element space are ordered in the same manner?

Thanks!

Details about how saving and setting up FE space in my code is implemented:

  • Saving the mesh:
    mesh.ngmesh.Save(‘folder/mesh.vol’)

  • Loading the mesh:
    mesh = ngsolve.Mesh(‘folder/mesh.vol’)

  • Saving the gridfunction (which has many components):
    gf_vec_array = gf.vec.FV().NumPy()
    np.save(filename, gf_vec_array)

  • Loading the gridfunction:
    array = np.load(filename)
    gridfunction.vec.FV().NumPy()[:] = array

  • Setting up FE space: (self.mesh is the mesh loaded in from the .vol-file, or the one generated at the start of the solution procedure and later saved)

U = ngsolve.H1(self.mesh, order=self.order)
G = ngsolve.H1(self.mesh, order=self.order, dirichlet=BOUNDARY_DICT[SEA])

list_of_spaces = [U for _ in range(2*self.M*(2*self.imax + 1))]
for _ in range(2*self.imax+1):
list_of_spaces.append(G)

X = ngsolve.FESpace(list_of_spaces)
self.femspace = X

Example what goes wrong:
When I tested why my loaded solution looked different than the solution I got when first solving the system, I found the following. For a fourth order basis, the first 312 first order basis functions were ordered in the exact same way. The 313th basis function was different.

By executing the following code and plotting right after solving and after loading the solution in a different file, I found the following figures.

empty_gridfunction = ngsolve.GridFunction(fespace)
empty_gridfunction.components[1].vec[313] = 1

After solving:

After loading:

the reason is that a refined (multilevel) mesh counts also edges on the coarser mesh, but Save/Load only keeps the finest mesh. Thus, the number of edges is different, and the ordering of coefficients is messed up.

The recommended way of serialising ngs-objects is pickling, it takes care of all that, see:

https://docu.ngsolve.org/latest/i-tutorials/appendix-pickling/pickling.html

Joachim

Thanks for the quick reply! I’ll try that and see if I can make it work

Unfortunately, pickling does not seem to work for this mesh. I’m getting the error “netgen.libngpy._meshing.NgException: Archive error: Polymorphic type netgen::DiscretePointsSeg<2> not registered for archive”, but I have no clue what that means. Do you know what happens here?

I don’t get the error when I load the mesh from the .vol-file and then try to pickle it. But I do get the error if I try pickling before mesh refinement.

Pickling the mesh also pickles the geometry, such that further geometry-adapted mesh refinement can be done.

It looks like pickling of SplineGeometry2D is only partially supported. Not sure if this will improve, since now the focus is (only) on OCC-Geometry.

You can remove the geometry from the mesh object, something like
mesh.ngmesh.SetGeometry(None)

Joachim

Thanks a lot for the suggestion, but it’s still not really working unfortunately. Right now, I set the geometry to None as you suggested, immediately after generating the mesh. After, I set up the weak forms as usual, and solve.

The solution is saved using

with open('testsol.pkl', 'wb') as file:
    pickle.dump(hydro.alpha_solution[0][1], file)

hydro.alpha_solution[0][1] is a gridfunction. It is immediately loaded again using

with open('testsol.pkl', 'rb') as file:
    alpha = pickle.load(file)

The result is the following:


I guess it’s still something with the ordering of the higher order basis functions, because the general spatial structure of the solution is the same (captured in the linear basis functions) but the higher order corrections are all messed up. I’m a bit at a loss on how to fix it now. I’m also not creating a new FEspace in this example; just plotting the gridfunction after loading it again.

Can you post a minimal example? I suspect some property of the fespace is not serialized correclty.

Hi Cristopher,

I was trying to make a minimal example outside of my model to keep things focused, but every attempt was successful. So here I have an example within my model code, but it runs in a second or two. The dependencies are numpy, matplotlib, sympy, pickle and dill (all accessible with pip).

The code to run is saveload_example.py.

The finite element space is constructed when initialising the Hydrodynamics object, and then in the private method _setup_fem_space. I don’t think there is any issue with the other files, at least not for this problem.

FYI: this problem is not urgent; it has to be fixed eventually, but there are plenty of other things to do for me within this code as well.

saveload_example.zip (165.5 KB)

Thank you!

Cas