ruizbaier/poroelasticmodelforacutemyocarditis

Finite element simulation of immune system dynamics coupled with large-strain poromechanics

poroelasticmodelforacutemyocarditis's Issues

The results do not coincide with the paper.

I have made modifications to the code in ex04 to replicate the example from Fig. 4 in the Barnafi21 paper, but unfortunately, it did not succeed. Numerically, our results are significantly smaller even though we used the same parameters.

Here is my updated code:

from fenics import *

parameters["form_compiler"]["representation"] = "uflacs"
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2

# parameter-dependent functions

normalize_vec = lambda u : u/sqrt(dot(u,u))
kappa  = lambda J,phi: kappa0*(1 + (1-phi0)**2/phi0**3*J*phi**3*(J-phi)**2)  # isotropic Kozeny-Carman

# time constants
t = Expression("t", t=0,degree=1)
dt = 1
Tfinal = 100
cont = 0

# ********* model constants ******* #

E = Constant(3.e4)
nu = Constant(0.2)
ell = Constant(0.0)

rhos = Constant(2.e-3)
rhof = Constant(1.e-3)
kappa0 = Constant(1.e-8)
alpha = Constant(1)
muf = Constant(1.e-3)
bb = Constant((0, 0))
I = Identity(2)

# neo-Hookean
mu_s = Constant(E/(2.*(1. + nu)))
lmbdas = Constant(E*nu/((1. + nu)*(1. - 2.*nu)))

# ********* Mesh and I/O ********* #
mesh = RectangleMesh(Point(0.0, 0.0), Point(8.0, 5.0), 100, 100)
bdry = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)


class Boundary_1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)
        

class Boundary_2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 8.0)


class Boundary_3(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0) 


class Boundary_4(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 5.0)


class Boundary_5(SubDomain):
    def inside(self, x, on_boundary):
       return on_boundary and near(x[1], 5.0) and (0.0 <= x[0] <= 1.0)  


boundary_1 = Boundary_1()
boundary_2 = Boundary_2()
boundary_3 = Boundary_3()
boundary_4 = Boundary_4()
boundary_5 = Boundary_5()


boundary_1.mark(bdry, 1)
boundary_2.mark(bdry, 2)
boundary_3.mark(bdry, 3)
boundary_4.mark(bdry, 4)
boundary_5.mark(bdry, 5)


fileO = File("outputs/boundary.pvd")
fileO << bdry

ds = Measure("ds", subdomain_data= bdry)  
nn = FacetNormal(mesh)

fileO_p = File("outputs/drainage1/p.pvd")
fileO_u = File("outputs/drainage1/u.pvd")
fileO_phi = File("outputs/drainage1/phi.pvd")

# ********* Finite dimensional spaces ********* #

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Bub = FiniteElement("Bubble", mesh.ufl_cell(), 4)
P1b = VectorElement(P1 + Bub)
P2vec = VectorElement("CG", mesh.ufl_cell(), 1)
Hh = FunctionSpace(mesh, MixedElement([P2vec, P1, P1]))

print("**************** Total Dofs = ", Hh.dim())

Sol = Function(Hh)
dSol = TrialFunction(Hh)
u, phi, p = split(Sol)
v, psi, q = TestFunctions(Hh)

# ******* Boundary conditions ********** #
u_D = project(Constant((0, 0)), Hh.sub(0).collapse())
bcU = DirichletBC(Hh.sub(0), u_D, bdry, 3)
bcU1 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 1)
bcU2 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 2)
bcP1 = DirichletBC(Hh.sub(2), 0, bdry, 4)
bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 5)
# bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 1)
# bcP3 = DirichletBC(Hh.sub(2), 0, bdry, 2)
# bcP4 = DirichletBC(Hh.sub(2), 0, bdry, 3)
bcH = [bcU, bcP1,bcP2, bcU1, bcU2]

# ******* Initial conditions ********** #
phi0 = Constant(0.33)
uold = project(Constant((0, 0)), Hh.sub(0).collapse())
phiold = project(phi0, Hh.sub(1).collapse())
pold = interpolate(Constant(0), Hh.sub(2).collapse())


# ********  Weak form ********** #

ndim = u.geometric_dimension()

F = I + grad(u)
J = det(F)
invF = inv(F)

Peff = mu_s*(F - invF.T) + lmbdas*ln(J)*invF.T
t_1 = 0.2*(lmbdas+2.0*mu_s)*sin(2*pi*t/15.0)

# Poromechanics
F1 = inner(Peff - alpha*p*J*invF.T, grad(v)) * dx \
    + psi*(J-1-phi+phi0)*dx \
    + rhof*J/dt*(phi-phiold)*q * dx \
    + rhof*dot(phi*J*invF*kappa(J, phi)/muf*invF.T*grad(p), grad(q)) * dx \
    - rhos*dot(bb, v) * dx \
    - dot(-J*t_1*invF.T*nn, v) * ds(5) \
    - rhof*J*ell*q*dx
    
    

Tang = derivative(F1, Sol, dSol)

# ********* Time loop ************* #

while (t.t < Tfinal):

    t.t += dt

    print("t=%.2f" % t.t)

    # ********* Solving ************* #
    
    solve(F1 == 0, Sol, J=Tang, bcs=bcH,
          solver_parameters={'newton_solver': {'linear_solver': 'mumps',
                                               'absolute_tolerance': 1.0e-7,
                                               'relative_tolerance': 1.0e-7}})
    uh, phih, ph = Sol.split()

    # ********* Writing ************* #

    fileO_p << ph
    fileO_u << uh
    fileO_phi << phih

    cont += 1

    # ********* Updating ************* #
    assign(uold, uh)
    assign(phiold, phih)
    assign(pold, ph)

# ************* End **************** #