Iterative POLAR elliptic operator with anisotropic TENSOR coefficient (cross terms). More...

#include <pops/core/foundation/types.hpp>
#include <pops/mesh/layout/box_array.hpp>
#include <pops/mesh/layout/distribution_mapping.hpp>
#include <pops/mesh/storage/fab2d.hpp>
#include <pops/mesh/execution/for_each.hpp>
#include <pops/mesh/geometry/geometry.hpp>
#include <pops/mesh/storage/mf_arith.hpp>
#include <pops/mesh/storage/multifab.hpp>
#include <pops/mesh/boundary/physical_bc.hpp>
#include <pops/parallel/comm.hpp>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <stdexcept>
#include <vector>

Include dependency graph for polar_tensor_operator.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	pops::detail::PolarApplyKernel
	L_int(phi) = div(A grad phi) in polar (apply). More...

struct	pops::detail::PolarJacobiApplyKernel
	out = (f - L0 phi) / \|diag\| – one Jacobi iteration (point-by-point relaxation) on the DIAGONAL polar stencil (cross terms excluded from the Jacobi splitting: they stay on the right-hand side via the residual). More...

struct	pops::detail::PolarCopyKernel
	Copy component 0 (dst <- src). More...

struct	pops::detail::PolarInvDiagKernel
	Computes idiag = 1 / diag of the diagonal polar stencil (for Jacobi). diag = polar_diag (< 0). More...

struct	pops::PolarKrylovResult
	Result of a polar BiCGStab solve: iterations, relative residual, convergence. More...

class	pops::PolarTensorKrylovSolver
	MATRIX-FREE BiCGStab Krylov solver for the FULL-tensor POLAR elliptic operator L_int(phi) = div(A grad phi), A = [[a_rr, a_rt], [a_tr, a_tt]] possibly NON symmetric. More...

Namespaces
namespace	pops

namespace	pops::detail

Concepts
concept	pops::PolarLinearSolver
	Contract of the iterative POLAR elliptic operators: same shape as PolarEllipticSolver (cf.

Enumerations
enum class	pops::PolarPrecond { pops::Jacobi , pops::RadialLine }
	Choice of the SIMPLE BiCGStab PRECONDITIONER (NO MG V-cycle – stagnation on polar 1/r^2, cf. More...

Functions
POPS_HD Real	pops::detail::polar_radial_div (const ConstArray4 &p, const ConstArray4 &arr, bool hrt, const ConstArray4 &art, int i, int j, Real ri, Real rfm, Real rfp, Real idr, Real idth)
	RADIAL FACE contribution to the diagonal + cross stencil, at (i, j).

POPS_HD Real	pops::detail::polar_azimuthal_div (const ConstArray4 &p, const ConstArray4 &att, bool htr, const ConstArray4 &atr, int i, int j, Real ri, Real idr, Real idth)
	AZIMUTHAL FACE contribution to the diagonal + cross stencil, at (i, j).

POPS_HD Real	pops::detail::polar_diag (const ConstArray4 &arr, const ConstArray4 &att, int i, int j, Real ri, Real rfm, Real rfp, Real idr, Real idth)
	Diagonal (coefficient of phi_{i,j}) of the diagonal POLAR stencil (radial + azimuthal), for the Jacobi preconditioner.

void	pops::apply_polar_tensor (const MultiFab &phi, const PolarGeometry &geom, MultiFab &lap, const MultiFab a_rr, const MultiFab a_tt, const MultiFab a_rt, const MultiFab a_tr)
	Applies L_int(phi) = div(A grad phi) in polar over the whole MultiFab.

Detailed Description

Iterative POLAR elliptic operator with anisotropic TENSOR coefficient (cross terms).

Path A stage 2a (foundational brick toward the polar Schur).

CONTEXT. To build the FULL STIFF Euler-Poisson system in POLAR geometry (a strongly magnetized regime, high omega_c), the condensed Schur stage (CondensedSchurSourceStepper, level 4 of docs/SCHUR_CONDENSATION_DESIGN.md) is needed on the polar side. The lock is the elliptic OPERATOR: the Schur condenses a FULL TENSOR operator A = I + c rho B^{-1} where B^{-1} is the Lorentz rotation, which injects CROSS terms a_rt / a_tr (and a theta-dependent coefficient as soon as rho or B_z varies in theta). The existing DIRECT PolarPoissonSolver (FFT-in-theta + tridiag-in-r, polar_poisson_solver.hpp) diagonalizes theta via the FFT: it REQUIRES a coefficient CONSTANT in theta and WITHOUT cross coupling. It is therefore structurally INCOMPATIBLE with the Schur operator. This operator is the ITERATIVE (matrix-free) counterpart that accepts the full tensor.

EQUATION. On a ring (r, theta), r in [r_min, r_max] > 0 (NO r=0 singularity), theta PERIODIC, we solve the operator in DIVERGENCE form with polar metric and tensor A in the ORTHONORMAL physical frame (e_r, e_theta): L(phi) = div(A grad phi), grad phi = (d_r phi, (1/r) d_theta phi), div F = (1/r) d_r(r F_r) + (1/r) d_theta(F_theta), A = [[a_rr, a_rt], [a_tr, a_tt]] (possibly NON symmetric: a_rt != a_tr). The flux is F = A grad phi: F_r = a_rr d_r phi + a_rt (1/r) d_theta phi F_theta = a_tr d_r phi + a_tt (1/r) d_theta phi ISOTROPIC case a_rr=a_tt=1, a_rt=a_tr=0: L = (1/r) d_r(r d_r phi) + (1/r^2) d_theta^2 phi (scalar polar Laplacian, cf. PolarPoissonSolver).

CONSERVATIVE FINITE-VOLUME DISCRETIZATION (order 2, like assemble_rhs_polar / PolarPoissonSolver). Cell (i, j) of volume weight r_i dr dtheta. We integrate div(F) over the cell and divide by the weight, which gives on each face the flux carried by its FV LENGTH:

RADIAL faces i+-1/2 (radius r_face, face weight r_face dtheta): DIAGONAL radial term: (1/r_i) [ r_{i+1/2} a_rr^{i+1/2} (phi_{i+1,j}-phi_{i,j})/dr
- r_{i-1/2} a_rr^{i-1/2} (phi_{i,j}-phi_{i-1,j})/dr ] / dr; CROSS radial term a_rt (1/r) d_theta phi at face i+-1/2 (d_theta averaged over 4 corners, 1/r_{face} local to the face);
AZIMUTHAL faces j+-1/2 (radius r_i, face weight dr): DIAGONAL azimuthal term: a_tt^{j+-1/2} (phi_{i,j+1}-2 phi+phi_{i,j-1}) / (r_i^2 dtheta^2); CROSS azimuthal term a_tr d_r phi at face j+-1/2 (d_r averaged over 4 corners). FACE coefficient = ARITHMETIC mean of the two adjacent centers (order 2; for the radial/azimuthal diagonal term of a smooth medium the arithmetic mean suffices; the cross term is not a normal flux -> arithmetic mean as well, cf. cartesian poisson_operator.hpp). The coefficients a_rr/a_tt/ a_rt/a_tr are at the CENTER of the cells (1 component each, ghosts filled by the caller).

METRIC VALIDITY: the radial faces use r_face(i) > 0 (ring, r_min > 0); no division by r=0. The stencil diagonal is the SUM of the face weights (for the Jacobi preconditioner); in the isotropic case it gives back EXACTLY the scalar polar stencil.

BOUNDARY CONDITIONS. theta PERIODIC (periodic ghosts, fill_ghosts with ylo/yhi=Periodic). r_min/r_max: PHYSICAL BC Dirichlet (value at the face, reflection ghost 2 v - interior) or homogeneous Neumann (Foextrap, ghost = interior). This is exactly the cartesian fill_ghosts applied in (r, theta) (index direction 0 = radial, 1 = azimuthal; cf. PolarGeometry convention). The Dirichlet data is folded into the residual (AFFINE operator) for the true residual r0, and LINEARIZED (offset c_bc subtracted) for the matvec in the loop (correction directions): same mechanics as the cartesian TensorKrylovSolver (krylov_solver.hpp).

SINGULAR OPERATOR (pure radial Neumann + periodic theta, no reaction): the CONSTANT is in the kernel of L_int; BiCGStab diverges without treatment. We FIX THE GAUGE by PROJECTION onto the subspace of ZERO FV MEAN (weighted mean r dr dtheta removed from the initial residual, the solution, and each preconditioned correction direction). This is the ITERATIVE counterpart of the mode-0 pinning of the direct PolarPoissonSolver, without perturbing the stencil. Detected automatically (both radial boundaries non-Dirichlet); at least ONE Dirichlet boundary => invertible operator, NO pinning, bit-identical Dirichlet path. The solution is then defined modulo a constant.

SOLVER: MATRIX-FREE BiCGStab Krylov (handles the NON symmetric cross term), PRECONDITIONED by a SIMPLE preconditioner (NO MG V-cycle: the scoping, polar_poisson_solver.hpp warning, reports that the MG V-cycle STAGNATES on the polar 1/r^2 anisotropy; we stick to Krylov + simple precond as requested). Two preconditioners (cf. enum PolarPrecond):

Jacobi (diagonal): the simplest, but the iteration count GROWS like 1/h^2 (ill-conditioned Laplacian). On a fine grid it PLATEAUS without converging (cf. test, block E). Sanity check / fallback.
RadialLine (DEFAULT): inverts EXACTLY the RADIAL tridiagonal of the diagonal block per theta line (Thomas), azimuthal diagonal lumped. Attacks the strong radial coupling + the 1/r^2 anisotropy; LOW iteration count with MODERATE growth (measured: isotropic ~ x2 per grid doubling, tensor ~ x2.4). Stays "simple" (no MG hierarchy, no coarse grid). HONEST LIMITATION: the CROSS term and the AZIMUTHAL coupling are NOT in the preconditioner (which only inverts the radial part) -> the iteration count of the tensor case grows faster than the isotropic one. This is acceptable for the foundational brick (SOLID convergence, no stagnation observed); an azimuthal preconditioner (theta line) or a robust MG-line would be a LATER refinement. BiCGStab has a FIXED memory footprint (r, rhat, p, v, s, t, + the preconditioned phat/shat) and does not store a growing Krylov basis.

SCOPE: MULTI-RANK MPI AND MULTI-BOX (several boxes PER RANK) by AZIMUTHAL splitting (theta only) for the RadialLine precond, or by FREE 2D TILING (r AND theta) for the Jacobi precond. The dot products (dot/L2 norm) are collective (all_reduce via mf_arith::dot) and local_size()==0-safe; they loop over ALL the local boxes (mf_arith). The matvec goes through fill_ghosts (MPI halo exchange + physical BC) which also fills the DIAGONAL CORNERS p(i+-1, j+-1) between adjacent boxes – read by the CROSS TERMS a_rt/a_tr of the 9-POINT stencil (fill_boundary enumerates the 8 shift directions, fill_physical_bc extends the radial BC into the theta halo: without this corner, the cross term would be WRONG at the box boundary, cf. test_polar_schur_multibox). The gauge projection project_mean accumulates sum/vol over ALL the local boxes then all_reduces (numerator + denominator). The RadialLine precond (Thomas in r) requires that a full RADIAL LINE (column i=0..nr-1 at fixed theta) fit in ONE box (otherwise the sequential sweep in r would cross a box boundary): we therefore require that EACH box cover the FULL RADIAL RANGE of the domain (splitting in theta only, BUT ANY number of theta boxes per rank). M^{-1} is then BLOCK-DIAGONAL per box (no communication along r); the tridiagonals are per LOCAL BOX (line_*_[li]). A splitting that cuts r UNDER RadialLine is REFUSED by an explicit guard (check_radial_columns, instead of a silent wrong result). The Jacobi fallback (PolarPrecond::Jacobi, per cell) has NO layout constraint: it accepts a full 2D TILING (cuts r AND theta), at the cost of an iteration count growing like 1/h^2. The theta splitting is legitimate anyway: theta is the PERIODIC direction and carries the unstable azimuthal mode; cutting r (few cells, strong radial coupling) would be counterproductive. SINGLE-RANK / SINGLE BOX: BIT-IDENTICAL path (all_reduce = identity in serial; the local_size() loop = the fab(0) loop when there is only one box covering the whole ring).

ADDITIVE: no existing path is touched. The DIRECT scalar PolarPoissonSolver stays UNTOUCHED (separate path); the CARTESIAN Schur stays BIT-IDENTICAL. This header is OPT-IN (stage 2b = wiring a polar SchurReconstructKernel + a coupler, out of scope here).

include/pops/numerics/elliptic/polar/polar_tensor_operator.hpp File Reference¶

Classes

Namespaces

Concepts

Enumerations

Functions

Detailed Description