System of PDEs does not work for specific geometry

1

Hello all,
I am solving the system of PDEs for my problem. To start, I have implemented the code for a simple square and everything looks good and I got desired results. However, when I am choosing my own geometry and I solve the exact problem, All of the values I get are zero. I did not apply any Dirichlet and Neumann boundary conditions and I have just source term.
this is the sample code for the square if needed I can upload my own geometry too. Could anybody tell me what I am doing wrong? I will be appreciated it.

This is the link for my geometry:


from dolfin import *
import numpy as np
from numpy import linalg as LA




Flag_fenics2017=False
if Flag_fenics2017:
    rank = MPI.rank(mpi_comm_world())
###Fenics 2018
else:
    rank = MPI.rank(MPI.comm_world)


set_log_level(30)


###########################################################
if (1):

    bmesh = UnitSquareMesh(5,5)
else:
    mesh_filename='whole_model.h5'
    bmesh = Mesh()
    if (rank == 0):
        print ('Loading mesh:', mesh_filename)
    hdf = HDF5File(bmesh.mpi_comm(), mesh_filename, "r")
    hdf.read(bmesh, "mesh",False)
WSS0=1
sci=1  #### rate of ldl influx
d_ldl=3e-3 ### LDL oxidation
d_l=2.4e-5   ### LDL diffusing out
C_ldl_in=9.44e-5 ## Influx of LDL into the leaflets in abscense of flow
tou_0=1   ### refrence WSS
k_l=1.2e-18    ### macrophage reaction with ox-LDL
d_m=1.15e-6    #### monoyte diffrentiation to macrophage
m_d=2.76e-6     #### monocyte apoptosis rate
m_r=6.6e-12     ### monocyte penetration rate at zero wss
C_m1=5.5e11     ## monocyte concentration at the lumen wall
alpha=5.76e-6   ## = M_r1/M_r2
M_r1=2.83e-11   ##### ox-LDL leading to foam cell   
M_r2=9.25e-24   ##### foam cell formatio8
k_f=1.7e-7      ### foam cell efferocytosis and defradation rate    



##### calcification phase constants 
V=27e-8   ### total aortic valve volume
lamda=1.28e-4       ##### latent TGF_B degradation rate
lamda_a=5.78e-3     #### active TGF_B activation rate by macrophage
a_F=5.2e-15         ### TGF_B activation rate by macrophages
k_prime=4.3e-6  ####    TGF to calcium
k_c=4.3e-6      ####   calcification rate
gamma=0.09          ####  calcification unit conversion
d_Ca=1              ####   calcification removal rate
strain_peak=0.2
### constant functions
def F_l_tou(wss):
    FL_tou=C_ldl_in/ ( 1 + wss/WSS0 )
    return FL_tou

fl=F_l_tou(0.5)
# sys.exit()

def F_m(wss):
    a=m_r/(1+wss/WSS0)*C_m1
    return a

def F_t(Cf):
    return (3.3e-7*Cf*V)/(Cf*V+2.84e4)


def h(CA):
    # epsil=(0.2*CA+30)/(CA+333.33)
    return 4.435e4*exp(6.404*CA)-4.435e4




time_step=24*3600*365
# tot_time=10000000


element_pde = FiniteElement('CG', bmesh.ufl_cell(),1)
# V_pde= FunctionSpace(bmesh, 'CG',1)
V_pde= FunctionSpace(bmesh, MixedElement([element_pde,element_pde,element_pde,element_pde,element_pde,element_pde,element_pde,element_pde]))
v_0p,v_1p,v_2p,v_3p,v_4p,v_5p,v_6p,v_7p= TestFunctions(V_pde)

up= Function(V_pde)

un = Function(V_pde)


up0,up1,up2,up3,up4,up5,up6,up7=split(up)

un0,un1,un2,un3,un4,un5,un6,un7=split(un)

C_F=up3+up4
k_p = Constant(time_step)
eps=Constant(0.1)

calcification_Agatston2 =  up7*( gamma )
strain = ( 0.2* calcification_Agatston2 + 30 ) /  (calcification_Agatston2 + 333.33)
wssq = 36 / (181 + calcification_Agatston2 );


F = ((up0 - un0) / k_p)*v_0p*dx + eps*dot(grad(up0), grad(v_0p))*dx - fl*v_0p*dx+d_ldl*up0*v_0p*dx+d_l*up0*v_0p*dx \
   +((up1 - un1) / k_p)*v_1p*dx + eps*dot(grad(up1), grad(v_1p))*dx - d_ldl*up0*v_1p*dx +k_l*up1*up2*v_1p*dx \
   +((up2 - un2) / k_p)*v_2p*dx + eps*dot(grad(up2), grad(v_2p))*dx - F_m(wssq)*up1*v_2p*dx + d_m*up2*v_2p*dx+m_d*up2*v_2p*dx\
   +((up3 - un3) / k_p)*v_3p*dx + eps*dot(grad(up3), grad(v_3p))*dx - d_m*up2*v_3p*dx + alpha*k_l*up1*up3*v_3p*dx \
   +((up4 - un4) / k_p)*v_4p*dx + eps*dot(grad(up4), grad(v_4p))*dx - alpha*k_l*up1*up3*v_4p*dx \
   +((up5 - un5) / k_p)*v_5p*dx + eps*dot(grad(up5), grad(v_5p))*dx - F_t(C_F)*v_5p*dx+lamda*up5*v_5p*dx+a_F*C_F*up5*v_5p*dx \
   +((up6 - un6) / k_p)*v_6p*dx + eps*dot(grad(up6), grad(v_6p))*dx - a_F*C_F*up5*v_6p*dx +lamda_a*up6*v_6p*dx \
   +((up7 - un7) / k_p)*v_7p*dx + eps*dot(grad(up7), grad(v_7p))*dx - gamma*k_prime*up6*(1+h(strain))*v_7p*dx

xdmffile_u = XDMFFile('pde_test.xdmf')


num_steps=101
t=0
count=0
for n in range(num_steps):
    _u_0,_u_1,_u_2,_u_3,_u_4,_u_5,_u_6,_u_7= up.split()

    _u_0.rename('1', 'Label')
    _u_1.rename('2', 'Label')
    _u_2.rename('3', 'Label')
    _u_3.rename('4', 'Label')
    _u_4.rename('5', 'Label')
    _u_5.rename('6', 'Label')
    _u_6.rename('7', 'Label')
    _u_7.rename('8', 'Label')

    if count%5==0:
        xdmffile_u.write(_u_0, t)
        xdmffile_u.write(_u_1, t)
        xdmffile_u.write(_u_2, t)
        xdmffile_u.write(_u_3, t)
        xdmffile_u.write(_u_4, t)
        xdmffile_u.write(_u_5, t)
        xdmffile_u.write(_u_6, t)
        xdmffile_u.write(_u_7, t)
    if rank ==0:
        print(count)

    solve(F==0,up)

    count+=1
    t += 1
    un.assign(up)
2

I haven’t inspected the PDE system in detail, but it may be that the exact solution is zero. This seems plausible when there are homogeneous Neumann boundary conditions (as implied by dropping the boundary term after integrating by parts) and zero volume sources. The initial condition in un also defaults to zero when initially defined, and doesn’t appear to be assigned any value before the time stepping loop.

3

My bad, I typed it wrong. I have a source term and I can get nonzero results for the simple square. but for the same problem in 3D geometry, I get zero values.

4

After some investigation, I discovered that this seems to be a scaling issue. Because you have a zero initial guess for up, the initial residual vector r_0 (prior to the first Newton iteration of the nonlinear solver) is:

r_0 = b - A x = b

where b is found by assembly of the rhs of your nonlinear form. On the first Newton iteration, it is equal to the integral of the constant volumetric source terms, and thus b (and r_0) are proportional to the volume of your mesh. For meshes of very small volume (such as your whole_model.h5; volume ~ 1.5e-6), the l_2 norm of r_0 is smaller than the default absolute_tolerance of the Newton solver, and so the solver exits without invoking the linear solver. Thus you obtain a 0 solution on every timestep. I illustrate this in the code below using a UnitCubeMesh and systematically decreasing the mesh volume using bmesh.scale(). As you can see, for decreasing scale aaa, the initial residual steadily decreases, until at aaa = 0.01 it becomes less than the default absolute_tolerance = 1e-10, and the Newton solver exits with

Newton solver finished in 0 iterations and 0 linear solver iterations.

To overcome this, you could try setting the absolute tolerance to 0, e.g.

solve(F == 0, up, solver_parameters={"newton_solver":{"absolute_tolerance":0}})

Since the solution (obtained for a UnitSquareMesh) appears to span 20 orders of magnitude between the largest (_u_3) and smallest (_u_5) components (see figure below), this system may be poorly scaled, but I don’t know enough about PDE scaling to offer any suggestions on how to improve it.

from dolfin import *
import numpy as np
from numpy import linalg as LA




Flag_fenics2017=False
if Flag_fenics2017:
    rank = MPI.rank(mpi_comm_world())
###Fenics 2018
else:
    rank = MPI.rank(MPI.comm_world)


set_log_level(20)
set_log_active(True)


###########################################################
for aaa in (1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01):
    bmesh = UnitCubeMesh(1,1,1)
    bmesh.scale(aaa)
    if rank == 0:
        print(f"aaa = {aaa}")

    WSS0=1
    sci=1  #### rate of ldl influx
    d_ldl=3e-3 ### LDL oxidation
    d_l=2.4e-5   ### LDL diffusing out
    C_ldl_in=9.44e-5 ## Influx of LDL into the leaflets in abscense of flow
    tou_0=1   ### refrence WSS
    k_l=1.2e-18    ### macrophage reaction with ox-LDL
    d_m=1.15e-6    #### monoyte diffrentiation to macrophage
    m_d=2.76e-6     #### monocyte apoptosis rate
    m_r=6.6e-12     ### monocyte penetration rate at zero wss
    C_m1=5.5e11     ## monocyte concentration at the lumen wall
    alpha=5.76e-6   ## = M_r1/M_r2
    M_r1=2.83e-11   ##### ox-LDL leading to foam cell   
    M_r2=9.25e-24   ##### foam cell formatio8
    k_f=1.7e-7      ### foam cell efferocytosis and defradation rate    



    ##### calcification phase constants 
    V=27e-8   ### total aortic valve volume
    lamda=1.28e-4       ##### latent TGF_B degradation rate
    lamda_a=5.78e-3     #### active TGF_B activation rate by macrophage
    a_F=5.2e-15         ### TGF_B activation rate by macrophages
    k_prime=4.3e-6  ####    TGF to calcium
    k_c=4.3e-6      ####   calcification rate
    gamma=0.09          ####  calcification unit conversion
    d_Ca=1              ####   calcification removal rate
    strain_peak=0.2
    ### constant functions
    def F_l_tou(wss):
        FL_tou=C_ldl_in/ ( 1 + wss/WSS0 )
        return FL_tou

    fl=F_l_tou(0.5)
    # sys.exit()

    def F_m(wss):
        a=m_r/(1+wss/WSS0)*C_m1
        return a

    def F_t(Cf):
        return (3.3e-7*Cf*V)/(Cf*V+2.84e4)


    def h(CA):
        # epsil=(0.2*CA+30)/(CA+333.33)
        return 4.435e4*exp(6.404*CA)-4.435e4




    time_step=24*3600*365
    # tot_time=10000000


    element_pde = FiniteElement('CG', bmesh.ufl_cell(),1)
    # V_pde= FunctionSpace(bmesh, 'CG',1)
    V_pde= FunctionSpace(bmesh, MixedElement([element_pde,element_pde,element_pde,element_pde,element_pde,element_pde,element_pde,element_pde]))
    v_0p,v_1p,v_2p,v_3p,v_4p,v_5p,v_6p,v_7p= TestFunctions(V_pde)

    up= Function(V_pde)

    un = Function(V_pde)


    up0,up1,up2,up3,up4,up5,up6,up7=split(up)

    un0,un1,un2,un3,un4,un5,un6,un7=split(un)

    C_F=up3+up4
    k_p = Constant(time_step)
    eps=Constant(0.1)

    calcification_Agatston2 =  up7*( gamma )
    strain = ( 0.2* calcification_Agatston2 + 30 ) /  (calcification_Agatston2 + 333.33)
    wssq = 36 / (181 + calcification_Agatston2 );


    F = ((up0 - un0) / k_p)*v_0p*dx + eps*dot(grad(up0), grad(v_0p))*dx - fl*v_0p*dx+d_ldl*up0*v_0p*dx+d_l*up0*v_0p*dx \
       +((up1 - un1) / k_p)*v_1p*dx + eps*dot(grad(up1), grad(v_1p))*dx - d_ldl*up0*v_1p*dx +k_l*up1*up2*v_1p*dx \
       +((up2 - un2) / k_p)*v_2p*dx + eps*dot(grad(up2), grad(v_2p))*dx - F_m(wssq)*up1*v_2p*dx + d_m*up2*v_2p*dx+m_d*up2*v_2p*dx\
       +((up3 - un3) / k_p)*v_3p*dx + eps*dot(grad(up3), grad(v_3p))*dx - d_m*up2*v_3p*dx + alpha*k_l*up1*up3*v_3p*dx \
       +((up4 - un4) / k_p)*v_4p*dx + eps*dot(grad(up4), grad(v_4p))*dx - alpha*k_l*up1*up3*v_4p*dx \
       +((up5 - un5) / k_p)*v_5p*dx + eps*dot(grad(up5), grad(v_5p))*dx - F_t(C_F)*v_5p*dx+lamda*up5*v_5p*dx+a_F*C_F*up5*v_5p*dx \
       +((up6 - un6) / k_p)*v_6p*dx + eps*dot(grad(up6), grad(v_6p))*dx - a_F*C_F*up5*v_6p*dx +lamda_a*up6*v_6p*dx \
       +((up7 - un7) / k_p)*v_7p*dx + eps*dot(grad(up7), grad(v_7p))*dx - gamma*k_prime*up6*(1+h(strain))*v_7p*dx

    xdmffile_u = XDMFFile('pde_test.xdmf')


    num_steps=1#101
    t=0
    count=0
    for n in range(num_steps):
        _u_0,_u_1,_u_2,_u_3,_u_4,_u_5,_u_6,_u_7= up.split()

        _u_0.rename('1', 'Label')
        _u_1.rename('2', 'Label')
        _u_2.rename('3', 'Label')
        _u_3.rename('4', 'Label')
        _u_4.rename('5', 'Label')
        _u_5.rename('6', 'Label')
        _u_6.rename('7', 'Label')
        _u_7.rename('8', 'Label')

        if count%5==0:
            xdmffile_u.write(_u_0, t)
            xdmffile_u.write(_u_1, t)
            xdmffile_u.write(_u_2, t)
            xdmffile_u.write(_u_3, t)
            xdmffile_u.write(_u_4, t)
            xdmffile_u.write(_u_5, t)
            xdmffile_u.write(_u_6, t)
            xdmffile_u.write(_u_7, t)
        if rank ==0:
            print(count)

        solve(F==0,up)

        count+=1
        t += 1
        un.assign(up)

Results_Unit_Square
Results_Unit_Square1075×763 74.3 KB

2 Likes
5

Wow, you saved me. I was trying everything and I could not find the issue. Thank you very much for your complete explanation!!
Bests,