Dear. Researchers
I hope you are doing well.
My current project is about modelling oxidation problem. To do so, I have developed a thermodynamically-consistent framework. Due to the complexity of the model and its full ingredients, I have decided to first simulate a simpler similar problem with the purpose of finding the optimal configuration of the solver. Accordingly, I have opted a paper published in Computational Materials Science in which FEM was used to solve reaction-diffusion problem within oxidation context. I implemented weak forms monolithically (listed in the paper) and applied associated boundary conditions as described in the paper. Upon, running the code, the solver becomes stagnated after a few iterations and I think something does not work properly. Despite dedicated efforts, I am unable to find the responsible bug in the code. I will be really appreciated for any assistance.
link to paper: Redirecting
in the following I present the full code:
# ============================================================
# Coupled chemo-mechanical oxidation model (Shen et al.)
# 2D plane strain, monolithic, dolfinx
# ============================================================
from mpi4py import MPI
import dolfinx
from dolfinx.mesh import create_rectangle
import ufl
import basix.ufl
from dolfinx.mesh import meshtags, locate_entities_boundary
from dolfinx.fem import functionspace, Function, locate_dofs_topological, Constant, locate_dofs_geometrical
from dolfinx.mesh import compute_midpoints
from dolfinx import geometry
import numpy as np
from dolfinx.io import XDMFFile
from dolfinx.cpp.mesh import CellType # not strictly required, but okay to import
from petsc4py import PETSc
from ufl.classes import *
from ufl.algorithms import *
from ufl import Measure
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
# ------------------------------------------------------------
# Geometry & mesh (single period, rectangle)
# ------------------------------------------------------------
BC_thickness = 100
Lx = 60 # wavelength
H_oxide = 1
Ly = BC_thickness + H_oxide # total height
nx,ny = 80, 300
domain = create_rectangle(
MPI.COMM_WORLD,
[[0.0, 0.0], [Lx, Ly]],
[nx, ny])
# ------------------------------------------------------------
# Material parameters (from Tables 1 & 2)
# ------------------------------------------------------------
E_BC = dolfinx.fem.Constant(domain, PETSc.ScalarType(110e9))
E_TGO = dolfinx.fem.Constant(domain,PETSc.ScalarType(320e9))
nu_BC = dolfinx.fem.Constant(domain, PETSc.ScalarType(0.33))
nu_TGO = dolfinx.fem.Constant(domain, PETSc.ScalarType(0.25))
beta = dolfinx.fem.Constant(domain, PETSc.ScalarType(0.05))
D = dolfinx.fem.Constant(domain, PETSc.ScalarType(2e-14))
eta = dolfinx.fem.Constant(domain, PETSc.ScalarType(1e-5))
L = dolfinx.fem.Constant(domain, PETSc.ScalarType(1.11e5))
Omega = dolfinx.fem.Constant(domain, PETSc.ScalarType(2.6e-5)) # molar volume
R = 8.314
T = dolfinx.fem.Constant(domain, PETSc.ScalarType(1273.0)) # 1000 C
dt = 1e-7 # 1 hour
# ------------------------------------------------------------
# Function spaces (monolithic)
# ------------------------------------------------------------
deg_u = 2
deg_s = 1
element_degree = 1
variant = basix.LagrangeVariant.equispaced
disp_el = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, variant, shape = (domain.topology.dim,))
cons_el = basix.ufl.element("Lagrange", domain.topology.cell_name(), element_degree, variant)
xi_el = basix.ufl.element("Lagrange", domain.topology.cell_name(), element_degree, variant)
trace_sigma = basix.ufl.element("Lagrange", domain.topology.cell_name(), element_degree, variant)
elements = basix.ufl.mixed_element([disp_el, cons_el, xi_el, trace_sigma])
W = dolfinx.fem.functionspace(domain, elements)
w = dolfinx.fem.Function(W) # current solution
w_old = dolfinx.fem.Function(W) # previous time step
dw = ufl.TrialFunction(W)
(v_u, v_c, v_n, v_tr_sig) = ufl.TestFunctions(W)
(u, c, n, tr_sig) = ufl.split(w)
(u_old, c_old, n_old,tr_sig_old) = ufl.split(w_old)
w.sub(0).x.array[:] = 0.0
w.sub(1).x.array[:] = 0.0
w.sub(2).x.array[:] = 0.0
w.sub(3).x.array[:] = 0.0
x = ufl.SpatialCoordinate(domain)
#ell = BC_thickness/200 #main
ell = BC_thickness/200
delta = np.sqrt(2) * ell
# max_raw = 0.5*(1 + np.tanh((ytop - y0)/delta))
#w.sub(2).interpolate(lambda x:0.5*(1 + np.tanh((x[1] - y0)/delta)) / max_raw)
w.sub(2).interpolate(lambda x: 0.5*(1 + np.tanh((x[1] - BC_thickness)/delta)))
w.x.scatter_forward()
# w.x.scatter_forward()
with XDMFFile(domain.comm, "phi_function_xi.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(w.sub(2))
# ------------------------------------------------------------
# Elasticity (plane strain)
# ------------------------------------------------------------
mu_BC = E_BC/(2*(1 + nu_BC))
mu_TGO = E_TGO/(2*(1 + nu_TGO))
lmbda_BC = E_BC * nu_BC / ((1.0 + nu_BC) * (1.0 - 2.0 * nu_BC))
lmbda_TGO = E_TGO * nu_TGO / ((1.0 + nu_TGO) * (1.0 - 2.0 * nu_TGO))
def eps(u):
return ufl.sym(ufl.grad(u))
def sigma(u_, c_, n_):
mu_eff = (1-n_)*mu_BC + n_*mu_TGO
lmbda_eff = (1-n_)*lmbda_BC + n_*lmbda_TGO
I = ufl.Identity(2)
chem_strain = beta * c_ * I
return lmbda_eff * ufl.tr(eps(u_) - chem_strain) * I + 2.0 * mu_eff * (eps(u_) - chem_strain)
def mean_stress(u, c):
return ufl.tr(sigma(u, c)) / 3.0
# ------------------------------------------------------------
# Weak forms
# ------------------------------------------------------------
# Mechanical equilibrium
F_mech = ufl.inner(sigma(u, c, n), eps(v_u)) * ufl.dx
# Diffusion flux with stress coupling
#coupling_coeff = dolfinx.fem.Constant(domain,PETSc.ScalarType(2e-23))
flux = -D * (ufl.grad(c) - (Omega / (R * T)) * c * ufl.grad(tr_sig))
F_diff = (
(c - c_old) / dt * v_c * ufl.dx
- ufl.inner(flux, ufl.grad(v_c)) * ufl.dx
+ L * eta *(1 - n) * c * v_c * ufl.dx
)
# Reaction (TGO evolution)
F_reac = (
(n - n_old) / dt * v_n * ufl.dx
- eta*(1 - n) * c * v_n * ufl.dx
)
F_sigma = tr_sig*v_tr_sig*ufl.dx - (1.0/3.0)*ufl.tr(sigma(u,c,n))*v_tr_sig*ufl.dx
F = F_mech + F_diff + F_reac + F_sigma
# ------------------------------------------------------------
# Boundary conditions
# ------------------------------------------------------------
# Bottom: v = 0
def bottom(x):
return np.isclose(x[1], 0.0)
u_y_D = PETSc.ScalarType(0.0)
bot_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim-1, bottom)
dofs_bc_uy = dolfinx.fem.locate_dofs_topological(W.sub(0).sub(1), domain.topology.dim-1, bot_facets)
bc_uy = dolfinx.fem.dirichletbc(u_y_D,dofs_bc_uy, W.sub(0).sub(1))
tags_bot_facets = np.full_like(bot_facets, 1)
bot_meshtags = meshtags(domain, domain.topology.dim-1, bot_facets, tags_bot_facets)
with XDMFFile(domain.comm, "shen_bot_factes_bc.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(bot_meshtags, domain.geometry)
def left(x):
return np.isclose(x[0], 0.0)
u_x_D = PETSc.ScalarType(0.0)
left_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim-1, left)
dofs_bc_ux = dolfinx.fem.locate_dofs_topological(W.sub(0).sub(0), domain.topology.dim-1, left_facets)
bc_ux = dolfinx.fem.dirichletbc(u_x_D,dofs_bc_ux, W.sub(0).sub(0))
tags_left_facets = np.full_like(left_facets, 1)
left_meshtags = meshtags(domain, domain.topology.dim-1, left_facets, tags_left_facets)
with XDMFFile(domain.comm, "shen_left_factes_bc.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(left_meshtags, domain.geometry)
# Oxygen concentration BCs
def top(x):
return np.isclose(x[1], Ly)
c_top_D = PETSc.ScalarType(1.55)
c_bot_D = PETSc.ScalarType(0.0)
t_sig_D = PETSc.ScalarType(0.0)
top_facets = dolfinx.mesh.locate_entities_boundary(domain,domain.topology.dim-1, top)
dofs_c_vox_top = dolfinx.fem.locate_dofs_topological(W.sub(1), domain.topology.dim-1,top_facets)
dofs_c_vox_bot = dolfinx.fem.locate_dofs_topological(W.sub(1), domain.topology.dim-1,bot_facets)
dofs_tr_top = dolfinx.fem.locate_dofs_topological(W.sub(3), domain.topology.dim-1,top_facets)
bc_cvox_top = dolfinx.fem.dirichletbc(c_top_D, dofs_c_vox_top, W.sub(1))
bc_cvox_bot = dolfinx.fem.dirichletbc(c_bot_D, dofs_c_vox_bot, W.sub(1))
bc_tr_sig = dolfinx.fem.dirichletbc(t_sig_D, dofs_tr_top, W.sub(3))
tags_top_facets = np.full_like(top_facets, 1)
top_meshtags = meshtags(domain, domain.topology.dim-1, top_facets, tags_top_facets)
with XDMFFile(domain.comm, "shen_top_factes_bc.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(top_meshtags, domain.geometry)
bcs = [bc_ux, bc_uy, bc_cvox_top, bc_cvox_bot, bc_tr_sig]
def corner(x):
return np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0)
corner_facets = locate_entities_boundary(domain, 0, corner)
dofs_corner = locate_dofs_topological(W.sub(0).sub(0), 0, corner_facets)
bc_corner = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corner, W.sub(0).sub(0))
bcs.append(bc_corner)
Gain = ufl.derivative(F, w, dw)
# ------------------------------------------------------------
# Nonlinear problem
# ------------------------------------------------------------
problem = NonlinearProblem(F, w, bcs = bcs , J = Gain)
solver = NewtonSolver(domain.comm, problem)
solver.atol = 1e-6
solver.rtol = 1e-5
solver.max_it = 50
#solver.line_search = "bt"
solver.relaxation_parameter = 0.8
solver.report = True
ksp = solver.krylov_solver
ksp.setType("gmres")
ksp.getPC().setType("lu")
opts = PETSc.Options()
prefix = ksp.getOptionsPrefix()
opts[f"{prefix}ksp_monitor"] = ""
opts[f"{prefix}ksp_converged_reason"] = ""
# Now apply options
ksp.setFromOptions()
w_old.x.array[:] = w.x.array
with XDMFFile(domain.comm, "visU.xdmf", "w") as xdmfu:
xdmfu.write_mesh(domain)
with XDMFFile(domain.comm, "visC.xdmf", "w") as xdmfc:
xdmfc.write_mesh(domain)
with XDMFFile(domain.comm, "visN.xdmf", "w") as xdmfn:
xdmfn.write_mesh(domain)
with XDMFFile(domain.comm, "visTr.xdmf", "w") as xdmftr:
xdmftr.write_mesh(domain)
# ------------------------------------------------------------
# Time loop
# ------------------------------------------------------------
t = 0.0
T_end = 10000*dt
step = 0
while t < T_end:
t += dt
step += 1
n_iter, converged = solver.solve(w)
assert converged
w.x.scatter_forward()
w_old.x.array[:] = w.x.array
u_out, c_out, n_out, tr_sig_out = w.split()
xdmfu.write_function(u_out,t)
xdmfc.write_function(c_out, t)
xdmfn.write_function(n_out, t)
xdmftr.write_function(tr_sig_out, t)
xdmfu.close()
xdmfc.close()
xdmfn.close()
xdmftr.close()