### Convergence issues in drift-diffusion equations

533

views

1

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

The residual was calculated as such

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 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.

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.

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.

> 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

### 1 Answer

0

File attached: ddtest_coupled_1D_mwe.py (9.58 KB)

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

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.

As a suggestion I would try to keep computing using quadrature, you claim converges, but a posteriori check the residuals using uflacs.