 # Nonlinear Poisson equation coupled with continuity equations

Hello,
starting from the example code given for the nonlinear Poisson equation by nate here, I wanted to extend this model to solve 2 continuity equations, (where the right hand sides depend on the dependent variable ne) along with the Poisson equation. However, it seems that I am doing something wrong, since all Newton iterations give me:
r (abs) = 0,000e+00 (tol = 1,000e-10) r (rel) = -nan (tol = 1,000e-09)
A minimal example of the code looks as follows:

from fenics import *

class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)

def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)

class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("ksp_monitor")
PETScOptions.set("pc_type", "ilu")

self.linear_solver().set_from_options()

# Create mesh
R = 1
nx = 100
mesh = IntervalMesh(nx, 0, R)

# Create spaces
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(mesh, ME_elem)
Phi_test, ni_test, ne_test = TestFunctions(ME)
u = Function(ME)
Phi, ni, ne = split(u)
u0 = Function(ME)
Phi0, ni0, ne0 = split(u0)

# Initial values
u = interpolate(Constant((0, 1, 1)), ME)
u0 = interpolate(Constant((0, 1, 1)), ME)

# Time stepping
t = 0.0
dt = 1
T = 2*dt

# Variational formulation of Poisson equation incl. b. c. --> Phi
f = ( ni - ne )

# Variational formulation of continuity equations for ions and electrons
# Time discretization done by implicit Euler scheme
R = ne
ui = -Dx(Phi,0) - Dx(ni,0)
ue = +Dx(Phi,0) - Dx(ne,0)
F += ( ni*ni_test - dt*ni*ui*Dx(ni_test, 0) - (ni0 + dt*R)*ni_test )*dx
F += ( ne*ne_test - dt*ne*ue*Dx(ne_test, 0) - (ne0 + dt*R)*ne_test )*dx

while (t < T):
J = derivative(F, u)
problem = Problem(J, F, [])
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
u0.vector()[:] = u.vector()
t += dt


Hi domi,

I know nothing about the physics, but I have to pieces of advice to give:

1. Your variational formulation considers no previous iteration, only the initial one. In addition, this ni=ni0 and ne=ne0, so I would expect 0 to be a solution of your problem. If this is the case, the errors you see are fine, as it find the solution on the first iteration and it is equal to the initial guess.
2. Following the previous comment, modify the initial conditions to something without a null derivative to activate the transport terms ni*Dx(ni), etc. This should make a difference.

Hope it helps,
Nico

1 Like

Hi Nico,

thank you for your answer, it helped me to fix the issue. Cheers,

domi

Please post your solution as it might help others with a similar issue:)

Hi dokken,
sure, no problem. In order to be of some (practical) use for other users, I extended the example to get a well posed (1d, cylindrically symmetric) PDE system that has a stationary state, i.e. I added another equation for the electron temperature and some boundary conditions to check that the solver is running correctly. The equation system that is now solved reads:
\partial_t n_i + \nabla \cdot n_i u_i = R_i
\partial_t n_e + \nabla \cdot n_e u_e = R_e
\partial \frac{3}{2} n_e T_e + \nabla \cdot Q_e = P_\mathrm{ext}
\nabla^2 \Phi = f,
with symmetry B.C. at 0 for all PDEs and Neumann B.C. for the continuity equations and the energy balance at R_m and a zero Dirichlet B.C. for the potential \Phi. u_i, u_e, Q_e, P_\mathrm{ext} and f are defined in the (almost) MWE, that reads:

from fenics import (NonlinearProblem, NewtonSolver, IntervalMesh, FiniteElement,
MixedElement, assemble, FunctionSpace, TestFunctions, Function,
interpolate, Expression, split, inner, grad, pi, dx, DirichletBC,
Constant, exp, ln, Dx, sqrt, ds, derivative, PETScKrylovSolver,
PETScFactory, near, PETScOptions, assign, File, plot)
import numpy as np
from scipy.constants import (epsilon_0, elementary_charge, m_p, m_e, k)
import matplotlib.pyplot as plt

class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)

def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)

class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "ilu")
self.linear_solver().set_from_options()

# Length 1D simulation domain (cylinrical sym. assumed)
R_m = 0.1

# Create mesh and define function space
nx = int(1e3) # cells
mesh = IntervalMesh(nx, 0, R_m)
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(mesh, ME_elem)

# Some parameters
dt = 1e-11
nm = 1e21 # Density neutrals 1/m**3

def boundary_right(x, on_boundary):
tol = 1E-14
return on_boundary and near(x, R_m, tol)

Phi_test, ni_test, ne_test, Te_test = TestFunctions(ME)

u  = Function(ME)
u0 = Function(ME)
# Initial values for Phi, ni, ne, Te
u  = interpolate(Constant((0, 1e15, 1e15, 1)), ME)
u0 = interpolate(Constant((0, 1e15, 1e15, 1)), ME)

Phi, ni, ne, Te = split(u)
Phi0, ni0, ne0, Te0 = split(u0)

x = Expression('x', degree = 1)

# Variational formulation of Poisson equation incl. b. c. --> Phi
f = -elementary_charge * ( ni - ne) / epsilon_0
bc_Poisson_right = DirichletBC(ME.sub(0), Constant(0), boundary_right)

# Variational formulation of continuity equation for ions incl. b. c. --> ni
X_Ioniz_H2_ndiss = 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)
R = X_Ioniz_H2_ndiss*nm*ne*1e5
Ti_K = 300
X2H2 = 1e-15
mui = elementary_charge/(2*m_p*(1/2*X2H2*nm))
ui = mui*(-Dx(Phi,0) - k/elementary_charge*Ti_K*Dx(ni,0)/ni)
F += ( ni*ni_test - dt*ni*ui*Dx(ni_test,0) - (ni0 + dt*R)*ni_test )*2*pi*x*dx
# Non-homogeneous Neumann BC at R_m
F += dt*ni_test*ni*ui*2*pi*x*ds

# Variational formulation of continuity equation for electrons incl. b. c. --> ne
X_eH2 = 1e-6*exp(-16.0325527377+0.412023668827*ln(Te)**1+-0.193551110298*ln(Te)**2+-0.0247678479005*ln(Te)**3+0.0105040269085*ln(Te)**4+-0.000152726360717*ln(Te)**5+-0.000538801463715*ln(Te)**6+9.70545683084e-05*ln(Te)**7+-5.10345827257e-06*ln(Te)**8)
nueleH2 = X_eH2*nm
mue = elementary_charge/(m_e*nueleH2)
ue = mue*(Dx(Phi,0) - Te*Dx(ne,0)/ne)
F += ( ne*ne_test - dt*ne*ue*Dx(ne_test,0) - (ne0 + dt*R)*ne_test )*2*pi*x*dx
# Thermal flux at the boundary
vthe = sqrt(8*elementary_charge*Te/(pi*m_e))
F += dt*ne_test*ne*0.25*vthe*2*pi*x*ds

# Variational formulation of power balance for electrons
kappae = 2.5*mue*ne*Te
qer = -kappae*elementary_charge*Dx(Te,0)
Qer = 2.5*elementary_charge*Te*ne*ue + qer
Pext = 1e5*exp(-1e6*(x-0.85*R_m)**4)
P  = -Pext
F += ( 1.5*ne*Te*Te_test - dt*Qer*Dx(Te_test,0) - (1.5*ne0*Te0 + dt*P)*Te_test )*2*pi*x*dx
F += dt*Te_test*2*elementary_charge*Te*ne*0.25*vthe*2*pi*x*ds

# Output file
#file = File("solution/output.pvd")

# Create VTK files for visualization output
vtkfile_Phi = File('solution/Phi.pvd')
vtkfile_ni = File('solution/ni.pvd')
vtkfile_ne = File('solution/ne.pvd')
vtkfile_Te = File('solution/Te.pvd')

# Time-stepping
t = 0.0
T = 1e-3
while (t < T):
t += dt
dt *= 1.2
u0.vector()[:] = u.vector()
J = derivative(F, u)
problem = Problem(J, F, [bc_Poisson_right])
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
_Phi, _ni, _ne, _Te = u.split()
vtkfile_Phi << (_Phi, t)
vtkfile_ni << (_ni, t)
vtkfile_ne << (_ne, t)
vtkfile_Te << (_Te, t)
#file << (u.split(), t)


For the time stepping, I chose the implicit euler scheme and I am able to ‘excite’ the system by the non-zero flux boundary conditions, i.e. I chose constant values for my initial guess instead of expressions as suggested by Nico.
I know from comsol, that for this system with the chosen parameters, a stationary state with n_{i,max} = n_{e,max} = 4.5e17, T_{e,max} = 4 and \Phi_{max} = 16 should be reached after roughly T=10^{-4}\,s.
However, the solver in the example above does not seem to reach a stationary state at all and n_i and n_e are somehow clamped to the initial value. Maybe this has to do with Nico’s first comment, that my code does not consider the previous iteration, but I don’t see why that should be the case, since u0 is updated by u in each time step with

u0.vector()[:] = u.vector().

Any hints on what I do wrong here are very welcome. Thank you in advance for your help, I really appreciate it.

domi

You should use u0.assign(u) instead of moving the vector. I haven’t looked through all the code as it is quite extensive, But this stands out to me.

Hi dokken,
I replaced u0.vector()[:] = u.vector() by u0.assign(u) in my code but this does not affect the solution.
cheers,
domi

If you move your assign until after the solve, you get a variation as far as I can see.
In the following code I also removed the expression for x with the SpatialCoordinate-function.

from fenics import (NonlinearProblem, NewtonSolver, IntervalMesh, FiniteElement,
MixedElement, assemble, FunctionSpace, TestFunctions, Function,
interpolate, Expression, split, inner, grad, pi, dx, DirichletBC,
Constant, exp, ln, Dx, sqrt, ds, derivative, PETScKrylovSolver,
PETScFactory, near, PETScOptions, assign, File, plot, SpatialCoordinate)
import numpy as np
from scipy.constants import (epsilon_0, elementary_charge, m_p, m_e, k)
import matplotlib.pyplot as plt

class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)

def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)

class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "ilu")
self.linear_solver().set_from_options()

# Length 1D simulation domain (cylinrical sym. assumed)
R_m = 0.1

# Create mesh and define function space
nx = int(1e3) # cells
mesh = IntervalMesh(nx, 0, R_m)
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(mesh, ME_elem)

# Some parameters
dt = 1e-11
nm = 1e21 # Density neutrals 1/m**3

def boundary_right(x, on_boundary):
tol = 1E-14
return on_boundary and near(x, R_m, tol)

Phi_test, ni_test, ne_test, Te_test = TestFunctions(ME)

u  = Function(ME)

# Initial values for Phi, ni, ne, Te
u  = interpolate(Constant((0, 1e15, 1e15, 1)), ME)
u0 = interpolate(Constant((0, 1e15, 1e15, 1)), ME)

Phi, ni, ne, Te = split(u)
Phi0, ni0, ne0, Te0 = split(u0)

x = SpatialCoordinate(mesh)

# Variational formulation of Poisson equation incl. b. c. --> Phi
f = -elementary_charge * ( ni - ne) / epsilon_0
bc_Poisson_right = DirichletBC(ME.sub(0), Constant(0), boundary_right)

# Variational formulation of continuity equation for ions incl. b. c. --> ni
X_Ioniz_H2_ndiss = 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)
R = X_Ioniz_H2_ndiss*nm*ne*1e5
Ti_K = 300
X2H2 = 1e-15
mui = elementary_charge/(2*m_p*(1/2*X2H2*nm))
ui = mui*(-Dx(Phi,0) - k/elementary_charge*Ti_K*Dx(ni,0)/ni)
F += ( ni*ni_test - dt*ni*ui*Dx(ni_test,0) - (ni0 + dt*R)*ni_test )*2*pi*x*dx
# Non-homogeneous Neumann BC at R_m
F += dt*ni_test*ni*ui*2*pi*x*ds

# Variational formulation of continuity equation for electrons incl. b. c. --> ne
X_eH2 = 1e-6*exp(-16.0325527377+0.412023668827*ln(Te)**1+-0.193551110298*ln(Te)**2+-0.0247678479005*ln(Te)**3+0.0105040269085*ln(Te)**4+-0.000152726360717*ln(Te)**5+-0.000538801463715*ln(Te)**6+9.70545683084e-05*ln(Te)**7+-5.10345827257e-06*ln(Te)**8)
nueleH2 = X_eH2*nm
mue = elementary_charge/(m_e*nueleH2)
ue = mue*(Dx(Phi,0) - Te*Dx(ne,0)/ne)
F += ( ne*ne_test - dt*ne*ue*Dx(ne_test,0) - (ne0 + dt*R)*ne_test )*2*pi*x*dx
# Thermal flux at the boundary
vthe = sqrt(8*elementary_charge*Te/(pi*m_e))
F += dt*ne_test*ne*0.25*vthe*2*pi*x*ds

# Variational formulation of power balance for electrons
kappae = 2.5*mue*ne*Te
qer = -kappae*elementary_charge*Dx(Te,0)
Qer = 2.5*elementary_charge*Te*ne*ue + qer
Pext = 1e5*exp(-1e6*(x-0.85*R_m)**4)
P  = -Pext
F += ( 1.5*ne*Te*Te_test - dt*Qer*Dx(Te_test,0) - (1.5*ne0*Te0 + dt*P)*Te_test )*2*pi*x*dx
F += dt*Te_test*2*elementary_charge*Te*ne*0.25*vthe*2*pi*x*ds

# Output file
#file = File("solution/output.pvd")

# Create VTK files for visualization output
vtkfile_Phi = File('solution/Phi.pvd')
vtkfile_ni = File('solution/ni.pvd')
vtkfile_ne = File('solution/ne.pvd')
vtkfile_Te = File('solution/Te.pvd')

# Time-stepping
t = 0.0
T = 1e-3
while (t < T):
t += dt
dt *= 1.2
J = derivative(F, u)
problem = Problem(J, F, [bc_Poisson_right])
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
u0.assign(u)
_Phi, _ni, _ne, _Te = u.split()
vtkfile_Phi << (_Phi, t)
vtkfile_ni << (_ni, t)
vtkfile_ne << (_ne, t)
vtkfile_Te << (_Te, t)
#file << (u.split(), t)


I also cannot see a significant change in the solution when doing the assign after the solve. Also the different definition of x does no change, i.e. n_e, n_i and T_e are clamped in volume to there initial values and I see no increase in time (as is expected from the comsol solution).
BTW: I messed up the rate, i.e. in order to obtain the steady state values stated above, the rate should be without the 1e5, i.e.

R = X_Ioniz_H2_ndiss*nm*ne


Sorry for that. Strangely, the solution does not seem to be affected by the actual rate R. If I set it e.g. to 0, the profiles do not change and I have to multiply with a huge number (e.g. 1e20), in order to get the n_e just a bit over the initial 1e15. Is it somehow possible that the rate R is not applied correctly?

I would say your problem is quite poorly scaled, as you solution is of magnitude 10^{15}, while you have parameters of order 10^{-6}.
Anyhow, your problem has so many variables as I cannot really do much more.
You need to start from making minimal examples, maybe just solving one of the equations above and check that everything behaves as expected.

I resolved the issue and have now a working model that produces the same stationary solution as comsol. The error I made was that dt was not defined as an expression and therefore not updated correctly in the variational form. (The problem is now defined outside the time loop.) Below is the code that works:

from petsc4py import PETSc
from fenics import (NonlinearProblem, NewtonSolver, IntervalMesh, FiniteElement,
MixedElement, assemble, FunctionSpace, TestFunctions, Function,
interpolate, Expression, split, inner, grad, pi, dx, DirichletBC,
Constant, exp, ln, Dx, sqrt, ds, derivative, PETScKrylovSolver,
PETScFactory, near, PETScOptions, assign, File, plot, SpatialCoordinate)
import numpy as np
import sys
from scipy.constants import (epsilon_0, elementary_charge, m_p, m_e, k)
import matplotlib.pyplot as plt

class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)

def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)

class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

#PETScOptions.set("ksp_rtol", 1.0e-3)
#PETScOptions.set("ksp_atol", 1.0e-3)
PETScOptions.set("ksp_type", "gmres")
#PETScOptions.set("ksp_monitor")
#PETScOptions.set("ksp_converged_reason")
PETScOptions.set("pc_type", "ilu")

self.linear_solver().set_from_options()
#print(PETSc.Options().getAll())

# Length 1D simulation domain (cylinrical sym. assumed)
R_m = 0.1

# Plot Rate coefficient
#Te_plt = np.linspace(0.1, 25)
#X_plt = 1e-6*np.exp(-35.3071744514+16.9033721249*np.log(Te_plt)**1+-7.81316801032*np.log(Te_plt)**2+2.47737696605*np.log(Te_plt)**3+-0.62396194653*np.log(Te_plt)**4+0.145044251711*np.log(Te_plt)**5+-0.0286523207942*np.log(Te_plt)**6+0.00352120706732*np.log(Te_plt)**7+-0.000181868246303*np.log(Te_plt)**8)
#plt.semilogy(Te_plt, X_plt)
#plt.xlim(0, 25)
#plt.ylim(1e-18, 1e-13)
#sys.exit(0)

# Plot power deposition profile
#r_ = np.linspace(0, R_m)
#Pext_ = 1e7*np.exp(-1e6*(r_-0.85*R_m)**4)
#plt.plot(r_, Pext_)
#sys.exit(0)

# Create mesh and define function space
nx = int(1e3) # cells
mesh = IntervalMesh(nx, 0, R_m)
#print(mesh.coordinates())
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
ME_elem = MixedElement([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
ME = FunctionSpace(mesh, ME_elem)

# Some parameters
nm = 1e21 # Density neutrals 1/m**3

def boundary_right(x, on_boundary):
tol = 1E-14
return on_boundary and near(x, R_m, tol)

Phi_test, ni_test, ne_test, Te_test = TestFunctions(ME)

u  = Function(ME)
u0 = Function(ME)
# Initial values for Phi, ni, ne, Te
u  = interpolate(Constant((0, 1e15, 1e15, 5)), ME)
u0 = interpolate(Constant((0, 1e15, 1e15, 5)), ME)

Phi, ni, ne, Te = split(u)
Phi0, ni0, ne0, Te0 = split(u0)

x = SpatialCoordinate(mesh)

# Time-stepping
t = 0.0
T = 1e-5
dt = Expression('dtvalue', dtvalue = 1e-10, degree=1)

# Variational formulation of Poisson equation incl. b. c. --> Phi
f = -elementary_charge * ( ni - ne) / epsilon_0
bc_Poisson_right = DirichletBC(ME.sub(0), Constant(0), boundary_right)

# Variational formulation of continuity equation for ions incl. b. c. --> ni
X_Ioniz_H2_ndiss = 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)
R = X_Ioniz_H2_ndiss*nm*ne
Ti_K = 300
X2H2 = 1e-15
mui = elementary_charge/(2*m_p*(0.5*X2H2*nm))
ui = mui*(-Dx(Phi,0) - k/elementary_charge*Ti_K*Dx(ni,0)/ni)
F += ( (ni-ni0)*ni_test - dt*ni*ui*Dx(ni_test,0) - dt*R*ni_test )*2*pi*x*dx
# Non-homogeneous Neumann BC at R_m
F += dt*ni_test*ni*ui*2*pi*x*ds

# Variational formulation of continuity equation for electrons incl. b. c. --> ne
X_eH2 = 1e-6*exp(-16.0325527377+0.412023668827*ln(Te)**1+-0.193551110298*ln(Te)**2+-0.0247678479005*ln(Te)**3+0.0105040269085*ln(Te)**4+-0.000152726360717*ln(Te)**5+-0.000538801463715*ln(Te)**6+9.70545683084e-05*ln(Te)**7+-5.10345827257e-06*ln(Te)**8)
nueleH2 = X_eH2*nm
mue = elementary_charge/(m_e*nueleH2)
ue = mue*(Dx(Phi,0) - Te*Dx(ne,0)/ne) #thermal diffusion not considered!
F += ( (ne-ne0)*ne_test - dt*ne*ue*Dx(ne_test,0) - dt*R*ne_test )*2*pi*x*dx
# Thermal flux at the boundary
vthe = sqrt(8*elementary_charge*Te/(pi*m_e))
F += dt*ne_test*ne*0.25*vthe*2*pi*x*ds

# Variational formulation of power balance for electrons
kappae = 2.5*mue*ne*Te
qer = -kappae*elementary_charge*Dx(Te,0)
Qer = 2.5*elementary_charge*Te*ne*ue + qer
Pext = 1e5*exp(-1e6*(x-0.85*R_m)**4)
P  = Pext
F += ( 1.5*elementary_charge*(ne*(Te-Te0)+Te*(ne-ne0))*Te_test - dt*Qer*Dx(Te_test,0) - dt*P*Te_test )*2*pi*x*dx
F += dt*Te_test*2.5*elementary_charge*Te*ne*0.25*vthe*2*pi*x*ds

J = derivative(F, u)
problem = Problem(J, F, [bc_Poisson_right])
custom_solver = CustomSolver()
custom_solver.parameters["relative_tolerance"] = 1e-3
custom_solver.parameters["maximum_iterations"] = 20

# Create VTK files for visualization output
vtkfile_Phi = File('solution_custom_solver/Phi.pvd')
vtkfile_ni = File('solution_custom_solver/ni.pvd')
vtkfile_ne = File('solution_custom_solver/ne.pvd')
vtkfile_Te = File('solution_custom_solver/Te.pvd')

while (t < T):
print('Time:' + str(t))
custom_solver.solve(problem, u.vector())
u0.assign(u)
_Phi, _ni, _ne, _Te = u.split()
vtkfile_Phi << (_Phi, t)
vtkfile_ni << (_ni, t)
vtkfile_ne << (_ne, t)
vtkfile_Te << (_Te, t)
t += dt.dtvalue
dt.dtvalue *= 1.05


There is one issue remaining: when I try to go to T = 1e-4 instead of 1e-5, then the Newton solver does not converge at 8.7e-5 (max. iterations exception), even though the system seems to be already in a steady state. I guess this issue comes from a too large time step (I increase by 5% in every step). I wonder whether there is a more intelligent way to adapt the time step or at least catch the exception and reduce the time step in this case. Any suggestion or hint for an example, where this has been done is therefore highly welcome. Thank you.

You can always add a try/except on the solve, see for instance https://www.w3schools.com/python/python_try_except.asp