CondensedSchurSourceStepper: Schur-condensed SOURCE STAGE (level 4 of docs/SCHUR_CONDENSATION_DESIGN.md) for the implicit source coupling potential / velocity / Lorentz. More...

#include <pops/core/foundation/types.hpp>
#include <pops/core/state/variables.hpp>
#include <pops/coupling/schur/core/schur_condensation.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/mg/geometric_mg.hpp>
#include <pops/numerics/elliptic/linear/krylov_solver.hpp>
#include <pops/numerics/linalg/lorentz_eliminator.hpp>
#include <stdexcept>

Include dependency graph for condensed_schur_source_stepper.hpp:

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

Go to the source code of this file.

Classes
struct	pops::detail::SchurReconstructKernel
	Reconstructs v^{n+theta} = B^{-1}(v^n - theta dt grad phi^{n+theta}) and writes mom = rho^n v^{n+theta} into the state. More...

class	pops::CondensedSchurSourceStepper
	Schur-condensed SOURCE STAGE, STANDALONE (transport frozen), GENERIC over any fluid block that exposes the Density / MomentumX / MomentumY roles (+ optional Energy). More...

Namespaces
namespace	pops

namespace	pops::detail

Detailed Description

CondensedSchurSourceStepper: Schur-condensed SOURCE STAGE (level 4 of docs/SCHUR_CONDENSATION_DESIGN.md) for the implicit source coupling potential / velocity / Lorentz.

This is a STANDALONE SOURCE stage (transport frozen): it does NOT replace the System::step path and does NOT hook into it (facade wiring = PR5). It COMPOSES the bricks already on master:

ElectrostaticLorentzCondensation (schur_condensation.hpp, #124): assembles A_op = I + c rho B^{-1} (eps_x/eps_y diag, a_xy/a_yx cross) plus the condensed RHS;
TensorKrylovSolver (krylov_solver.hpp, #122): matrix-free BiCGStab preconditioned by the GeometricMG V-cycle on the SYMMETRIC part (the operator is NON self-adjoint as soon as B_z != 0);
LorentzEliminator (lorentz_eliminator.hpp, #118): closed-form B^{-1} for the reconstruction.

SEQUENCE (docs/SCHUR_CONDENSATION_DESIGN.md section 5, level 4), once per source stage:

ASSEMBLE A_op = I + c rho^n B^{-1} (c = theta^2 dt^2 alpha) and the condensed RHS rhs_schur = -Lap phi^n - theta dt alpha div(rho^n B^{-1} v^n) [#124]
SOLVE L_schur(phi^{n+theta}) = rhs_schur, L_schur = -div(A_op grad phi) [#122]
RECONSTRUCT v^{n+theta} = B^{-1}(v^n - theta dt grad phi^{n+theta}), mom = rho^n v [#118]
ENERGY (if Energy role present): see CHOICE below
EXTRAPOLATE U^{n+theta} -> U^{n+1}: see CHOICE below
FILL the ghosts (fill_boundary, MPI halos) and publish phi^{n+1} + grad phi^{n+1}.

SOLVE SIGN CONVENTION (read before any modification). The builder #124 writes the stage in the form L_schur(phi) = rhs_schur with L_schur(phi) = -div(A_op grad phi) (cf. the header of schur_condensation.hpp). The TensorKrylovSolver, on the other hand, solves L_int(phi) = rhs_kry with the poisson_operator convention L_int = div(A grad phi) - kappa phi (kappa = 0 here). So L_schur = -L_int, and solving L_schur(phi) = rhs_schur amounts to solving L_int(phi) = -rhs_schur. We thus pass rhs_kry = -rhs_schur to the solver (rhs() field negated). Safeguard: c = 0 and B_z = 0 -> A = I, the solve becomes Lap phi^{n+theta} = Lap phi^n -> phi^{n+theta} = phi^n (up to a constant), and the reconstruction degenerates into the explicit electrostatic push v^{n+theta} = v^n - theta dt grad phi^n (case B = 0). The builder already tests this safeguard on the assembly side (test C2).

EXTRAPOLATION CHOICE (step 5), DOCUMENTED. theta-stage -> n+1 by LINEAR EXTRAPOLATION: U^{n+1} = U^n + (1/theta) (U^{n+theta} - U^n). This is the form cited by docs/SCHUR_CONDENSATION_DESIGN.md (section 5, level 4, step 5). At theta = 1 (pure implicit) the stage lands DIRECTLY at n+1 and the extrapolation is the identity (factor 1/theta = 1). At theta = 1/2 (Crank-Nicolson) the factor is 2: the half-step state is extrapolated to the full step. Applied to phi and v (and thus mom = rho v); rho is FROZEN in the source (all its transport is in the hyperbolic stage), so rho^{n+1} = rho^n.

ENERGY CHOICE (step 4), DOCUMENTED. The source d_t v = -grad phi + v x Omega does WORK on the fluid via the electrostatic force -grad phi (the Lorentz rotation v x Omega does no work: |B v|^2 = |v|^2 up to a det dilation; in the absence of dissipation the kinetic energy changes only through the field). With rho frozen, the INTERNAL energy is conserved in the source; only the KINETIC energy varies. We thus update the total energy by the kinetic-energy INCREMENT: E^{n+1} = E^n + (1/2) rho^n ( |v^{n+1}|^2 - |v^n|^2 ). Consistent with a total energy E = E_internal + (1/2) rho |v|^2 where only the kinetic term moves under the source. Applied ONLY if the Energy role is present; otherwise the energy is ignored (isothermal model / no energy equation), like the builder #124 which already ignores Energy.

DEVICE / MPI (docs/GPU_RUNTIME_PORT.md). All kernels of this stage are NAMED FUNCTORS device-clean (no extended lambda first-instantiated cross-TU, nvcc limit #64/#97). The MultiFab buffers (operator A_op, RHS, phi^n, grad phi, v^n) are ALLOCATED ONCE at construction and REUSED on each step() call (the #124 design notes the cost of a reallocation per step). All loops iterate over local_size(): a rank WITHOUT a box (local_size()==0, MPI partition) executes no kernel and never derefs fab(0). The Krylov solve is COLLECTIVE (dot/all_reduce over all ranks, including empty ones): no deadlock.

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

Classes

Namespaces

Detailed Description