### Convergence issues in drift-diffusion equations

533
views
1
9 months ago by
Hello all,

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:

$\begin{gather} 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)) \end{gather}$

$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)

Thanks,
Chaffra
Community: FEniCS Project
I would believe more uflacs when dealing with numbers of type 1e-20, see https://bitbucket.org/fenics-project/ffc/issues/127 and also https://bitbucket.org/fenics-project/ffc/issues/139.

As a suggestion I would try to keep computing using quadrature, you claim converges, but a posteriori check the residuals using uflacs.
written 9 months ago by Jan Blechta
I seem to be seeing the converse effect: quadrature converges and uflacs does not. It would be nice to find out what's being dropped by quadrature in this case but it seems to give me a physical answer. As anyone tried running the example to see if non convergence is reproducible?

If you change the variable n_i from 1e-3 to > 1e10 both uflacs and quadrature converge but not for small values of n_i like <1e-3.
written 9 months ago by Chaffra Affouda
I am still having problem with this. Can anyone else try to reproduce? My worry is that when the quadrature back-end is deprecated some of my examples will stop working. Uflacs tend to diverge or stop converging with a high residual. For the quadrature backend the residual goes to zero as expected.
written 9 months ago by Chaffra Affouda
We can't be looking for a bug in uflacs based on 300 lines of code application.

I suggested that you try to compute the residuals which you claim tend to zero by uflacs. It is possible that they don't go to zero because it is known that quadrature has troubles with small numbers. See the issues in my comment above.
written 9 months ago by Jan Blechta
In fact this

> If you change the variable n_i from 1e-3 to > 1e10 both uflacs and quadrature converge but not for small values of n_i like <1e-3.

is consistent with the issues above and just supports a hypothesis that uflacs is right and quadrature is wrong. The fact that you code does not converge with uflacs might not be a bug in uflacs but a property of your algorithm.
written 9 months ago by Jan Blechta
What do you mean by compute the residual? The residual is computed by the NewtonSolver. Am I missing something? Are you saying it's possible my residual is not zero based on the forms I have written? In that case the quadrature representation would be giving me some rounding errors that would somehow lead to convergence?
written 9 months ago by Chaffra Affouda
Yes. Assemble the residuals which you claim tend to zero by uflacs.
written 9 months ago by Jan Blechta
Sorry for my confusion. So I should assemble the residual with the solution given by quadrature (since that's the one that seems to converge) using the uflacs representation to see if my residual is indeed close zero.
written 9 months ago by Chaffra Affouda
Exactly................
written 9 months ago by Jan Blechta

I did what you suggested and both quadrature and uflacs give the same residual. I updated the example in the question for reference. Something is off either in my script or uflacs. The residual for the 3 coupled non-linear Poisson's equations is ~4e-2 for the bias used, which I think is ok. But the problem remains, uflacs does not converge. Maybe it's the way I am writing the forms, I don't know...

The residual was calculated as such

F_residual = abs(div(eps*eps0*grad(u*VT)) + rho*q) #7.7e-7
F_residual += abs(-div(mun*n*grad(un*VT)) - R) #1.7e-2
F_residual += abs(-div(mup*p*grad(up*VT)) + R)​ #2.58E-2
written 9 months ago by Chaffra Affouda
Convergence of Newton solver is established on weak residuals, i.e., assemble(F) module boundary conditions. I am sure these strong residuals are relevant - they feature second derivative of function which does not have integrable second derivatives.
written 9 months ago by Jan Blechta
do you mean to say, "assemble(F) modulo boundary conditions?"  If so, do you mean that in this context the boundary conditions do not matter?  I would feel smart if I have interpreted you correctly.
written 9 months ago by pf4d
No. Boundary conditions matter. See implementation of NewtonSolver for computation of residuals.
written 9 months ago by Jan Blechta
Right, I knew this, I was only trying to decipher what you mean by "assemble(F) module boundary conditions"
written 9 months ago by pf4d

0
9 months ago by
File attached: ddtest_coupled_1D_mwe.py (9.58 KB)

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?