Dear @c_katke,
There are quite a few issues with your code.
First of all:
- 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.
- 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)