Hello Everyone,
I’m trying to solve a Reaction-Diffusion-Migration Equation, where I have two concentrations a mobile concentration c_f and an immobile one c_b.
I’m using a MixedElement for the definition of my FunctionSpace. My problem is, that after the first time step to solution does not converge, because of that I was wondering whether I have a mistake in my solving loop.
Here is a snippet of the iterations for the first and 2. time-step:
And Here is a minimal example of my code:
from fenics import *
from decimal import Decimal
lam = 1*10**(-2) # Forward-Reactionrateconstant
mu = lam*10**(-3) # Backward-Reactionrateconstant
D_0 = float(Decimal(1*10**(-14))) # Diffusionskoeff. m^2/s
D = Constant(D_0)
c_ini_f = Constant(6.53) # Anfangs-Konz. der freien Konz. mol/m^3 -> Umrechnung von /m^3 in /mm^3
cb_max = Constant(42.5)
# Konstanten für Migrationsterm
z_f = 1; z_f = Constant(z_f) # Ladung = 3+
z_b = -1; z_b = Constant(z_b)
eps_r = 57 #für Muskel aus rahman2020
eps_0 = 8.854189*10**(-12) # As/Vm
eps = eps_0*eps_r; eps = Constant(eps)
R = 8.314; R = Constant(R) # Gaskonstant
Temp = 293; Temp=Constant(Temp) # Temperature in Kelvin
F = 9.648*10**4; F= Constant(F) # Faraday-Konstante
##############################################################################
numEle = 100
z_min = 0; z_max = 10*10**(-6)
mesh = IntervalMesh (numEle, z_min, z_max )
# Zeitvariablen
T = 5*60 # final time: Anzahl min in Sekunden
num_steps = T*10 # number of time steps
dt = T / num_steps # time step size
##############################################################################
# Definition der Funktionen #
##############################################################################
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 2) # für gebundene Konzentration
P2 = FiniteElement('Lagrange', mesh.ufl_cell(), 2) # für freie Konzentration
P3 = FiniteElement('Lagrange', mesh.ufl_cell(), 1) # für elektrostatisches Potential psi
# Make a mixed space
TH=MixedElement([P1,P2,P3])
V = FunctionSpace(mesh, TH)
(dc_b,dc_f, d_psi) = TestFunctions(V)
dc = TrialFunction(V)
c = Function(V)
(c_b, c_f, psi) = split(c) # psi = elektrische Potential
time_c = Function(V) # Hilfsgröße zum splitten
c0_b, c0_f, psi_0 = split(time_c)
def boundary(x, on_boundary):
return on_boundary
class oben(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],z_max,tol)
class unten(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],z_min,tol)
top = oben()
bottom = unten()
bc_f_t = DirichletBC(V.sub(1), Constant((c_ini_f)),top)
bc_psi_b = DirichletBC(V.sub(2), Constant(0), bottom)
bc_psi_t = DirichletBC(V.sub(2), Constant(0), top)
bc_psi_all = DirichletBC(V.sub(2), Constant(0), boundary)
bc_ges = [bc_f_t,bc_psi_t]
Diff_eqn = dc_f*(c_f-c0_f)/dt *dx + D*inner(nabla_grad(c_f)+z_f*c_f*F/(R*Temp)*nabla_grad(psi),nabla_grad(dc_f))*dx +\
dc_f*lam*c_f*(cb_max-c_b)*dx - dc_f*mu*c_b*dx
El_Pot_eqn= inner(nabla_grad(d_psi),nabla_grad(psi))*dx -\
d_psi*F/eps * (z_f*c_f + z_b*(cb_max-c_b))*dx
Reac_egn = dc_b*(c_b-c0_b)/dt *dx - dc_b*lam*c_f*(cb_max-c_b)*dx + dc_b*mu*c_b*dx
eqn = Diff_eqn + El_Pot_eqn + Reac_egn
Jac = derivative(eqn, c, dc)
t = 0 # Zeit
##############################################################################
# Lösungsschleife #
##############################################################################
for n in range(num_steps):
# Update current time
t += dt
solve(eqn == 0, c, bc_ges, J=Jac, solver_parameters={"newton_solver":
{"relative_tolerance": 1e-6,
"absolute_tolerance": 1e-15}})
(c_b,c_f,psi) = c.split(True)
# assign(time_c, [c_b,c_f,psi])
assign(time_c.sub(0), c_b); assign(time_c.sub(1), c_f)
If it helps: here are the equations iam using:
The concentration c_f is bound at specific sites. The concentration of this sites is modeled with the term: (max. amount of c_b - c_b).
The reaction is of the type: c_f + (max. amount of c_b - c_b) -> c_b
Thanks in advance!