Adios4dolfinx misreads written function

I/O
1

Hi,

I am trying to solve wave equation using a mesh generated with gmsh. When adios reads the exported function, it mismatches with the solution. Here is MWE;

from ufl import FiniteElement, Measure, inner, grad, TrialFunction, TestFunction, dx
from dolfinx.fem import form, Function, FunctionSpace
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.io.gmshio import model_to_mesh
from slepc4py import SLEPc
from petsc4py import PETSc
from pathlib import Path
from mpi4py import MPI
import adios4dolfinx
import numpy as np
import dolfinx
import gmsh
import sys
import os

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add(__name__)

cube = gmsh.model.occ.addBox(0,0,0, 1, 1, 1, tag=1)

gmsh.model.occ.synchronize()

gmsh.option.setNumber('Mesh.MeshSizeMin', 1)

gmsh.model.mesh.generate(3)

boundaries_entities = [i[1] for i in gmsh.model.getEntities(2)]
for tag, surf in enumerate(boundaries_entities):
    gmsh.model.addPhysicalGroup(2,[surf],tag)

gmsh.model.addPhysicalGroup(3,[1],1)

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

path = os.path.dirname(os.path.abspath(__file__))

mesh, subdomains, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)

degree = 2

v_cg = FiniteElement("CG", mesh.ufl_cell(), degree)
V = FunctionSpace(mesh, v_cg)

c = 340

u, v = TrialFunction(V), TestFunction(V)
A = assemble_matrix(form(-c**2* inner(grad(u), grad(v))*dx))
A.assemble()
C = assemble_matrix(form(inner(u , v) * dx))
C.assemble()

target = PETSc.ScalarType(2*np.pi*180)

E = SLEPc.EPS().create(MPI.COMM_WORLD)
C = - C
E.setOperators(A, C)
st = E.getST()
st.setType('sinvert')
E.setTarget(target**2)
E.setFromOptions()
E.solve()

A = E.getOperators()[0]
eig = E.getEigenvalue(0)
print("eig: ", np.sqrt(eig))
p_dir = Function(V)
p_dir_i = Function(V)
E.getEigenvector(0, p_dir.vector, p_dir_i.vector)

path_out = Path("adios_dir/checkpoint.bp")

# write it
adios4dolfinx.write_mesh(mesh, path_out, engine="BP4")
adios4dolfinx.write_function(p_dir, path_out, engine="BP4")

# read it back
mesh_import = adios4dolfinx.read_mesh(MPI.COMM_WORLD, path_out, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
v_cg = FiniteElement("CG", mesh_import.ufl_cell(), degree)
V_import = FunctionSpace(mesh_import, v_cg)
p_dir_import = Function(V_import)
adios4dolfinx.read_function(p_dir_import, path_out, engine="BP4")

np.testing.assert_allclose(p_dir.x.array, p_dir_import.x.array, atol=1e-14)

gives;

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=1e-14

Mismatched elements: 41 / 63 (65.1%)
Max absolute difference: 0.34638213
Max relative difference: 2.94448321e+15
 x: array([ 8.114930e-16-7.606003e-15j, -7.413879e-02+1.565202e-01j,
       -7.519621e-02+1.587526e-01j, -7.413879e-02+1.565202e-01j,
       -7.005317e-02+1.478948e-01j, -6.627179e-02+1.399116e-01j,...
 y: array([ 8.114930e-16-7.606003e-15j, -7.413879e-02+1.565202e-01j,
       -7.519621e-02+1.587526e-01j, -7.413879e-02+1.565202e-01j,
       -7.005317e-02+1.478948e-01j, -6.627179e-02+1.399116e-01j,...

I could not figure out what is wrong here. I am using dolfinx v0.7.3 and adios4dolfinx v0.7.2.

2

I feel like we have discussed this before. Here you assume that adios4dolfin reads in the mesh in the same way as it was created. Please use “original_checkpoint” i.e. write_function_on_input_mesh for this.
See: Checkpoint on input mesh — ADIOS2Wrappers
for a detailed description. This means that you actually do not need to store the mesh, but you can read the data directly in on your mesh object and V.

1 Like