Direct POLAR Poisson solver on an annulus (PolarGeometry), Phase 2a. 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/geometry/geometry.hpp>
#include <pops/mesh/storage/multifab.hpp>
#include <pops/mesh/boundary/physical_bc.hpp>
#include <pops/numerics/elliptic/poisson/poisson_fft.hpp>
#include <pops/parallel/comm.hpp>
#include <cmath>
#include <complex>
#include <concepts>
#include <stdexcept>
#include <vector>

Include dependency graph for polar_poisson_solver.hpp:

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

Go to the source code of this file.

Classes
class	pops::PolarPoissonSolver

Namespaces
namespace	pops

Concepts
concept	pops::PolarEllipticSolver

Detailed Description

Direct POLAR Poisson solver on an annulus (PolarGeometry), Phase 2a.

Solves the polar Poisson equation on the ANNULAR mesh (r, theta) of the annular polar grid (Phase 1, PolarGeometry): (1/r) d_r(r d_r phi) + (1/r^2) d_theta^2 phi = f with theta PERIODIC (the annulus covers [0, 2pi)) and a PHYSICAL BC at r_min / r_max (Dirichlet or homogeneous Neumann). The annulus EXCLUDES r = 0 (r_min > 0): NO coordinate singularity, unlike the full disk.

METHOD (FFT-in-theta + tridiagonal-in-r, robust, direct – NO multigrid). Scoping reported that the MG V-cycle can STAGNATE on the 1/r^2 operator in polar. The classic robust method: since theta is periodic with constant coefficient, an FFT in theta DECOUPLES the azimuthal modes m (the 1/r^2 d_theta^2 stencil diagonalizes under the DFT). For each mode m we then solve a 1D radial ODE (1/r) d_r(r d_r phi_m) - (lambda_theta(m)/r^2) phi_m = f_m by a tridiagonal solve (Thomas) in r. This is EXACT per mode, robust (no iterative convergence issue) and cheap (O(ntheta log ntheta) in theta + O(nr) per mode).

The 1D FFT (fft1d / dft1d_direct) is REUSED verbatim from poisson_fft.hpp (radix-2 on powers of 2, fallback to direct O(n^2) DFT otherwise -> any ntheta accepted).

SPECTRAL AZIMUTHAL EIGENVALUE: the Fourier basis diagonalizes d_theta^2 EXACTLY. The mode of DFT index m (m in [0, ntheta)) corresponds to the signed wavenumber k(m) = m if m <= ntheta/2, otherwise m - ntheta (aliasing); d_theta^2 e^{i k theta} = -k^2 e^{i k theta}, so the eigenvalue is -k(m)^2 (and NOT the 2-point stencil (2cos-2)/dtheta^2, which is only an O(dtheta^2) approximation of -k^2). Using -k(m)^2 makes the theta direction SPECTRAL (exact for band-limited data, i.e. azimuthal content carried by a small number of modes). The azimuthal term at cell (i, m) is therefore (-k(m)^2 / r_i^2) phi_hat(i, m), DIAGONAL in m (1/r_i^2 local to row i).

RADIAL DISCRETIZATION (conservative finite volumes, order 2, like assemble_rhs_polar): (1/r_i) [ r_{i+1/2}(phi_{i+1}-phi_i)/dr - r_{i-1/2}(phi_i - phi_{i-1})/dr ] / dr that is, per mode m, the tridiagonal system (in phi_hat(., m)): sub-diag a_i = (1/r_i) r_{i-1/2} / dr^2 super-diag c_i = (1/r_i) r_{i+1/2} / dr^2 diag b_i = -(a_i + c_i) - k(m)^2 / r_i^2 r_i = r_cell(i), r_{i+/-1/2} = r_face(i+1)/r_face(i) (all > 0 on the annulus).

BOUNDARY CONDITIONS IN r (Dirichlet or homogeneous Neumann, via BCRec.xlo / .xhi):

Dirichlet (value v at the face): reflection ghost phi_{-1} = 2 v - phi_0. The coefficient a_0 then folds into the diagonal (b_0 -= a_0) and 2 a_0 v moves to the right-hand side. Same at r_max with c_{nr-1}, v = xhi_val.
Homogeneous Neumann (Foextrap, zero gradient at the face): ghost phi_{-1} = phi_0 -> a_0 folds into b_0 += a_0 (the flux at the wall is zero). Same c_{nr-1} at r_max. theta stays PERIODIC (handled by the FFT, no azimuthal ghost).

MODE m = 0 + TWO Neumann boundaries: the pure radial operator has a kernel (the constant), so the tridiagonal is SINGULAR. We fix the gauge by pinning phi_hat(0, 0) = 0 (phi of zero radial mean for the constant mode), as the FFT solver pins the k = 0 mode. With at least one Dirichlet boundary the operator is invertible for every m (no gauge to fix).

EllipticSolver CONTRACT (rhs()/phi()/solve()/residual()/geom()): MODEL adapted to polar. geom() returns the PolarGeometry (not Geometry) – this solver is NOT yet wired into System.step (Phase 2b); it only composes as a standalone polar elliptic solver. ADDITIVE: no Cartesian path (geometric_mg, PoissonFFTSolver) is touched.

SCOPE: single-rank, single box covering the annulus (like PoissonFFTSolver). The FFT-in-theta + tridiag-in-r requires the full theta row AND the full radial column on one rank; the distributed case would require a parallel transpose (out of scope for Phase 2a). HARD guard (active in Release) if n_ranks() > 1 or ba.size() != 1, thrown on ALL ranks (no deadlock); solve()/residual() are local_size()==0-safe.

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

Classes

Namespaces

Concepts

Detailed Description