Strange behaviour with write_checkpoint / read_checkpoint

Iactor · March 3, 2020, 11:09am

Hi,
I get a strange behaviour when trying to reload functions from xdmf files.
I have a Function u that is my solution of a linear variational problem and I save it in this way (it has two components representing real and imaginary part)
‘’’
u1, u2 = u.split(deepcopy = False)
with dol.XDMFFile(‘Re.xdmf’) as outfile:
outfile.write_checkpoint(u1, ‘re’)
with dol.XDMFFile(‘Im.xdmf’) as outfile:
outfile.write_checkpoint(u2, ‘im’)
‘’’

I then calculate a series of functionals using

‘’’
M = [(dol.inner(epsilon2(par_sim[‘materials’][i+1], par_sim[‘omega’]) * dol.grad(u[0]), dol.grad(u[0])) * dx(i+1)
+ dol.inner(epsilon2(par_sim[‘materials’][i+1], par_sim[‘omega’]) * dol.grad(u[1]), dol.grad(u[1])) * dx(i+1) ) for i in range(N_subdomains)]
for i in range(N_subdomains):
print (dol.assemble(M[i]))
‘’’
and everything works fine.

However if I reload u from the xdmf files using
‘’’
u = dol.Function(W)
u1 = dol.Function(W_scalar)
u2 = dol.Function(W_scalar)
assigner = dol.FunctionAssigner(W, [W_scalar, W_scalar])

with dol.XDMFFile(‘Re.xdmf’) as infile:
infile.read_checkpoint(u1, ‘re’)
with dol.XDMFFile(‘Im.xdmf’) as infile:
infile.read_checkpoint(u2, ‘im’)
assigner.assign(u, [u1,u2])
‘’’

and I calculate again the functionals I get a nan value just for the first functional.
All the remaining values seem to be the same.
Any idea of what can be wrong

dokken · March 3, 2020, 11:28am

Could you please supply a minimal example that can reproduce this behavior?
Also, is there a particular reason for using deepcopy=False on the first line?

Iactor · March 3, 2020, 12:19pm

Hi,
the code for u is quite long but is essentially solving the complex Poisson equation in a domain with different subdomains.

i get the same behaviour with the simpler functional

M2 = [(dol.inner( dol.grad(u[0]), dol.grad(u[0])) * dx(i+1)
+ dol.inner( dol.grad(u[1]), dol.grad(u[1])) * dx(i+1) ) for i in range(N_subdomains)]

that is essentially the energy norm calculated on each subdomain.
Again the first subdomain (dx(1)) returns a finite result when computed on u as calculated but returns a nan if I use the function loaded from file.

I checked the integration measure by computing the volumes of each subdomain (which I know analytically) and everything is OK.
I also removed deepcopy = False but it has no effect.
Best.
I

dokken · March 3, 2020, 12:34pm

Without a reproducable example it is really tricky to help you. Instead of solving a complex pde, could you just project some functions to a real and imaginary part and see if you have the same issue?

Iactor · March 3, 2020, 2:21pm

Hi, I am working on isolating the problem
I switched to a simpler real formulation of my problem (now the FunctionSpace W is just P1 elements over my mesh)

 u = dol.Function(W)
 problem = dol.LinearVariationalProblem(a, L, u, bc)
 solver = dol.LinearVariationalSolver(problem)
 solver.solve()

I then saved the solution u using

with dol.XDMFFile('potential.xdmf') as outfile:
         outfile.write_checkpoint(u, 'phi')

and reloaded it into another Function taken from the same Functionspace using

v = dol.Function(W)
with dol.XDMFFile( 'potential.xdmf') as infile:
            infile.read_checkpoint(v, 'phi')

I then checked the vertex values of the two functions with

index = np.argwhere(u.compute_vertex_values() != v.compute_vertex_values())
print(index)
print(u.compute_vertex_values()[index])
print(v.compute_vertex_values()[index])

and got

[[18423]]
[[-368.76013452]]
[[nan]]

a single vertex (out of a mesh of 37790 vertices) where the two functions are different, u evaluating to a finite value and v to nan.

Of course this single value is changing the value of the functional calculated over the domain where it is located to nan.

This looks like a bug in the input-output of data.
I will try to find a minimal example for which the problem happens (I tried with simpler problems and everything was ok).
Best.
I

Iactor · March 3, 2020, 3:29pm

Hi again.
stripping a bit of unessential details from my code I get this example.
I am using
dolfin version ‘2019.1.0’
meshio version ‘3.3.1’

import dolfin as dol
import numpy as np
import meshio

class Delta_Function(dol.UserExpression):
    '''Class for mollified Dirac delta function'''
    def __init__(self, x_0, eps):
        '''parameters 
           x_0 : 1D array like. Center
           eps : Float. Standard deviation'''
        self.x_0 = x_0
        self.eps2 = eps**2
        self.dim = len(x_0)
        self.norm = (np.sqrt(2. *np.pi * self.eps2))**self.dim
        dol.UserExpression.__init__(self)
    def eval(self, values, x):
        r2 = np.sum([(x[i]-self.x_0[i])**2 for i in range(self.dim)])
        values[0] = np.exp(-0.5 * r2 /self.eps2)/self.norm
    def value_shape(self):
        return ()
    

def convert_to_XDMF(filename):
    '''reads a .msh file and writes mesh, subdomains, and boundary subdomains 
    to XDMF files'''
    #check if input file i in msh
    if filename.split(sep = '.')[-1] != 'msh':
        raise TypeError('.msh file required')
    msh = meshio.read(filename)
    #write mesh xdmf file
    meshio.write(''.join(filename.split(sep ='.')[:-1]) + '_mesh.xdmf', 
                meshio.Mesh(points = msh.points, 
                            cells = {'tetra': msh.cells['tetra']}))
    #write boundary subdomains xdmf file
    meshio.write(''.join(filename.split(sep ='.')[:-1]) + '_boundary_subdomains.xdmf', 
                meshio.Mesh(points = msh.points, 
                            cells = {'triangle': msh.cells['triangle']},
                            cell_data = {'triangle': {'name_to_read' : msh.cell_data['triangle']['gmsh:physical']}}))
    #write subdomains xdmf file
    meshio.write(''.join(filename.split(sep ='.')[:-1]) + '_subdomains.xdmf', 
                meshio.Mesh(points = msh.points, 
                            cells = {'tetra' : msh.cells['tetra']},
                            cell_data = {'tetra' : {'name_to_read' : msh.cell_data['tetra']['gmsh:physical']}}))

def load_XDMF_files(filename):
    '''load XDMF files'''
    mesh = dol.Mesh()
    with dol.XDMFFile(filename) as infile:
        infile.read(mesh)
    #create subdomain cell function
    mvc = dol.MeshValueCollection('size_t', mesh, 3)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    subdomain_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    #create boundary subdomain facet function
    mvc = dol.MeshValueCollection('size_t', mesh, 2) 
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_boundary_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    return mesh, subdomain_marker, boundary_marker

convert_to_XDMF('cavity.msh')

mesh, subdomain_marker, boundary_marker = load_XDMF_files('cavity_mesh.xdmf')   #define measures
dx = dol.Measure('dx', domain = mesh, subdomain_data = subdomain_marker)
ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)

element = dol.FiniteElement('P', mesh.ufl_cell(), degree = 1)
W = dol.FunctionSpace(mesh, element)
u = dol.TrialFunction(W)
v = dol.TestFunction(W)
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
L = 4. * np.pi * Delta_Function([0.0, 0.0, 0.52], 0.05) * dol.inner(dol.grad(v), dol.Constant((0,0,1))) * dx 
bc = [dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 1),
          dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 2)]

u = dol.Function(W)
problem = dol.LinearVariationalProblem(a, L, u, bc)
solver = dol.LinearVariationalSolver(problem)
solver.solve()

with dol.XDMFFile('potential.xdmf') as outfile:
            outfile.write_checkpoint(u, 'phi')
v = dol.Function(W)
with dol.XDMFFile('potential.xdmf') as infile:
        infile.read_checkpoint(v, 'phi')

index = np.argwhere(u.compute_vertex_values() != v.compute_vertex_values())
print(index)
print(u.compute_vertex_values()[index])
print(v.compute_vertex_values()[index])
print(v.compute_vertex_values().shape)
print(v.compute_vertex_values().shape)

The .msh file is generated using gmsh 4.5.2 for linux and the following .geo file.
Opening it with the gui looks ok and no error is raised anywhere in the code.

SetFactory("OpenCASCADE");

lc = 0.1;
lc_dipole = 0.01;
lc_cavity = 0.004;

//substrate 
//lower face
Point(1) = { -1.0 , -1.0, 0.0, lc };
Point(2) = { 1.0 , -1.0, 0.0, lc };
Point(3) = { 1.0 , 1.0, 0.0, lc };
Point(4) = { -1.0 , 1.0, 0.0, lc };
//upper face
Point(5) = { -1.0 , -1.0, 0.285, lc };
Point(6) = { 1.0 , -1.0, 0.285, lc };
Point(7) = { 1.0 , 1.0, 0.285, lc };
Point(8) = { -1.0 , 1.0, 0.285, lc };

//lower face
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
//vertical lines
Line(5) = {1,5};
Line(6) = {2,6};
Line(7) = {3,7};
Line(8) = {4,8};
//upper face
Line(9) = {5,6};
Line(10) = {6,7};
Line(11) = {7,8};
Line(12) = {8,5};

//lower face
Curve Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};
//vertical faces
Curve Loop(2) = {1,6,-9,-5};
Plane Surface(2) = {2};
Curve Loop(3) = {2,7,-10,-6};
Plane Surface(3) = {3};
Curve Loop(4) = {3,8,-11,-7};
Plane Surface(4) = {4};
Curve Loop(5) = {4,5,-12,-8};
Plane Surface(5) = {5};
//upper face
Curve Loop(6) = {9,10,11,12};
Plane Surface(6) = {6};

//volume
Surface Loop(1) = {1,2,3,4,5,6};
Volume(1) = {1};
//aux points
Point(60) = {-0.4, -0.4, 0.285-0.2, lc};
Point(61) = {0.4, -0.4, 0.285-0.2, lc};
Point(62) = {0.4, 0.4, 0.285-0.2, lc};
Point(63) = {-0.4, 0.4, 0.285-0.2, lc};
Point(64) = {-0.4, -0.4, 0.285, lc};
Point(65) = {0.4, -0.4, 0.285, lc};
Point(66) = {0.4, 0.4, 0.285, lc};
Point(67) = {-0.4, 0.4, 0.285, lc};
Point{60} In Volume{1};
Point{61} In Volume{1};
Point{62} In Volume{1};
Point{63} In Volume{1};
Point{64} In Surface{6};
Point{65} In Surface{6};
Point{66} In Surface{6};
Point{67} In Surface{6};

// cavity
//lower face
Point(9) = { -0.1 , -0.1, 0.285, lc_cavity };
Point(10) = { 0.1 , -0.1, 0.285, lc_cavity };
Point(11) = { 0.1 , 0.1, 0.285, lc_cavity };
Point(12) = { -0.1 , 0.1, 0.285, lc_cavity };
//upper face
Point(13) = { -0.1 , -0.1, 0.295, lc_cavity };
Point(14) = { 0.1 , -0.1, 0.295, lc_cavity };
Point(15) = { 0.1 , 0.1, 0.295, lc_cavity };
Point(16) = { -0.1 , 0.1, 0.295, lc_cavity };

//lower face
Line(13) = {9,10};
Line(14) = {10,11};
Line(15) = {11,12};
Line(16) = {12,9};
//vertical faces
Line(17) = {9,13};
Line(18) = {10,14};
Line(19) = {11,15};
Line(20) = {12,16};
//upper face
Line(21) = {13,14};
Line(22) = {14,15};
Line(23) = {15,16};
Line(24) = {16,13};

//lower face
Curve Loop(7) = {13,14,15,16};
Plane Surface(7) = {7};
//vertical faces
Curve Loop(8) = {13,18,-21,-17};
Plane Surface(8) = {8};
Curve Loop(9) = {14,19,-22,-18};
Plane Surface(9) = {9};
Curve Loop(10) = {15,20,-23,-19};
Plane Surface(10) = {10};
Curve Loop(11) = {16,17,-24,-20};
Plane Surface(11) = {11};
//upper face
Curve Loop(12) = {21,22,23,24};
Plane Surface(12) = {12};
//volume
Surface Loop(2) = {7,8,9,10,11,12};
Volume(2) = {2};

//hbn 
//lower face
Point(17) = { -1.0 , -1.0, 0.295, lc};
Point(18) = { 1.0 , -1.0, 0.295, lc};
Point(19) = { 1.0 , 1.0, 0.295, lc};
Point(20) = { -1.0 , 1.0, 0.295, lc};
//upper face
Point(21) = { -1.0 , -1.0, 0.32, lc};
Point(22) = { 1.0 , -1.0, 0.32, lc};
Point(23) = { 1.0 , 1.0, 0.32, lc};
Point(24) = { -1.0 , 1.0, 0.32, lc};

//lower face
Line(25) = {17,18};
Line(26) = {18,19};
Line(27) = {19,20};
Line(28) = {20,17};
//vertical lines
Line(29) = {17,21};
Line(30) = {18,22};
Line(31) = {19,23};
Line(32) = {20,24};
//upper face
Line(33) = {21,22};
Line(34) = {22,23};
Line(35) = {23,24};
Line(36) = {24,21};

//lower face
Curve Loop(13) = {25,26,27,28};
Plane Surface(13) = {13};
//vertical faces
Curve Loop(14) = {25,30,-33,-29};
Plane Surface(14) = {14};
Curve Loop(15) = {26,31,-34,-30};
Plane Surface(15) = {15};
Curve Loop(16) = {27,32,-35,-31};
Plane Surface(16) = {16};
Curve Loop(17) = {28,29,-36,-32};
Plane Surface(17) = {17};
//upper face
Curve Loop(18) = {33,34,35,36};
Plane Surface(18) = {18};

//volume
Surface Loop(3) = {13,14,15,16,17,18};
Volume(3) = {3};

//aux points
Point(70) = {-0.4, -0.4, 0.295, lc};
Point(71) = {0.4, -0.4, 0.295, lc};
Point(72) = {0.4, 0.4, 0.295, lc};
Point(73) = {-0.4, 0.4, 0.295, lc};

Point(74) = { -0.1 , -0.1, 0.32, lc_cavity };
Point(75) = { 0.1 , -0.1, 0.32, lc_cavity };
Point(76) = { 0.1 , 0.1, 0.32, lc_cavity };
Point(77) = { -0.1 , 0.1, 0.32, lc_cavity };

Point{70} In Surface{13};
Point{71} In Surface{13};
Point{72} In Surface{13};
Point{73} In Surface{13};
Point{74} In Surface{18};
Point{75} In Surface{18};
Point{76} In Surface{18};
Point{77} In Surface{18};

//air
//upper face
Point(25) = { -1.0 , -1.0, 1.12, lc };
Point(26) = { 1.0 , -1.0, 1.12, lc };
Point(27) = { 1.0 , 1.0, 1.12, lc };
Point(28) = { -1.0 , 1.0, 1.12, lc };
//vertical lines
Line(37) = {21,25};
Line(38) = {22,26};
Line(39) = {23,27};
Line(40) = {24,28};
//upper face
Line(41) = {25,26};
Line(42) = {26,27};
Line(43) = {27,28};
Line(44) = {28,25};

//vertical faces
Curve Loop(19) = {33,38,-41,-37};
Plane Surface(19) = {19};
Curve Loop(20) = {34,39,-42,-38};
Plane Surface(20) = {20};
Curve Loop(21) = {35,40,-43,-39};
Plane Surface(21) = {21};
Curve Loop(22) = {36,37,-44,-40};
Plane Surface(22) = {22};
//upper face
Curve Loop(23) = {41,42,43,44};
Plane Surface(23) = {23};
//volume
Surface Loop(4) = {18,19,20,21,22,23};
Volume(4) = {4};

//auxiliary points
Point(50) = {0.0, 0.0, 0.52, lc_dipole};
Point(51) = {-0.4, -0.4, 0.52+0.2, lc};
Point(52) = {0.4, -0.4, 0.52+0.2, lc};
Point(53) = {0.4, 0.4, 0.52+0.2, lc};
Point(54) = {-0.4, 0.4, 0.52+0.2, lc};
Point{50} In Volume{4};
Point{51} In Volume{4};
Point{52} In Volume{4};
Point{53} In Volume{4};
Point{54} In Volume{4};

v() = BooleanFragments{ Volume{1}; Delete;}{ Volume{2,3,4}; Delete; };

Physical Volume('substrate', 1) = {#v()-3};
Physical Volume('cavity', 2) = {#v()-2};
Physical Volume('hBN', 3) = {#v()-1};
Physical Volume('air', 4) = {#v()};

Physical Surface('backgate', 1) = {24};
Physical Surface('patterned_gate', 2) = {30, 8, 9, 10, 11, 29};

Thanks to whoever can help.
Iacopo

dokken · March 3, 2020, 6:10pm

I cannot reproduce your issue. I separated the mesh generation to a separate file, calling the mesh i load in mesh_mesh.xdmf.

import dolfin as dol
import numpy as np
import meshio

class Delta_Function(dol.UserExpression):
    '''Class for mollified Dirac delta function'''
    def __init__(self, x_0, eps):
        '''parameters
           x_0 : 1D array like. Center
           eps : Float. Standard deviation'''
        self.x_0 = x_0
        self.eps2 = eps**2
        self.dim = len(x_0)
        self.norm = (np.sqrt(2. *np.pi * self.eps2))**self.dim
        dol.UserExpression.__init__(self)
    def eval(self, values, x):
        r2 = np.sum([(x[i]-self.x_0[i])**2 for i in range(self.dim)])
        values[0] = np.exp(-0.5 * r2 /self.eps2)/self.norm
    def value_shape(self):
        return ()


def load_XDMF_files(filename):
    '''load XDMF files'''
    mesh = dol.Mesh()
    with dol.XDMFFile(filename) as infile:
        infile.read(mesh)
    #create subdomain cell function
    mvc = dol.MeshValueCollection('size_t', mesh, 3)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    subdomain_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    #create boundary subdomain facet function
    mvc = dol.MeshValueCollection('size_t', mesh, 2)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_boundary_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    return mesh, subdomain_marker, boundary_marker


mesh, subdomain_marker, boundary_marker = load_XDMF_files('mesh_mesh.xdmf')   #define measures
dx = dol.Measure('dx', domain = mesh, subdomain_data = subdomain_marker)
ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)

element = dol.FiniteElement('P', mesh.ufl_cell(), degree = 1)
W = dol.FunctionSpace(mesh, element)
u = dol.TrialFunction(W)
v = dol.TestFunction(W)
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
L = 4. * np.pi * Delta_Function([0.0, 0.0, 0.52], 0.05) * dol.inner(dol.grad(v), dol.Constant((0,0,1))) * dx
bc = [dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 1),
          dol.DirichletBC(W, dol.Constant(0.), boundary_marker, 2)]

u = dol.Function(W)
problem = dol.LinearVariationalProblem(a, L, u, bc)
solver = dol.LinearVariationalSolver(problem)
solver.solve()

with dol.XDMFFile('potential.xdmf') as outfile:
            outfile.write_checkpoint(u, 'phi')
v = dol.Function(W)
with dol.XDMFFile('potential.xdmf') as infile:
        infile.read_checkpoint(v, 'phi')

index = np.argwhere(u.compute_vertex_values() != v.compute_vertex_values())
print(index)
print(u.compute_vertex_values()[index])
print(v.compute_vertex_values()[index])
print(v.compute_vertex_values().shape)
print(v.compute_vertex_values().shape)

gives me

[]
[]
[]
(50129,)
(50129,)

So are you sure you are running the latest version of dolfin?

Iactor · March 3, 2020, 7:54pm

Very strange.
You are getting a different number of vertices for the function, mine had 37790.
Which version of gmsh are you using?
Maybe I can put my .msh file somewhere and attach a link.
I

dokken · March 3, 2020, 8:44pm

Im used gmsh 4.4.1 for this example

Iactor · March 4, 2020, 8:43am

Ok, I tried to mesh with gmsh 4.4.1 (linux 64).
I also ran gmsh -check and everything is ok (as was ok with the other version).
Did you give any special option to gmsh? I used just gmsh -3.
The result is different in terms of number of nodes of the final mesh (but still different from yours).

The output of my code is (I added a command to print the version of the modules I am importing and corrected the last instruction that was meant to compare the number of nodes of u and v and was printing twice the number of nodes of u)

dolfin version = 2019.1.0
numpy version = 1.18.1
meshio version = 3.3.1
/usr/lib/python3/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
[[16702]]
[[0.]]
[[nan]]
(47432,)
(47432,)

I see only this deprecation warning (that appears during the conversion to XDMF) but everything else looks ok.

Shall I file a bug report?
Thanks.
I

dokken · March 4, 2020, 10:01am

I was able to reproduce the problem now. But it seems to be mesh related, as there are nan-values in the original u-vector. I’ve reduced your example to the following:

import dolfin as dol
import numpy as np

def load_XDMF_files(filename):
    '''load XDMF files'''
    mesh = dol.Mesh()
    with dol.XDMFFile(filename) as infile:
        infile.read(mesh)
    #create subdomain cell function
    mvc = dol.MeshValueCollection('size_t', mesh, 3)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    subdomain_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    #create boundary subdomain facet function
    mvc = dol.MeshValueCollection('size_t', mesh, 2)
    with dol.XDMFFile(''.join(filename.split(sep ='_')[:-1]) + '_boundary_subdomains.xdmf') as infile:
        infile.read(mvc, 'name_to_read')
    boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    return mesh, subdomain_marker, boundary_marker


mesh, subdomain_marker, boundary_marker = load_XDMF_files('mesh_mesh.xdmf')   #define measures
# mesh = dol.UnitSquareMesh(10,10)
element = dol.FiniteElement('P', mesh.ufl_cell(), degree = 1)
W = dol.FunctionSpace(mesh, element)

x = dol.SpatialCoordinate(mesh)
u = dol.project(x[0]*x[1], W)

with dol.XDMFFile('potential.xdmf') as outfile:
    outfile.write_checkpoint(u, 'phi',0.0, dol.XDMFFile.Encoding.HDF5, False)
v = dol.Function(W)
with dol.XDMFFile('potential.xdmf') as infile:
    infile.read_checkpoint(v, 'phi', 0)

u_array = u.vector().get_local()
v_array = v.vector().get_local()
print(np.argwhere(np.isnan(u_array)),np.argwhere(np.isnan(v_array)))
# assert np.allclose(u_array, v_array)

The print yields:

[[47429]
 [47430]
 [47431]] [[44291]]

Thus there is something with the original mesh that causes issues in checkpointing.

dokken · March 4, 2020, 10:05am

By changing the resolution in your mesh file:

lc = 0.1;
lc_dipole = 0.02;
lc_cavity = 0.008;

or

lc = 0.1;
lc_dipole = 0.01;
lc_cavity = 0.0055;

the code i posted above runs

Iactor · March 4, 2020, 11:00am

Trying to sum up
-the original u vector u.vector().get_local() has three nans at positions [[47429], [47430], [47431]] (i get the same on my pc with gmsh 4.4.1)
-the loaded function v has instead only one nan at [[44291]]
-beside this the two vectors are equal

-the vertex values of u (u.compute_vertex_values() ) do not have nans
-the vertex values of v have one nan at position [[16702]]
-again this is the only difference between the two vectors

-the problem disappears when changing the mesh
-no warning or error is raised apart from a deprecation warning in converting the mesh

Do you have any idea of what is going on?
I am quite lost.
Thanks a lot for the support.
I

dokken · March 4, 2020, 11:07am

Loading the mesh does not do quality checks, so no errors will be raised. I guess not the load_checkpoint and compute_vertex_values does something buggy when encountering a nan-value. As this boils down to a meshing specific question, I’m not sure if I can help much more. Then you would need to reproduce this behavior on a small mesh. (As inspecting the current mesh by hand is impossible).

Iactor · March 4, 2020, 11:23am

Yes, what I am really surprised about is the LinearSolver (or in your case the project operator) that runs without problems but returns a nan in the computed solution (in its vector, not in the vertex values).
Actually this may be the difference between the original and the loaded solution.
Probably save_checkpoint and load_checkpoint just work with the dofs and not with vertex_values (althought for P1 elements are the same apart from reordering) that are computed only when the function is called.

Topic		Replies	Views
XDMF write_checkpoint questions I/O	3	2320	April 3, 2023
Loading xdmf data back in I/O	21	8025	February 15, 2022
Write and read XDMF I/O	2	1132	March 22, 2020
Does XDMF read_checkpoint avaliable in DOLFINX I/O io	5	172	July 2, 2024
Read_checkpoint written with older dolfin version I/O	0	389	June 22, 2021

Strange behaviour with write_checkpoint / read_checkpoint

Related topics