Hello Everyone,
Let’s consider transport of the acid in the gel and the supernatant domain attached to each other. This transport is characterized by different diffusivity of the acid in the gel domain and in the supernatant domain. For simplicity we can ignore the advection term in the advection-diffusion equation. Hence in the gel domain -
boundary condition at the interface-
These boundary conditions are dynamic. These equations were simulated on a 2-d mesh.
The code -
#importing required libraries
from fenics import *
from ufl import tanh
from ufl import exp
Tf = 0.00002 #final time
dt = Constant(1e-7) #initial timestep
#constants
Hdiff_gel = Constant(10000.)
Hdiff_sup = Constant(1000000.)
#defining the mesh
mesh = RectangleMesh(Point(0,0),Point(10,1),40,20)
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
P2 = VectorElement('P', mesh.ufl_cell(), 2)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh, element)
V1 = FunctionSpace(mesh,'P',1)
v1, v2 = TestFunctions(V)
u = Function(V)
un = Function(V)
dv = TrialFunction(V)
#u1 = Free acid volume fraction in the gel
#u2 = Free acid volume fraction in the supernatant domain
#defining initial conditions
class InitialCondition(UserExpression):
def eval(self, values, x):
values[0] = 0.
values[1] = (1/(2*1.7724))*exp(-(x[0]-5)*(x[0]-5)/4)*0.006*0.5*(1-tanh(20*(x[1]-0.5)))
def value_shape(self):
return(2,)
u0 = InitialCondition()
un = interpolate(u0,V)
u1, u2 = split(u); #splitting u
un1, un2 = split(un);
#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)
#boundary values in the weak statement
n = FacetNormal(mesh)
flux_acid = Constant(0.)
flux_acid_sup = -Hdiff_gel*dot(grad(u1),n)
#weak statement of the equations
F1 = (u1 - un1)*v1*dx + dt*(Hdiff_gel*(inner(grad(u1),grad(v1)))*dx -flux_acid*v1*ds(1) - flux_acid*v1*ds(0) -flux_acid*v1*ds(2))
F2 = (u2 - un2)*v2*dx + dt*(Hdiff_sup*inner(grad(u2),grad(v2))*dx - flux_acid*v2*ds(0) - flux_acid*v2*ds(1) - flux_acid*v2*ds(2) - flux_acid_sup*v2*ds(3))
F = F1 + F2
#creating a VTK file for visualization output
vtkfile_u1 = File('boundary_check/acid_gel.pvd')
vtkfile_u2 = File('boundary_check/acid_sup.pvd')
t = 0
J = derivative(F,u,dv)
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)
#continuity in the acid concentration as DirichletBC
acid_cont_1 = Function(V1)
class Acid_cont(UserExpression):
def eval(self,values,x):
values[:] = acid_cont_1(x)
acid_cont = Acid_cont()
bc1 = DirichletBC(V.sub(0),acid_cont,bxU)
# Define solver
problem = NonlinearVariationalProblem(F,u,bc1,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'
projector_h = Projector(u2, acid_cont_1, solver_type="cg")
#time integration
while t < Tf:
_un1,_un2 = un.split() # saving files
vtkfile_u1 << (_un1,t)
vtkfile_u2 << (_un2,t)
t += float(dt)
projector_h.project()
solver.solve()
un.assign(u)
For a random point on the gel supernatant domain boundary, I plotted the \phi_{gel} and \phi_{sup} for first few hundred time steps and they should be exactly equal as implemented by DirichletBC but there is a slight mismatch. The graph for \phi_{gel} is shifted slightly to the right. I believe that this might be happening because solver is using \phi_{sup} from the previous timestep in the DirichletBC. How can I get around this issue? I’m running this code in parallel on HPC cluster.
Again thank you so much for your time and the help!