Hello,
I have a problem with the implementation of an implicit adaptive time stepping controller. As a first test, I wanted to advance the solution by 1/2 dt, then by dt, and subsequently calculate the relative L2_error resulting from the solutions at these two “trial steps”. The optimal time step is then calculated as dt *= 1e-8 / L2_error. If I run this model it reaches the final time T = 1e-3 s after only 7 steps (starting from t = 1e-11 s). However, if I use a constant time step (e.g. by simply removing dt *= 1e-8 / L2_error) I see that the solution changes quite heavily on small time steps, i.e. the low L2_error and the corresponding large optimum time steps are somehow wrong. I attached the following working example:
from cmath import pi
import numpy as np
from ufl import (SpatialCoordinate, FiniteElement, MixedElement, split, TestFunctions, ds, dx, grad, inner, ln, exp)
from dolfinx import fem, mesh, io
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.fem import FunctionSpace, Function, petsc
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from scipy.constants import elementary_charge as e
from scipy.constants import (epsilon_0, m_p, m_e, k)
# General parameters
t = 0.0
T = 1e-3
dt = 1e-11
n_n = 1e21 # density neutrals in units of m**-3
T_n = 300 # temperature neutrals in units of K
T_i = 300 # temperature ions in units of K
m_i = 40 * m_p
# Initial conditions
def V_0(x):
values = np.zeros((1, x.shape[1]))
values[0] = 0.0
return values
def ni_0(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1e15
return values
def ne_0(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1e15
return values
def Ee_0(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1.5 * 1e15 * e * 3
return values
def generate_mesh(x_left, x_right):
# uniform mesh
NumberOfCells = int(1e3)
h = (x_right - x_left) / NumberOfCells # cell size for equidistant mesh
msh = mesh.create_interval(MPI.COMM_WORLD, NumberOfCells, points=(x_left, x_right))
return msh
# domain boundaries
x_left = -0.1
x_right = +0.1
msh = generate_mesh(x_left, x_right)
x = SpatialCoordinate(msh)
CG1_elem = FiniteElement("CG", msh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(msh, ME_elem)
V_test, ni_test, ne_test, Ee_test = TestFunctions(ME)
# Function at time t_i
u = Function(ME)
# Function at time t_(i-1)
u_ = Function(ME)
# Apply initial conditions
u.sub(0).interpolate(V_0)
u.sub(1).interpolate(ni_0)
u.sub(2).interpolate(ne_0)
u.sub(3).interpolate(Ee_0)
u_.sub(0).interpolate(V_0)
u_.sub(1).interpolate(ni_0)
u_.sub(2).interpolate(ne_0)
u_.sub(3).interpolate(Ee_0)
V, ni, ne, Ee = split(u)
V_, ni_, ne_, Ee_ = split(u_)
Te = 2*Ee / (3*ne*e)
Te_ = 2*Ee_ / (3*ne_*e)
# Left boundary
facet_l_bdry = mesh.locate_entities_boundary(msh, dim=0, marker=lambda x: np.isclose(x[0], x_left))
dof_l_bdry = fem.locate_dofs_topological(V=ME.sub(0), entity_dim=0, entities=facet_l_bdry)
# Right boundary
facet_r_bdry = mesh.locate_entities_boundary(msh, dim=0, marker=lambda x: np.isclose(x[0], x_right))
dof_r_bdry = fem.locate_dofs_topological(V=ME.sub(0), entity_dim=0, entities=facet_r_bdry)
# rate coefficients for elastic collisions
X_en = 1e-6*exp(-15.70627755836476+-0.00436884406897324*ln(Te)**1+-0.1135732340628103*ln(Te)**2+-0.0011873155516964292*ln(Te)**3)
X_in = 1e-15
# rate coefficients for ionization
X_ion = 10*1e-6*exp(-35.3071744514+16.9033721249*ln(Te)**1+-7.81316801032*ln(Te)**2+2.47737696605*ln(Te)**3+-0.62396194653*ln(Te)**4+0.145044251711*ln(Te)**5+-0.0286523207942*ln(Te)**6+0.00352120706732*ln(Te)**7+-0.000181868246303*ln(Te)**8)
# Weak expression of Poisson equation
RHS_PoissonEq = -e * ( ni - ne ) / epsilon_0
F = ( -inner(grad(V), grad(V_test)) - V_test * RHS_PoissonEq ) * dx
# Boundary conditions for Poisson equation
PoissonEq_dirichlet_bc_l = fem.dirichletbc(value=ScalarType(0), dofs=dof_l_bdry, V=ME.sub(0))
PoissonEq_dirichlet_bc_r = fem.dirichletbc(value=ScalarType(0), dofs=dof_r_bdry, V=ME.sub(0))
# Weak expression of particle transport equation for ni
Ri = X_ion*n_n*ne
mui = e/(m_i*(0.5*X_in*n_n))
ui = mui*(-V.dx(0) - k/e*T_i*ni.dx(0)/ni)
F += ni_test*(ni-ni_) * dx
F += -dt* ni_test*Ri * dx
F += -dt* ni_test.dx(0)*(ni*ui) * dx
# Non-homogeneous Neumann BCs
F += dt* ni_test*(ni*abs(ui)) * ds
# Weak expression of particle transport equation for ne
Re = X_ion*n_n*ne
nuelen = X_en * n_n
mue = e/(m_e*nuelen)
ue = mue*(+V.dx(0) - Te*ne.dx(0)/ne - Te.dx(0))
F += ne_test*(ne-ne_) * dx
F += -dt* ne_test*Re * dx
F += -dt* ne_test.dx(0)*(ne*ue) * dx
# Non-homogeneous Neumann BCs
vthe = (8*e*Te/pi/m_e)**0.5
F += dt* ne_test*(0.25*ne*vthe) * ds
# Weak expression of energy transport equation for electrons
mobe = e/(m_e*X_en*n_n)
kappae = 2.5*mobe*ne*Te
qex = -kappae*e*Te.dx(0)
Qex = 2.5*e*Te*ne*ue + qex
elastic_losses = ne*m_e/m_p * ( 1.5*X_en*n_n*(e*Te-k*T_n) )
Pext = 1e7
F += Ee_test*(Ee-Ee_) * dx
F += -dt* Ee_test*Pext * dx
F += -dt* Ee_test*(-e*inner(ue,-V)) * dx
F += -dt* Ee_test.dx(0)*Qex * dx
# # Non-homogeneous Neumann BCs
F += dt* Ee_test*(1/3*Ee*vthe) * ds
with io.XDMFFile(MPI.COMM_WORLD, "output/output.xdmf", "w") as file:
file.write_mesh(msh)
while (t < T):
# calculate with dt/2
dt *= 0.5
problem = petsc.NonlinearProblem(F, u, bcs=[PoissonEq_dirichlet_bc_l, PoissonEq_dirichlet_bc_r])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.rtol = 1e-6
n, converged = solver.solve(u)
assert(converged)
est_sol_dt_half = u.x.array.copy()
# recalculate using dt
dt *= 2
problem = petsc.NonlinearProblem(F, u, bcs=[PoissonEq_dirichlet_bc_l, PoissonEq_dirichlet_bc_r])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(u)
assert(converged)
est_sol_dt = u.x.array.copy()
# determine optimal time step
L2_error = np.linalg.norm(est_sol_dt - est_sol_dt_half) / np.linalg.norm(est_sol_dt_half)
dt *= 1e-8 / L2_error
# # calculate solution again using optimal time step
u.x.array[:] = u_.x.array
problem = petsc.NonlinearProblem(F, u, bcs=[PoissonEq_dirichlet_bc_l, PoissonEq_dirichlet_bc_r])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(u)
assert(converged)
u_.x.array[:] = u.x.array
t += dt
#print(t)
print(L2_error)
V_plot = u.sub(0).collapse()
file.write_function(V_plot, t)
Any help is highly appreciated.
