Issue with scaling up the simulation while running in parallel

1

Hello all,

The set of equations I’m simulating describe the time evolution of the 2 dimensional hydrogel gel subjected to the osmotic stress due to presence of the chemical ions inside. On the top of the gel, there exists a supernatant domain and chemical ions can move in or out of the gel due to diffusion as well as advection. The fluid flow in the supernatant domain (governed by Stokes equation) is coupled to the fluid flow in the gel domain(governed by Darcy’s flow) by flux continuity at the gel-supernatant domain boundary. The supernatant domain on top of the gel is flipped at the gel-supernatant boundary, such that the whole system can be simulated with one rectangular mesh describing both gel domain and supernatant domain.

I run the code in parallel on HPC cluster(with mpirun -np 2 python3 hydrogel_swelling.py). The code works fine for the mesh size 10x10 and 20x20 but for higher mesh sizes it hangs. I could still see the processes running but I think it fails to solve the equations even in the first iteration. Initially I thought it could be memory issue but I never got out of memory error and there seem to be plenty of memory remaining on the cluster node (only 5 GB utilized out of 256 GB) when checked with htop command.

I would really appreciate any insight or help. I have attached the code below along with the output before it hangs

code -

#importing required libraries

from fenics import *


from ufl import tanh
from ufl import nabla_div
from ufl import exp
from time import perf_counter
from datetime import date

t1_start = perf_counter()

Tf = 1   #final time
dt = Constant(1e-7) #initial timestep

#constants
lambda_ = Constant(1.)
mu = Constant(1.)
epsilon_p = Constant(0.1)
chi = Constant(50.)
gamma = Constant(100.)
xi = Constant(1.)
zeta = Constant(100.)
Cudiff = Constant(10.)
Hdiff = Constant(10000.)
r = Constant(5.*pow(10., 5.))
phi0 = Constant(0.006)
phiac = Constant(0.003)
phiac_hcl = Constant(0.006)
phiac_hcl2 = Constant(0.006)
delta = Constant(2.)
alpha = Constant(1000.)
beta = Constant(1e-9)
ccn = Constant(-1.)
ccn_1 = Constant(0.)
ccn_2 = Constant(0.)
ccn_3 = Constant(0.)
ccn_4 = Constant(0.)
ccr = Constant(1.)


#defining the mesh
mesh = RectangleMesh(Point(0,0),Point(10,1),50,50)
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
P2 = VectorElement('P', mesh.ufl_cell(), 2) 
element = MixedElement([P2,P1,P1,P1,P1,P1,P1,P1,P2,P1])
V = FunctionSpace(mesh, element)
V1 = FunctionSpace(mesh,'P',1)
V2 = VectorFunctionSpace(mesh,'P',2)

v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 = TestFunctions(V)

u = Function(V)
un = Function(V)
un_1 = Function(V)
un_2 = Function(V)
un_3 = Function(V)
un_4 = Function(V)
du = TrialFunction(V)
dv = TrialFunction(V)
d = u.geometric_dimension()

#u1 = deformation
#u2 = pore pressure
#u3 = Free copper volume fraction in the gel
#u4 = Free acid volume fraction in the gel
#u5 = Bound copper volume fraction in the gel
#u6 = Bound acid volume fraction in the gel
#u7 = Free copper volume fraction in the supernatant domain
#u8 = Free acid volume fraction in the supernatant domain
#u9 = Flow velocity in the supernatant domain
#u10 = Pressure in the supernatant domain


#defining initial conditions


class InitialCondition(UserExpression):
    def eval(self, values, x):
        values[0] = 0.
        values[1] = -x[1]
        values[2] = 0.
        values[3] = 0.
        values[4] = 0.
        values[5] = phiac
        values[6] = 0.
        values[7] = 0.
        values[8] = (1/(2*1.7724))*exp(-(x[0]-5)*(x[0]-5)/4)*phiac_hcl2*0.5*(1-tanh(20*(x[1]-0.5)))   
        values[9] = 0.
        values[10] = 0.
        values[11] = 0.
    def value_shape(self):
        return(12,)
        
u0 = InitialCondition()
un = interpolate(u0,V)

u1, u2, u3, u4, u5, u6, u7, u8, u9, u10 = split(u); #splitting u
un1, un2, un3, un4, un5, un6, un7, un8, un9, un10 = split(un);
un_11, un_12, un_13, un_14, un_15, un_16, un_17, un_18, un_19, un_110 = split(un_1);
un_21, un_22, un_23, un_24, un_25, un_26, un_27, un_28, un_29, un_210 = split(un_2);
un_31, un_32, un_33, un_34, un_35, un_36, un_37, un_38, un_39, un_310 = split(un_3);
un_41, un_42, un_43, un_44, un_45, un_46, un_47, un_48, un_49, un_410 = split(un_4);



#defining bdf functions
def BDFparameters(nn, ccn, ccn_1, ccn_2, ccn_3, ccn_4, ccr):
    if nn==0:
        ccn.assign(Constant(-1.)); ccr.assign(Constant(1.)); ccn_1.assign(Constant(0.))
        ccn_2.assign(Constant(0.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))                
    elif nn==1:
        ccn.assign(Constant(-4./3.)); ccr.assign(Constant(2./3.)); ccn_1.assign(Constant(1./3.))
        ccn_2.assign(Constant(0.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))
    elif nn==2:
        ccn.assign(Constant(-18./11.)); ccr.assign(Constant(6./11.)); ccn_1.assign(Constant(9./11.))
        ccn_2.assign(Constant(-2./11.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))
    elif nn==3:
        ccn.assign(Constant(-48./25.)); ccr.assign(Constant(12./25.)); ccn_1.assign(Constant(36./25.))
        ccn_2.assign(Constant(-16./25.)); ccn_3.assign(Constant(3./25.)); ccn_4.assign(Constant(0.))
    else:
        ccn.assign(Constant(-300./137.)); ccr.assign(Constant(60./137.)); ccn_1.assign(Constant(300./137.))
        ccn_2.assign(Constant(-200./137.)); ccn_3.assign(Constant(75./137.)); ccn_4.assign(Constant(-12./137.))

def BDFvariables(nn, u, un, un_1, un_2, un_3, un_4):
    if nn==0:
        un.assign(u)    
    elif nn==1:
        un_1.assign(un); un.assign(u)  
    elif nn==2:
        un_2.assign(un_1); un_1.assign(un); un.assign(u)  
    elif nn==3:
        un_3.assign(un_2); un_2.assign(un_1); un_1.assign(un); un.assign(u)  
    else:
        un_4.assign(un_3); un_3.assign(un_2); un_2.assign(un_1); un_1.assign(un); un.assign(u)
        

def an1(uud):
    return 1./2. *(1+tanh(pow(10,3)*uud))

        
#defining boundaries

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)

class BoundaryL(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)
bxL = BoundaryL()
bxL.mark(boundary_markers, 0) 
        
class BoundaryR(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 10, DOLFIN_EPS)
bxR = BoundaryR()
bxR.mark(boundary_markers, 2)
        
class BoundaryD(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0, DOLFIN_EPS)
bxD = BoundaryD()
bxD.mark(boundary_markers, 1)
        
class BoundaryU(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1, DOLFIN_EPS)
bxU = BoundaryU()
bxU.mark(boundary_markers, 3)
        
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)


#defining stress and strain

def epsilon(u1):
    return 0.5*(nabla_grad(u1) + nabla_grad(u1).T)
    
def sigma(u1):
    return (lambda_*epsilon_p*nabla_div(u1) - u2 + chi*u6 + gamma*u5)*Identity(d) + 2*mu*epsilon_p*epsilon(u1)
    
    
#boundary values in the weak statement
T = Constant((0,0))
f = Constant((0,0))
n = FacetNormal(mesh)
g0 = Constant(0.)
g_p = Constant(0.)
g3 = -u9[1] - alpha*u3.dx(1).dx(1)
flux_copper = Constant(0.)
flux_acid = Constant(0.)
#boundary constants for pressure eq

flux_copper_sup = -xi*Cudiff/(delta)*dot(grad(u3),n)                                                                                     
flux_acid_sup = -xi*Hdiff/(delta)*dot(grad(u4),n)


        
#weak statement of the equations
                                                                                                                                
F1 = inner(sigma(u1), epsilon(v1))*dx-dot(T,v1)*ds(0) - dot(T,v1)*ds(2) - dot(T,v1)*ds(3)

F2 = epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(v2))*dx - dt*(ccr*g0*v2*ds(0) + ccr*g0*v2*ds(2) + ccr*g0*v2*ds(1) + ccr*g3*v2*ds(3) + ccr*inner(grad(u2),grad(v2))*dx - ccr*g_p*v2*ds(0) -ccr*g_p*v2*ds(1) - ccr*g_p*v2*ds(2))

F3 = (u3 + ccn*un3 + ccn_1*un_13 + ccn_2*un_23 + ccn_3*un_33 + ccn_4*un_43)*v3*dx + epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(u3))*v3*dx + dt*(-ccr*inner(grad(u2),grad(u3))*v3*dx +  ccr*xi*Cudiff*(inner(grad(u3),grad(v3)))*dx- ccr*r*u4*u5*v3*dx + ccr*r*u3*(phi0-2.*u5-u6)*an1(phi0-2.*u5-u6)*v3*dx  - ccr*flux_copper*v3*ds(1) - ccr*flux_copper*v3*ds(0) - ccr*flux_copper*v3*ds(2))

F4 = (u4 + ccn*un4 + ccn_1*un_14 + ccn_2*un_24 + ccn_3*un_34 + ccn_4*un_44)*v4*dx + epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(u4))*v4*dx + dt*(-ccr*inner(grad(u2),grad(u4))*v4*dx + ccr*xi*Hdiff*(inner(grad(u4),grad(v4)))*dx + ccr*r*u4*(phiac_hcl-u6)*v4*dx-ccr*flux_acid*v4*ds(1) - ccr*flux_acid*v4*ds(0) -ccr*flux_acid*v4*ds(2))

F5 = (u5 + ccn*un5 + ccn_1*un_15 + ccn_2*un_25 + ccn_3*un_35 + ccn_4*un_45)*v5*dx + dt*(ccr*r*u4*u5*v5*dx - ccr*r*u3*(phi0-2.*u5-u6)*an1(phi0-2.*u5-u6)*v5*dx)

F6 = (u6 + ccn*un6 + ccn_1*un_16 + ccn_2*un_26 + ccn_3*un_36 + ccn_4*un_46)*v6*dx - dt*(ccr*u4*(phiac_hcl-u6)*v6*dx)

F7 = (u7 + ccn*un7 + ccn_1*un_17 + ccn_2*un_27 + ccn_3*un_37 + ccn_4*un_47)*v7*dx + dt*(ccr*(zeta*Cudiff/delta*delta)*inner(grad(u7),grad(v7))*dx + ccr*inner(u9,grad(u7))*v7*dx - ccr*flux_copper*v7*ds(0) - ccr*flux_copper*v7*ds(1) - ccr*flux_copper*v7*ds(2) - ccr*flux_copper_sup*v7*ds(3))

F8 = (u8 + ccn*un8 + ccn_1*un_18 + ccn_2*un_28 + ccn_3*un_38 + ccn_4*un_48)*v8*dx + dt*(ccr*(zeta*Hdiff/delta*delta)*inner(grad(u8),grad(v8))*dx + ccr*inner(u9,grad(u8))*v8*dx - ccr*flux_acid*v8*ds(0) - ccr*flux_acid*v8*ds(1) - ccr*flux_acid*v8*ds(2) - ccr*flux_acid_sup*v8*ds(3))

F9 = beta*inner(grad(u9),grad(v9))*dx + u10*div(v9)*dx + div(u9)*v10*dx - inner(f,v9)*dx 

F = F1 + F2 + F3 + F4 + F5 + F6 + F7 + F8 + F9

#creating a VTK file for visualization output
vtkfile_u1 = File('./deform.pvd')
vtkfile_u2 = File('./pressure.pvd')
vtkfile_u3 = File('./phi_o.pvd')
vtkfile_u4 = File('./phi_o_hcl.pvd')
vtkfile_u5 = File('./phi_b.pvd')
vtkfile_u6 = File('./phi_b_hcl.pvd')
vtkfile_u7 = File('./phi_a.pvd')
vtkfile_u8 = File('./phi_a_hcl.pvd')
vtkfile_u9 = File('./sol_vel.pvd')
vtkfile_u10 = File('./sol_pressure.pvd')
                                                                                                                                      

#Boundary conditions
noflux = Constant((0.,0.)) 
bc1 = DirichletBC(V.sub(0), Constant((0.,0.)),bxD)
bc5 = DirichletBC(V.sub(8), Constant((0.,0.)),bxR)
bc6 = DirichletBC(V.sub(8), Constant((0.,0.)),bxL)
bc7 = DirichletBC(V.sub(9), Constant(0.),bxD)

t = 0
nn = 0

J = derivative(F,u,dv)

reflection = as_tensor([[1,0],
                        [0,-1]])


#time integration

while t < Tf:

   
    
    _un1,_un2,_un3,_un4,_un5,_un6,_un7,_un8,_un9,_un10 = un.split()   # saving files
        
    vtkfile_u1 << (_un1,t)
    vtkfile_u2 << (_un2,t)
    vtkfile_u3 << (_un3,t)
    vtkfile_u4 << (_un4,t)
    vtkfile_u5 << (_un5,t)
    vtkfile_u6 << (_un6,t)
    vtkfile_u7 << (_un7,t)
    vtkfile_u8 << (_un8,t)
    vtkfile_u9 << (_un9,t)
    vtkfile_u10 << (_un10,t)
    
    t  += float(dt)
    
    
    #continuity in the copper concentration as DirichletBC
    class Copper_cont(UserExpression):
        def eval(self,values,x):
            values[:] = u[7](x)
            
        
                
    copper_cont = Copper_cont()
    
    bc2 = DirichletBC(V.sub(2),copper_cont,bxU)
    
    #continuity in the acid concentration as DirichletBC
    class Acid_cont(UserExpression):
        def eval(self,values,x):
            values[:] = u[8](x)
    acid_cont = Acid_cont()  
              
    bc3 = DirichletBC(V.sub(3),acid_cont,bxU)
    
    #Osmotic pressure 
    u3_der = project(u[3].dx(1),V1,solver_type="cg")
    
    class Osmotic_p(UserExpression):
        def eval(self,values,x):
            values[:] = -alpha*u3_der(x)
            
                
    
    osmotic_p = Osmotic_p()
    
    bc4 = DirichletBC(V.sub(1),osmotic_p,bxU)
    
    
    #velocity continuity at the gel boundary
    vel_flip = reflection*(-grad(u2) + epsilon_p*((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41)/(dt*ccr)))
    
    rel_vel = project(vel_flip ,V2,solver_type="cg")
    class Solvent_vel(UserExpression):
        def eval(self,values,x):
            values[:] = rel_vel(x)
        def value_shape(self):
            return(2,)
            
    solvent_vel = Solvent_vel()
    bc8 = DirichletBC(V.sub(8), solvent_vel, bxU)
    
   
    bcs = [bc1,bc2,bc3,bc4,bc5,bc6,bc7,bc8]
    
    BDFparameters(nn, ccn, ccn_1, ccn_2, ccn_3, ccn_4, ccr)
    
    problem = NonlinearVariationalProblem(F,u,bcs,J=J)
    solver  = NonlinearVariationalSolver(problem)
    
    
    prm = solver.parameters
    prm['nonlinear_solver']='newton'
    prm["newton_solver"]["absolute_tolerance"]= 1e-7   
    prm["newton_solver"]["relative_tolerance"] = 1e-6 
    prm["newton_solver"]["maximum_iterations"] = 25
    prm['newton_solver']['relaxation_parameter'] = 1.0
    prm["newton_solver"]["error_on_nonconvergence"]=True
    prm['newton_solver']['linear_solver'] = 'superlu_dist'   
        
    
    solver.solve()
    
    nn += 1
        
    BDFvariables(nn, u, un, un_1, un_2, un_3, un_4)
            
    
    
t1_end = perf_counter()    

print(f"Total runtime in seconds =", t1_end-t1_start)

The output -

Process 0: Solving nonlinear variational problem.
Process 1: Solving nonlinear variational problem.
/home/chinobangs/.conda/envs/fenicsproject/lib/python3.9/site-packages/ufl/exproperators.py:336: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if arg in ("+", "-"):
2

Dear @c_katke,
There are quite a few issues with your code.
First of all:

  1. You should never define all variables inside your time loop. The code you have posted does 10 million time steps, and thus you are initializing your solver 10 million times.
  2. You should define all variables for a form outside a temporal loop. By constantly changing the form, you add extreme amounts of overhead to the code.

I’ve attempted to fix a few of the issue in the code at the end of this post. However, as you have not supplied a minimal code, I am reluctant to give you any further assistance. Please see the guidelines: Read before posting: How do I get my question answered?

#importing required libraries

from fenics import *


from ufl import tanh
from ufl import nabla_div
from ufl import exp
from time import perf_counter
from datetime import date

t1_start = perf_counter()

Tf = 1   #final time
dt = Constant(1e-7) #initial timestep

#constants
lambda_ = Constant(1.)
mu = Constant(1.)
epsilon_p = Constant(0.1)
chi = Constant(50.)
gamma = Constant(100.)
xi = Constant(1.)
zeta = Constant(100.)
Cudiff = Constant(10.)
Hdiff = Constant(10000.)
r = Constant(5.*pow(10., 5.))
phi0 = Constant(0.006)
phiac = Constant(0.003)
phiac_hcl = Constant(0.006)
phiac_hcl2 = Constant(0.006)
delta = Constant(2.)
alpha = Constant(1000.)
beta = Constant(1e-9)
ccn = Constant(-1.)
ccn_1 = Constant(0.)
ccn_2 = Constant(0.)
ccn_3 = Constant(0.)
ccn_4 = Constant(0.)
ccr = Constant(1.)


#defining the mesh
mesh = RectangleMesh(Point(0,0),Point(10,1),20,20)
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
P2 = VectorElement('P', mesh.ufl_cell(), 2) 
element = MixedElement([P2,P1,P1,P1,P1,P1,P1,P1,P2,P1])
V = FunctionSpace(mesh, element)
V1 = FunctionSpace(mesh,'P',1)
V2 = VectorFunctionSpace(mesh,'P',2)

v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 = TestFunctions(V)

u = Function(V)
un = Function(V)
un_1 = Function(V)
un_2 = Function(V)
un_3 = Function(V)
un_4 = Function(V)
du = TrialFunction(V)
dv = TrialFunction(V)
d = u.geometric_dimension()

#u1 = deformation
#u2 = pore pressure
#u3 = Free copper volume fraction in the gel
#u4 = Free acid volume fraction in the gel
#u5 = Bound copper volume fraction in the gel
#u6 = Bound acid volume fraction in the gel
#u7 = Free copper volume fraction in the supernatant domain
#u8 = Free acid volume fraction in the supernatant domain
#u9 = Flow velocity in the supernatant domain
#u10 = Pressure in the supernatant domain


#defining initial conditions


class InitialCondition(UserExpression):
    def eval(self, values, x):
        values[0] = 0.
        values[1] = -x[1]
        values[2] = 0.
        values[3] = 0.
        values[4] = 0.
        values[5] = phiac
        values[6] = 0.
        values[7] = 0.
        values[8] = (1/(2*1.7724))*exp(-(x[0]-5)*(x[0]-5)/4)*phiac_hcl2*0.5*(1-tanh(20*(x[1]-0.5)))   
        values[9] = 0.
        values[10] = 0.
        values[11] = 0.
    def value_shape(self):
        return(12,)
        
u0 = InitialCondition()
un = interpolate(u0,V)

u1, u2, u3, u4, u5, u6, u7, u8, u9, u10 = split(u); #splitting u
un1, un2, un3, un4, un5, un6, un7, un8, un9, un10 = split(un);
un_11, un_12, un_13, un_14, un_15, un_16, un_17, un_18, un_19, un_110 = split(un_1);
un_21, un_22, un_23, un_24, un_25, un_26, un_27, un_28, un_29, un_210 = split(un_2);
un_31, un_32, un_33, un_34, un_35, un_36, un_37, un_38, un_39, un_310 = split(un_3);
un_41, un_42, un_43, un_44, un_45, un_46, un_47, un_48, un_49, un_410 = split(un_4);



#defining bdf functions
def BDFparameters(nn, ccn, ccn_1, ccn_2, ccn_3, ccn_4, ccr):
    if nn==0:
        ccn.assign(Constant(-1.)); ccr.assign(Constant(1.)); ccn_1.assign(Constant(0.))
        ccn_2.assign(Constant(0.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))                
    elif nn==1:
        ccn.assign(Constant(-4./3.)); ccr.assign(Constant(2./3.)); ccn_1.assign(Constant(1./3.))
        ccn_2.assign(Constant(0.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))
    elif nn==2:
        ccn.assign(Constant(-18./11.)); ccr.assign(Constant(6./11.)); ccn_1.assign(Constant(9./11.))
        ccn_2.assign(Constant(-2./11.)); ccn_3.assign(Constant(0.)); ccn_4.assign(Constant(0.))
    elif nn==3:
        ccn.assign(Constant(-48./25.)); ccr.assign(Constant(12./25.)); ccn_1.assign(Constant(36./25.))
        ccn_2.assign(Constant(-16./25.)); ccn_3.assign(Constant(3./25.)); ccn_4.assign(Constant(0.))
    else:
        ccn.assign(Constant(-300./137.)); ccr.assign(Constant(60./137.)); ccn_1.assign(Constant(300./137.))
        ccn_2.assign(Constant(-200./137.)); ccn_3.assign(Constant(75./137.)); ccn_4.assign(Constant(-12./137.))

def BDFvariables(nn, u, un, un_1, un_2, un_3, un_4):
    if nn==0:
        un.assign(u)    
    elif nn==1:
        un_1.assign(un); un.assign(u)  
    elif nn==2:
        un_2.assign(un_1); un_1.assign(un); un.assign(u)  
    elif nn==3:
        un_3.assign(un_2); un_2.assign(un_1); un_1.assign(un); un.assign(u)  
    else:
        un_4.assign(un_3); un_3.assign(un_2); un_2.assign(un_1); un_1.assign(un); un.assign(u)
        

def an1(uud):
    return 1./2. *(1+tanh(pow(10,3)*uud))

        
#defining boundaries

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)

class BoundaryL(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)
bxL = BoundaryL()
bxL.mark(boundary_markers, 0) 
        
class BoundaryR(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 10, DOLFIN_EPS)
bxR = BoundaryR()
bxR.mark(boundary_markers, 2)
        
class BoundaryD(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0, DOLFIN_EPS)
bxD = BoundaryD()
bxD.mark(boundary_markers, 1)
        
class BoundaryU(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1, DOLFIN_EPS)
bxU = BoundaryU()
bxU.mark(boundary_markers, 3)
        
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)


#defining stress and strain

def epsilon(u1):
    return 0.5*(nabla_grad(u1) + nabla_grad(u1).T)
    
def sigma(u1):
    return (lambda_*epsilon_p*nabla_div(u1) - u2 + chi*u6 + gamma*u5)*Identity(d) + 2*mu*epsilon_p*epsilon(u1)
    
    
#boundary values in the weak statement
T = Constant((0,0))
f = Constant((0,0))
n = FacetNormal(mesh)
g0 = Constant(0.)
g_p = Constant(0.)
g3 = -u9[1] - alpha*u3.dx(1).dx(1)
flux_copper = Constant(0.)
flux_acid = Constant(0.)
#boundary constants for pressure eq

flux_copper_sup = -xi*Cudiff/(delta)*dot(grad(u3),n)                                                                                     
flux_acid_sup = -xi*Hdiff/(delta)*dot(grad(u4),n)


        
#weak statement of the equations
                                                                                                                                
F1 = inner(sigma(u1), epsilon(v1))*dx-dot(T,v1)*ds(0) - dot(T,v1)*ds(2) - dot(T,v1)*ds(3)

F2 = epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(v2))*dx - dt*(ccr*g0*v2*ds(0) + ccr*g0*v2*ds(2) + ccr*g0*v2*ds(1) + ccr*g3*v2*ds(3) + ccr*inner(grad(u2),grad(v2))*dx - ccr*g_p*v2*ds(0) -ccr*g_p*v2*ds(1) - ccr*g_p*v2*ds(2))

F3 = (u3 + ccn*un3 + ccn_1*un_13 + ccn_2*un_23 + ccn_3*un_33 + ccn_4*un_43)*v3*dx + epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(u3))*v3*dx + dt*(-ccr*inner(grad(u2),grad(u3))*v3*dx +  ccr*xi*Cudiff*(inner(grad(u3),grad(v3)))*dx- ccr*r*u4*u5*v3*dx + ccr*r*u3*(phi0-2.*u5-u6)*an1(phi0-2.*u5-u6)*v3*dx  - ccr*flux_copper*v3*ds(1) - ccr*flux_copper*v3*ds(0) - ccr*flux_copper*v3*ds(2))

F4 = (u4 + ccn*un4 + ccn_1*un_14 + ccn_2*un_24 + ccn_3*un_34 + ccn_4*un_44)*v4*dx + epsilon_p*inner((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41),grad(u4))*v4*dx + dt*(-ccr*inner(grad(u2),grad(u4))*v4*dx + ccr*xi*Hdiff*(inner(grad(u4),grad(v4)))*dx + ccr*r*u4*(phiac_hcl-u6)*v4*dx-ccr*flux_acid*v4*ds(1) - ccr*flux_acid*v4*ds(0) -ccr*flux_acid*v4*ds(2))

F5 = (u5 + ccn*un5 + ccn_1*un_15 + ccn_2*un_25 + ccn_3*un_35 + ccn_4*un_45)*v5*dx + dt*(ccr*r*u4*u5*v5*dx - ccr*r*u3*(phi0-2.*u5-u6)*an1(phi0-2.*u5-u6)*v5*dx)

F6 = (u6 + ccn*un6 + ccn_1*un_16 + ccn_2*un_26 + ccn_3*un_36 + ccn_4*un_46)*v6*dx - dt*(ccr*u4*(phiac_hcl-u6)*v6*dx)

F7 = (u7 + ccn*un7 + ccn_1*un_17 + ccn_2*un_27 + ccn_3*un_37 + ccn_4*un_47)*v7*dx + dt*(ccr*(zeta*Cudiff/delta*delta)*inner(grad(u7),grad(v7))*dx + ccr*inner(u9,grad(u7))*v7*dx - ccr*flux_copper*v7*ds(0) - ccr*flux_copper*v7*ds(1) - ccr*flux_copper*v7*ds(2) - ccr*flux_copper_sup*v7*ds(3))

F8 = (u8 + ccn*un8 + ccn_1*un_18 + ccn_2*un_28 + ccn_3*un_38 + ccn_4*un_48)*v8*dx + dt*(ccr*(zeta*Hdiff/delta*delta)*inner(grad(u8),grad(v8))*dx + ccr*inner(u9,grad(u8))*v8*dx - ccr*flux_acid*v8*ds(0) - ccr*flux_acid*v8*ds(1) - ccr*flux_acid*v8*ds(2) - ccr*flux_acid_sup*v8*ds(3))

F9 = beta*inner(grad(u9),grad(v9))*dx + u10*div(v9)*dx + div(u9)*v10*dx - inner(f,v9)*dx 

F = F1 + F2 + F3 + F4 + F5 + F6 + F7 + F8 + F9

#creating a VTK file for visualization output
vtkfile_u1 = File('./deform.pvd')
vtkfile_u2 = File('./pressure.pvd')
vtkfile_u3 = File('./phi_o.pvd')
vtkfile_u4 = File('./phi_o_hcl.pvd')
vtkfile_u5 = File('./phi_b.pvd')
vtkfile_u6 = File('./phi_b_hcl.pvd')
vtkfile_u7 = File('./phi_a.pvd')
vtkfile_u8 = File('./phi_a_hcl.pvd')
vtkfile_u9 = File('./sol_vel.pvd')
vtkfile_u10 = File('./sol_pressure.pvd')
                                                                                                                                      

#Boundary conditions
noflux = Constant((0.,0.)) 
bc1 = DirichletBC(V.sub(0), Constant((0.,0.)),bxD)
bc5 = DirichletBC(V.sub(8), Constant((0.,0.)),bxR)
bc6 = DirichletBC(V.sub(8), Constant((0.,0.)),bxL)
bc7 = DirichletBC(V.sub(9), Constant(0.),bxD)

t = 0
nn = 0

J = derivative(F,u,dv)

reflection = as_tensor([[1,0],
                        [0,-1]])



#continuity in the copper concentration as DirichletBC
class Copper_cont(UserExpression):
    def eval(self,values,x):
        values[:] = u[7](x)
        
    
            
copper_cont = Copper_cont()

bc2 = DirichletBC(V.sub(2),copper_cont,bxU)

#continuity in the acid concentration as DirichletBC
class Acid_cont(UserExpression):
    def eval(self,values,x):
        values[:] = u[8](x)
acid_cont = Acid_cont()  
            
bc3 = DirichletBC(V.sub(3),acid_cont,bxU)

#Osmotic pressure 
u3_der = project(u[3].dx(1),V1,solver_type="cg")

class Osmotic_p(UserExpression):
    def eval(self,values,x):
        values[:] = -alpha*u3_der(x)
        
            

osmotic_p = Osmotic_p()

bc4 = DirichletBC(V.sub(1),osmotic_p,bxU)


class Projector():
    """
    Projector class that reuses matrices and the left hand side vector for every solve.
    """
    def __init__(self, expr, uh , solver_type):
        self.solver_type = solver_type
        u = TrialFunction(uh.function_space())
        v = TestFunction(uh.function_space())
        a = inner(u, v) * dx
        self.L = inner(expr, v) * dx        
        self.A = assemble(a)
        self.b = assemble(self.L)
        self.uh = uh
    def project(self):
        assemble(self.L, tensor=self.b)
        solve(self.A, self.uh.vector(), self.b, self.solver_type)

rel_vel  = Function(V2)

class Solvent_vel(UserExpression):
    def eval(self,values,x):
        values[:] = rel_vel(x)
    def value_shape(self):
        return(2,)
        
solvent_vel = Solvent_vel()
bc8 = DirichletBC(V.sub(8), solvent_vel, bxU)


bcs = [bc1,bc2,bc3,bc4,bc5,bc6,bc7,bc8]

# Define solver
problem = NonlinearVariationalProblem(F,u,bcs,J=J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver']='newton'
prm["newton_solver"]["absolute_tolerance"]= 1e-7   
prm["newton_solver"]["relative_tolerance"] = 1e-6 
prm["newton_solver"]["maximum_iterations"] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
prm["newton_solver"]["error_on_nonconvergence"]=True
# prm['newton_solver']['linear_solver'] = 'superlu_dist'   



BDFparameters(nn, ccn, ccn_1, ccn_2, ccn_3, ccn_4, ccr)


vel_flip = reflection*(-grad(u2) + epsilon_p*((u1 + ccn*un1 + ccn_1*un_11 + ccn_2*un_21 + ccn_3*un_31 + ccn_4*un_41)/(dt*ccr)))
projector = Projector(vel_flip, rel_vel, solver_type="cg")


#time integration
while t < Tf:
    
    _un1,_un2,_un3,_un4,_un5,_un6,_un7,_un8,_un9,_un10 = un.split()   # saving files
    vtkfile_u1 << (_un1,t)
    vtkfile_u2 << (_un2,t)
    vtkfile_u3 << (_un3,t)
    vtkfile_u4 << (_un4,t)
    vtkfile_u5 << (_un5,t)
    vtkfile_u6 << (_un6,t)
    vtkfile_u7 << (_un7,t)
    vtkfile_u8 << (_un8,t)
    vtkfile_u9 << (_un9,t)
    vtkfile_u10 << (_un10,t)
    
    t  += float(dt)
    
    
 
    
    
    #velocity continuity at the gel boundary
    
    # Update BC     
    projector.project()

    solver.solve()
    
    nn += 1
        
    BDFvariables(nn, u, un, un_1, un_2, un_3, un_4)
            
    
    
t1_end = perf_counter()    

print(f"Total runtime in seconds =", t1_end-t1_start)
3

Dear @dokken,

Thank you so much for your help and I apologize for not providing a minimal code. The correction you had suggested were helpful and I was able to scale up the simulations to (400x200) mesh size. After implementation of the new code, I did notice a small error in the boundary condition. I have opened a new topic for the same as I think it might not be relevant to the title of this topic.

Thanks again,
Chinmay