Issue with scaling up the simulation while running in parallel

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)