Reading xdmf file from FEniCSx with python is not working

[psj1866]

Hello,
I am using Paraview (version 5.10.1) to visualize the xdmf file created from the FEniCSx code.
Actually, I can visualize my code well if I just load the file by using GUI :

However, I am trying to use Python shell in the Paraview to load the same file but it doesn’t work : I tried XdmfReaderT, ExodusIIReader, Xdmf3ReaderS but none of them succeeded.
Below is one of the codes that I tried

reader = Xdmf3ReaderS(FileName = "~/Data_Terzaghi/p_10.xdmf")

where after execution, I only load the empty file with no data.
And here I attach the codes for file generation :

# Solve system & output results
# ----------------------------------------------------------------
# time stepping
time = np.linspace(t_i, t_f, Nsteps+1)

# output file
#xdmf = XDMFFile(domain.comm, "terzaghi.xdmf", "w")
#xdmf.write_mesh(domain)
with io.XDMFFile(MPI.COMM_WORLD, "u_"+str(n)+".xdmf", "w") as xdmf1:
    xdmf1.write_mesh(domain)
with io.XDMFFile(MPI.COMM_WORLD, "p_"+str(n)+".xdmf", "w") as xdmf2:
    xdmf2.write_mesh(domain)

# Writing initial condition
t = t_i
u, p = xt.split()
u.name = "u"
p.name = "p"

xdmf1.write_function(u, t_i)
xdmf2.write_function(p, t_i)
   
#Solution
for (i, dt) in enumerate(np.diff(time)):

    t = time[i+1]
    print('-----------------------------------------')
    print('>> t =', t, '[sec]')
    
    #Variational formulation
    Res = inner(grad(eta), sigma_e(ut)) * dx - inner(pt, div(eta)) * dx - dot(eta, trac) * ds(4)\
        + inner(psi, div((ut - u_n)/dt)) * dx - inner(grad(psi), w_Darcy(pt)) * dx

    Jac = derivative(Res, xt) # Jacobian
    Prob   = NonlinearProblem(Res, xt, BC, Jac)
    solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, Prob)

    # Set nonlinear solver parameters
    solver.rtol = 1e-7 #relative_tolerance
    solver.atol = 1e-7 #absolute_tolerance 
    solver.max_it = 100 #maximum_iterations
    solver.report = True
    
    # solve system
    solver.solve(xt)
    #xt.x.scatter_forward()   

    # update storage variable
    x_n.vector.array[:] = xt.vector.array[:]
    
    # data for visualization
    u, p = xt.split()
    u.name = "u"
    p.name = "p"

    xdmf1.write_function(u, t)
    xdmf2.write_function(p, t)
    
xdmf1.close()
xdmf2.close()

Thanks for reading!

[Mamadou]

Hi dear Fenicsx users,

I am using fenicsx to approximate the solution of my problem. Following these tutorials, I am able now to define the mesh, the Dirichlet boundary and initials conditions, the variational formulation of my problem, and solve the problem (i.e., we will use the following naming convention: m (mathematically m^{n+1} is known as a trial function in the variational form. m_ is the most recently computed approximation (m^{n+1} available as a Function object), m_n is m^{n} , and the same convention goes for T(T^{n+1}), T_n, T^n.

Here I attach the code:

import os

import matplotlib.pyplot as plt

import numpy as np

import ufl

from dolfinx.io import VTXWriter





import math
from math import sqrt
# For plotting
import matplotlib.pyplot as plt






from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)

from dolfinx import log, plot
#from dolfinx.plot import vtk_mesh
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, TrialFunction, derivative, dot
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType

import math
from math import sqrt
from mpi4py import MPI
from petsc4py import PETSc

try:
    import pyvista as pv
    import pyvistaqt as pvqt
    have_pyvista = True
    if pv.OFF_SCREEN:
        pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
    print("pyvista and pyvistaqt are required to visualise the solution")
    have_pyvista = False

# Save all logging to file
log.set_output_file("log.txt")

import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

#lx = 1.0  #length of x-axis


# parmeters
r_0 = 0
r_b = 1  # the coordinate of right boundary


# Define mesh            #nums =  [8, 16, 32,64,128,255] number of cells
nx, ny = 15, 15      
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([r_0, r_0]), np.array([r_b, r_b])],
                               [nx, ny], mesh.CellType.quadrilateral)





#
###################################################
#### Parameters : time and load and material ######
###################################################
# 

am = fem.Constant(msh, ScalarType(3.e-06 ))  
k = fem.Constant(msh, ScalarType(0.12))  
c = fem.Constant(msh, ScalarType(1284.))    
rho_0 = fem.Constant(msh, ScalarType(500.))      
R = fem.Constant(msh, ScalarType(2400.))    
alpha = k/(c*rho_0)
beta = R/(c*rho_0)
epsilon = fem.Constant(msh, ScalarType(0.1))  
nu = fem.Constant(msh, ScalarType(1.e-02))        




## Time parametrization
t         = 0                # Start time
Tf        = 10.            # End time
num_steps = 100    # Number of time steps
dt = (Tf-t)/num_steps # Time step size

 

#
###################################################
########## Mixed Space of resolution ##############
###################################################
# 

# Define ME function space
# Lagrange" is type of the finite element, and 1 is the degree
m1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
T1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)

ME = FunctionSpace(msh, m1 * T1)

#Q = FunctionSpace(msh, m1) #ME.sub(0)
#V = FunctionSpace(msh, T1) #ME.sub(1)


# Create the initial and solution functions of space
u = Function(ME)  # current solution
u_n = Function(ME)  # solution from previous converged step (i.e., initial solution)
#du = TrialFunction(ME)


#We create four python functions for determining the facets to apply boundary conditions to




def right(x):
    return np.isclose(x[0], r_b)  #right (x = 0.1)


def top(x):
    return np.isclose(x[1], r_b)  #top (y = 0.1)




# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_right, ME.sub(0))


boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_top, ME.sub(0))




boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_right)
bc_T_right = dirichletbc(PETSc.ScalarType(60.), boundary_dofs_right, ME.sub(1))


boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_top)
bc_T_top = dirichletbc(PETSc.ScalarType(60.), boundary_dofs_top, ME.sub(1))



bcs = [bc_m_right, bc_m_top, bc_T_right, bc_T_top]  #dirichlet boundary conditions




# 
# 
###################################################
############## Initial Values #####################
###################################################
# 

Mi = 0.30 #Initial moisture content
# 
M_n_, M_n_to_ME = ME.sub(0).collapse()                             
FM_n_ = Function(M_n_)
with FM_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Mi)) 
     
# # Re assign in u_n
u_n.x.array[M_n_to_ME] = FM_n_.x.array
u_n.x.scatter_forward()

Ti = 10. #Initial temperature content
T_n_, T_n_to_ME = ME.sub(1).collapse()
FT_n_ = Function(T_n_)                                         
with FT_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Ti)) 
     
# # Re assign in u_n
u_n.x.array[T_n_to_ME] = FT_n_.x.array
u_n.x.scatter_forward()
#  


# Identify the unknowns from the function

# Split mixed functions
m, T = ufl.split(u)
m_n, T_n = ufl.split(u_n)



# Set up the test functions (q, v)
q, v = ufl.TestFunctions(ME)


#  
###################################################
############## Variationnal form ##################
###################################################
# 


#Variational formulation of the system
F0 = inner(m, q) * dx - inner(m_n, q) * dx + dt*am * inner(grad(m), grad(q)) * dx + dt*am*nu * inner(grad(T), grad(q)) * dx
#F1 = c*rho_0* inner(T, v) * dx - c*rho_0* inner(T_n, v) * dx + dt*k * inner(grad(T), grad(v)) * dx - R*epsilon * inner(m, v) * dx + R*epsilon * inner(m_n, v) * dx
F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt*alpha * inner(grad(T), grad(v)) * dx - beta*epsilon * inner(m, v) * dx + beta*epsilon * inner(m_n, v) * dx
F = F0 + F1



# Non linear problem definition
du  = TrialFunction(ME)
#Compute the jacobian
Jac = derivative(F, u, du)
problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)

#Solve variational problem
# Create nonlinear problem and Newton solver
#problem = NonlinearProblem(F, u, bcs = bcs)     #Replace
#problem = NonlinearProblem(F, u, bcs = bcs)  
  
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"      #residual    #incremental
# Absolute tolerance
solver.atol = 5e-10
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12
#solver.max_it = 70 
solver.report = True
#solver.max_it = 100        #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 2


info = solver.solve(u)
print(info)

I am trying to use Paraview to visualize my results but it doesn’t work (i.e., empty file with no data).
Here is the code to create the loop over time.

# 
###################################################
################ Pre-processing ###################
###################################################
# 
# Create an output xdmf file to store the values
# xdmf = XDMFFile(mesh.comm, "Column.xdmf", "w")
# xdmf.write_mesh(mesh)

# We prepare output files for the velocity and pressure data, and write the mesh and initial conditions to file
# Create an output xdmf file to store the values
file = XDMFFile(msh.comm, "problem_rectangle_mesh_4th_march2024_alhaatim.xdmf", "w")
file.write_mesh(msh)

#with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", [u_n.sub(0).collapse()], "w") as xdmf1:
with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", "w") as xdmf1:
    xdmf1.write_mesh(msh)
with io.XDMFFile(MPI.COMM_WORLD, "temperature_4th_March2024.xdmf", "w") as xdmf2:
    xdmf2.write_mesh(msh)



# 
###################################################
################ Processing #######################
###################################################
# time steps to evaluate the pressure in space:
#n0, n1, n2 = 200,400,800
# 
t = 0
#L2_p = np.zeros(num_steps, dtype=PETSc.ScalarType)
for n in range(num_steps):
	t += dt
	try:
		num_its, converged = solver.solve(u)
	except:
		if MPI.COMM_WORLD.rank == 0:
			print("*************") 
			print("Solver failed")
			print("*************") 
			pass
	u.x.scatter_forward()
      
	# Update Value
	u_n.x.array[:] = u.x.array
	u_n.x.scatter_forward()
	__m, __T = u.split()
	# 
	# Export the results
	__m.name = "Moisture"
	__T.name = "Temperature"
      
	xdmf1.write_function(__m,t)
	xdmf2.write_function(__T,t)
	# 
    

# 
file.close()
#

xdmf1.close()
xdmf2.close() 

Thank you for your time and your patience.

[dokken]

Your code is way too complex for most people for have time to parse it. Please try to reduce the problem to something smaller, as there is so much information, and multiple xdmf files in this code.
Please also check what is written in each of the xdmf files (they are text files and can be opened in a text editor), and post the relevant (non-working files) here.

[Mamadou]

Thank you Dr. Dokken for your reply.

I am trying to approximate the following system of PDEs (with m and T are the unknowns functions).

\dfrac{\partial {m}}{\partial t} \ = \ a_m \Delta m + a_m \nu \Delta T \ , \hspace{1cm} in \ \Omega \times [0, tf] \hspace{4cm} \textbf{(1)}

c \rho_{0} \dfrac{\partial {T}}{\partial t} \ = k \Delta T + R \epsilon \dfrac{\partial {m}}{\partial t} \ , \hspace{1cm} in \ \Omega \times [0, tf] \hspace{4.2cm} \textbf{(2)}

with initial conditions m(x, y, 0) = 0.30 and T(x, y, 0) = 10.
The boundary boundary conditions are:
m(x=1, y, t) = 0.12, m(x, y=1, t), T(x=1, y, t) = 60 and T(x, y=1, t) = 60.

To do so, I am following this tutorial Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial, using a mixed space element as in Solving Stokes equation with NewtonSolver - #4 by dokken and https://fenicsproject.org/olddocs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html.

Here is the implementation of the mesh, mixed element, boundary conditions and initial conditions, variational formulation. Using a NonProblem solver, it converges in two iterations.

import os

import matplotlib.pyplot as plt

import numpy as np

import ufl

from dolfinx.io import VTXWriter





import math
from math import sqrt
# For plotting
import matplotlib.pyplot as plt






from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)

from dolfinx import log, plot
#from dolfinx.plot import vtk_mesh
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, TrialFunction, derivative, dot
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType

import math
from math import sqrt
from mpi4py import MPI
from petsc4py import PETSc

try:
    import pyvista as pv
    import pyvistaqt as pvqt
    have_pyvista = True
    if pv.OFF_SCREEN:
        pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
    print("pyvista and pyvistaqt are required to visualise the solution")
    have_pyvista = False

# Save all logging to file
log.set_output_file("log.txt")

import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

#lx = 1.0  #length of x-axis


# parmeters
r_0 = 0
r_b = 1  # the coordinate of right boundary


# Define mesh            #nums =  [8, 16, 32,64,128,255] number of cells
nx, ny = 15, 15      
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([r_0, r_0]), np.array([r_b, r_b])],
                               [nx, ny], mesh.CellType.quadrilateral)





#
###################################################
#### Parameters : time and load and material ######
###################################################
# 

am = fem.Constant(msh, ScalarType(3.e-06 ))  
k = fem.Constant(msh, ScalarType(0.12))  
c = fem.Constant(msh, ScalarType(1284.))    
rho_0 = fem.Constant(msh, ScalarType(500.))      
R = fem.Constant(msh, ScalarType(2400.))    
alpha = k/(c*rho_0)
beta = R/(c*rho_0)
epsilon = fem.Constant(msh, ScalarType(0.1))  
nu = fem.Constant(msh, ScalarType(1.e-02))        




## Time parametrization
t         = 0                # Start time
Tf        = 10.            # End time
num_steps = 100    # Number of time steps
dt = (Tf-t)/num_steps # Time step size

 

#
###################################################
########## Mixed Space of resolution ##############
###################################################
# 

# Define ME function space
# Lagrange" is type of the finite element, and 1 is the degree
m1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
T1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)

ME = FunctionSpace(msh, m1 * T1)

#Q = FunctionSpace(msh, m1) #ME.sub(0)
#V = FunctionSpace(msh, T1) #ME.sub(1)


# Create the initial and solution functions of space
u = Function(ME)  # current solution
u_n = Function(ME)  # solution from previous converged step (i.e., initial solution)
#du = TrialFunction(ME)


#We create four python functions for determining the facets to apply boundary conditions to




def right(x):
    return np.isclose(x[0], r_b)  #right (x = 0.1)


def top(x):
    return np.isclose(x[1], r_b)  #top (y = 0.1)




# 
###################################################
######## Dirichlet boundary conditions #############
###################################################
# 

boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_right, ME.sub(0))


boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_top, ME.sub(0))




boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_right)
bc_T_right = dirichletbc(PETSc.ScalarType(60.), boundary_dofs_right, ME.sub(1))


boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_top)
bc_T_top = dirichletbc(PETSc.ScalarType(60.), boundary_dofs_top, ME.sub(1))



bcs = [bc_m_right, bc_m_top, bc_T_right, bc_T_top]  #dirichlet boundary conditions




# 
# 
###################################################
############## Initial Values #####################
###################################################
# 

Mi = 0.30 #Initial moisture content
# 
M_n_, M_n_to_ME = ME.sub(0).collapse()                             
FM_n_ = Function(M_n_)
with FM_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Mi)) 
     
# # Re assign in u_n
u_n.x.array[M_n_to_ME] = FM_n_.x.array
u_n.x.scatter_forward()

Ti = 10. #Initial temperature content
T_n_, T_n_to_ME = ME.sub(1).collapse()
FT_n_ = Function(T_n_)                                         
with FT_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Ti)) 
     
# # Re assign in u_n
u_n.x.array[T_n_to_ME] = FT_n_.x.array
u_n.x.scatter_forward()
#  


# Identify the unknowns from the function

# Split mixed functions
m, T = ufl.split(u)
m_n, T_n = ufl.split(u_n)



# Set up the test functions (q, v)
q, v = ufl.TestFunctions(ME)


#  
###################################################
############## Variationnal form ##################
###################################################
# 


#Variational formulation of the system
F0 = inner(m, q) * dx - inner(m_n, q) * dx + dt*am * inner(grad(m), grad(q)) * dx + dt*am*nu * inner(grad(T), grad(q)) * dx
#F1 = c*rho_0* inner(T, v) * dx - c*rho_0* inner(T_n, v) * dx + dt*k * inner(grad(T), grad(v)) * dx - R*epsilon * inner(m, v) * dx + R*epsilon * inner(m_n, v) * dx
F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt*alpha * inner(grad(T), grad(v)) * dx - beta*epsilon * inner(m, v) * dx + beta*epsilon * inner(m_n, v) * dx
F = F0 + F1



# Non linear problem definition
du  = TrialFunction(ME)
#Compute the jacobian
Jac = derivative(F, u, du)
problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)

#Solve variational problem
# Create nonlinear problem and Newton solver
#problem = NonlinearProblem(F, u, bcs = bcs)     #Replace
#problem = NonlinearProblem(F, u, bcs = bcs)  
  
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"      #residual    #incremental
# Absolute tolerance
solver.atol = 5e-10
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12
#solver.max_it = 70 
solver.report = True
#solver.max_it = 100        #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 2


info = solver.solve(u)
print(info)

The next step is to create a loop over time. I built three xdmf files with xdmf1 and xdmf2 containing the solutions m and T.

# 
###################################################
################ Pre-processing ###################
###################################################
# 
# Create an output xdmf file to store the values
# xdmf = XDMFFile(mesh.comm, "Column.xdmf", "w")
# xdmf.write_mesh(mesh)

# We prepare output files for the velocity and pressure data, and write the mesh and initial conditions to file
# Create an output xdmf file to store the values
file = XDMFFile(msh.comm, "problem_rectangle_mesh_4th_march2024_alhaatim.xdmf", "w")
file.write_mesh(msh)

#with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", [u_n.sub(0).collapse()], "w") as xdmf1:
with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", "w") as xdmf1:
    xdmf1.write_mesh(msh)
with io.XDMFFile(MPI.COMM_WORLD, "temperature_4th_March2024.xdmf", "w") as xdmf2:
    xdmf2.write_mesh(msh)



# 
###################################################
################ Processing #######################
###################################################
# time steps to evaluate the pressure in space:
#n0, n1, n2 = 200,400,800
# 
t = 0
#L2_p = np.zeros(num_steps, dtype=PETSc.ScalarType)
for n in range(num_steps):
	t += dt
	try:
		num_its, converged = solver.solve(u)
	except:
		if MPI.COMM_WORLD.rank == 0:
			print("*************") 
			print("Solver failed")
			print("*************") 
			pass
	u.x.scatter_forward()
      
	# Update Value
	u_n.x.array[:] = u.x.array
	u_n.x.scatter_forward()
	__m, __T = u.split()
	# 
	# Export the results
	__m.name = "Moisture"
	__T.name = "Temperature"
      
	xdmf1.write_function(__m,t)
	xdmf2.write_function(__T,t)
	# 
    

# 
file.close()
#

xdmf1.close()
xdmf2.close()

The above doesn’t return an error, but when I try to use Paraview to visualise my results, it doesn’t work (i.e. empty file with no data).

Thank you for your time and your patience.

[Mamadou]

@Mamadou

###################################################
############## Initial Values #####################
###################################################
# 

Mi = 0.30 #Initial moisture content
# 
M_n_, M_n_to_ME = ME.sub(0).collapse()                             
FM_n_ = Function(M_n_)
with FM_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Mi)) 
     
# # Re assign in u_n
u_n.x.array[M_n_to_ME] = FM_n_.x.array
u_n.x.scatter_forward()

Ti = 10. #Initial temperature content
T_n_, T_n_to_ME = ME.sub(1).collapse()
FT_n_ = Function(T_n_)                                         
with FT_n_.vector.localForm() as initial_local:
	initial_local.set(ScalarType(Ti)) 
     
# # Re assign in u_n
u_n.x.array[T_n_to_ME] = FT_n_.x.array
u_n.x.scatter_forward()
#  

Please let me know if my initial conditions are well implemented (i.e., m(x, y, 0) = 0.30 and T(x, y, 0) = 10.)

[dokken]

@Mamadou

file = XDMFFile(msh.comm, "problem_rectangle_mesh_4th_march2024_alhaatim.xdmf", "w")
file.write_mesh(msh)

#with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", [u_n.sub(0).collapse()], "w") as xdmf1:
with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", "w") as xdmf1:
    xdmf1.write_mesh(msh)
with io.XDMFFile(MPI.COMM_WORLD, "temperature_4th_March2024.xdmf", "w") as xdmf2:
    xdmf2.write_mesh(msh)

You are closing these files before trying to use them

@Mamadou

xdmf1.write_function(__m,t)
	xdmf2.write_function(__T,t)

please initiate xdmf1 and xdmf2 in the same way as the variable file:

@Mamadou

file = XDMFFile(msh.comm, "problem_rectangle_mesh_4th_march2024_alhaatim.xdmf", "w")

Note that your example is still very far away from a minimal example and the likelihood of anyone having time to parse all of it is minimal. Please make an effort in minimizing the example to handle the one thing you are struggling with, not a full application code.

[Mamadou]

Hello everyone,

Thank you Prof Dokken for your time and patience.
I have found a lot of discussions on Saved XDMF file immediately crashes Paraview when opened and [Reading xdmf file from FEniCSx with python is not working] (Reading xdmf file from FEniCSx with python is not working) to name a few.

Following this discussion https://fenicsproject.discourse.group/t/paraview-crash-when-exporting-xdmf-or-vtx-files/13270/5 I am able to visualise my results.

The problem was the choice of XDMF reader. I tested XDMF3 READER instead of XDMF3 READER (Saved XDMF file immediately crashes Paraview when opened).

Thank you again for your help.