TIFF: Gyrofluid Turbulence in Full-f and Full-k

A model and code (TIFF) for isothermal gyrofluid computation of quasi-two-dimensional interchange and drift wave turbulence in magnetized plasmas with arbitrary fluctuation amplitudes (full-f) and arbitrary polarization wavelengths (full-k) is introduced. The model reduces to the paradigmatic Hasegawa-Wakatani model in the limits of small turbulence amplitudes (delta-f), cold ions (without finite Larmor radius effects), and homogeneous magnetic field. Several solvers are compared for the generalized Poisson problem, that is intrinsic to the full-f gyrofluid (and gyrokinetic) polarization equation, and a novel implementation based on a dynamically corrected Fourier method is proposed. The code serves as a reference case for further development of three-dimensional full-f full-k models and solvers, and for fundamental exploration of large amplitude turbulence in the edge of magnetized plasmas.


I. INTRODUCTION: GYROFLUID MODELS AND POLARIZATION
Turbulence in magnetized plasmas, which is generally driven by ubiquitous gradients of density and temperature, is a subject of considerable interest and importance in fusion energy research with toroidal magnetic high-temperature plasma confinement experiments such as tokamaks and stellarators [1]. Models for instabilities and turbulence in magnetized plasmas are based on gyrokinetic, gyrofluid or fluid drift descriptions [2,3].
Gyrofluid theory as a fluid-like model for the low-frequency dynamics of gyrocenter densities in magnetized plasmas including finite Larmor radius (FLR) effects has been pioneered in the late 1980s [4]. In the early 1990s collisionless Landau closure has been applied to derive three-dimensional gyrofluid models [5] for computations in slab [6] and toroidal geometry [7].
The approach was further generalized to electromagnetic models [8,9]. Energetically consistent 6-moment electromagnetic gyrofluid equations with FLR effects have subsequently been systematically derived from the corresponding gyrokinetic theory [10,11].
Gyrofluid theory has its place in the hierarchy of magnetized plasma models somewhere in between gyrokinetic and drift-reduced fluid theory [9][10][11][12][13][14][15]. The quality of a model of course not only depends on the rank within some straight hierarchy, but also on multiple secondary modelling assumptions that are introduced in any practical application and numerical simulation. Gyrofluids have the advantage over fluid theories to consistently model some kinetic effects such as finite Larmor radii or Landau damping. They do not depend on (Braginskii type) collisional closures, but collisional effects can be amended. Compared to gyrokinetic modelling of a 5-dimensional distribution function, the 3-dimensional gyrofluid moment models can be computationally much more efficient.
The development of full-f gyrokinetic codes started a bit more than a decade ago (see, for example, refs. [15,25,[28][29][30][31][32][33][34][35][36][37]), and often involves long-term development in large research groups or collaborations. All these (perpetually evolving) codes employ different levels of secondary modelling assumptions and approximations, such as for geometry (full torus vs. flux tube; circular vs. X-point; core vs. edge/SOL), for polarization (linear vs. nonlinear), or for collision operators, electromagnetic effects, and (SOL and core) boundary conditions. Full-f gyrofluid models are directly related to gyrokinetics, as they are usually obtained by taking appropriate fluid moments of the full-f gyrokinetic equations [18]. Therefore the difference between delta-f and full-f gyrofluid models involves the same choice of assumptions on the level of the gyrokinetic Lagrangian action, which in principle requires a consistent decision between either small perturbations but applicability for arbitary wave lengths, or arbitrary perturbation and flow amplitudes but restriction to long wave lengths [14,16,17].
For simulations of full-f gyrofluid models, over the last years two separate code implementations were developed alongside, which are both based on the same or similar model sets [18,38], but use largely different numerical approaches: The modular open source code suite FELTOR ("Full-F ELectromagnetic code in TORoidal geometry") [39,40] has been primarily developed and maintained by Wiesenberger and Held et al. [41][42][43][44]; and the code family TOEFL is being primarily developed by Kendl and includes the 2d code branch "TIFF", which is reported herein. The acronym TOEFL denotes "TOkamak Edge FLuid", and TIFF is "TOEFL In Full-f and Full-k". In its first (drift-Alfvén fluid in C/C++) implementation, TOEFL was largely tantamount to DALF [45] but in another language.
The 3d delta-f gyrofluid version T3P in the TOEFL set had been designed to be comparable with the GEM3 model and code by Scott [46] and Ribeiro [47], and was applied in refs. [48][49][50][51], and as 2d delta-f gyrofluid model reduction in refs. [50,52,53]. Full-f low-k 2d and 3d gyrofluid versions of TOEFL in the previous long wavelength approximation implementation, with otherwise similar numerics as employed here for TIFF, have already been applied in refs. [54][55][56]. These previous delta-f and full-f low-k simulations can serve as test cases for the implementation of the present full-f full-k model in its respective parameter limits. The dual FELTOR vs. TOEFL code development strategy allows cross-verification, but they for example also have optimized applicability for different geometries. Whereas FELTOR can be applied to full global 3d torus geometry, TOEFL is presently designed for locally field-aligned 3d flux-tube type toroidal simulation geometry.
What ever gyrofluid moment sets (such as thermal or isothermal, 2d or 3d) are treated in the codes, all gyrocenter densities are evolved and coupled to obtain the electric potential ϕ(x, t) via the polarization equation, which is isomorphic to its full-f gyrokinetic equivalent, and is nothing but the gyrocenter density (N s (x, t)) formulation of quasi-neutrality, summing up all species (index s) charge and polarization densities.
Until recently, both full-f gyrofluid code families have treated the usual, consistent longwavelength form of the full-f polarization equation [16][17][18]: This can be re-cast into a generalized 2d Poisson equation: where ε(x) ≡ s m s N s (x)/B 2 , and σ(x) ≡ − s q s G 1s N s (x). Here m s is the mass and q s = Z s e the charge of plasma species s with gyrocenter density N s . The plasma species include electrons and at least one but often several ion species. In this present work only one main ion species (index s = i) in addition to the electrons (index s = e) is considered.
In the case of a pure e-i plasma, the polarization of the electrons can be neglected because m e ≪ m i , so that ε(x) ≈ m i N i (x)/B 2 . The gyrocenter densities N s in σ are affected by gyroaveraging, which is denoted by the operator G 1s . This can be expressed as G 1 (b s ) = G 1/2 0 (b s ) in wavenumber space with the gyro-screening operator G 0s ≡ I 0 (b s )e −bs for b s = (ρ s k ⊥ ) 2 [14].
This form (derived from velocity integration of the gyrokinetic pendant including Bessel functions) makes use of the modified Bessel function of the first kind I 0 . The Larmor radius ρ s = √ T s m s /(eB) of the particle species (s), with temperature T s and mass m s , normalizes the k x and k y wavenumber components in the 2d plane locally perpendicular to the magnetic field B = Be z . The gyro-averaging operators can be efficiently approximated [14] by their The full dynamical nonlinearity ∇ · ε(x, y, t) ∇ ⊥ ϕ in the long-wavelength polarization equation is here retained. Most full-f gyrokinetic implementations so far approximate this term by using either only radial variations of the static background ε 0 (x), or completely linearise it to ∼ ε 0 ∇ 2 ⊥ ϕ and so neglect spatial and temporal variations in the polarization. However, the general delta-f form of the polarization equation for small perturbationsÑ s on a reference background density N 0 , differs from this linearised long-wave length polarization by an additional FLR contribution proportional to the Larmor radius ρ s , and only agrees in lowest order Taylor expansion in b = (ρ s k ⊥ ) 2 [14]: In any case consistency throughout the equations has to be ensured, for example in full-f models by keeping E-cross-B energy terms in the generalised potential [14,18]. Computation of delta-f gyrofluid turbulence with the exact model in comparison to the delta-f long-wavelength approximation such as in eq. (5) show that the differences can be rather pronounced. This has motivated the development of a consistent arbitrary wavelength full-f gyrofluid polarization model [58]. A first implementation of this "full-f full-k" model in the code FELTOR and application to simulation of interchange driven "blob" perturbations in a magnetized plasma is presented by Held and Wiesenberger in ref. [59]. The results therein clearly show the relevance of arbitrary wave length polarization for interchange drift modes.
In the present work, an independent code implementation (TIFF) of the "full-f full-k" model and application to 2d drift instabilities and turbulence is given. In the respective interchange blob mode limit, the results are cross-verified with the recent FELTOR code results. Different solvers for the underlying generalized Poisson problem are compared. An efficient solver based on a dynamically corrected Fourier method is proposed and tested.

II. FULL-F FULL-K 2D GYROFLUID TURBULENCE MODEL
The arbitrary wave length full-f model derived in ref. [58] and first applied by Held and Wiesenberger in ref. [59], is in the following implemented in a form which also includes the full-f formulation of the gyrofluid generalization [43,44] of the quasi-2d (modified) Hasegawa-Wakatani (HW) drift wave turbulence model [60,61].

A. Gyrocenter density equations
The set of isothermal full-f full-k gyrofluid equations is based on dynamical evolution equations of the gyrocenter densities for each species s in the general form of a continuity equation: The velocities U s = U E + U B + U ∥ include the gyro-center E-cross-B drifts , and parallel velocities U ∥ . In contrast to the corresponding fluid continuity equation for particle densities, the gyrofluid formulation for gyrocenter densities does not contain polarization drifts, whose effects are covered by the relation of the gyrocenter densities within the polarization equation for the electric potential ϕ. The resulting set of equations can be seen as a variation of the vorticity-streamfunction formulation for a 2D Euler fluid model. The electric potential here takes the role of a streamfunction for the advecting E-cross-B velocity.
The gyrofluid potentials ϕ s = Γ 1s ϕ+Ψ s in the E-cross-B drift U E include the gyro-average part, and the consistent full-f full-k form [58] of the polarization contribution through the E-cross-B energy as Ψ s = (m/2qB 2 )|∇ ⊥ √ Γ 0 ϕ| 2 . For electrons ϕ e ≈ ϕ can be used because of the small mass ratio and associated small Larmor radius in comparison to ions. The potential ϕ is retrieved from solution of the polarization equation.

B. Polarization equation
The consistent arbitrary wavelength full-f polarisation equation has been derived in ref. [58] and is used here in isothermal (constant gyroradius) form as in ref. [59]: where the polarization densities are given as Here, in the (isothermal) constant gyroradius approximation the Γ 1 and √ Γ 0 operators are self-adjoint, and may for example be evaluated efficiently in k space. The general form for ρ s = ρ s (x) is given in ref. [58].
For the present implementation in the code TIFF only one ion species (e.g. Deuterium) is treated, and the electron polarization is neglected. The full-f full-k polarization equation then can be re-written as: where σ(x) ≡ − s Z s Γ 1s N s (x), with Γ 0 (here for ions only) and Γ 1s given in second order accurate Padé approximation [58,59], as given above, and ε i (x) = m i N i (x)/B 2 . In the constant gyroradius (isothermal) case the √ Γ 0 operators commute with the ∇ operators, so that also the arbitrary wave length polarization eq. (9) can again be re-cast into the usual form of a generalised 2d Poisson equation as ∇ · ε i ∇ ⊥ ϕ G = σ G . This allows to re-use common solvers (see Appendix) for this type of problem to obtain ϕ G and thus ϕ from known ε i and σ. For this contribution is here referred to as "curvature" although in a strict sense a quasi-2d model with straight magnetic field lines B = B(x)e z only has a gradient-B effect; however, in a toroidal system both curvature and gradient-B contributions occur in combination, and can be treated similarly within joint 3d expressions in (gyro-) fluid models if T ∥ ≈ T ⊥ .) The 2d advective drift operators (e z × ∇ϕ s ) · ∇f ≡ [ϕ s , f ] are here expressed in Poisson . The 2d gyrofluid potential field thus has the meaning of a stream function for the turbulent E-cross-B flows.
The divergence of the gradient-B flux gives the diamagnetic curvature term ∇·(N s U B ) = (N s T s /q s B)[ln B, N s ]. In contrast to delta-f models, these terms are not linearized and the full gyrocenter densities are in this form retained as multipliers to the Poisson brackets.

D. Parallel closure
In the present quasi-2d model the parallel velocity U ∥ contribution can be approximated by means of the Hasegawa-Wakatani closure [60]. From the full-f electron parallel momentum equation [18] in the quasi-stationary limit a relation in the form of a generalized Ohm's law is obtained as eη ∥ J ∥ = T e ∇ ∥ ln N e − e∇ ∥ ϕ. With J ∥ ≈ −eN e U ∥e , by assumption of a Spitzer resistivity η ∥ = 0.51m e ν e /(n e e 2 ), and applying that the electron gyrocenter density N e ≈ n e can be approximated well by the electron particle density n e , an expression for the term ∇ · (N e U ∥e ) ≡ −Λ ce in eq. (5) can be obtained [43]. The ion velocity contribution in the parallel response can be neglected because of the high ion inertia compared to electrons.
For this purpose a full-f non-adiabatic coupling parameter α ≡ T e k 2 ∥ /(η ∥ e 2 n 0 ω 0 ) = n e T e k 2 ∥ /(0.51m e ν e n 0 ω 0 ) can be defined for a selected parallel wavenumber k ∥ with ω 0 = eB/m i . The electron collision frequency ν e is in principle proportional to n e , inversely to T 3/2 e , and to the (in general density and temperature dependent) Coulomb logarithm. The usual approximation of a constant Coulomb logarithm is applied here, and in the present isothermal model only the density dependence ν e (n e ) ∼ n e ∼ N e needs to be discussed [43]: in the classical delta-f fluid HW model the collision frequency is assumed constant, so that α is a free constant parameter, whereas for the present full-f model the dependence is kept as α = N e α 0 . The final non-adiabatic "ordinary" HW drive term for electrons in the full-f form [43] is Λ c = αn 0 ω 0 [(eϕ/T e ) − ln(N/⟨N ⟩)]. The angled brackets denote a zonal average, which in the present 2d geometry amounts to averaging in the y direction.

E. Normalization to dimensionless form
The preceding evaluation of perpendicular and parallel fluxes gives the (still dimensional) density equations (5) alternatively as where the nonadiabatic coupling term Λ cs only is contributing for electrons (Λ ci ≡ 0).
Time t in the partial time derivative is normalized with respect to L ⊥ /c 0 , where c 0 = T e /m i is the thermal speed, and L ⊥ is a "typical" perpendicular length scale. For local pressure gradient driven systems this is usually set equal to the background gradient length, here L ⊥ ≡ L n = |∂ x ln n 0 (x)| −1 , so that temporal normalization relates to the diamagnetic drift frequency ω * = c 0 /L n . For global gradient driven systems often the minor torus radius L ⊥ ≡ a is rather used. In case of model systems with absent background gradient, such as for interchange driven "blob" setups, L ⊥ ≡ ρ 0 ≡ √ T e m i /(eB) is chosen as the drift scale, and time normalization is related to the ion gyration frequency ω 0 = c 0 /ρ 0 . These choices can be set by specifying δ ≡ ρ 0 /L ⊥ as a free input parameter.
The perpendicular spatial derivatives in x and y are always normalized with respect to ρ 0 , so that∂ t ≡ (L ⊥ /c 0 )∂ t , and {·, ·} ≡ ρ 2 0 [·, ·]. Densities and the magnetic field are normalized to reference quantitiesN = N/N 0 andB = B/B 0 , respectively, and the electric potential aŝ ϕ ≡ eϕ/T e . Temperature ratios are defined as τ s ≡ T s /(T e Z s ).
Dividing eq. (10) by N s and multiplying with L ⊥ /c 0 gives: The normalized dissipative coupling termΛ as a free parameter, appears only in the electron equation, and Λ ci ≡ 0. The "modified HW" model for toroidally more concordant zonal flow treatment requires to useΛ ce =α[(φ − ⟨φ⟩) − (ln N e − ⟨ln N e ⟩)] instead [43]. In these forms both are directly consistent with their respective common delta-f limits (ln N e →Ñ e /N 0 ).
In the present quasi-2d setup it is assumed that the magnetic field has locally only a taken as a constant parameter. This is also consistent with the assumption of a constant gyroradius, throughout the computational domain which is supposed to have a radial extension L x ≪ R much smaller than the major torus radius R. For global simulations across the whole torus cross section this condition would need to be relaxed. The x direction is here chosen to be local radially outwards on the outboard midplane side of a torus (compare Fig. 1), where κ > 0, andB(x) ≈ 1 − κx ≈ 1 is weakly decreasing with x. This givesΛ Bs ≡ K(φ s ) + τ s K(lnN s ). The curvature strength can be evaluated as κ ≈ 2ρ 0 /R. The temperature ratio is always τ e = −1 for electrons, and a free parameter for ions in the order of τ i ∼ +1. The ion temperature ratio τ i thus controls the FLR effects, and in addition here also the ion diamagnetic contribution to gradient-B (interchange) drive for inhomogeneous magnetic fields (κ ̸ = 0).
The formulation in terms of lnN as the dynamical variable for the time evolution of densities ensures thatN is always positive definite, which is a requirement for the solution

F. Normalized polarization
In a quasi-neutral magnetized plasma under fluid drift ordering, the divergences of the above perpendicular and parallel fluxes are balanced by the divergence of the polarization drift. In gyrofluid and gyrokinetic models this is taken into account by enforcing quasineutrality s q s n s ≡ 0, and replacing the particle densities n s herein in terms of their respective gyrocenter and polarization densities, with the need to determine a consistent electric potential as a result.

G. Selection of edge and scrape-off-layer model scenarios
The usual applications for quasi-2d simulations of nonlinear drift dynamics in magnetized plasmas are, for example, fundamental studies on either (a) gyrofluid turbulence and zonal flows with FLR effects, or on (b) FLR effects on interchange dynamics of warm "blob" perturbation propagation. Both cases are commonly treated in separate studies, where for (a) Λ Bs ≡ 0, and instead for (b) Λ c ≡ 0 is set respectively, for better separability and understanding of basic underlying mechanisms.
"Blob" transport is regarded as most relevant in the scrape-off-layer (SOL) of fusion plasmas. Studies with single or few seeded blobs are important to reveal fundamental mechanisms, but in tokamaks blobs are presumed to be generated rather "randomly" around the separatrix, so that turbulent drift wave vortices and zonal flow effects from the closed field line (CFL) edge region will play a large role for consistent studies of "blobby" (intermittent) SOL transport. Several 2d fluid codes (for example as in refs. [62][63][64]) take into account both drift wave and interchange effects and different background conditions by assigning two regions, that are defined by a "separatrix" location x s within the same 2d computational domain.
An additional effect that should be taken into account in the SOL region is sheath coupling of the open field lines with limiter or divertor material walls. This effect intrinsically includes parallel kinetics along the magnetic field direction. Gyrokinetic, gyrofluid and fluid models for sheath coupling conditions may be approximated under assumption of Bohm conditions.
Any perpendicular 2d fluid-type approximation will likely miss relevant physics here, but can in principle be included also in the present 2d full-f gyrofluid model. For test purposes presently a (partly inconsistent) sheath instability "toy model" is optionally included in the TIFF code with additional coupling terms added to the right hand side of eq. (11), which needs further improvement before any application. For example, a 2d delta-f form of the gyrofluid sheath coupling terms was introduced by Ribeiro [47], as andΛ Si → γ DÑe , with sheath coupling parameter γ D as defined in ref. [47], and . This model add-on in principle will allow to examine simplified edge-SOL coupled turbulence and flow dynamics including turbulent generation of SOL blobs, once a consistent full-f full-k sheath coupling term is derived and implemented in the future.
For this reason, the present paper only discusses simulations without sheath coupling.
All coupling termsΛ on the right hand side of eq. (11) are then selectively only applied in the respective regions of interest, for example formally by multiplication with Heaviside type step functions λ(x s ) for a given relative separatrix position 0 ≤ (x s /L x ) ≤ 1.
Application cases and systematic physics studies will be presented elsewhere. Here the focus is on introduction of the code TIFF and its presently underlying model as a reference for later applications and further developments, such as toroidal geometry and inclusion of thermal and electromagnetic dynamics in a (field-aligned) 3d extension. A main aspect here is also on introduction and testing of an efficient dynamically corrected Fourier solver for the generalized Poisson problem in the polarization.

III. NUMERICAL SOLUTION ALGORITHM
The normalized equations solved in the TIFF code are: The general procedure for solution of this set of equations is as follows: (#1) SpecifyN e (x) andN i (x) on an equidistant grid in a 2d rectangular (x, y) domain, either as initial condition, or subsequently updated in each time step by eqs. (12)(13).
The constant gyroradius assumption allows to evaluate all gyro-averaging operations efficiently in k space, hereN Gi (k) =N i (k)/(1 + τ ik 2 /2). In the TIFF code presently the 2d discrete Fourier transform from the FFTW3 library [65] is used for transformations between k and x space representations.
(#4) Prepare the input functions to eq.
(#5) Obtainφ G from eq. (14) with one of the solvers discussed in the Appendix.

The time step update (#9) first requires evaluation of the advective Poisson brackets
and all coupling termsΛ. The brackets {·, ·} are here presently solved with the energy and enstrophy conserving (but not shock capturing) fourth order Arakawa scheme [66,67].
InΛ Bs the curvature operators K(f ) = κ∂ y f are evaluated by (fourth order) centered finite differencing over the (periodic) y direction. Evaluation ofΛ Ss andΛ ce are straightforward.Λ ce includes calculation of zonal averages ⟨f ⟩(x) = (1/n y ) j f i,j , which on a 2d rectangular local grid simply requires summation over all n y grid points j of the y direction.
Time step updating of f ≡ lnN s in the form ∂ t f = F here uses a fourth order accurate three-step Adams-Bashforth method with Karniadakis weights [68]: where c 0 = 18/11, c 1 = 9/11, c 2 = 2/11 and c F = 6/11. The (normalized) time step size ∆t has to be small enough for CFL stability. For further numerical stabilization of the otherwise explicit scheme an artificial sub-grid type viscosity term Λ ν is added. For this a hyperviscosityΛ ν = −ν 4∇ 4 lnN s is here applied. The coefficient 0 ≤ ν 4 ≪ 1 is heuristically chosen to prevent grid instability at smallest scales, as a sink for the nonlinear direct vorticity cascade, and to (slightly) damp out Gibbs type noise that can appear at under-resolved strong gradients due to the not shock capturing nature of the Arakawa discretization. The value of ν 4 (h) needs to be adapted when spatial grid resolution h is changed to achieve optimum results. If desired, for example for comparability to other code implementations, an ordinary "physical" viscous diffusion termΛ µ = +µ(∇ 2N s )/N s could be added to the right hand side of eqs. (12)(13). For usual application of the gyrofluid model to hot fusion edge plasmas the actual viscosity µ in general would be too small to be resolved efficiently by direct numerical simulation.
For solution of the generalized Poisson problem as in step (#5) of the outlined algorithm, presently three methods are implemented in TIFF as described in more detail in the Appendix: An iterative preconditioned conjugate gradient (PCG) solver, an iterative red-black successive over-relaxation (SOR) solver, and a novel dynamically corrected Fourier (DCF) solver. All solvers make use of the results for ϕ (t−1) and ϕ (t−2) of the previous time steps, either for extrapolated initialisation of the iterative schemes, or for the (therefore denoted "dynamical") correction of the approximate Fourier method. The DCF method leads to predictable run times of the generalized Poisson problem, that only depend on the grid size but not on other simulation parameters, whereas the number of iterations and run time can vary significantly for the PCG and SOR solvers for any specified accuracy, in particular when densities (and thusε i (x)) are strongly inhomogeneous.

IV. INITIAL AND BOUNDARY CONDITIONS
Initially the gyrocenter densityN i (x, y) andN e (x, y) fields have to be specified. For full-f initial background profiles a radially exponential decline is set, each withN 0 (x) = Here L p = x s denotes either the width of the pedestal region if a separatrix at x s < L x is applied, or the width L x of the whole x domain if no SOL region is treated. The initial lnN s , profiles (and in its limit the delta-fÑ profiles) are thus linear. On top of this initial background, either single perturbations, such as a Gaussian density blob, or a pseudo-random "bath" of modes, is added with a given amplitude. The initial ion gyrocenter density is either set equal to the electron density, or a "vorticity free" 1N e is used so thatσ = 0. Boundary conditions are applied in several instances: on density profiles in order to maintain a background gradient if required, for the intrisic boundary value problem of the Poisson equation, for solution of the gyro operators in k space, and on finite difference operators in x space. The "poloidal" y direction is assumed periodic. The box length L y in units of ρ 0 has to be chosen much larger than typical perpendicular correlation lengths for turbulence simulations, or much larger than blob or vortex scales for simulation of such structures, in order to minimize self-interaction effects across the periodic y boundary.

A. Density profile boundary conditions
Density profiles could be maintained by specifying source termsΛ Q (x) around the inner (x = 0) radial boundary and sink terms around the outer (x = L x ) boundary.
The presently considered application scenarios are for example a full-f gyrofluid generalization of HW turbulence simulations, or seeded blob simulations on a (usually) constant background density. For comparability to common delta-f implementations of these scenarios it is adequate to maintain an average density profile by prescribing fixed boundary densities N (x = 0) ≡N L andN (x = L x ) ≡N R , while the full densities may still self-consistently evolve in between. Seeded blob simulations could for example useN L =N R = 1.
"Classical" fluid or delta-f gyrofluid HW turbulence simulations usually decouple a constant background gradient in the advective derivative as δ −1 {ϕ,Ñ + N 0 (x)} = δ −1 {ϕ,Ñ } + g∂ y ϕ with a gradient parameter g = δ −1 (ρ s /L n ). For the "diamagnetic drift frequency" normalization of time this amounts to g = 1, but g could also be kept as a free parameter [61]. The restriction in delta-f fluid or gyrofluid simulations to small fluctuations that are B. Parametric transition to delta-f limit within full-f equations Specification of the density boundary in this way allows for a consistent treatment and verification of the delta-f limit of small perturbation amplitudes on large relative background densities within the full-f code. For gradient driven turbulence this can be achieved by reducing all of the background variation, drift scale, and initial perturbation amplitudes by the same small factor ϵ. For example, settingN L = 1 + ϵ · 0.5 andN R = 1 − ϵ · 0.5 for the same box size of L x = 64 gives δ ≈ ϵ · 0.015. When initial (for example blob) perturbation amplitudes ϵ · ∆n are in the full-f model chosen smaller by the same factor (for example ϵ = 1/100) as in a corresponding delta-f code (with initial amplitude ∆n), then this corresponds to the respective delta-f setup and enables a direct comparison (and code cross-verfication) of this limit when ϵ → 0.

C. Radial boundary density condition
It is desirable to reduce density fluctuations and maintain zero vorticity and/or flows on narrow inner and outer radial boundary layers. The boundary layer region of "width" The parameters a L and a R are set to 1 for boundary gradient driven cases (as described above), or can be respectively set to 0 if source and/or sink termsΛ Q (x), or a free outflow condition on the outer boundary, are activated.
The vorticity around the radial boundaries can be approximately set to zero, when equivalently the right hand sideσ =N e (x)−Z i Γ 1iNi (x) of the polarization equation, corresponding to the first order polarization density with FLR effects, is set to zero.
For this purpose it is most convenient to first define the boundary values of the gyroaveraged ion gyrocenter densityN Gi and then compute the consistent electron and ion gyrocenter densities in each time step (#3 in the algorithm), by re-settingN Gi  Existence of a unique solution to the generalized Poisson problem ∇ · P = ∇ · ε∇ϕ = σ requires fulfillment of a compatibility boundary condition. Taking the domain integral over the (x, y) area S, one has S dx ∇·P = S dx σ, and therefore δS dl P·n S = S dx σ. This is ensured by the above vorticity-free and flow-free (n S ·∇ϕ = ∂ x ϕ = u y ≡ 0) conditions on the boundary δS in x (and periodicity in y). The condition correponds to global conservation of the polarization charge density.

D. Mirror padding for Fourier transforms
The algorithm involves four evaluations of gyro-operators (in steps #2, 3,4,6), which is in the present isothermal model achieved with high accuracy in k space by applying (FFTW3 library) Fourier transforms. Actually the evaluation ofN N Gi in Padé approximation for the vorticity free boundary condition (step #3) could alternatively be achieved by (considerably faster but less accurate) finite differencing in x space. Although more costly, this evaluation is here for consistency also done in k space.
Standard Fourier transformation requires periodic boundary conditions, else discontinuities would introduce Gibbs noise artefacts. The radial "physical" domain 0 ≤x ≤ 1 however usually includes density profiles withN s (x = 0) ̸ =N s (x = 1). For this reason, ghost domains in an extendedx direction are introduced for the (x, y) field arrays before solution by Fourier transforms to k space, and also before solution of the polarization equation.
Doubling the (initially "quarter-wave" physical) domain to 0 ≤x ≤ 2, by copying the initial array f (x,ŷ) into the region 1 <x ≤ 2 and defining symmetrically mirrored It is as usual favourable to use power-of-two numbers of grid points in thex andŷ directions, so that the FFTW3 library will make use of Fast Fourier Transform (FFT) algorithms. Any other grid size can be chosen, but FFTW3 then automatically uses somewhat (depending on the grid size) slower Discrete Fourier Transform (DFT) methods.

V. PARALLELIZATION AND REPRODUCIBILITY
The present development and production run platforms for the TIFF code are multi-core shared memory office workstations, so that multi-core and/or multi-thread parallelization is simply achieved by OpenMP (OMP) parallelization of grid array sized double loops, and use of the respective OMP library of FFTW3 for the Fourier transforms. The frequent transistions between single and parallel regions, required by the above algorithm, prevent good scaling. Efficient speed-up is usually achieved (depending on the input parameters and on the hardware) when between 8 and 32 parallel threads are used.
Bitwise reproducibility of subsequent executions on the same system can be achieved when the "ESTIMATE" flag for the transform planner routine of the FFTW3 library is used, which is desirable for verification and testing purposes. The "MEASURE" flag is faster in execution (which would be desirable for long production runs with for example millions of The 2d dynamical field arraysN e (i, j),N i (i, j),φ(i, j),ω(i, j),σ(i, j), andÑ e (i, j) = N e (i, j) −N 0 (i) are regularly written out (e.g. to a RAM disc), and can be viewed (already during run time) with any 2d visualization software (for which presently a simple gnuplot script is used). In addition, 1d cross sections of several arrays f (i, j 0 ) and f (i 0 , j), and averaged radial profiles ⟨f ⟩ j (i) are written.
Spectra ⟨f (k x )⟩ j and ⟨f (k y )⟩ i are computed at each diagnostic output step for several quantities, such as kinetic energy, enstrophy, and density and potential power spectra.
Energetic and transport quantities are recorded as time traces of (x, y) domain averages.
The (normalized) E-cross-B advective electron particle transport is obtained as Q n (t) = (1/S) dxN e (x, y)∂ yφ (x, y), where S = L x Ly, and the y derivative is provided by simple centered finite differencing (as also in the following diagnostics).
The full-f global thermal free energy is given by (compare ref. [59]): The total energy is an ideal conserved quantity, and can be used as a diagnostic for saturation of a turbulent state.
Further diagnostics can be included, for example output and computation of difference norms for solver testing against a constructed solution, or center-of-mass calculation of an interchange driven "blob". The equivalently normalized delta-f set of equations is: The ( form of this delta-f polarization is optionally obtained by setting the term τ ik 2 ≡ 0 herein. The gyrofluid ion potential (step #7) is simplyφ i = Γ 1iφ .
For diagnostics, the delta-f total thermal free energy is evaluated as The dynamical context is important, because the underlying "Teague method", introduced in different context, for approximate solution of a generalized Poisson-type problem, is in itself not very accurate. Here however an efficient dynamical (or recursive) correction based on the solution from the previous time step is newly applied with the "Teague method", which is shown to improve its accuray by two orders of magnitude, and by that can achieve similar accuracy as iterative solvers with a "reasonable" (affordable) number of iterations in a turbulence code.

A. Teague's original method
In optics, the "transport of intensity equation" (TIE) is an approximate relation for the intensity and the wave phase of a coherent beam in an optical field [69,70]. The underlying mathematical form of the TIE is basically equivalent to the 2d generalized Poisson problem: where I is the known (measured) planar intensity distribution at a distance z, k the wave number, and ψ the sought phase of the wave. Various solution methods have been used on the TIE in the field of applied optics [70], but a now widely used approximate Fourier method has been suggested by Teague in 1983 [69] by introducing what is now commonly known as Teague's auxiliary function.
In the following this original "Teague's method" is described in terms of the notation introduced in the previous sections in context of the gyrofluid polarization equation (and not in the TIE notations used in optics). Starting with an auxiliary 2d scalar function p(x) is introduced that is supposed to satisfy This reduces the generalized Poisson problem to an ordinary 2d Poisson equation The right hand side contains the (by now) known quantities ε(x) and p(x), and can be evaluated by standard fourth order 2d centered finite differencing in x and y. Eq. (22) can therefore again be solved (e.g. in k space) to obtain ϕ(x). When Fourier transforms are used, the method involves four evaluations (two forward and two backward) of the 2d transform F, which can be formally expressed (compare ref. [71] for the TIE version) as: Rather, the complete Helmholtz decomposition of the vector field P(x) is: with a scalar potential p(x) and a vector potential H(x). This general form provides the same solution path for p(x) like above, as still ∇ · P = ∇ 2 p = σ holds. The next step, generalizing the result of eq. (22), now gives: with ∇·(1/ε)∇×H = (∇×H)·∇(1/ε) = {(1/ε), η}, where the 2d Poisson bracket notation (defined above) is used, and η is the z component of (unknown) H(x, y) = η(x, y) e z .
To find a constraint on η, one can apply the curl operator on both sides, instead of the divergence operator ∇ · ∇ϕ = ∇ · ..., correspondingly as in eq. (25). This gives the condition This determining relation for the unknown η however again has the form of a generalized Poisson equation, which sends us back to start...
Another constraint can be generated by taking the curl on P, which gives the relation ∇ × (ε∇ϕ) = ∇ × ∇ × H. This can be rephrased as: This first shows that a sufficient condition for the accuracy of Teague's method would be that {ϕ, ε} = ∇ϕ × ∇ε ≡ 0, which holds if the isocontours of ϕ and ε were aligned everywhere [72]. The relation can also be used to derive constraints for the error norm [72].
If applied to gyrofluid simulations, the error is surely not negligible as {ϕ, ε} ∼ {ϕ, N i } for constant magnetic field, and thus directly related to the advective ion nonlinearity.

B. Iterative and dynamical corrections
Iterative methods have been suggested to correct the error by the original Teague approximation. In ref. [74] a Picard-type iteration is applied, that uses an initial solution ϕ 0 by Teague's method (without need to refer to the vector potential) to compute an according source term σ 0 ≡ ∇ · ε∇ϕ 0 , then use ∆σ = σ 0 − σ in another turn of application of Teague's method to obtain a correction ∆ϕ, and repeat until ∆ϕ is smaller than a specified error.
The set of equations (25) and (27)  has to be carried out in every of very many time steps in a dynamical turbulence simulation.
It will be shown that already after one or two iterations a high accuracy is reached. This motivates another possible non-iterative correction to Teague's method in the context of a dynamical simulation with time evolution of ϕ(t) in the presence of small time steps. In simulations of fully developed turbulence, and in particular when an explicit finite difference scheme is used like in the present code, the time step ∆t, as it appears in eq. (15), is small, and so accordingly is the difference between successive solutions ϕ (t) and ϕ (t−1) .
The idea is to therefore evaluate with an "old" solution from previous time steps in a dynamical simulation, and with this approximation compute In contrast to iterative schemes applied within each time step, this reduces the additional expense for the correction to only one (for example FFT based) inversion of the conventional of the previous time step. In the present code an extrapolation is applied as: with a free estimator factor a ∈ (0, 1). Because this extension of Teague's method uses previous results of a time evolving simulation, and evaluates the multiple occuring inversions of the Poisson problem with Fourier solvers in k space, it is in the following abbreviated as the "DCF" (dynamically corrected Fourier) method. The accuracy however will here depend on the size of the (small) dynamical time step.

C. Unit test of generalized Poisson solvers
The specific implementation of the PCG, SOR and DCF solvers as used in the TIFF code is in detail described in the Appendix. The code can be run with a unit testing option, by setting a flag in the input parameter file, that initialises analytical constructed functions ε c and σ c and only calls the specified Poisson solver once, so that the numerical solutions for ϕ can be directly compared with the analytical function ϕ c . The main purpose of this test is to determine the applicability and accuracy of the novel dynamically corrected Fourier (DCF) solver in comparison to the established SOR and PCG methods.
The constructed and numerical solutions are shown in Fig. 2 (left) as a function of x at y 0 = (L y /2 − 5). The differences ∆ϕ(x, y) between constructed and numerical solutions can be visualised and compared by their global L 2 norms as a function of resolution or iterations.
In Fig. 2 (right) the global error err = ||∆ϕ|| 2 is generally large for g = 1 when a background gradient in ε c (x) (solid lines) is present and discontinuities in the derivatives occur at the x boundaries. For comparison, also results for the above constructed solution but with g = 0 (constant background, periodic ε c (x), dashed lines) is shown, with reduced more "ideal" errors. In practice, the g = 0 case would for example correspond to a "seeded blob" simulation scenario, whereas the g = 1 case would be relevant to gradient driven drift wave turbulence scenarios.    Fig. 2), obtained with the "recursively corrected Fourier" (RCF) method after 4 iterations.
The second order accurate SOR scheme (orange lines in Fig. 2) has much slower convergence and does for this resolution never reach the accuracy of the (fourth order) RCF and PCG schemes. Usually several hundred SOR iterations are required for acceptable results x (i.e., a slope of -2 in the log-log plot). The "ideal" RCF scheme (bold black line) here denotes a correction term calculated once from the constructed solution ϕ c instead of recursively, which also shows a fourth order (N −4 x ) slope, as expected. The gradually thicker black lines (from top to bottom) show, topmost, the error of the original "Teague method" without correction, which is for usual resolutions always in the range of 2·10 −3 for g = 0 (left figure). The RCF(1) scheme uses one recursion of the solver in the unit test problem, and already gives an improvement of the error by around two orders of magnitude. The RCF(2) scheme with two recursions follows the "ideal" dependence until around N x = 256 and for higher resolution saturates at an error of around 2 · 10 −7 . The gradually thicker blue lines (from top to bottom) show the error by the PCG scheme after 4, 5 and 50 iterations. The n iter = 5 case for PCG has for high resolution an error in between the results from the RCF(1) and RCF(2) schemes.
The right Fig. 3 shows the corresponding results for g = 1 in the presence of a background gradient, which in the present implementation of the solvers in TIFF overall reduces the achievable accuracy. The general behaviour of the different solvers is similar to the g = 0 results, but no "ideal" scaling of the error (-4 power law) with resolution is obtained any more. The RCF(1) and PCG(4) schemes also show similar errors for medium to large resolution.
The dynamically corrected Fourier (DCF) scheme needs three evaluations of the standard Poisson problem, which is here presently achieved by a fast Fourier solver. In a dynamical simulation the error correction then is calculated from the previous time step solution. The PCG scheme needs one evaluation of a standard Poisson inversion (e.g. by Fourier solver) per iteration step, so that PCG (n = 3 − 4) has approximately the same computational cost as a DCF evaluation (which in addition needs another call of a Poisson bracket evaluation).
The necessary number of iterations in PCG or SOR to reach a given accuracy can strongly depend on the complexity of the problem, basically determined by the degree of variability of ε(x, y, t). The computational expense of the DCF method only depends on the size of the grid. The DCF scheme therefore presents itself as a viable alternative method (for moderate accuray) with predictable run times.

D. Cross-verification of full-f full-k blob simulation
The unit testing of the PCG, DCF and SOR solvers above has been applied on the "pure" generalized Poisson equation ∇ · ε∇ϕ = σ. The implementation in the TIFF code for solution of the full-f full-k gyrofluid polarization eq. (9) requires for τ i ̸ = 0 two additional applications of (here Fourier based) solvers for the √ Γ 0 operator, and further (also Fourier based) evaluations of the gyro-averaging operator Γ 1 N i in σ, and of Γ 1 ϕ in the ion gyrocenter density advection.
As a further test, an overall cross-verification of the code is intended by running a seeded "blob" simulation for parameters that correspond to a recent study with the full-f full-k gyrofluid isothermal version of the FELTOR code, described by Held and Wiesenberger in ref. [59] and shown in Figures 4-7   differ between the code results, because in ref. [59] FELTOR uses a discontinuous Galerkin method with 300 grid cells and 5 polynomial coefficients, which would correspond to 1500 grid points for the TIFF solvers. Overall the agreement between the codes and between the different solvers used in TIFF can be regarded as satisfactory. In particular, the new DCF method has been shown to achieve strong agreement with the other solvers, and in particular with the also fourth order PCG scheme.
The results shown in Fig. 5 have been obtained on a Threadripper PRO 5975WX Linux workstation with 256 GB RAM, running on 32 threads of the CPU. The run time was 216 min with the SOR solver, also 216 min for the PCG solver, 67 min for the DCF solver, and 60 min for the uncorrected Teague solver. The times for the zero vorticity runs (Fig. 4) were 212 min for the SOR solver, 116 min for the PCG solver, and 68 min for the DCF The morphology of 2d drift wave turbulence in the full-f gyrofluid OHW model is illustrated in Fig. 6 as a snapshot in the saturated turbulent phase, for parameters described below. The left frame shows the electric potential ϕ(x, y), which acts as a streamfunction for the advecting turbulent E-cross-B velocity. The right frame shows the vorticity Ω(x, y) = ∇ 2 ϕ with the characteristic thin vorticity sheaths induced by FLR spin-up [53].  The new DCF scheme appears to be both sufficiently accurate and efficient to be considered for further use. It should be noted that the 2d FFT evaluation is a memory bound application and profits from high available bandwidth. The relative performance between the solvers can therefore differ between hardware systems. These relative energy and mass change rates as a function of time are shown in Fig. 9 for the same turbulence simulation parameters as above, obtained with the DCF scheme.
A linear regression for the values in 500 ≤ t ≤ 5000 gives a relative tendency ⟨R E ⟩(t) ∼ −4 · 10 −6 − 7 · 10 −9 t for the energy, and ⟨R M ⟩(t) ∼ −8 · 10 −8 − 2 · 10 −10 t for the particle number. Both constitute very small loss rates and can be regarded as sufficiently good conservation property of the numerical scheme. The standard deviations of the nonlinear growth rate fluctuations are s(R E ) = 4 · 10 −3 and s(R M ) = 2 · 10 −4 .
C. Full-f model in small amplitude limit vs. delta-f model The full-f full-k model should agree in the limit of small amplitudes with the deltaf model, as discussed in sections IV B and VII. To test this we apply largely the same parameters as above to OHW simulations with the DCF solver:α = 0.2, τ i = 1, κ = 0,

X. CONCLUSIONS AND OUTLOOK
An isothermal quasi-2d gyrofluid model and a code (TIFF) for arbitrary amplitude (fullf) and arbitrary wavelength (full-k) drift instabilities and turbulence in magnetized plasmas has been introduced. A major aspect was on testing of a here newly suggested "dynamically corrected Fourier" (DCF) solver for the generalized Poisson equation, which is based on the original approximate Teague method. The generalized (a.k.a. "variable") Poisson problem appears in the solution of the gyrofluid (similar to the gyrokinetic) polarization equation, which couples the evolving gyrocenter plasma densities with the electric potential. The present approach retains the complete FLR effects and spatial and dynamical variations of the polarization density. The new DCF solver was shown to be a viable, sufficiently accurate and efficient method for dynamical turbulence simulations coupled to the generalized Poisson problem, which motivates its further use in future extensions of the code.
In its recursive "RCF" form (with for example four iteration) as a high accuracy extension of Teague's method the scheme could also find applications in the field of optics as a standalone application for solution of the transport of intensity equation (TIE).
The main purpose of the present work is as reference and test case for the TIFF model and code, whereas detailed physics studies with the present code, such as for example on FLR effects on zonal flows, or on properties of turbulently generated "blobby" (intermittent) transport, will be discussed elsewhere.
The extension of the 2d TIFF code to a magnetic field-aligned flux-tube like 3d geometry can directly make use of the here tested polarization solvers, as these always only act in the locally perpendicular 2d drift plane. Future developments of such a corresponding 3d full-f full-k gyrofluid edge turbulence code will include the addition of temperature evolution equations (and through this applicability on temperature gradient driven modes and use of Landau damping mechanisms), and generalization to electromagnetic drift-Alfvén dynamics.

Funding
This work was supported by the Austrian Science Fund (FWF) project P33369.

Data availability
The code TIFF is openly available at: https://git.uibk.ac.at/c7441036/tiff Appendix A: PCG scheme Conjugate gradient methods for the solution of linear systems are introduced in numerical textbooks such as by LeVeque [75]. Here a pre-conditioned conjugate gradient (PCG) scheme, with a pre-conditioner and algorithm suggested by Fisicaro et al. [77], is applied.
A generalized Poisson operator A ≡ ∇ · ε∇ is defined so that a solution of the system Aϕ = σ is sought. The residual of an approximate solution ϕ (n) , say obtained in step n during an iteration, is r (n) = σ − Aϕ (n) . Within the context of a dynamical simulation with small time steps ∆t, an initial guess can be obtained by extrapolating the converged solutions from previous time steps towards with a free estimation factor a ∈ (0, 1), and ϕ (0) (t = 0) = 0 in the first time step. From this the initial residual r (0) = σ − Aϕ (0) for the iteration is obtained.
The PCG scheme basically searches for the minimum of a (quadratic) residual function f (ϕ) = (1/2)(ϕ, Aϕ) − (ϕ, σ), which is given from a 2d Hessian Taylor expansion, by eval- denotes the inner product of two matrices A and B. A pre-conditioning matrix q (0) is evaluated once (in each dynamical time step) and remains constant throughout the PCG iteration.
The algorithm allows to efficiently re-use several already calculated terms and products.
A classical conjugate gradient (CG) solver without preconditioning of the residual in step (1) would set P = 1 and v (n) = r (n) in the algorithm. The idea is that a preconditioned residual v = P −1 r = P −1 σ − P −1 Aϕ is minimized in the same way as the original r.
Application of a suitable preconditioner can speed up convergence significantly, but needs to be chosen carefully. Here the preconditioner suggested in ref. [77] for the generalized Poisson problem is applied as P() ≡ √ ε∇ 2 √ ε(), and accordingly q (0) ≡ √ ε∇ 2 √ ε above. (4) is greatly simplified by this preconditioner [77]. For the inversion in step (1) here needs to be evaluated including the application of a standard fast Poisson solver, which is here presently achieved by an FFT solver in k space.
The additional expense of calling a standard Poisson (∇ 2 u = f ) solver at each iteration step is paid off by the in general rapid convergence of this PCG scheme. The default choice for the standard Poisson solver in the TIFF code is by FFT, but optionally also a (generally slower) SOR scheme is available. The order of this PCG solver is determined by the order of evaluation of the Laplacians (in q (0) , or in a non-FFT standard Poisson solver), which is here achieved in fourth order accuracy.
The algorithm above could also be converted into a (preconditioned) steppest descent scheme by setting β (n) = 0, which more serves didactical than practical purposes because of its generally slower convergence.
Starting with an initial guess ϕ 0 (x), the relation (32) can be iterated, with grid values of ϕ (n) i,j applied on the r.h.s. in order compute the updated ϕ (n+1) i,j on the left, until the residual error norm ||R (n+1) || with R (n+1) = ϕ (n+1) − ϕ (n) is smaller than a specified limit.
The grid array is swept in odd ("red") -even ("black") order, as in eq. (32) the values of ϕ (n) i,j for even indices i or j depend only on odd indexed grid values of ϕ (n−1) , and vice versa. The over-relaxation parameter ω is determined by Chebyshev acceleration according to ω 0 odd = 1 and ω 0 even = 1/(1 − r 2 /2) as initial values, and further ω = 1/(1 − r 2 ω/4), updated in each half-sweep. Here the spectral radius is calculated as r = [cos(π/n x ) + cos(π/n y )]/2 for a homogeneous equidistant grid with grid point numbers n x and n y . Each "red" and "black" sweep through the 2d checkerboard-like grid is loop parallelized with OpenMP, respectively.
Convergence can be significantly accelerated by an initial guess for ϕ 0 (x) based on the solution from previous time step(s) within a dynamical simulation for small time steps ∆t, by extrapolating (in similar spirit as for the over-relaxation factor above) ϕ 0 ≡ ϕ(t − 1) + a · (ϕ(t − 1) − ϕ(t − 2)).
with a free prediction factor a ∈ (0, 1). The first time step uses ϕ 0 = 0, in the second step a = 0, and in later times a factor between 0.5 and 1.0 has been found to be most efficient.
The rate of convergence depends on the degree of (non)uniformity of ε(x). In classical electrostatics the permittivity is usually only weakly varying within one medium, but may have discontinuities between neighbouring media, and is often set in complicated geometries, so that mostly rather finite element schemes instead of finite difference schemes are employed.
In the present application on the gyrofluid polarization equation in a small local section of an edge plasma, the grid can be simply chosen as rectangular, but the polarization density can be strongly inhomogeneous for large amplitude turbulent fluctuations.
The present SOR implementation is only second order accurate, but straightforward to implement and can serve as a reference for the fourth order PCG and DCF schemes. For evaluation of the standard Poisson problem (∇ 2 u = f ) also a fourth order accurate SOR scheme is available in the TIFF code, which uses discretization by a Collatz Mehrstellenverfahren (9pt stencil for the Laplacian on u, and a 5pt stencil correction for f ).

Appendix C: DCF scheme
The here introduced "dynamically corrected Fourier" (DCF) approach for the generalized Poisson problem adds a correction term, computed from the result of the previous time step within a dynamical simulation, to "Teague's method" (compare main text above). Teague's approximate method [69] had originally been devised for solution of the TIE (transport of intensity equation) in optics [70], and can be efficiently evaluated in Fourier space [71].
It is here innovatively applied on solution of the electric potential ϕ(x) from the gyrofluid (or gyrokinetic) polarization equation, and within a dynamical context. For small time step sizes ∆t this allows to re-use solutions of ϕ old from past times, instead of a (more expensive) iterative error correction. An extrapolated prediction with a free estimation parameter a ∈ (0, 1) based on two previous solutions is used for the corrector, as ϕ old ≡ ϕ (t−1) + a · (ϕ (t−1) − ϕ (t−2) ).
The generalized Poisson equation ∇ · ε∇ϕ = σ is formulated by means of a polarization density vector field P ≡ ε∇ϕ ≡ ∇p + ∇ × H in terms of a scalar potential p and a vector potential H = η e z . Input quantities are the known 2d fields σ(x) and ε(x).
By application of the divergence on ∇ · P = σ the scalar field p is first obtained from: The inversion of the Laplacian in the standard Poisson problem is here obtained in k space by Fourier transforms, using the FFTW3 library, as p = −F −1 k −2 Fσ.
Further, an approximation for A is obtained (see main text above) from ∇ × P as: This estimated quantity is used in the extended, dynamically corrected relation: Like above, the inverted Laplacians in eqs. (37) and (38) can be efficiently solved in Fourier space. The term ∇ · (1/ε)∇p is evaluated by standard (fourth order) centered finite differences, and the Poisson bracket can be computed either also by simple fourth order centered differences, or by re-using the Arakawa scheme (introduced for the advective terms in the dynamical gyrocenter density equations).
For a solver unit test (with a constructed solution) the set of equations (37) and (38) can also be applied recursively (instead of within a dynamical context) by re-setting ϕ → ϕ old directly after each iteration step, starting with η o = 0. In the main text above this is refered to as "RCF" for recursively corrected Fourier scheme.