Reading xdmf file from FEniCSx with python is not working

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!