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.
Copyright 2018 B.D.Dudson, M. Loiten, J. Omotani
Contact: Ben Dudson, email@example.com
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:
In fact we allow for the slightly more general form
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 fromwhere phiCur is phi of the current iteration [and DC(f) is the constant-in-z component of f]
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]
Invert \(phi\) to find the voricity usingwhere phiNext is the newly obtained \(phi\)
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.
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) )
EAbsLInf > atol
ERelLInf > rtol
EAbsLInf > EAbsLInf(previous step)
If yesRestart iteration
underrelax_factor *= 0.9
- increase curCount and start from step 1
phiCur = phiNext
If number of iteration is above maxit, throw exception
Stop: Function returns phiNext
Stop: Function returns phiNext