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

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()))

Dokken, Thank you really much for your answer.

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))

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

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

I/O

dokken September 30, 2020, 9:58am 2

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

Topic		Replies	Views
Strange behaviour with write_checkpoint / read_checkpoint I/O	14	1623	March 4, 2020
Access function in time dependent XDMFFile Uncategorized	14	4389	October 26, 2019
Loading xdmf data back in I/O	21	8700	February 15, 2022
Saving an xdmf file - dolfinx v 0.7.0 dolfinx io	6	1823	February 26, 2024
Issues with Checkpointing in Mixed Element Function Spaces (FEniCS Legacy) General legacy-dolfin , checkpointing	2	160	September 25, 2024

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

Related topics