I am so grateful for your response. I have modified the solver and now it runs but the convergence is largely deteriorated compared with the Fenics-version of the code. please note that I do not change the specification and structure of Fenics counterpart while transferring the code into Fenicsx. Besides, the types of linear solver and line search method is MUMPS and “bt” for both cases. However, after running the code I have received converged_reason equal -8 which implies convergence to local minimum of weak form. During the debugging process, I have recognized that Dirichlet boundary conditions are not properly enforced. I applied the procedure described in the forum but it was unsuccessful. I will put the selected part of the code. Any assistance will be appreciated.
‘’'python
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
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
import os
from pathlib import Path
from mpi4py import MPI
from dolfinx import default_real_type, log, plot
from petsc4py.PETSc import ScalarType
-----------------------------
Time parameters
-----------------------------
num_steps = 72000
t = 0.0
t_end = 3600
dt = t_end/num_steps
BC_thickness = 100e-6
Oxide_thickness = 1e-6
Lx = 60e-6
Ly = BC_thickness + Oxide_thickness
refine_length = BC_thickness - 10e-6
-----------------------------
Parameters
-----------------------------
x_bar = Ly
BC_bar = BC_thickness/x_bar
Oxide_bar = Oxide_thickness/x_bar
Lx_bar = Lx/x_bar
Ly_bar = Ly/x_bar
refine_length_bar = refine_length/x_bar
print(refine_length_bar)
c_characteristic = 0.32e6
-----------------------------
Geometry
-----------------------------
nx, ny = 300,1000
domain = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([Lx_bar,Ly_bar])], [nx,ny])
tdim = domain.topology.dim
def refine_cells(x, tol=1e-16):
return x[1] >= refine_length_bar + tol
cells2 = dolfinx.mesh.locate_entities(domain,tdim,refine_cells)
print(“cells shape”,np.shape(cells2))
domain.topology.create_connectivity(tdim,tdim-1)
domain_topo = domain.topology
associated_edges = dolfinx.mesh.compute_incident_entities(domain_topo,cells2,2,1)
refined_domain,,= dolfinx.mesh.refine(domain,associated_edges)
with XDMFFile(MPI.COMM_WORLD,“refined_domain_scaled.xdmf”,“w”) as xdmf:
xdmf.write_mesh(refined_domain)
element_degree = 1
variant = basix.LagrangeVariant.equispaced
c_el = basix.ufl.element(“Lagrange”, refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)
n_el = basix.ufl.element(“Lagrange”, refined_domain.topology.cell_name(), element_degree, variant, dtype = default_real_type)
total_elements = basix.ufl.mixed_element([c_el,n_el])
ME = dolfinx.fem.functionspace(refined_domain, total_elements)
w = dolfinx.fem.Function(ME)
w.name = “WW”
w_old = dolfinx.fem.Function(ME)
dw = ufl.TrialFunction(ME)
c_new,n_new = ufl.split(w)
c_old,n_old = ufl.split(w_old)
w.sub(0).x.array[:] = 0.0
w.sub(1).x.array[:] = 0.0
x= ufl.SpatialCoordinate(refined_domain)
w.sub(1).interpolate(lambda x: np.where(x[1] >= BC_bar, 1.0, 0.0))
def top(x):
return np.isclose(x[1], Ly_bar)
def bottom(x):
return np.isclose(x[1], 0.0)
c_space, _ = ME.sub(0).collapse()
c_bc_top = dolfinx.fem.Function(c_space)
c_bc_top.x.array[:] = 1.55
top_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, tdim-1, top)
top_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(0), c_space), tdim-1, top_facets)
c_space_bot, _ = ME.sub(0).collapse()
c_bc_bot = dolfinx.fem.Function(c_space)
c_bc_bot.x.array[:] = 0.0
bottom_facets = dolfinx.mesh.locate_entities_boundary(refined_domain, tdim-1, bottom)
bot_dofs = dolfinx.fem.locate_dofs_topological((ME.sub(0), c_space), tdim-1, bottom_facets)
bc_c_top = dolfinx.fem.dirichletbc(c_bc_top, top_dofs,c_space)
bc_c_bot = dolfinx.fem.dirichletbc(c_bc_bot,bot_dofs,c_space)
bcs_problem = [bc_c_top, bc_c_bot]
‘’’