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:
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)
