Question about XDMFFile in dolfin

I/O
1

Hello,

I am using dolfin and would like to restart a simulation from a previous one. I am storing the data of my previous simulation using XDMFFile class. And I use the write(u, t) method to write my solution on the xdmf file where u is my solution (class Function) and t is a float (the simulation time). So I launch my simulation, and get the result in my xdmf file as expected. Then I restart the simulation, and write on the same file. Say my ending time was t = 10.0, and my simulation time step is dt = 0.1 then the next time step is t = 10.1. So when I restart the simulation from the previous data file, my script automatically save the next solution u on t = 10.1. So I would expect that my xdmf file contains the solution from t = 0 to t = 10.0 and also the solution at t = 10.1. But it only contains the solution at t = 10.1, and all the previous time steps are erased. I didn’t managed to understand how this class works, so I created this post. Could someone tell me how to happen data in XDMFFile without erasing previous data ?

Thanks and regards,
TOP1

2

A minimal working example would help us debug the issue. Are you sure you set the optional argument in XDMFFile.write_checkpoint, append=True?

3 
from dolfin import *
from mpi4py import MPI
import meshio

export_path = "/home/your_path/"

mesh_comm = MPI.COMM_WORLD

from solver_gym_fenics import create_mesh

mesh = Mesh()

# The mesh generated with gmsh's python API
mesh_file = "cylinder_mesh.msh"

# xdmf file that will contain the mesh
xdmf_file = "mesh.xdmf"

# xdmf file that will contain the 1D elements mesh (facets)
facet_file = "facet_mesh.xdmf"

# read mesh from msh file
msh = meshio.read(mesh_file)

# create the 1D element mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write(facet_file, line_mesh)

# create the xdmf mesh
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write(xdmf_file, triangle_mesh)

# writing xdmf file into Mesh object
with XDMFFile(xdmf_file) as infile:
	infile.read(mesh)

# retrieving the facets
mvc1 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(facet_file) as infile:
	infile.read(mvc1, "name_to_read")
facet = cpp.mesh.MeshFunctionSizet(mesh, mvc1)




velocity_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
pressure_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

W = FunctionSpace(mesh, velocity_element * pressure_element)

w = Function(W)
wall_marker, inlet_marker, outlet_marker = range(1, 4)
obstacle_marker = 4


solution_file = XDMFFile(mesh_comm, export_path + "mwe.xdmf")
restart = False
if not restart: 
	i = 0
	bc = DirichletBC(W.sub(0), Constant((1, -1)), facet, wall_marker)
	solution_file.write(mesh)
else:
	i = 1
	bc = DirichletBC(W.sub(0), Constant((-1, 1)), facet, wall_marker)

bc.apply(w.vector())




solution_file.write_checkpoint(w.sub(0), "mixed_function", i, XDMFFile.Encoding.HDF5, True)
# solution_file.write(w, i)
solution_file.close()

Here you can find my mwe. I launch it, it creates the mesh file and a solution file. But there is an error in paraview, and I can’t know why. I hope you understood what I wanted to do in my post, so that you can explain me what I am doing wrong with the XDMFFile class.

This mwe needs a cylinder_mesh file, which is created using this script:

############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy

from env_data import r, c_x, c_y, gdim


from solver_gym_fenics import printr0

# Calculating the time that mesh generation will take
tic = time.perf_counter()

gmsh.initialize()

c_x, c_y = c_x[0], c_y[0]

gdim = 2

gmsh.model.add("cylinder")

# target mesh size for any point
lc = 0.05
# lc = 0.1

############################################
#---------- Geometry creation -------------#
############################################
## 0 dim : Points definition
# Rectangle delimiting the flow
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(2.2, 0, 0, lc, 2)
gmsh.model.geo.addPoint(2.2, 0.41, 0, lc, 3)
gmsh.model.geo.addPoint(0, 0.41, 0, lc, 4)

# Points to create the circle for the cylinder (4 points + center)
gmsh.model.geo.addPoint(c_x, c_y, 0, lc, 5) # center
gmsh.model.geo.addPoint(c_x + r, c_y, 0, lc, 6)
gmsh.model.geo.addPoint(c_x, c_y + r, 0, lc, 7)
gmsh.model.geo.addPoint(c_x - r, c_y, 0, lc, 8)
gmsh.model.geo.addPoint(c_x, c_y - r, 0, lc, 9)

# Points to create the polygon (pentagon in this case) that will include the cylinder and the wake
# In that pentagon, the mesh will be refined
gmsh.model.geo.addPoint(c_x - 2*r, c_y, 0, lc, 10)
gmsh.model.geo.addPoint(c_x, c_y + 2*r, 0, lc, 11)
gmsh.model.geo.addPoint(2.2, c_y + 2*r, 0, lc, 12)
gmsh.model.geo.addPoint(2.2, c_y - 2*r, 0, lc, 13)
gmsh.model.geo.addPoint(c_x, c_y - 2*r, 0, lc, 14)

## 1 dim: Lines and Curves
# Lines for the rectangle
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 13, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addLine(12 , 3, 5)

# Lines for the pentagone
gmsh.model.geo.addLine(10, 11, 6)
gmsh.model.geo.addLine(11, 12, 7)
gmsh.model.geo.addLine(12, 13, 8)
gmsh.model.geo.addLine(13, 14, 9)
gmsh.model.geo.addLine(14, 10, 10)

# Creation of the circle for the cylinder (4 CircleArcs joining the 4 points previously defined)
gmsh.model.geo.addCircleArc(6, 5, 7, tag=11)
gmsh.model.geo.addCircleArc(7, 5, 8, tag=12)
gmsh.model.geo.addCircleArc(8, 5, 9, tag=13)
gmsh.model.geo.addCircleArc(9, 5, 6, tag=14)


## 2 dim
# Creation of the plane surface: rectangle - pentagone
gmsh.model.geo.addCurveLoop([1, 2, 9, 10, 6, 7, 5, 3, 4], 1)
gmsh.model.geo.addCurveLoop(list(range(6, 11)), 2)
gmsh.model.geo.addPlaneSurface([1], 1)


############################################
#------ Parameters & Physical groups ------#
############################################
# Creation of a curve for the polygon
curved_loop_list = [2]

# Creation of curves loops for the cylinder
gmsh.model.geo.addCurveLoop([11, 12, 13, 14], 3)
curved_loop_list.append(3)

# Creation of the plane surface: polygon - cylinder
gmsh.model.geo.addPlaneSurface(curved_loop_list, 2)

gmsh.model.geo.synchronize()

## Fields
# Distance field to calculate distance from the three cylinders
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "CurvesList", list(range(11, 15)))
gmsh.model.mesh.field.setNumber(1, "Sampling", 100)

# Threshold field to refine the mesh around the three cylinders
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "SizeMin", lc / 9)
gmsh.model.mesh.field.setNumber(2, "SizeMax", lc)
gmsh.model.mesh.field.setNumber(2, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(2, "DistMax", r)

# Constant mesh size inside polygon
gmsh.model.mesh.field.add("Constant", 3)
gmsh.model.mesh.field.setNumber(3, "VIn", lc/6)
gmsh.model.mesh.field.setNumber(3, "VOut", lc)
gmsh.model.mesh.field.setNumbers(3, "SurfacesList", [2])

# Distance field to calculate distance from the outlet
gmsh.model.mesh.field.add("Distance", 4)
gmsh.model.mesh.field.setNumbers(4, "CurvesList", [2, 8, 5])
gmsh.model.mesh.field.setNumber(4, "Sampling", 100)

# Threshold field to refine the mesh around the three cylinders
gmsh.model.mesh.field.add("Threshold", 5)
gmsh.model.mesh.field.setNumber(5, "InField", 4)
gmsh.model.mesh.field.setNumber(5, "SizeMin", lc/2)
gmsh.model.mesh.field.setNumber(5, "SizeMax", lc)
gmsh.model.mesh.field.setNumber(5, "DistMin", r/2)
gmsh.model.mesh.field.setNumber(5, "DistMax", r)

# Final field as in t10.py, take the minimum mesh size for all Threshold fields.
gmsh.model.mesh.field.add("Min", 6)
gmsh.model.mesh.field.setNumbers(6, "FieldsList", [2, 3, 5])
gmsh.model.mesh.field.setAsBackgroundMesh(6)

## Create Physical Groups
gmsh.model.addPhysicalGroup(1, [1, 3], name="walls") # walls: up and down sides of the rectangle
gmsh.model.addPhysicalGroup(1, [4], name="inlet") # inlet: left side of the rectangle
gmsh.model.addPhysicalGroup(1, [2, 8, 5], name="outlet") # outlet: right side of the rectangle
gmsh.model.addPhysicalGroup(1, list(range(11, 15)), name="cylinder") # obstacel
gmsh.model.addPhysicalGroup(2, [1, 2], name="control_volume", tag=5) # Reunion of the two plane surfaces: where fluid will flow

gmsh.model.geo.synchronize()

# printing how long the mesh creation took
toc = time.perf_counter()
print("Mesh generation took {} secondes !".format(toc-tic))


############################################
#----- Mesh Saving & Post-processing ------#
############################################
# Generate and write to file
gmsh.model.mesh.generate(gdim)
gmsh.write("cylinder_mesh.msh")
gmsh.write("cylinder_mesh.vtk")

# Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

# retrieving points of the fluid domain (our physical group in gmsh)
points = gmsh.model.mesh.getNodesForPhysicalGroup(2, 5)[1]

# getting the coordinates of the mesh points
coordinates = []
for i in range(0, len(points), 3):
	coordinates.append(list(points[i:i+3]))
two_D_coord = []
for coord in coordinates:
	two_D_coord.append(coord[:-1])

# saving the coordinates a file
with open('mesh_data.npy', 'wb') as f:
    numpy.save(f, two_D_coord)

So what I would like to do:

  1. create the mesh using the second script.
  2. launch the first script with restart = False
  3. launch the first script again with restart = True.

And I would like to have a solution_file.xdmf containing 2 steps of the solution, one created when restart = False, and the other created when restart = True.

4

Is my mwe clearly explained ?

5

It seems you’re always appending to the file when you’re using XDMF.write_checkpoint. Does this mean you delete /home/your_math/mwe.xdmf before you run your script each time? If you’re appending with the same time index, then this can cause issues trying to find relevant data for reading by dolfin or paraview with degenerate indices.

6

No I do not delete it. Here is what I do exactly:

  1. create the mesh “cylinder_mesh.msh” using second script (gmsh python’s API)

  2. launch the first script with restart = False. This script should create a file mwe.xdmf containing my solution at step i=0.

  3. then I relaunch the script but with restart = True. I would like it to add the solution at step i = 1 to the file mwe.xdmf. So I would mwe.xdmf to contain my solution at steps i = 0 and i = 1.

7

Could you help me with this issue ? Have you run my mwe ? I really need to do this simple thing (write solution in a already existing xdmf file), and I can’t find any example or even documentation to help me.