PolarCondensedSchurSourceStepper: SOURCE STAGE condensed by Schur in POLAR geometry (ring (r, theta)). More...

#include <pops/core/foundation/types.hpp>
#include <pops/core/state/variables.hpp>
#include <pops/coupling/schur/core/schur_source_kernels.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/numerics/elliptic/polar/polar_tensor_operator.hpp>
#include <pops/numerics/linalg/lorentz_eliminator.hpp>
#include <pops/parallel/comm.hpp>
#include <stdexcept>

Include dependency graph for polar_condensed_schur_source_stepper.hpp:

Go to the source code of this file.

Classes
struct	pops::detail::PolarSchurOperatorCoeffKernel
	POLAR condensed tensor coefficients A = I + c rho B^{-1} at cell centers. More...

struct	pops::detail::PolarSchurExplicitFluxKernel
	EXPLICIT flux F = rho B^{-1} v^n = B^{-1}(mr, mtheta) at the center (physical components (e_r, e_theta)). More...

struct	pops::detail::PolarSchurRhsAssembleKernel
	rhs_polar(i,j) = lap_polar(i,j) (= Lap_polar phi^n) + g * div_polar F, second-order centered POLAR divergence of a vector field F = (F_r, F_theta) at the center (ghosts filled). More...

struct	pops::detail::PolarSchurReconstructKernel
	Reconstructs v^{n+theta} = B^{-1}(v^n - theta dt grad_polar phi^{n+theta}) and writes mom = rho^n v^{n+theta}. More...

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

class	pops::PolarCondensedSchurSourceStepper
	SOURCE STAGE condensed by Schur in POLAR geometry, STANDALONE (transport frozen), GENERIC over any polar fluid block that exposes the roles Density / MomentumX (radial) / MomentumY (azimuthal) (+ optional Energy). More...

Namespaces
namespace	pops

namespace	pops::detail

Detailed Description

PolarCondensedSchurSourceStepper: SOURCE STAGE condensed by Schur in POLAR geometry (ring (r, theta)).

POLAR counterpart of the Cartesian CondensedSchurSourceStepper (#126, condensed_schur_source_stepper.hpp), for the implicit source coupling potential / velocity / Lorentz of a STIFF magnetized fluid on a ring mesh. It is a STANDALONE SOURCE stage (transport frozen). WIRED into the facade since Path A step 2c: System::set_source_stage builds it when the geometry is polar and SystemStepper::run_source_stage invokes it after transport (cf. python/system.cpp) – the old "does not get wired in" mention was stale (audit 2026-06).

SEPARATE, ADDITIVE PATH (Path A step 2b). The CARTESIAN Schur stays BIT-IDENTICAL (condensed_schur_source_stepper.hpp UNTOUCHED); this header touches no existing path. It COMPOSES the bricks already on master:

PolarTensorKrylovSolver (polar_tensor_operator.hpp, #210): matrix-free BiCGStab for L_int(phi) = div(A grad phi) in POLAR METRIC, FULL non-symmetric tensor A = [[a_rr, a_rt], [a_tr, a_tt]] (RadialLine preconditioner);
LorentzEliminator (lorentz_eliminator.hpp, #118): B^{-1} closed form (2x2 rotation-dilation).

WHY THE SAME 2x2 B^{-1} AS IN CARTESIAN. The polar velocity (v_r, v_theta) is expressed in the LOCAL ORTHONORMAL frame (e_r, e_theta). The magnetic field is carried by z (B = B_z z_hat, orthogonal to the (r, theta) plane). The Lorentz force v x B is a ROTATION in the plane, INDEPENDENT of the orientation of the local orthonormal frame: (v x B)_r = +B_z v_theta, (v x B)_theta = -B_z v_r, exactly the Cartesian form (v x B)_x = +B_z v_y, (v x B)_y = -B_z v_x. Thus the SAME LorentzEliminator applies to (v_r, v_theta), and the condensed tensor A = I + c rho B^{-1} (c = theta^2 dt^2 alpha) has the same entries, simply RE-LABELED (x -> r, y -> theta): a_rr = 1 + c rho binv_11 a_rt = c rho binv_12 a_tr = c rho binv_21 a_tt = 1 + c rho binv_22 with B^{-1} = (1/det)[[1, w], [-w, 1]], w = theta dt B_z, det = 1 + w^2. This is EXACTLY the FULL non-symmetric tensor that PolarTensorKrylovSolver (#210) accepts: the cross term a_rt/a_tr IS the Lorentz rotation. The POLAR METRIC (1/r, 1/r^2) is carried by the solver operator; the stepper only assembles the cell-centered coefficients and the right-hand side.

SEQUENCE (once per source stage), polar counterpart of docs/SCHUR_CONDENSATION_DESIGN.md level 4:

ASSEMBLE A = I + c rho^n B^{-1} (4 cell-centered coefficient fields) and the condensed right-hand side rhs_polar = Lap_polar phi^n + theta dt alpha div_polar(rho^n B^{-1} v^n);
SOLVE L_int(phi) = div(A grad phi) = rhs_polar via PolarTensorKrylovSolver;
RECONSTRUCT v^{n+theta} = B^{-1}(v^n - theta dt grad_polar phi^{n+theta}), mom = rho^n v;
ENERGY (if Energy role): E^{n+1} = E^n + (1/2) rho (|v^{n+1}|^2 - |v^n|^2);
EXTRAPOLATE U^{n+theta} -> U^{n+1}: f^{n+1} = f^n + (1/theta)(f^{n+theta} - f^n);
FILL the ghosts (theta periodic, r physical) before returning.

SOLVE SIGN CONVENTION (read before any modification). The Cartesian builder (#124) writes the stage L_schur(phi) = rhs_schur with L_schur(phi) = -div(A grad phi), and rhs_schur = -Lap phi^n - theta dt alpha div(rho B^{-1} v^n). The PolarTensorKrylovSolver (#210) solves L_int(phi) = +div(A grad phi) = rhs (POSITIVE sign, cf. test_polar_tensor_elliptic_mms: the RHS there equals directly f = div(A grad phi_exact)). We have L_schur = -L_int, so solving L_schur(phi) = rhs_schur amounts to solving L_int(phi) = -rhs_schur. The POLAR right-hand side fed to the solver is therefore -rhs_schur: rhs_polar = -rhs_schur = +Lap_polar phi^n + theta dt alpha div_polar(rho B^{-1} v^n). Safeguard: c = 0 and B_z = 0 -> A = I, the solve becomes Lap_polar phi^{n+theta} = Lap_polar phi^n -> phi^{n+theta} = phi^n (up to the gauge), and the reconstruction degenerates to the explicit electrostatic push v^{n+theta} = v^n - theta dt grad_polar phi^n (case B = 0).

RIGHT-HAND SIDE DISCRETIZATION (consistent with the polar operator #210 and assemble_rhs_polar).

Lap_polar phi^n: applies apply_polar_tensor with A = I (coefficients at 1) = (1/r) d_r(r d_r phi)
- (1/r^2) d_theta^2 phi, conservative radial FV stencil + 2-point azimuthal FD (the SAME as the solve operator, guaranteeing the exact degeneracy at c=0, B_z=0).
div_polar(F), F = rho B^{-1} v^n = B^{-1}(mr, mtheta) at the center: centered FV polar divergence div F (i,j) = (1/r_i) (r_{i+1/2} - r_{i-1/2} spread by centered difference) ... discretized by div F = (1/r_i)(d_r(r F_r))_center + (1/r_i) (d_theta F_theta)_center, that is div F = (1/r_i) [ (r_{i+1} F_r(i+1) - r_{i-1} F_r(i-1)) / (2 dr) + (F_th(i,j+1) - F_th(i,j-1)) / (2 dtheta) ]. Centered, second order (polar counterpart of the Cartesian centered divergence of builder #124).

POLAR GRADIENT (reconstruction). grad_polar phi = (d_r phi, (1/r) d_theta phi), CENTERED differences: (grad phi)_r = (phi(i+1,j) - phi(i-1,j)) / (2 dr) (grad phi)_theta = (phi(i,j+1) - phi(i,j-1)) / (2 r_i dtheta). Same centered differences as the RHS divergence (term-by-term consistency at solve precision, as the Cartesian SchurReconstructKernel uses a centered gradient).

SCOPE: MULTI-RANK MPI AND MULTI-BOX (several boxes PER RANK) by AZIMUTHAL split (theta only), like PolarTensorKrylovSolver (#210). The body of step() already iterates over local_size() everywhere and goes through fill_ghosts (MPI halo exchange + physical BC, DIAGONAL CORNERS included for the cross terms of the 9-point stencil) for all fields (bz_/a_rr_/a_tt_/a_rt_/a_tr_/fr_/ft_/phi/state) : it is structurally distributed and multi-box. The elliptic solve (PolarTensorKrylovSolver) supports multi-rank AND multi-box under the theta-split constraint (each box covers the full radial range for the RadialLine preconditioner; the solver safeguard check_radial_columns enforces it – Jacobi fallback without constraint for a 2D tiling that cuts r). Single-rank / single box: BIT-IDENTICAL path. Multi-box validation: test_polar_schur_multibox (multi-box vs mono-box parity bit-close, cross terms AND 2D Jacobi tiling). Device: all kernels of the path are NAMED device-clean functors (recipe #93); CUDA-device validation = ci-full.

DEVICE. All kernels are NAMED device-clean FUNCTORS (recipe #93: no extended lambda first-instantiated cross-TU, nvcc limit #64/#97). The buffers are ALLOCATED ONCE at construction and REUSED on each step(). The Krylov solver is created per call (fixed buffers).

include/pops/coupling/schur/source/polar_condensed_schur_source_stepper.hpp File Reference¶

Classes

Namespaces

Detailed Description