Convergence issues in drift-diffusion equations
9 months ago by
I use quadrature representation with FFC to solve non-linear
(drift-diffusion) problems in DOLFIN. With the new UFLACS representation
I have noticed that though my program runs faster it sometimes does not
converge in situations where the quadrature representation converges.
I am attaching a simple example describing the situation where a problem would converge with
quadrature but not with uflacs. My suspicion is that the float precision
in uflacs or some other parameters need to be tweaked. I would appreciate any suggestion.
I am trying to solve the following coupled system of non-linear Poisson's equations in dolfin:
V_t\nabla \cdot E - \rho = 0\\
V_t \nabla \cdot J_n - R = 0 \\
V_t \nabla \cdot J_p + R = 0 \\
V_t = kT/q \\
n = n_i \exp(u-u_n) \\
p = n_i \exp(u_p-u) \\
\rho = q(p+d-n) \\
E = -\epsilon \nabla u \\
J_n = -\mu_n n \nabla u_n \\
J_p = -\mu_p p \nabla u_p \\
R = B(np-n_i^2) + (np-n_i^2)/(\tau_n(p+n_i) + \tau_p(n+n_i))
\(u, u_n, u_p\) are the unknowns and all the non-defined parameters are physical constants.
File attached: ddtest_coupled_1D.py (8.56 KB)
Community: FEniCS Project
9 months ago by
I think I am making progress. To compute the residual correctly I am now using
Q = FunctionSpace(mesh,'CG',4) parameters['form_compiler']['quadrature_degree'] 8
This should be correct because the source terms have the form \(np\) so if \(n,p\) live in CG2 then \(np\) is in CG4 (not sure if this is the correct reasoning).
I no longer see convergence problem with uflacs. It was also important to scale the 3 equations so that they have the same units for the NewtonSolver residual vector. However one problem remains. The residual is noisy when the the solution is almost flat. I am attaching the new script to see if it makes sense. The noise always occurs on the side where I am applying the non-zero boundary condition (right side here). Clearly the fact that \(J_n+J_p=constant\) is not strongly enforced. Can this noise be removed somehow or does anyone know a way to enforce current conservation in these equations?
Please login to add an answer/comment or follow this question.