Problem on storing and reading sub-solution (as the norm difference between the initial and read solution is not near 0)

I/O
1

Hi everyone,

I have a problem in post-processing my solution (maybe trivial) that I cannot solve.

So, I want to save and load a FEM “sub-solution” extracting from a spliting like u below

w = Function(W) // where W is a mixed space (velocity and pressure for Stokes)

solve(a == L, w, bc, solver_parameters={'linear_solver':'mumps'});

(u, p) = w.split()

I don’t have any compilation error (in my weak formulation and the saving/reading) but when I try to use the “erronorm” function on the difference u^{compute}-u^{read}, I don’t have a null L_2-norm:

 ############################
 # save and read "u" from the w splitting
 ############################
  V1 = VectorFunctionSpace(mesh, "CG", 2)

 # Save solution
output_file = HDF5File(mesh.mpi_comm(), "/home/ConvergenceMaillage/maillage/u.h5", "w")
output_file.write(u11, "solution")
output_file.close()

# Load solution
 U = Function(V1)
 input_file = HDF5File(mesh.mpi_comm(),  "/home/ConvergenceMaillage/maillage/u.h5", "r")
 input_file.read(U, "solution")
 input_file.close()

E =errornorm(U,u11, norm_type='L2', degree_rise=4)

This should be at most equal to an epsilon because they are the same solution but I get somehting like 1.687.
Note that I have ||u^{compute}-u^{compute}||_{L_2}=0 (same for u^{read}).

Do I need to interpolate on a special space (as I am on a mixed-space W and only want the solution of W.sub(0)=V=velocity component) ? Is this problem coming from my method to save and read the solution ?

Thank you in advance for your answer.

Adel

2

You should use write_checkpoint and read_checkpoint: Loading xdmf data back in

3

Dokken,

Thank you really much for your answer.

4

Dokken,

I’m really sorry to bother you with a well documented issue (cf. all the other thread) but I dont understant why this is not working in my case.

I did like you said : using checkpoint (.read and .write).

Can you please correct me if you see any error on the post-processing :

from __future__ import print_function
from dolfin import *
from fenics import *
from dolfin.cpp.mesh import *
from mshr import*
import sys
import meshio
import matplotlib.pyplot as plt
import dolfin as do
import numpy as np
from numpy import linalg as la
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp, Measure,\
DirichletBC, FunctionSpace, Constant, TrialFunction, \
TestFunction, dot, grad, dx, Function, solve

mesh = Mesh()

with XDMFFile(sys.argv[1]) as infile:
infile.read(mesh)
File("Dolfin_circle_mesh.pvd").write(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile(sys.argv[2]) as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("Dolfin_circle_facets.pvd").write(mf)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

tol = 1E-14

boundaries = MeshFunction("size_t", mesh, 
mesh.topology().dim()-1)

ds = Measure("ds")(subdomain_data=boundaries)

ds_G = Measure("ds", domain=mesh, subdomain_data=mf, 
subdomain_id=100001)
ds_D = Measure("ds", domain=mesh, subdomain_data=mf, 
subdomain_id=100003)
ds_B = Measure("ds", domain=mesh, subdomain_data=mf, 
subdomain_id=100002)
ds_H = Measure("ds", domain=mesh, subdomain_data=mf, 
subdomain_id=100004)
ds_O = Measure("ds", domain=mesh, subdomain_data=mf, 
subdomain_id=100000)

n1 = FacetNormal(mesh)
f1 = Constant((0, 0))
mu1 = 0.001
mu = Constant(mu1)
rho = Constant(1000)

bc_os1 = DirichletBC(W.sub(0), (0.0,0.0), mf, 100000)

bc_noflow1 = DirichletBC(W.sub(0).sub(1), 0.0, mf, 100004)
bc_2noflow1 = DirichletBC(W.sub(0).sub(1), 0.0, mf, 100002)
bcp_outflow1 = DirichletBC(W.sub(1), 0.0, mf, 100003)

bc11 = [bc_2noflow1,bc_noflow1, bc_os1,bcp_outflow1]

a11 = mu1*inner(grad(u),grad(v))*dx(mesh)\
- div(v)*p*dx(mesh) - q*div(u)*dx(mesh)
L11 = -dot(1*v,n1)*ds_G

w11 = Function(W)
solve(a11 == L11, w11, bc11, solver_parameters= 
{'linear_solver':'mumps'});
(u11, p11) = w11.split()

V1 = VectorFunctionSpace(mesh, "CG", 2)

with XDMFFile(MPI.comm_world, "/home/amoreno3/ConvergenceMaillage/maillage/vitesseDokken7.xdmf") as outfile:
 outfile.write_checkpoint(u11, "u11", 0, append=False)

f1 = Function(V1)

with XDMFFile(MPI.comm_world, "/home/amoreno3/ConvergenceMaillage/maillage/vitesseDokken7.xdmf") as infile:
infile.read_checkpoint(f1, "u11")

E =errornorm(f1,u11, norm_type='L2', degree_rise=4)

print(E)
print(len(u11))
print(len(f1))
5

For further issues, please make a minimal working code example (an example I can run on my computer and reproduce the issue). Please make sure to remove all unecessary imports. You are import several things multiple times in your code, (from dolfin import *, from fenics import *, from dolfin import Mesh …) all does the same thing. You also import dolfin as do.

Please use the example below as a guideline on minimal working code examples:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10,10)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)
w = Function(W)

# Create functionspaces for each of the elements
V1 = FunctionSpace(mesh, V)
Q1 = FunctionSpace(mesh, Q)
f = FunctionAssigner(W, [V1,Q1])

# Create a non-zero velocity profile
u_vel = project(Expression(("x[0]", "x[1]"),degree=2), V1)

# Assign velocity profile to w
print(w.vector().norm("l2"))
f.assign(w, [u_vel, Function(Q1)])
print(w.vector().norm("l2"))

u, p = w.split()
# Save checkpoint of u
with XDMFFile(MPI.comm_world, "vitesseDokken7.xdmf") as outfile:
    outfile.write_checkpoint(u, "u", 0, append=False)

# Load checkpoint of u 
f1 = Function(V1)
# Removing this loading will make the test fail
with XDMFFile(MPI.comm_world, "vitesseDokken7.xdmf") as infile:
    infile.read_checkpoint(f1, "u")

# Compare loaded function with original function
E =errornorm(f1,u, norm_type='L2', degree_rise=4) 
u_org = w.split(deepcopy=True)[0].vector().get_local()
assert(np.allclose(u_org, f1.vector().get_local()))
6

Thanks for the demo. That’s work now.
Have a good day