Why the solution is not appearing

Rene_Thierry_Djoumes · May 4, 2024, 10:27am

Hello everyone,

I really need your help. I’ve managed to solve a problem with FEniCSx, but I’m having trouble displaying the solution paraview. It seems that the solution is correctly present, as you can see in the attached image.

However, to my surprise, the only result displayed is the one corresponding to the last solution written in the code (ui). Below, you’ll find a part of the code where I split the solution Ksol and save it in xdmf format.

I don’t understand why the solution is marked as (partial). I’m not sure if I’m saving the solutions incorrectly.

Thank you in advance for your assistance.

# import pyhton and Fenicsx dependecy
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import geometry,fem, io
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem import FunctionSpace, Function, Constant, Expression
from dolfinx.fem.petsc import LinearProblem, assemble, NonlinearProblem, assemble_matrix, create_vector, assemble_vector
from dolfinx.nls.petsc import NewtonSolver
import dolfinx.fem.petsc
import ufl 
from ufl import dx, grad, inner, ds, Measure, TrialFunctions, TestFunctions, SpatialCoordinate, lhs, rhs, split, lt, gt

........

#import the gmsh generated by gmsh: it should be in gmsh format

msh, cell, facet_tags = gmshio.read_from_msh("GI.msh",MPI.COMM_WORLD, gdim=3)

# time intergration
t = 0.0; dt = 0.1; Tfinal=100

freqSave = 10; inc = 0

stim_t0 = dt; stim_t2 = 5000; stim_dur = 10; stim_amp = 1

# ******* Create mesh and define function spaces **********

P1 = ufl.FiniteElement("CG", msh.ufl_cell(), 1)
P2v = ufl.VectorElement("CG", msh.ufl_cell(), 2)
FiberSpace = FunctionSpace(msh,("CG",1))

Mh = FunctionSpace(msh, P1)
Nh = FunctionSpace(msh, ufl.MixedElement([P1,P1,P1,P1]))

phih = Function(Mh)

#print(" **************** Electr Dof = ", Nh.dim())

# Define Test and Trial function
(ul,vl,ui,vi)   = TrialFunctions(Nh)
(w1,w2,m1,m2)  = TestFunctions(Nh)

Ksol = Function(Nh)

#ul,vl,ui,vi = ufl.split(Ksol)
#Ksol.name = "Electrophysiology"

Ksol0 = Function(Nh)
ul_old,vl_old,ui_old,vi_old = ufl.split(Ksol)



# ********* Electrical parameters ********* #


# ******** Save solution in Xdmf format ***

fileO = io.XDMFFile(msh.comm, "outE/ThermoElect.xdmf", "w")
fileO.write_mesh(msh)

# *********  Forcing terms and electrical stimulation ******* #


waveS0 = Function(Mh)
waveS0.interpolate(lambda x: stim_amp*(np.less(-6,x[0]))*(np.less(x[0],6))*(np.less(-6,x[1]))*(np.less(x[1],6))*(np.less(x[2],5)))



# Define the expression using FEniCSx

'''def WaveS0(amp, x):
    return amp * (lt(-6, x[1]) and lt(-6, x[0]) and lt(x[2], 5))

# Create an instance of the expression
waveS0 = WaveS0(stim_amp,x)'''




# ***** Linear weak form **** 

KL = ......

# initial condition

#Ksol.x.array[:] = 0.0

# assign previous solution

#Ksol0.x.array[:] = Ksol.x.array

bilinear_form = fem.form(KL)
LHS = assemble_matrix(bilinear_form)
LHS.assemble()

solver = PETSc.KSP().create(msh.comm)
solver.setOperators(LHS)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

while (t <= Tfinal):

    print("t=%.2f" % t)
    
    KR = .......
    
    #RHS = rhs(KR)
    #LHS = lhs(KL)
    #LHS1 = assemble_matrix(KL)
    linear_form = fem.form(KR)
    RHS = create_vector(linear_form)
    

    # Update the right hand side reusing the initial vector
    with RHS.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(RHS, linear_form)
    # solver optimisation

    # Solve linear problem
    solver.solve(RHS, Ksol.vector)
    ul,vl,ui,vi = Ksol.split()
    Ksol.x.scatter_forward()
  

    Ksol0.x.array[:] = Ksol.x.array
    ul.name = "ul"; ui.name = "ui"; vl.name = "vl"; vi.name = "vi"
    #fileO.write_function(Ksol, t)
    if (inc % freqSave == 0):
        fileO.write_function(ul, t)
        fileO.write_function(vl, t)
        fileO.write_function(vi, t)
        fileO.write_function(ui, t)
    t += dt;inc += 1
fileO.close()

francesco-ballarin · May 4, 2024, 12:12pm

This is really a paraview question, in future you may want to consider looking on their website for documentation, or their forum for questions.

In the second row of icons there is a “Time” dropdown. Change the time there to see solutions at different time steps.

Rene_Thierry_Djoumes · May 4, 2024, 1:02pm

Thank you for your reply.

I really hesitated between the discourse group of Paraview and FEniCS, but I have already obtained this result with FEniCS and the solutions are appearing in Paraview as you can see on the figure. I am in the process of rewriting my code from FEniCS to FEniCSx.

The problem isn’t ‘Time’. My problem is that all the solutions are displayed with solid color (no color appearing). I’ve tried applying some filters in Paraview and it doesn’t work.
and also when I try ‘plot selection over Time’ only one solution is present in the panel while I’ve save four solutions as you can see in the code.
As I’m new to FEniCSx, I thought maybe I was saving the solutions in the wrong way.

Thank you

francesco-ballarin · May 4, 2024, 1:08pm

In that case look at the available filters, especially those that contain the word “block”. I think one of them is called “Extract blocks”, and you should try that.

If that still doesn’t work, produce a minimal (but complete) code, because I won’t reply any further if all we can see are paraview screenshots.

Rene_Thierry_Djoumes · May 4, 2024, 3:08pm

Sorry to bother you again.
Here is the very minimal code (simple cas in 2D).
Thanks in advance for your help

# import pyhton and Fenicsx dependecy
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import geometry,fem, io, mesh
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem import FunctionSpace, Function, Constant, Expression
from dolfinx.fem.petsc import LinearProblem, assemble, NonlinearProblem, assemble_matrix, create_vector, assemble_vector
from dolfinx.nls.petsc import NewtonSolver
import dolfinx.fem.petsc
import ufl 
from ufl import dx, grad, inner, ds, Measure, TrialFunctions, TestFunctions, SpatialCoordinate, lhs, rhs, split, lt, gt

# important Functions for the electrophysiology Problem
def Fl(ul,vl):
    return kl*ul*(ul-al)*(1-ul) 

def Gl(ul,vl):
    return eps_l*(gam_l*(ul-bet_l)-vl)


#import the gmsh generated by gmsh: it should be in gmsh format
#msh, cell, facet_tags = gmshio.read_from_msh("GI.msh",MPI.COMM_WORLD, gdim=3)

# Define mesh
nx, ny = 50, 50
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
                               [nx, ny], mesh.CellType.triangle)

# time intergration
t = 0.0; dt = 0.1; Tfinal=100

freqSave = 10; inc = 0

stim_t0 = dt; stim_t2 = 5000; stim_dur = 10; stim_amp = 1

# ******* Create mesh and define function spaces **********

P1 = ufl.FiniteElement("CG", msh.ufl_cell(), 1)
P2v = ufl.VectorElement("CG", msh.ufl_cell(), 2)
FiberSpace = FunctionSpace(msh,("CG",1))

Mh = FunctionSpace(msh, P1)
Nh = FunctionSpace(msh, ufl.MixedElement([P1,P1]))

phih = Function(Mh)

#print(" **************** Electr Dof = ", Nh.dim())

# Define Test and Trial function
(ul,vl)   = TrialFunctions(Nh)
(w1,w2)  = TestFunctions(Nh)

Ksol = Function(Nh)


Ksol0 = Function(Nh)
ul_old,vl_old= ufl.split(Ksol)



'''ul_old = Function(Mh); vl_old = Function(Mh)
ui_old = Function(Mh); vi_old = Function(Mh)
ul_old.name = "ul_old"; vl_old.name = "vl_old"; ui_old.name = "ui_old"; vi_old.name = "vi_old"'''
#ul_old.interpolate(initial_condition); vl_old.interpolate(initial_condition)
#ui_old.interpolate(initial_condition); vi_old.interpolate(initial_condition)

# ********* Electrical parameters ********* #

diffScale =Constant(msh, 1.0e-2)
#***** LM Layer***
kl = Constant(msh, 10.0)
bet_l = Constant(msh, 0.0)
eps_l = Constant(msh, 0.15)
Dli = Constant(msh, 0.3)
al = Constant(msh, 0.06)
gam_l = Constant(msh, 8.0)
apl_l = Constant(msh, 1.0)
Dl = Constant(msh, 0.4)


# Save solution in Xdmf format
fileO = io.XDMFFile(msh.comm, "outE/ThermoElect.xdmf", "w")
fileO.write_mesh(msh)

# *********  Forcing terms and electrical stimulation ******* #
#waveS0=(lambda x: stim_amp*(-6<x[1]<6 and -6<x[0]<6 and x[2]<5))
waveS0 = Function(Mh)

waveS0.interpolate(lambda x: stim_amp*(np.less(-2,x[0]))*(np.less(x[0],-1.5))*(np.less(-2,x[1]))*(np.less(x[1],-1.5)))

#waveS0 = Expression("1*(-6<x[1]<6 && -6<x[0]<6 && x[2]<5 )")
# Define the expression using FEniCSx


def Istim(t):
    if (stim_t0 <= t and t <= stim_t0 + stim_dur):
        return waveS0
    #if (stim_t2 <= t and t <= stim_t2 + stim_dur):
        #return waveS0
    else:
        return Constant(msh, 0.0)
def Istim1(t):
    if (stim_t0 <= t and t <= stim_t0 + stim_dur):
        return waveS0
    #if (stim_t2 <= t and t <= stim_t2 + stim_dur):
        #return waveS0
    else:
        return Constant(msh, 0.0)

# ***** Linear weak form **** 

KL = ul/dt*w1*dx + inner(Dl*grad(ul),grad(w1))*dx\
     + vl/dt * w2 * dx\
     

# initial condition

#Ksol.x.array[:] = 0.0

# assign previous solution

#Ksol0.x.array[:] = Ksol.x.array

bilinear_form = fem.form(KL)
LHS = assemble_matrix(bilinear_form)
LHS.assemble()

solver = PETSc.KSP().create(msh.comm)
solver.setOperators(LHS)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

while (t <= Tfinal):

    print("t=%.2f" % t)
    
    KR = ul_old/dt*w1*dx + Fl(ul_old,vl_old)*w1*dx \
         + vl_old/dt * w2 * dx + Gl(ul_old,vl_old) * w2 *dx \
         + Istim(t)*w1*dx 
    
    
    #RHS = rhs(KR)
    #LHS = lhs(KL)
    #LHS1 = assemble_matrix(KL)
    linear_form = fem.form(KR)
    RHS = create_vector(linear_form)
    

    # Update the right hand side reusing the initial vector
    with RHS.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(RHS, linear_form)
    # solver optimisation

    # Solve linear problem
    solver.solve(RHS, Ksol.vector)
    ul,vl = Ksol.split()
    
    Ksol.x.scatter_forward()
    #ul.interpolate(ul)
    Ksol0.x.array[:] = Ksol.x.array
    ul.name = "ul"; vl.name = "vl"
    #fileO.write_function(Ksol, t)
    if (inc % freqSave == 0):
        fileO.write_function(ul, t)
        fileO.write_function(vl, t)
    t += dt;inc += 1
fileO.close()

francesco-ballarin · May 4, 2024, 3:44pm

Great. For the sake of future posts, you may consider simplifying the code even further when asked to post a minimal example. For instance here it was enough to interpolate a known expression into the FE space and then export the result, I didn’t need to see the actual problem you were solving.

Going back to the topic of the post: when opening the xdmf file in paraview, you are asked to choose a reader type. Try to pick Xdmf3ReaderT (instead of the one with trailing S, which is the default). With that you’ll see the correct plot for ul and vl.

Rene_Thierry_Djoumes · May 4, 2024, 4:13pm

Thank you very mush for your help. it works perfectly.
In future I’ll consider reduced again the minimal code.

Topic		Replies	Views
Reading xdmf file from FEniCSx with python is not working I/O	6	1168	March 8, 2024
Error in paraview visualization of xdmf-file with multiple solution variables dolfinx	2	1086	October 21, 2021
Converting legacy FEniCS code to FEniCSx causes issue General	19	495	May 29, 2024
Writing XDMF file leads to inf & nan values dolfinx	4	304	March 6, 2024
Read xdmf file to do pyvista visualizations I/O io , dolfinx	4	379	October 31, 2024

Why the solution is not appearing

Related topics