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[0], 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[0]', degree = 1)
# Variational formulation of Poisson equation incl. b. c. --> Phi
f = -elementary_charge * ( ni - ne) / epsilon_0
F = (-inner( grad(Phi), grad(Phi_test)) - Phi_test*f )*2*pi*x*dx
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
