File naulin_laplace.hxx#

Iterative solver to handle non-constant-in-z coefficients

Scheme suggested by Volker Naulin: solve

\[\begin{split}\begin{eqnarray} \nabla^2(\phi[i+1]) + 1/DC(C_1 D)\nabla_\perp(DC(C_2))\nabla_\perp(\phi[i+1]) + DC(A/D)\phi[i+1] \\ = rhs(\phi[i]) + 1/DC(C_1 D)\nabla_\perp(DC(C_2))\nabla_\perp(\phi[i]) + DC(A/D)\phi[i] \end{eqnarray}\end{split}\]

using standard FFT-based solver, iterating to include other terms by evaluating them on rhs using \(\phi\) from previous iteration. DC part (i.e. Field2D part) of \(C_1 D\), \(C_2\) and \(A,D\) is kept in the FFT inversion to improve convergence by including as much as possible in the direct solve and so that all Neumann boundary conditions can be used at least when \(DC(A/D)!=0\)

Explanation of the procedure

A way to invert the equation \(\Omega^D = \nabla\cdot(n\nabla_\perp \phi)\) invented by Naulin, V. In an orthogonal system, we have that:

\[\begin{split}\begin{eqnarray} \Omega^D &=& \nabla\cdot(n\nabla_\perp \phi)\ \ &=& n \nabla_\perp^2 \phi + \nabla n\cdot\nabla_\perp \phi\\ &=& n \Omega + \nabla n\cdot\nabla_\perp \phi\\ &=& n \Omega + \nabla_\perp n\cdot\nabla_\perp \phi \end{eqnarray}\end{split}\]

Rearranging gives

\[\begin{eqnarray} \Omega &=& \frac{\Omega^D}{n} - \nabla_\perp \ln(n)\cdot\nabla_\perp \phi\ \ \nabla_\perp^2 \phi &=& \frac{\Omega^D}{n} - \nabla_\perp \ln(n)\cdot\nabla_\perp \phi \end{eqnarray}\]

In fact we allow for the slightly more general form

\[\begin{eqnarray} \nabla_\perp^2 \phi + <\frac{A}{D}>\phi &=& rhs/D - \frac{1}{D\,C1} \nabla_\perp C2\cdot\nabla_\perp \phi - (\frac{A}{D} - <\frac{A}{D}>)\phi \end{eqnarray}\]

The iteration can be under-relaxed to help it converge. Amount of under-relaxation is set by the parameter ‘underrelax_factor’. 0<underrelax_factor<=1, with underrelax_factor=1 corresponding to no under-relaxation. The amount of under-relaxation is temporarily increased if the iteration starts diverging, the starting value uof underrelax_factor can be set with the initial_underrelax_factor option.

The iteration now works as follows:

Get the vorticity from

vort = (vortD/n) - grad_perp(ln_n)*grad_perp(phiCur)
[Delp2(phiNext) + 1/DC(C2*D)*grad_perp(DC(C2))*grad_perp(phiNext) + DC(A/D)*phiNext
 = b(phiCur)
 = (rhs/D) - (1/C1/D*grad_perp(C2)*grad_perp(phiCur) - 1/DC(C2*D)*grad_perp(DC(C2))*grad_perp(phiCur)) - (A/D - DC(A/D))*phiCur]

where phiCur is phi of the current iteration [and DC(f) is the constant-in-z component of f]

Invert \(phi\) to find the voricity using

phiNext = invert_laplace_perp(vort)
[set Acoef of laplace_perp solver to DC(A/D)
 and C1coef of laplace_perp solver to DC(C1*D)
 and C2coef of laplace_perp solver to DC(C2)
 then phiNext = invert_laplace_perp(underrelax_factor*b(phiCur) - (1-underrelax_factor)*b(phiPrev))]
 where b(phiPrev) is the previous rhs value, which (up to rounding errors) is
 the same as the lhs of the direct solver applied to phiCur.

where phiNext is the newly obtained \(phi\)

Calculate the error at phi=phiNext

error3D = Delp2(phiNext) + 1/C1*grad_perp(C2)*grad_perp(phiNext) + A/D*phiNext - rhs/D
        = b(phiCur) - b(phiNext)
as b(phiCur) = Delp2(phiNext) + 1/DC(C2*D)*grad_perp(DC(C2))*grad_perp(phiNext) + DC(A/D)*phiNext
up to rounding errors

Calculate the infinity norms of the error

EAbsLInf = max(error3D)
ERelLInf = EAbsLInf/sqrt( max((rhs/D)^2) )

Check whether
```
EAbsLInf > atol
```
- If yes
  - Check whether
    ERelLInf > rtol
  - If yes
    - Check whether
      EAbsLInf > EAbsLInf(previous step)
      - If yes
        underrelax_factor *= 0.9
        
        Restart iteration
      - If no
        
        Set
        phiCur = phiNext
        
        increase curCount and start from step 1
        
        If number of iteration is above maxit, throw exception
  - If no
    - Stop: Function returns phiNext
- if no
  - Stop: Function returns phiNext

class LaplaceNaulin : public Laplacian #

Solves the 2D Laplacian equation.

Public Functions

LaplaceNaulin(Options *opt = NULL, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr, Solver *solver = nullptr)#

~LaplaceNaulin() = default#

inline virtual void setCoefA(const Field2D &val) override#: Set coefficients for inversion. Re-builds matrices if necessary.

inline virtual void setCoefA(const Field3D &val) override#

inline virtual void setCoefC(const Field2D &val) override#

inline virtual void setCoefC(const Field3D &val) override#

inline virtual void setCoefC1(const Field3D &val) override#

inline virtual void setCoefC1(const Field2D &val) override#

inline virtual void setCoefC2(const Field3D &val) override#

inline virtual void setCoefC2(const Field2D &val) override#

inline virtual void setCoefD(const Field3D &val) override#

inline virtual void setCoefD(const Field2D &val) override#

inline virtual void setCoefEx(const Field2D &val) override#

inline virtual void setCoefEz(const Field2D &val) override#

inline virtual bool uses3DCoefs() const override#: Does this solver use Field3D coefficients (true) or only their DC component (false)

inline virtual FieldPerp solve(const FieldPerp &b) override#

inline virtual FieldPerp solve(const FieldPerp &b, const FieldPerp &x0) override#

virtual Field3D solve(const Field3D &b, const Field3D &x0) override#

inline virtual Field3D solve(const Field3D &b) override#

inline virtual void setGlobalFlags(int f) override#

inline virtual void setInnerBoundaryFlags(int f) override#

inline virtual void setOuterBoundaryFlags(int f) override#

inline BoutReal getMeanIterations() const#

inline void resetMeanIterations()#

void outputVars(Options &output_options, const std::string &time_dimension) const override#

void setCoefA(const Field2D &val) = 0: Set coefficients for inversion. Re-builds matrices if necessary.

inline void setCoefA(const Field3D &val)

inline void setCoefA(BoutReal r)#

void setCoefC(const Field2D &val) = 0

inline void setCoefC(const Field3D &val)

inline void setCoefC(BoutReal r)#

inline void setCoefC1(const Field2D &val)

inline void setCoefC1(const Field3D &val)

inline void setCoefC1(BoutReal r)#

inline void setCoefC2(const Field2D &val)

inline void setCoefC2(const Field3D &val)

inline void setCoefC2(BoutReal r)#

void setCoefD(const Field2D &val) = 0

inline void setCoefD(const Field3D &val)

inline void setCoefD(BoutReal r)#

void setCoefEx(const Field2D &val) = 0

inline void setCoefEx(const Field3D &val)#

inline void setCoefEx(BoutReal r)#

void setCoefEz(const Field2D &val) = 0

inline void setCoefEz(const Field3D &val)#

inline void setCoefEz(BoutReal r)#

FieldPerp solve(const FieldPerp &b) = 0

Field3D solve(const Field3D &b)

Field2D solve(const Field2D &b)#

inline FieldPerp solve(const FieldPerp &b, const FieldPerp &x0)

Field3D solve(const Field3D &b, const Field3D &x0)

Field2D solve(const Field2D &b, const Field2D &x0)#

Private Functions

LaplaceNaulin(const LaplaceNaulin&)#

LaplaceNaulin &operator=(const LaplaceNaulin&)#

void copy_x_boundaries(Field3D &x, const Field3D &x0, Mesh *mesh)#: Copy the boundary guard cells from the input ‘initial guess’ x0 into x. These may be used to set non-zero-value boundary conditions

Private Members

Field3D Acoef#

Field3D C1coef#

Field3D C2coef#

Field3D Dcoef#

std::unique_ptr<Laplacian> delp2solver = {nullptr}#: Laplacian solver used to solve the equation with constant-in-z coefficients.

BoutReal rtol#: Solver tolerances.

BoutReal atol#

BoutReal rtol_accept#: Accept these tolerances if number of iterations exceeds maxits.

BoutReal atol_accept#

int maxits#: Maximum number of iterations.

BoutReal initial_underrelax_factor = {1.}#: Initial choice for under-relaxation factor, should be greater than 0 and less than or equal to 1. Value of 1 means no underrelaxation

BoutReal underrelax_threshold = {1.5}#: Factor by which the error may increase since the previous iteration before triggering a sub-cycle with decreasing underrelax_factor to try and prevent divergence.

BoutReal underrelax_decrease_factor = {0.9}#: Factor by which to decrease underrelax_factor at each stage of the sub-loop triggered if underrelax_threshold is crossed.

int underrelax_decrease_maxits = {10}#: Maximum number of iterations in the decreasing-underrelax_factor subcycle before trying to continue the main iteration loop.

BoutReal underrelax_recovery = {1.1}#: If underrelax_factor has been decreased, increase it again by this factor at the end of a successful iteration to try and speed up subsequent convergence.

BoutReal naulinsolver_mean_its#: Mean number of iterations taken by the solver.

BoutReal naulinsolver_mean_underrelax_counts = {0.}#: Mean number of times the underrelaxation factor is reduced.

int ncalls#: Counter for the number of times the solver has been called.