File naulin_laplace.cxx#

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

Scheme suggested by Volker Naulin: solve Delp2(phi[i+1]) + 1/DC(C1*D)*Grad_perp(DC(C2))*Grad_perp(phi[i+1]) + DC(A/D)*phi[i+1] = rhs(phi[i]) + 1/DC(C1*D)*Grad_perp(DC(C2))*Grad_perp(phi[i]) + DC(A/D)*phi[i] 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 C1*D, C2 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.

CHANGELOG

Copyright 2018 B.D.Dudson, M. Loiten, J. Omotani

Contact: Ben Dudson, benjamin.dudson@york.ac.uk

This file is part of BOUT++.

BOUT++ is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

BOUT++ is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with BOUT++. If not, see http://www.gnu.org/licenses/.

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{split}\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}\end{split}\]

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:

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

  2. 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\)

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

  4. Calculate the infinity norms of the error

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

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