I follow the tutorial about the heat equation but when I open the paraview and select the xdmf file , paraview crashes after I hit apply. I choose the XdmfReaderT and the version of paraview is 5.11.0.The code is the same as the tutorial and is implemented with vscode on WSL with Ubuntu22-04 while paraview is installed on windows.I don’t know the reason.Thanks in advance!
What tutorial? Could you specify if you are using legacy dolfin or dolfinx?
Could you Also have a look at the tips in
and subsequent posts for debugging
Hello dokken, I follow the Fenicsx tutorial using dolfinx.The code corresponding with xdmf is as follows
xdmf=io.XDMFFile(domain.comm,"Heat_eq/diffusion.xdmf","w")
xdmf.write_mesh(domain)
uh=fem.Function(V)
uh.name="uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh,t)
The xdmf file is as follows where I choose a small time step with a small mesh
<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0" xmlns:xi="http://www.w3.org/2001/XInclude">
<Domain>
<Grid Name="mesh" GridType="Uniform">
<Topology TopologyType="Triangle" NumberOfElements="200" NodesPerElement="3">
<DataItem Dimensions="200 3" NumberType="Int" Format="HDF">diffusion.h5:/Mesh/mesh/topology</DataItem>
</Topology>
<Geometry GeometryType="XY">
<DataItem Dimensions="121 2" Format="HDF">diffusion.h5:/Mesh/mesh/geometry</DataItem>
</Geometry>
</Grid>
<Grid Name="uh" GridType="Collection" CollectionType="Temporal">
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="0" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="121 1" Format="HDF">diffusion.h5:/Function/uh/0</DataItem>
</Attribute>
</Grid>
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="1" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="121 1" Format="HDF">diffusion.h5:/Function/uh/1</DataItem>
</Attribute>
</Grid>
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="2" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="121 1" Format="HDF">diffusion.h5:/Function/uh/2</DataItem>
</Attribute>
</Grid>
</Grid>
</Domain>
</Xdmf>
I cannot spot anything directly wrong with the xdmf file. Could you please post a complete code that can reproduce the xdmf file.
Also, could you try using ASCII encoding and see if that works (to check if the issue is with HDF5 or not).
Please post the result of the ascii code here as well.
(call xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w", encoding=io.XDMFFile.Encoding.ASCII))
instead of the original way of opening the file.
Thanks a lot.Now it can be viewed with paraview.
The original code is as follows
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh,fem,io,plot
t=0
T=2.0
num_steps=2
dt=T/num_steps
nx,ny=5,5
domain=mesh.create_rectangle(MPI.COMM_WORLD,[np.array([-2,-2]),np.array([2,2])],
[nx,ny],mesh.CellType.triangle)
V=fem.FunctionSpace(domain,("CG",1))
def initial_condition(x,a=5):
return np.exp(-a*(x[0]**2+x[1]**2))
u_n=fem.Function(V)
u_n.name="u_n"
u_n.interpolate(initial_condition)
fdim=domain.topology.dim-1
boundary_facets=mesh.locate_entities_boundary(
domain,fdim,lambda x:np.full(x.shape[1],True,dtype=bool)
)
bc=fem.dirichletbc(PETSc.ScalarType(0),
fem.locate_dofs_topological(V,fdim,boundary_facets),
V)
xdmf=io.XDMFFile(domain.comm,"Heat_eq/diffusion.xdmf","w",encoding=io.XDMFFile.Encoding.ASCII)
xdmf.write_mesh(domain)
uh=fem.Function(V)
uh.name="uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh,t)
import ufl
u,v=ufl.TrialFunction(V),ufl.TestFunction(V)
f=fem.Constant(domain,PETSc.ScalarType(0))
a=u*v*ufl.dx+dt*ufl.dot(ufl.grad(u),ufl.grad(v))*ufl.dx
L=(u_n+dt*f)*v*ufl.dx
bilinear_form=fem.form(a)
linear_form=fem.form(L)
A=fem.petsc.assemble_matrix(bilinear_form,bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(linear_form)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
import pyvista
grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(V))
def plot_function(t, uh):
"""
Create a figure of the concentration uh warped visualized in 3D at timet step t.
"""
p = pyvista.Plotter()
# Update point values on pyvista grid
grid.point_data[f"u({t})"] = uh.x.array.real
# Warp mesh by point values
warped = grid.warp_by_scalar(f"u({t})", factor=1.5)
# Add mesh to plotter and visualize in notebook or save as figure
actor = p.add_mesh(warped)
if not pyvista.OFF_SCREEN:
p.show()
else:
pyvista.start_xvfb()
figure_as_array = p.screenshot(f"diffusion_{t:.2f}.png")
# Clear plotter for next plot
p.remove_actor(actor)
plot_function(0, uh)
for i in range(num_steps):
t+=dt
with b.localForm() as loc_b:
loc_b.set(0)
fem.petsc.assemble_vector(b,linear_form)
fem.petsc.apply_lifting(b,[bilinear_form],[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b,[bc])
solver.solve(b, uh.vector)
uh.x.scatter_forward()
# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array
# Write solution to file
xdmf.write_function(uh, t)
xdmf.close()
And now I change the code
xdmf=io.XDMFFile(domain.comm,"Heat_eq/diffusion.xdmf","w")
to
xdmf=io.XDMFFile(domain.comm,"Heat_eq/diffusion.xdmf","w",encoding=io.XDMFFile.Encoding.ASCII)
and now it can be opened in paraview.
The ASCII code of the xdmf file is as follows
<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0" xmlns:xi="http://www.w3.org/2001/XInclude">
<Domain>
<Grid Name="mesh" GridType="Uniform">
<Topology TopologyType="Triangle" NumberOfElements="50" NodesPerElement="3">
<DataItem Dimensions="50 3" NumberType="Int" Format="XML">0 1 2
0 3 2
4 0 3
3 2 5
4 6 3
3 7 5
8 4 6
6 3 7
7 5 9
8 10 6
6 11 7
7 12 9
13 8 10
10 6 11
11 7 12
12 9 14
13 15 10
10 16 11
11 17 12
12 18 14
19 13 15
15 10 16
16 11 17
17 12 18
18 14 20
19 21 15
15 22 16
16 23 17
17 24 18
18 25 20
21 15 22
22 16 23
23 17 24
24 18 25
21 26 22
22 27 23
23 28 24
24 29 25
26 22 27
27 23 28
28 24 29
26 30 27
27 31 28
28 32 29
30 27 31
31 28 32
30 33 31
31 34 32
33 31 34
33 35 34
</DataItem>
</Topology>
<Geometry GeometryType="XY">
<DataItem Dimensions="36 2" Format="XML">1.2 -2
2 -2
2 -1.2
1.2 -1.2
0.4000000000000004 -2
2 -0.3999999999999999
0.4000000000000004 -1.2
1.2 -0.3999999999999999
-0.3999999999999999 -2
2 0.4000000000000004
-0.3999999999999999 -1.2
0.4000000000000004 -0.3999999999999999
1.2 0.4000000000000004
-1.2 -2
2 1.2
-1.2 -1.2
-0.3999999999999999 -0.3999999999999999
0.4000000000000004 0.4000000000000004
1.2 1.2
-2 -2
2 2
-2 -1.2
-1.2 -0.3999999999999999
-0.3999999999999999 0.4000000000000004
0.4000000000000004 1.2
1.2 2
-2 -0.3999999999999999
-1.2 0.4000000000000004
-0.3999999999999999 1.2
0.4000000000000004 2
-2 0.4000000000000004
-1.2 1.2
-0.3999999999999999 2
-2 1.2
-1.2 2
-2 2
</DataItem>
</Geometry>
</Grid>
<Grid Name="uh" GridType="Collection" CollectionType="Temporal">
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="0" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="36 1" Format="XML">1.538828043396807e-12
4.248354255291589e-18
1.538828043396818e-12
5.573903692694586e-07
9.261360220567754e-10
9.261360220567754e-10
0.0003354626279025113
0.0003354626279025113
9.261360220567754e-10
9.261360219822879e-10
0.0003354626279025121
0.2018965179946552
0.0003354626279024659
1.538828043396818e-12
1.538828043273042e-12
5.573903692694606e-07
0.2018965179946555
0.2018965179946549
5.573903692693832e-07
4.248354255291589e-18
4.248168606834276e-18
1.538828043396818e-12
0.0003354626279025121
0.2018965179946552
0.0003354626279025107
1.538827931665216e-12
9.261360220567754e-10
0.0003354626279025113
0.0003354626279025113
9.261360220381567e-10
9.261360220567754e-10
5.573903692694586e-07
9.261360220567478e-10
1.538828043396807e-12
1.538828043396824e-12
4.248354255291589e-18
</DataItem>
</Attribute>
</Grid>
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="1" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="36 1" Format="XML">0
0
0
0.00923558672835428
0
0
0.02105375362093878
0.02105375362093877
0
-5.270820257352938e-18
0.02373766414695158
0.05276348257012871
0.02373766414695155
0
-2.714207384890636e-18
0.01222370336719974
0.05513789206935638
0.05513789206935632
0.01222370336719971
0
-4.071311077335953e-18
0
0.02373766414695158
0.05276348257012871
0.02373766414695156
-7.906230386029409e-18
0
0.02105375362093878
0.02105375362093877
-1.168718101237574e-18
0
0.009235586728354282
-5.126780515870722e-19
0
0
0
</DataItem>
</Attribute>
</Grid>
<Grid Name="uh" GridType="Uniform">
<xi:include xpointer="xpointer(/Xdmf/Domain/Grid[@GridType='Uniform'][1]/*[self::Topology or self::Geometry])" />
<Time Value="2" />
<Attribute Name="uh" AttributeType="Scalar" Center="Node">
<DataItem Dimensions="36 1" Format="XML">0
0
0
0.00599294266644615
0
0
0.01092696699634184
0.01092696699634183
0
-2.579723095934986e-18
0.01161803997357187
0.01973539960832139
0.01161803997357186
0
-1.557562897379698e-18
0.007014639684245328
0.02021004520438805
0.02021004520438805
0.007014639684245318
0
-2.336344346069547e-18
0
0.01161803997357187
0.01973539960832139
0.01161803997357187
-3.86958464390248e-18
0
0.01092696699634184
0.01092696699634184
-6.065685174328948e-19
0
0.005992942666446152
-3.326751466773499e-19
0
0
0
</DataItem>
</Attribute>
</Grid>
</Grid>
</Domain>
</Xdmf>
You would need to inspect the h5 files (for instance with h5py), as I do not get any different xdmf file than you using hdf5 as encoding.
One thing to note is that in the code you shared, there are a different number of elements in the mesh than in the xdmf file you shared in Paraview crashes when opening xdmf file - #3 by Discoverer
Appreciate for your reply.I am not familiar with h5 files and with hdf5 as encoding it still cannot be opened in paraview. Maybe something wrong with h5 file?
The first xdmf file I share is with nx=10 while the later one is with nx =5.They can all be opened in paraview using ASCII encoding.
To inspect an h5 file, you need to use a tool that can decode the binary, and you can look at the raw data. See: Quick Start Guide — h5py 3.8.0 documentation
Thank you so much.It seems that the h5 file is okay. I use the h5py to inspect the raw data (nx=5,ny=5 ,totally 50 elements)
using
f=h5py.File('diffusion.h5','r')
and f has /Function/uh/0,/Function/uh/1,/Function/uh/2 and /Mesh/mesh/geometry and /Mesh/mesh/topology.
I compare the raw data with the ASCII encoding and they are the same.So h5 encoding seems right.Can there be any other problems? 
As you have not stated what error message Paraview throws when crashing, I cannot really help any further (as I cannot reproduce the behavior with Paraview 5.11 on my linux laptop).
I really thank you for your kind reply.It just crashes after a few seconds without throwing any error. So I really dont know what happened.