Dear. Researchers
I am trying to model 2D oxidation problem using reaction-diffusion approach. To achieve this goal, I have utilized fully implicit backward Euler for time-discretization as well as physical parameters like diffusion coefficient (D), reaction kinetics parameter(eta) and molarity. In this problem we have two time scales one for diffusion and another for reaction. Based on the parameters, the reaction possesses fast kinetics compared with diffusion. Accordingly, I choose time-step lower than reaction time-scale to resole this process. Upon, running the code, newton solvers executes two steps successfully but afterwards it diverges. I will be really appreciate if someone assist me to resolve the issue. I will put the code in the following:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# -----------------------------
# Time parameters
# -----------------------------
num_steps = 10000
t = 0.0
t_end = 100
dt = t_end/num_steps
dt_var = Constant(dt)
BC_thickness = 295e-6
Oxide_thickness = 5e-6
Lx = 60e-6
Ly = 300e-6
# -----------------------------
# Geometry
# -----------------------------
nx, ny = 300,1000
mesh = RectangleMesh(Point(0,0), Point(Lx,Ly),nx,ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(),1)
M = FunctionSpace(mesh,MixedElement([P1,P1]))
w = Function(M)
w_old = Function(M)
dw = TrialFunction(M)
c_new,n_new = split(w)
c_old, n_old = split(w_old)
v_c,v_n = TestFunctions(M)
# -----------------------------
# Parameters
# -----------------------------
D = Constant(2e-13)
eta = Constant(1e-5)
Molarity = Constant(1.11e5)
class K(UserExpression):
def eval(self,values,x):
ell = BC_thickness/100
delta = np.sqrt(2)*ell
values[0] = 0.0
if x[1] >= BC_thickness:
values[1] = 1.0
else:
values[1] = 0.0
# delta = 5*mesh.hmin()
# values[1] = 0.5*(1 + np.tanh((x[1] - BC_thickness)/delta))
def value_shape(self):
return(2,)
u_init = K(mesh)
w.interpolate(u_init)
file_1 = File("c_initial2D.pvd")
file_2 = File("n_initial2D.pvd")
file_1 << w.split(deepcopy= True)[0]
file_2 << w.split(deepcopy= True)[1]
w_old.assign(w)
def top(x,on_boundary):
return on_boundary and near(x[1], Ly)
def bot(x,on_boundary):
return on_boundary and near(x[1], 0.0)
c_bot = Constant(0.0)
c_top = Constant(1.55)
bc_c1 = DirichletBC(M.sub(0), c_top,top)
bc_c2 = DirichletBC(M.sub(0), c_bot,bot)
bcs_problem = [bc_c1,bc_c2]
# Form_c = ((c_new - c_old)/dt_var * v_c * dx
# + D*inner(grad(c_new), grad(v_c))*dx
# + Molarity*((n_new - n_old)/dt_var)*v_c*dx)
Form_c = ((c_new - c_old)/dt_var * v_c * dx
+ D*inner(grad(c_new), grad(v_c))*dx
+ Molarity*eta*(1-n_new)*c_new*v_c*dx)
Form_n = ((n_new - n_old)/dt_var * v_n * dx
- eta*c_new*(1 - n_new)*v_n*dx)
Form_total = Form_c + Form_n
Gain = derivative(Form_total, w,dw)
# -----------------------------
problem = NonlinearVariationalProblem(Form_total, w, bcs_problem,J = Gain)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters["newton_solver"]
prm["absolute_tolerance"] = 1e-6
prm["relative_tolerance"] = 1e-5
prm["maximum_iterations"] = 50
prm["linear_solver"] = "mumps"
prm["relaxation_parameter"] = 0.5
prm["report"] = True
file_c = File("results2D1/c_fully_implicit2D1.pvd")
file_n = File("results2D1/n_fully_implicit2D1.pvd")
# -----------------------------
# Time loop
# -----------------------------
front_positions = []
time_values = []
for n in range(num_steps):
t += dt
# solve(Form_total == 0, w , bcs_problem, J = Gain)
solver.solve()
c_sol, n_sol = w.split(deepcopy= True)
file_c << (c_sol,t)
file_n << (n_sol,t)
w_old.assign(w)
