test_momentum_convective_candidates.c File Reference

A4a focused convective-candidate study (states A-C only). More...

#include "test_support.h"
#include "momentumsolvers.h"
#include "rhs.h"
#include "setup.h"
#include "Boundaries.h"
#include <petscblaslapack.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

Include dependency graph for test_momentum_convective_candidates.c:

Go to the source code of this file.

Data Structures
struct	DofMap
struct	SeamDiagnostics
struct	GlobalVecStats
struct	StableCFLResult

Macros
#define	HAVE_I(ii)   InGhostRange((ii), info.gxs, info.gxm)
#define	HAVE_J(jj)   InGhostRange((jj), info.gys, info.gym)
#define	HAVE_K(kk)   InGhostRange((kk), info.gzs, info.gzm)
#define	HAVE_I(ii)   InGhostRange((ii), info.gxs, info.gxm)
#define	HAVE_J(jj)   InGhostRange((jj), info.gys, info.gym)
#define	HAVE_K(kk)   InGhostRange((kk), info.gzs, info.gzm)
#define	ROWSUM(cmp)

Enumerations
enum	StableCFLStatus { STABLE_CFL_NONE , STABLE_CFL_FINITE , STABLE_CFL_EXCEEDS_SCAN }
enum	CandState {   STATE_A , STATE_A_X , STATE_A_Y , STATE_B ,   STATE_C }
enum	PMetric { METRIC_RHO , METRIC_NORM }
enum	ResidualStage { STAGE_UCAT = 0 , STAGE_CART_RESID = 1 , STAGE_RCT = 2 , STAGE_FINAL = 3 }

Functions
static PetscInt	PeriodicRepCount (PetscInt npts)
	Returns the number of independent periodic representatives in one direction.
static PetscInt	DofMapExpectedCount (DMDALocalInfo info)
	Counts all independent component-staggered representatives used by ComputeRHS.
static PetscErrorCode	DofMapBuild (UserCtx *user, DofMap *map)
	Builds the serial periodic independent face-DOF map used by dense Jacobians.
static PetscErrorCode	DofMapBuildOwned (UserCtx *user, DofMap *map)
	Builds the rank-owned periodic independent face-DOF map for MPI checks.
static PetscErrorCode	DofMapDestroy (DofMap *map)
	Releases storage owned by an active-DOF map.
static PetscReal	CmpGet (Cmpnts c, PetscInt comp)
static PetscErrorCode	EvalConvResidual (UserCtx *user, Vec Ucont_in, Vec Rhs, const DofMap *map, PetscReal *Ract)
static PetscErrorCode	PerturbDof (UserCtx *user, Vec Ucont, const DofMap *map, PetscInt m, PetscReal delta)
static PetscErrorCode	GetDof (UserCtx *user, Vec Ucont, const DofMap *map, PetscInt m, PetscReal *val)
	Reads one active contravariant component from a global vector.
static PetscReal	PeriodicCellAngle (PetscInt idx, PetscInt npts)
	Returns the cell-centered periodic angle using duplicated endpoint planes.
static PetscReal	PeriodicFaceAngle (PetscInt idx, PetscInt npts)
	Returns a face-representative periodic angle for component-staggered Ucont.
static Cmpnts	AnalyticVelocity (CandState st, PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
	Evaluates one of the three analytic Cartesian candidate states.
static Cmpnts	DirectUcontVelocity (CandState st, PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
	Evaluates the declared direct component-staggered State B Ucont field.
static PetscReal	CmpDiffInf (Cmpnts a, Cmpnts b)
	Computes the componentwise infinity norm of the difference between two vectors.
static PetscErrorCode	ComputeDeclaredSeamMismatch (CandState st, DMDALocalInfo info, PetscReal seam[3])
	Computes analytic periodic seam mismatches for each coordinate direction.
static PetscBool	InGhostRange (PetscInt idx, PetscInt lo, PetscInt n)
	Reports whether a global index is present in a rank's local ghosted range.
static PetscErrorCode	ComputeLocalDuplicateMismatch (UserCtx *user, Vec local, PetscReal seam[3])
	Computes duplicate-plane mismatch in a local vector view.
static PetscErrorCode	ComputeLocalOuterGhostMismatch (UserCtx *user, Vec local, PetscReal seam[3])
	Computes outer periodic ghost mismatch for local Ucat.
static PetscErrorCode	ConfigureCandidateFixture (SimCtx *simCtx, UserCtx *user)
	Configures the minimal context for centered inviscid periodic convection tests.
static PetscErrorCode	SetUcatField (UserCtx *user, CandState st)
static PetscErrorCode	SetDirectUcontField (UserCtx *user, CandState st)
	Sets the direct State B component-staggered Ucont field on owned entries.
static PetscErrorCode	ComputeMaxDiscreteDivergence (UserCtx *user, PetscReal *maxdiv)
	Computes max discrete divergence of the current local Ucont field.
static PetscErrorCode	BuildBaseState (UserCtx *user, CandState st, Vec Ubase, PetscReal *repeat_inf, PetscReal *maxdiv, SeamDiagnostics *seam)
static PetscErrorCode	ComputeMaxGradientContribution (UserCtx *user, PetscReal *gradmax)
	Computes the global maximum Cartesian velocity-gradient row-sum used by Candidate D.
static PetscErrorCode	DenseSpectralRadius (const PetscReal *A, PetscInt n, PetscReal *rho, PetscReal *maxRealPart)
static PetscErrorCode	DenseMaxRealRightEigenpair (const PetscReal *A, PetscInt n, PetscReal *lamr, PetscReal *lami, PetscReal *vr_out, PetscReal *vi_out)
	Extracts the right eigenpair whose eigenvalue has largest real part.
static PetscErrorCode	DenseRKPolynomialSpectralRadius (const PetscReal *J, PetscInt n, PetscReal dtau, PetscReal *rho)
	Computes the spectral radius of the RK polynomial by applying it to eig(J).
static PetscErrorCode	DenseSigmaMax (const PetscReal *A, PetscInt n, PetscReal *smax, PetscReal *v1)
static PetscReal	DenseNonNormality (const PetscReal *A, PetscInt n)
static PetscReal	DenseSkewnessDefect (const PetscReal *A, PetscInt n)
	Computes the normalized Frobenius defect from skew symmetry.
static PetscErrorCode	PrintSpectrumSummary (const char *name, const PetscReal *J, PetscInt n)
	Prints eigenvalue and norm summary for one dense operator.
static void	DenseMatVec (const PetscReal *J, const PetscReal *x, PetscReal *y, PetscInt n)
	Applies a dense column-major matrix to an active-space vector.
static void	FillAuditDirection (const DofMap *map, PetscInt kind, PetscReal *x)
	Fills one deterministic active-space vector used by dense/matrix-free checks.
static PetscErrorCode	AddActiveVector (UserCtx *user, Vec U, const DofMap *map, const PetscReal *x, PetscReal scale)
	Adds a dense active-space vector into a global contravariant vector.
static PetscErrorCode	DenseShiftIdentity (const PetscReal *A, PetscInt n, PetscReal shift, PetscReal *B)
	Copies a dense matrix and adds a scalar shift to its diagonal.
static void	MatMul (const PetscReal *A, const PetscReal *B, PetscReal *out, PetscInt n)
static PetscErrorCode	RKPolynomial (const PetscReal *J, PetscReal dtau, PetscInt n, PetscReal *P)
static PetscErrorCode	AmplificationMetric (const PetscReal *J, PetscInt n, PetscReal dtau, PMetric which, PetscReal *metric)
	Evaluates either spectral-radius or 2-norm amplification for one pseudo-time step.
static PetscErrorCode	StableCFL (const PetscReal *J, PetscInt n, PetscReal lam, PMetric which, StableCFLResult *result)
static const char *	StableCFLStatusText (StableCFLResult r)
	Returns human-readable text for a stable-CFL search result.
static PetscErrorCode	PrintStableCFLLine (const char *candidate, StableCFLResult eig, StableCFLResult norm)
	Prints one candidate's eigenvalue and norm stable-CFL statuses.
static void	ApplyP (const PetscReal *P, const PetscReal *x, PetscReal *out, PetscInt n)
static PetscErrorCode	PrintFrozenAmplificationTable (const char *title, const PetscReal *J, PetscInt n, const PetscReal lams[3], const char *cn[3])
	Prints frozen RK amplification tables for the supplied operator and candidates.
static PetscReal	VecNorm2Array (const PetscReal *x, PetscInt n)
	Computes the Euclidean norm of a dense vector.
static PetscReal	VecNormInfArray (const PetscReal *x, PetscInt n)
	Computes the infinity norm of a dense vector.
static PetscReal	DenseFrobenius (const PetscReal *A, PetscInt n)
	Computes the Frobenius norm of a dense column-major matrix.
static PetscBool	DofTouchesPeriodicRepresentative (const DofMap *map, PetscInt m, DMDALocalInfo info)
	True if an active representative lies on a plane adjacent to periodic duplicates.
static PetscErrorCode	PrintPeriodicSpaceAudit (UserCtx *user, const DofMap *map)
	Prints the component-wise periodic storage count actually used by the active map.
static PetscErrorCode	PrintBlockAndSymmetricLocalization (UserCtx *user, const DofMap *map, const PetscReal *J)
	Prints 3x3 component block norms and localized symmetric rows.
static void	FillStaggeredFourierBasis (const DofMap *map, PetscInt wx, PetscInt wy, PetscInt wz, PetscReal *Q)
	Fills the six real basis vectors for one staggered same-wavevector subspace.
static PetscErrorCode	StaggeredFourierSymbol (const DofMap *map, const PetscReal *J, PetscInt wx, PetscInt wy, PetscInt wz, PetscReal A6[36], PetscReal *leak)
	Builds the 6x6 projected real symbol and leakage for one wavevector.
static PetscErrorCode	SymbolEigenSummary (const PetscReal A6[36], PetscReal *maxre, PetscReal wr_out[6], PetscReal wi_out[6])
	Eigenvalue summary for a 6x6 real symbol.
static PetscErrorCode	PrintFourierChecks (const DofMap *map, const PetscReal *J)
	Projects production J on full staggered six-dimensional Fourier subspaces.
static PetscErrorCode	TraceStateAResidualStages (UserCtx *user, Vec Ubase, Vec Rhs, const DofMap *map, const PetscReal *J)
	Trace observable production residual stages for one deterministic perturbation.
static PetscReal	DenseRelativeDiff (const PetscReal *A, const PetscReal *B, PetscInt n, PetscReal denom_ref)
	Computes a normalized Frobenius difference between two dense matrices.
static PetscErrorCode	ExtractActiveVector (UserCtx *user, Vec U, const DofMap *map, PetscReal *x)
	Extracts active-space entries from a global contravariant vector.
static PetscReal	DofWeight (const DofMap *map, PetscInt m)
	Returns a deterministic checksum weight for an active DOF.
static PetscErrorCode	ActiveStats (const DofMap *map, const PetscReal *x, GlobalVecStats *stats)
	Computes global active-vector norms and checksum.
static PetscErrorCode	FillDeterministicDirection (UserCtx *user, Vec V, const DofMap *map, PetscReal *x)
	Fills a globally normalized deterministic active-space perturbation direction.
static PetscErrorCode	FourStage (UserCtx *user, Vec U0full, PetscReal dtau, Vec Rhs, const DofMap *map, PetscReal *Rscratch, Vec Uwork, Vec Uout)
static PetscErrorCode	SetAnchoredStage (UserCtx *user, Vec U0full, PetscReal scale, const DofMap *map, const PetscReal *Ract, Vec Ustage)
	Forms one anchored RK stage state from the base state and active residual.
static PetscErrorCode	BuildStageStates (UserCtx *user, Vec U0full, PetscReal dtau, Vec Rhs, const DofMap *map, PetscReal *Rscratch, Vec Y1, Vec Y2, Vec Y3)
	Builds the first three anchored RK stage states for a base vector.
static PetscErrorCode	BuildFDJacobian (UserCtx *user, Vec Ucenter, PetscReal epsrel, Vec Rhs, const DofMap *map, PetscReal *Rp, PetscReal *Rm, Vec Uwork, PetscReal *J)
	Builds a centered finite-difference Jacobian of the production convective residual.
static PetscErrorCode	CheckDenseJacobianAction (UserCtx *user, Vec Ubase, Vec Rhs, Vec Uwork, const DofMap *map, const PetscReal *J, PetscReal epsrel)
	Verifies that the assembled dense Jacobian has the same action as production FD Jv.
static PetscErrorCode	MatrixFreeJv (UserCtx *user, Vec Ubase, Vec Rhs, const DofMap *map, const PetscReal *v, PetscReal eps, PetscReal *jv)
	Matrix-free production Jacobian-vector product on the active space.
static PetscReal	DenseAdditivityError (const PetscReal *A, const PetscReal *B, const PetscReal *C, PetscInt n)
	Normalized Frobenius norm of A - (B+C).
static PetscReal	DenseCommutatorNorm (const PetscReal *A, const PetscReal *B, PetscInt n)
	Computes C = A*B - B*A and returns its normalized Frobenius norm.
static PetscErrorCode	RunStateADirectionalMechanismAudit (const PetscReal *Jxy, UserCtx *user_xy, Vec Ubase_xy, Vec Rhs_xy, const DofMap *map_xy, PetscReal epsrel)
	Builds Jx, Jy, Jxy and prints additivity/commutator/eigenpair diagnostics.
static PetscErrorCode	ExtractLocalActiveVector (UserCtx *user, Vec local, const DofMap *map, PetscReal *x)
	Extract active-space entries from a local vector.
static PetscErrorCode	MapCartesianResidualToRct (UserCtx *user, Vec Rc, Vec Rct)
	Mirrors ComputeRHS's Cartesian residual -> contravariant local Rct mapping.
static PetscErrorCode	AverageRctToFinalActive (UserCtx *user, Vec Rct, const DofMap *map, PetscReal *out)
	Mirrors ComputeRHS's final Rct face averaging and active cleanup for P=0/no body force.
static PetscErrorCode	EvalResidualStage (UserCtx *user, Vec Ucont_in, Vec Rhs, const DofMap *map, ResidualStage stage, PetscReal *out)
	Evaluates one observable/mirrored residual-path stage for a Ucont input.
static PetscErrorCode	BuildStageJacobian (UserCtx *user, Vec Ucenter, PetscReal epsrel, Vec Rhs, const DofMap *map, ResidualStage stage, PetscReal *Rp, PetscReal *Rm, Vec Uwork, PetscReal *J)
	Builds a finite-difference Jacobian for one residual-path stage.
static PetscErrorCode	RunStateAResidualPathIsolation (UserCtx *user, Vec Ubase, Vec Rhs, Vec Uwork, const DofMap *map, PetscReal epsrel)
	Prints the stage where positive-real spectrum first appears.
static PetscErrorCode	BuildStageTangent (const PetscReal *J0, const PetscReal *J1, const PetscReal *J2, const PetscReal *J3, PetscReal dtau, PetscInt n, PetscReal *T4)
	Builds the exact four-stage tangent from stage-dependent Jacobians.
static PetscErrorCode	BuildPhiJacobian (UserCtx *user, Vec Ucenter, PetscReal dtau, Vec Rhs, const DofMap *map, PetscReal *Rscratch, Vec Upert, Vec Ustage, Vec PhiP, Vec PhiM, PetscReal *xp, PetscReal *xm, PetscReal epsrel, PetscReal *JPhi)
	Builds a finite-difference Jacobian of the complete nonlinear four-stage map.
static PetscErrorCode	RunRKTangentDiagnostics (UserCtx *user, CandState st, Vec Ubase, Vec Rhs, const DofMap *map, PetscReal epsrel, const PetscReal lams[3], const char *cn[3], const PetscReal *J0)
	Runs stage-dependent RK tangent and direct nonlinear perturbation diagnostics.
static PetscErrorCode	RunState (CandState st, const char *name)
static PetscErrorCode	TestStateA (void)
	Runs the State A candidate harness.
static PetscErrorCode	TestStateB (void)
	Runs the State B candidate harness.
static PetscErrorCode	TestStateC (void)
	Runs the State C candidate harness.
static PetscErrorCode	RunStateAGridAuditOne (PetscInt N)
	Runs one State A grid-size audit case.
static PetscErrorCode	RunStateASplitOne (CandState st, const char *label)
	Runs one State A transport-direction split.
static PetscErrorCode	TestStateAGridAudit (void)
	Runs the State A grid-dependence and active-space audit.
static PetscErrorCode	WriteStateADecompReference (GlobalVecStats r0, GlobalVecStats jv, GlobalVecStats phi, const MomStabilityReport *rep)
	Writes the one-rank State A matrix-free decomposition reference.
static PetscErrorCode	ReadAndCompareStateADecompReference (GlobalVecStats r0, GlobalVecStats jv, GlobalVecStats phi, const MomStabilityReport *rep)
	Compares a distributed State A matrix-free check against the one-rank reference.
static PetscErrorCode	StateADecompDirectionalCheck (UserCtx *user, Vec Ubase, Vec Rhs, PetscReal lcC, const MomStabilityReport *rep)
	Runs State A matrix-free residual, Jv, and four-stage MPI decomposition checks.
static PetscErrorCode	RunDecompBaseline (CandState st, const char *name, PetscReal expB, PetscReal expC, PetscReal expD)
	Runs one decomp baseline and compares scalar estimates with regenerated references.
static PetscErrorCode	TestDecompA (void)
	Runs the State A decomposition baseline.
static PetscErrorCode	TestDecompB (void)
	Runs the State B decomposition baseline.
static PetscErrorCode	TestDecompC (void)
	Runs the State C decomposition baseline.
static PetscErrorCode	ConfigureMPIReferenceOptions (void)
	Reads optional paired-run MPI reference path and token.
int	main (int argc, char **argv)
	PETSc entry point for the focused convective-candidate harness.

Variables
static char	g_ref_path [PETSC_MAX_PATH_LEN] = ""
static char	g_ref_token [128] = ""
static PetscBool	g_ref_path_set = PETSC_FALSE
static PetscBool	g_ref_token_set = PETSC_FALSE

Detailed Description

A4a focused convective-candidate study (states A-C only).

Builds the finite-difference Jacobian of the ACTUAL production convection-only residual (ComputeRHS with inviscid, P=0, centered, periodic, single block) on a tiny periodic Cartesian grid, then compares the B/C/D estimator candidates against rho(J), sigma_max(J), the exact 4-stage RK matrix polynomial P(z)=1+z+z^2/2+z^3/6+z^4/24, and a direct anchored 4-stage perturbation cross-check. Shadow-only: changes no production default.

States: A uniform divergence-free; B nonzero discrete divergence; C divergence-free shear.

Definition in file test_momentum_convective_candidates.c.

---

Data Structure Documentation

DofMap

struct DofMap

Definition at line 27 of file test_momentum_convective_candidates.c.

Data Fields
PetscInt	n
PetscInt	expected_n
PetscInt *	comp
PetscInt *	ci
PetscInt *	cj
PetscInt *	ck

SeamDiagnostics

struct SeamDiagnostics

Definition at line 28 of file test_momentum_convective_candidates.c.

Data Fields
PetscReal	declared[3]
PetscReal	ucat_global[3]
PetscReal	ucat_ghost[3]
PetscReal	ucont_global[3]

GlobalVecStats

struct GlobalVecStats

Definition at line 34 of file test_momentum_convective_candidates.c.

Data Fields
PetscReal	n2
PetscReal	ninf
PetscReal	checksum

StableCFLResult

struct StableCFLResult

Definition at line 48 of file test_momentum_convective_candidates.c.

Data Fields
StableCFLStatus	status
PetscReal	cfl

Macro Definition Documentation

HAVE_I [1/2]

#define HAVE_I

(

ii

)

InGhostRange((ii), info.gxs, info.gxm)

HAVE_J [1/2]

#define HAVE_J

(

jj

)

InGhostRange((jj), info.gys, info.gym)

HAVE_K [1/2]

#define HAVE_K

(

kk

)

InGhostRange((kk), info.gzs, info.gzm)

HAVE_I [2/2]

#define HAVE_I

(

ii

)

InGhostRange((ii), info.gxs, info.gxm)

HAVE_J [2/2]

#define HAVE_J

(

jj

)

InGhostRange((jj), info.gys, info.gym)

HAVE_K [2/2]

#define HAVE_K

(

kk

)

InGhostRange((kk), info.gzs, info.gzm)

ROWSUM

#define ROWSUM

(

cmp

)

Value:

    ( \
    PetscAbsReal(Ajc*(C.x*duc.cmp + E.x*due.cmp + Z.x*duz.cmp)) + \
    PetscAbsReal(Ajc*(C.y*duc.cmp + E.y*due.cmp + Z.y*duz.cmp)) + \
    PetscAbsReal(Ajc*(C.z*duc.cmp + E.z*due.cmp + Z.z*duz.cmp)) )

Enumeration Type Documentation

StableCFLStatus

enum StableCFLStatus

Enumerator
STABLE_CFL_NONE
STABLE_CFL_FINITE
STABLE_CFL_EXCEEDS_SCAN

Definition at line 41 of file test_momentum_convective_candidates.c.

             {
    STABLE_CFL_NONE,
    STABLE_CFL_FINITE,
 
    STABLE_CFL_EXCEEDS_SCAN
} StableCFLStatus;

CandState

enum CandState

Enumerator
STATE_A
STATE_A_X
STATE_A_Y
STATE_B
STATE_C

Definition at line 179 of file test_momentum_convective_candidates.c.

{ STATE_A, STATE_A_X, STATE_A_Y, STATE_B, STATE_C } CandState;

PMetric

enum PMetric

Enumerator
METRIC_RHO
METRIC_NORM

Definition at line 801 of file test_momentum_convective_candidates.c.

{ METRIC_RHO, METRIC_NORM } PMetric;

ResidualStage

enum ResidualStage

Enumerator
STAGE_UCAT
STAGE_CART_RESID
STAGE_RCT
STAGE_FINAL

Definition at line 1650 of file test_momentum_convective_candidates.c.

             {
    STAGE_UCAT = 0,
    STAGE_CART_RESID = 1,
    STAGE_RCT = 2,
    STAGE_FINAL = 3
} ResidualStage;

Function Documentation

PeriodicRepCount()

static PetscInt PeriodicRepCount ( PetscInt  npts )

inlinestatic

Returns the number of independent periodic representatives in one direction.

Definition at line 56 of file test_momentum_convective_candidates.c.

{ return npts - 2; }

Here is the caller graph for this function:

DofMapExpectedCount()

static PetscInt DofMapExpectedCount ( DMDALocalInfo  info )

inlinestatic

Counts all independent component-staggered representatives used by ComputeRHS.

Definition at line 61 of file test_momentum_convective_candidates.c.

{
    return 3 * PeriodicRepCount(info.mx) * PeriodicRepCount(info.my) * PeriodicRepCount(info.mz);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DofMapBuild()

static PetscErrorCode DofMapBuild ( UserCtx *  user, DofMap *  map  )

static

Builds the serial periodic independent face-DOF map used by dense Jacobians.

Production periodic synchronization copies global plane 0 from mx-2 and plane mx-1 from 1 (and analogously in y/z), so representatives 1..m-2 contain each independent component-staggered Ucont face DOF exactly once. Perturbations only touch these reps; EvalConvResidual() then calls SynchronizePeriodicStaggeredFields() to update duplicates.

Definition at line 74 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    const PetscInt lxs = 1, lxe = mx-1, lys = 1, lye = my-1, lzs = 1, lze = mz-1;
    PetscInt cnt = 0;
    PetscFunctionBeginUser;
    map->expected_n = DofMapExpectedCount(info);
    map->n = (lxe-lxs)*(lye-lys)*(lze-lzs)*3;
    PetscCall(PetscMalloc4(map->n, &map->comp, map->n, &map->ci, map->n, &map->cj, map->n, &map->ck));
    for (PetscInt k = lzs; k < lze; k++)
        for (PetscInt j = lys; j < lye; j++)
            for (PetscInt i = lxs; i < lxe; i++)
                for (PetscInt c = 0; c < 3; c++) {
                    map->comp[cnt] = c; map->ci[cnt] = i; map->cj[cnt] = j; map->ck[cnt] = k; cnt++;
                }
    PetscCheck(map->n == map->expected_n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ,
               "periodic independent DOF count mismatch: got %" PetscInt_FMT ", expected %" PetscInt_FMT,
               map->n, map->expected_n);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DofMapBuildOwned()

static PetscErrorCode DofMapBuildOwned ( UserCtx *  user, DofMap *  map  )

static

Builds the rank-owned periodic independent face-DOF map for MPI checks.

Definition at line 99 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    const PetscInt xs = info.xs, xe = info.xs + info.xm;
    const PetscInt ys = info.ys, ye = info.ys + info.ym;
    const PetscInt zs = info.zs, ze = info.zs + info.zm;
    const PetscInt lxs = PetscMax(xs, 1), lxe = PetscMin(xe, mx-1);
    const PetscInt lys = PetscMax(ys, 1), lye = PetscMin(ye, my-1);
    const PetscInt lzs = PetscMax(zs, 1), lze = PetscMin(ze, mz-1);
    PetscInt cnt = 0;
    PetscFunctionBeginUser;
    map->expected_n = DofMapExpectedCount(info);
    map->n = PetscMax(0,lxe-lxs)*PetscMax(0,lye-lys)*PetscMax(0,lze-lzs)*3;
    PetscCall(PetscMalloc4(map->n, &map->comp, map->n, &map->ci, map->n, &map->cj, map->n, &map->ck));
    for (PetscInt k = lzs; k < lze; k++)
        for (PetscInt j = lys; j < lye; j++)
            for (PetscInt i = lxs; i < lxe; i++)
                for (PetscInt c = 0; c < 3; c++) {
                    map->comp[cnt] = c; map->ci[cnt] = i; map->cj[cnt] = j; map->ck[cnt] = k; cnt++;
                }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DofMapDestroy()

static PetscErrorCode DofMapDestroy ( DofMap *  map )

static

Releases storage owned by an active-DOF map.

Definition at line 126 of file test_momentum_convective_candidates.c.

{ PetscFunctionBeginUser; PetscCall(PetscFree4(map->comp, map->ci, map->cj, map->ck)); PetscFunctionReturn(0); }

Here is the caller graph for this function:

CmpGet()

static PetscReal CmpGet ( Cmpnts  c, PetscInt  comp  )

inlinestatic

Definition at line 130 of file test_momentum_convective_candidates.c.

{ const PetscReal *p = (const PetscReal*)&c; return p[comp]; }

Here is the caller graph for this function:

EvalConvResidual()

static PetscErrorCode EvalConvResidual ( UserCtx *  user, Vec  Ucont_in, Vec  Rhs, const DofMap *  map, PetscReal *  Ract  )

static

Definition at line 135 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***r;
    PetscFunctionBeginUser;
    PetscCall(VecCopy(Ucont_in, user->Ucont));
    {
        const char *fld[] = {"Ucont"};
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fld));   /* global->local + periodic */
    }
    PetscCall(ComputeRHS(user, Rhs));                                  /* Contra2Cart + Convection + mapping */
    PetscCall(DMDAVecGetArrayRead(user->fda, Rhs, &r));
    for (PetscInt m = 0; m < map->n; m++)
        Ract[m] = CmpGet(r[map->ck[m]][map->cj[m]][map->ci[m]], map->comp[m]);
    PetscCall(DMDAVecRestoreArrayRead(user->fda, Rhs, &r));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PerturbDof()

static PetscErrorCode PerturbDof ( UserCtx *  user, Vec  Ucont, const DofMap *  map, PetscInt  m, PetscReal  delta  )

static

Definition at line 153 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***a;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArray(user->fda, Ucont, &a));
    { PetscReal *p = (PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]]; p[map->comp[m]] += delta; }
    PetscCall(DMDAVecRestoreArray(user->fda, Ucont, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

GetDof()

static PetscErrorCode GetDof ( UserCtx *  user, Vec  Ucont, const DofMap *  map, PetscInt  m, PetscReal *  val  )

static

Reads one active contravariant component from a global vector.

Definition at line 166 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***a;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArray(user->fda, Ucont, &a));
    { const PetscReal *p = (const PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]]; *val = p[map->comp[m]]; }
    PetscCall(DMDAVecRestoreArray(user->fda, Ucont, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PeriodicCellAngle()

static PetscReal PeriodicCellAngle ( PetscInt  idx, PetscInt  npts  )

inlinestatic

Returns the cell-centered periodic angle using duplicated endpoint planes.

Definition at line 184 of file test_momentum_convective_candidates.c.

{
    const PetscInt nuniq = npts - 1;
    const PetscInt ip = (idx == npts - 1) ? 0 : idx;
    return 2.0*PETSC_PI*((PetscReal)ip)/((PetscReal)nuniq);
}

Here is the caller graph for this function:

PeriodicFaceAngle()

static PetscReal PeriodicFaceAngle ( PetscInt  idx, PetscInt  npts  )

inlinestatic

Returns a face-representative periodic angle for component-staggered Ucont.

Definition at line 194 of file test_momentum_convective_candidates.c.

{
    const PetscInt nuniq = PeriodicRepCount(npts);
    PetscInt ip = idx - 1;
    if (idx == 0) ip = nuniq - 1;
    else if (idx == npts - 1) ip = 0;
    return 2.0*PETSC_PI*((PetscReal)ip)/((PetscReal)nuniq);
}

Here is the call graph for this function:

Here is the caller graph for this function:

AnalyticVelocity()

static Cmpnts AnalyticVelocity ( CandState  st, PetscInt  i, PetscInt  j, PetscInt  k, PetscInt  mx, PetscInt  my, PetscInt  mz  )

inlinestatic

Evaluates one of the three analytic Cartesian candidate states.

Definition at line 206 of file test_momentum_convective_candidates.c.

{
    Cmpnts v;
    const PetscReal x = PeriodicCellAngle(i, mx);
    const PetscReal y = PeriodicCellAngle(j, my);
    (void)k; (void)mz;
    if (st == STATE_A) {
        v.x = 0.7; v.y = -0.4; v.z = 0.0;
    } else if (st == STATE_A_X) {
        v.x = 0.7; v.y = 0.0; v.z = 0.0;
    } else if (st == STATE_A_Y) {
        v.x = 0.0; v.y = -0.4; v.z = 0.0;
    } else if (st == STATE_B) {
        (void)x; v.x = 0.0; v.y = 0.0; v.z = 0.0;
    } else {
        v.x = 1.0 + 0.5*PetscSinReal(y); v.y = 0.0; v.z = 0.0;
    }
    return v;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DirectUcontVelocity()

static Cmpnts DirectUcontVelocity ( CandState  st, PetscInt  i, PetscInt  j, PetscInt  k, PetscInt  mx, PetscInt  my, PetscInt  mz  )

inlinestatic

Evaluates the declared direct component-staggered State B Ucont field.

Definition at line 230 of file test_momentum_convective_candidates.c.

{
    Cmpnts v = {0.0, 0.0, 0.0};
    (void)j; (void)k; (void)my; (void)mz;
    if (st == STATE_B) v.x = PetscSinReal(PeriodicFaceAngle(i, mx));
    return v;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CmpDiffInf()

static PetscReal CmpDiffInf ( Cmpnts  a, Cmpnts  b  )

inlinestatic

Computes the componentwise infinity norm of the difference between two vectors.

Definition at line 242 of file test_momentum_convective_candidates.c.

{
    return PetscMax(PetscAbsReal(a.x-b.x), PetscMax(PetscAbsReal(a.y-b.y), PetscAbsReal(a.z-b.z)));
}

Here is the caller graph for this function:

ComputeDeclaredSeamMismatch()

static PetscErrorCode ComputeDeclaredSeamMismatch ( CandState  st, DMDALocalInfo  info, PetscReal  seam[3]  )

static

Computes analytic periodic seam mismatches for each coordinate direction.

Definition at line 250 of file test_momentum_convective_candidates.c.

{
    PetscReal loc[3] = {0.0, 0.0, 0.0}, glo[3];
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    if (st == STATE_B) {
        for (PetscInt k = 1; k < mz-1; k++)
            for (PetscInt j = 1; j < my-1; j++) {
                loc[0] = PetscMax(loc[0], CmpDiffInf(DirectUcontVelocity(st, 0, j, k, mx, my, mz),
                                                     DirectUcontVelocity(st, mx-2, j, k, mx, my, mz)));
                loc[0] = PetscMax(loc[0], CmpDiffInf(DirectUcontVelocity(st, mx-1, j, k, mx, my, mz),
                                                     DirectUcontVelocity(st, 1, j, k, mx, my, mz)));
            }
        for (PetscInt k = 1; k < mz-1; k++)
            for (PetscInt i = 1; i < mx-1; i++) {
                loc[1] = PetscMax(loc[1], CmpDiffInf(DirectUcontVelocity(st, i, 0, k, mx, my, mz),
                                                     DirectUcontVelocity(st, i, my-2, k, mx, my, mz)));
                loc[1] = PetscMax(loc[1], CmpDiffInf(DirectUcontVelocity(st, i, my-1, k, mx, my, mz),
                                                     DirectUcontVelocity(st, i, 1, k, mx, my, mz)));
            }
        for (PetscInt j = 1; j < my-1; j++)
            for (PetscInt i = 1; i < mx-1; i++) {
                loc[2] = PetscMax(loc[2], CmpDiffInf(DirectUcontVelocity(st, i, j, 0, mx, my, mz),
                                                     DirectUcontVelocity(st, i, j, mz-2, mx, my, mz)));
                loc[2] = PetscMax(loc[2], CmpDiffInf(DirectUcontVelocity(st, i, j, mz-1, mx, my, mz),
                                                     DirectUcontVelocity(st, i, j, 1, mx, my, mz)));
            }
    } else {
        for (PetscInt k = 0; k < mz; k++)
            for (PetscInt j = 0; j < my; j++)
                loc[0] = PetscMax(loc[0], CmpDiffInf(AnalyticVelocity(st, 0, j, k, mx, my, mz),
                                                     AnalyticVelocity(st, mx-1, j, k, mx, my, mz)));
        for (PetscInt k = 0; k < mz; k++)
            for (PetscInt i = 0; i < mx; i++)
                loc[1] = PetscMax(loc[1], CmpDiffInf(AnalyticVelocity(st, i, 0, k, mx, my, mz),
                                                     AnalyticVelocity(st, i, my-1, k, mx, my, mz)));
        for (PetscInt j = 0; j < my; j++)
            for (PetscInt i = 0; i < mx; i++)
                loc[2] = PetscMax(loc[2], CmpDiffInf(AnalyticVelocity(st, i, j, 0, mx, my, mz),
                                                     AnalyticVelocity(st, i, j, mz-1, mx, my, mz)));
    }
    PetscCallMPI(MPI_Allreduce(loc, glo, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    seam[0] = glo[0]; seam[1] = glo[1]; seam[2] = glo[2];
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InGhostRange()

static PetscBool InGhostRange ( PetscInt  idx, PetscInt  lo, PetscInt  n  )

inlinestatic

Reports whether a global index is present in a rank's local ghosted range.

Definition at line 299 of file test_momentum_convective_candidates.c.

{ return (PetscBool)(idx >= lo && idx < lo + n); }

ComputeLocalDuplicateMismatch()

static PetscErrorCode ComputeLocalDuplicateMismatch ( UserCtx *  user, Vec  local, PetscReal  seam[3]  )

static

Computes duplicate-plane mismatch in a local vector view.

Definition at line 305 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***a;
    PetscReal loc[3] = {0.0, 0.0, 0.0}, glo[3];
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->fda, local, &a));
#define HAVE_I(ii) InGhostRange((ii), info.gxs, info.gxm)
#define HAVE_J(jj) InGhostRange((jj), info.gys, info.gym)
#define HAVE_K(kk) InGhostRange((kk), info.gzs, info.gzm)
    if (HAVE_I(0) && HAVE_I(mx-2))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
                loc[0] = PetscMax(loc[0], CmpDiffInf(a[k][j][0], a[k][j][mx-2]));
    if (HAVE_I(mx-1) && HAVE_I(1))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
                loc[0] = PetscMax(loc[0], CmpDiffInf(a[k][j][mx-1], a[k][j][1]));
    if (HAVE_J(0) && HAVE_J(my-2))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[1] = PetscMax(loc[1], CmpDiffInf(a[k][0][i], a[k][my-2][i]));
    if (HAVE_J(my-1) && HAVE_J(1))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[1] = PetscMax(loc[1], CmpDiffInf(a[k][my-1][i], a[k][1][i]));
    if (HAVE_K(0) && HAVE_K(mz-2))
        for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[2] = PetscMax(loc[2], CmpDiffInf(a[0][j][i], a[mz-2][j][i]));
    if (HAVE_K(mz-1) && HAVE_K(1))
        for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[2] = PetscMax(loc[2], CmpDiffInf(a[mz-1][j][i], a[1][j][i]));
#undef HAVE_I
#undef HAVE_J
#undef HAVE_K
    PetscCall(DMDAVecRestoreArrayRead(user->fda, local, &a));
    PetscCallMPI(MPI_Allreduce(loc, glo, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    seam[0] = glo[0]; seam[1] = glo[1]; seam[2] = glo[2];
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeLocalOuterGhostMismatch()

static PetscErrorCode ComputeLocalOuterGhostMismatch ( UserCtx *  user, Vec  local, PetscReal  seam[3]  )

static

Computes outer periodic ghost mismatch for local Ucat.

Definition at line 352 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***a;
    PetscReal loc[3] = {0.0, 0.0, 0.0}, glo[3];
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->fda, local, &a));
#define HAVE_I(ii) InGhostRange((ii), info.gxs, info.gxm)
#define HAVE_J(jj) InGhostRange((jj), info.gys, info.gym)
#define HAVE_K(kk) InGhostRange((kk), info.gzs, info.gzm)
    if (HAVE_I(-1) && HAVE_I(1))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
                loc[0] = PetscMax(loc[0], CmpDiffInf(a[k][j][-1], a[k][j][1]));
    if (HAVE_I(mx) && HAVE_I(mx-2))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
                loc[0] = PetscMax(loc[0], CmpDiffInf(a[k][j][mx], a[k][j][mx-2]));
    if (HAVE_J(-1) && HAVE_J(1))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[1] = PetscMax(loc[1], CmpDiffInf(a[k][-1][i], a[k][1][i]));
    if (HAVE_J(my) && HAVE_J(my-2))
        for (PetscInt k = 1; k < mz-1; k++) if (HAVE_K(k))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[1] = PetscMax(loc[1], CmpDiffInf(a[k][my][i], a[k][my-2][i]));
    if (HAVE_K(-1) && HAVE_K(1))
        for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[2] = PetscMax(loc[2], CmpDiffInf(a[-1][j][i], a[1][j][i]));
    if (HAVE_K(mz) && HAVE_K(mz-2))
        for (PetscInt j = 1; j < my-1; j++) if (HAVE_J(j))
            for (PetscInt i = 1; i < mx-1; i++) if (HAVE_I(i))
                loc[2] = PetscMax(loc[2], CmpDiffInf(a[mz][j][i], a[mz-2][j][i]));
#undef HAVE_I
#undef HAVE_J
#undef HAVE_K
    PetscCall(DMDAVecRestoreArrayRead(user->fda, local, &a));
    PetscCallMPI(MPI_Allreduce(loc, glo, 3, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    seam[0] = glo[0]; seam[1] = glo[1]; seam[2] = glo[2];
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ConfigureCandidateFixture()

static PetscErrorCode ConfigureCandidateFixture ( SimCtx *  simCtx, UserCtx *  user  )

static

Configures the minimal context for centered inviscid periodic convection tests.

Definition at line 399 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    for (int f = 0; f < 6; f++) user->boundary_faces[f].mathematical_type = PERIODIC;
    simCtx->dt = 0.1; simCtx->step = 5; simCtx->StartStep = 0;     /* BDF2 -> a0=1.5 */
    simCtx->ren = 1.0e6; simCtx->invicid = 1; simCtx->les = 0; simCtx->rans = 0;
    simCtx->central = 1; simCtx->clark = 0; simCtx->TwoD = 0; simCtx->block_number = 1;
    simCtx->bulkVelocityCorrection = 0.0; simCtx->moveframe = 0; simCtx->rotateframe = 0;
    if (!user->lNu_t) PetscCall(DMCreateLocalVector(user->da, &user->lNu_t));
    PetscCall(VecSet(user->P, 0.0)); PetscCall(UpdateLocalGhosts(user, "P"));
    PetscCall(VecSet(user->lNvert, 0.0)); PetscCall(VecSet(user->Nvert, 0.0));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetUcatField()

static PetscErrorCode SetUcatField ( UserCtx *  user, CandState  st  )

static

Definition at line 414 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***u;
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &u));
    /* only owned region for a GLOBAL vec */
    for (PetscInt k = info.zs; k < info.zs+info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys+info.ym; j++)
            for (PetscInt i = info.xs; i < info.xs+info.xm; i++) {
                u[k][j][i] = AnalyticVelocity(st, i, j, k, mx, my, mz);
            }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &u));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetDirectUcontField()

static PetscErrorCode SetDirectUcontField ( UserCtx *  user, CandState  st  )

static

Sets the direct State B component-staggered Ucont field on owned entries.

Definition at line 434 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***u;
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    PetscCall(VecSet(user->Ucont, 0.0));
    PetscCall(DMDAVecGetArray(user->fda, user->Ucont, &u));
    for (PetscInt k = info.zs; k < info.zs+info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys+info.ym; j++)
            for (PetscInt i = info.xs; i < info.xs+info.xm; i++)
                u[k][j][i] = DirectUcontVelocity(st, i, j, k, mx, my, mz);
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucont, &u));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMaxDiscreteDivergence()

static PetscErrorCode ComputeMaxDiscreteDivergence ( UserCtx *  user, PetscReal *  maxdiv  )

static

Computes max discrete divergence of the current local Ucont field.

Definition at line 453 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    Cmpnts ***uc;
    PetscReal dv = 0.0, dv_global;
    PetscFunctionBeginUser;
    PetscCall(UpdateLocalGhosts(user, "Ucont"));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lUcont, &uc));
    for (PetscInt k = info.zs; k < info.zs+info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys+info.ym; j++)
            for (PetscInt i = info.xs; i < info.xs+info.xm; i++) {
                if (i<1||i>mx-2||j<1||j>my-2||k<1||k>mz-2) continue;
                const PetscReal d = (uc[k][j][i].x - uc[k][j][i-1].x)
                                  + (uc[k][j][i].y - uc[k][j-1][i].y)
                                  + (uc[k][j][i].z - uc[k-1][j][i].z);
                dv = PetscMax(dv, PetscAbsReal(d));
            }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lUcont, &uc));
    PetscCallMPI(MPI_Allreduce(&dv, &dv_global, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    *maxdiv = dv_global;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildBaseState()

static PetscErrorCode BuildBaseState ( UserCtx *  user, CandState  st, Vec  Ubase, PetscReal *  repeat_inf, PetscReal *  maxdiv, SeamDiagnostics *  seam  )

static

Definition at line 479 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscReal err = 0.0, err_global;
    Vec target;
    Cmpnts ***ur, ***ut;
    PetscFunctionBeginUser;
 
    PetscCall(VecDuplicate(user->Ucat, &target));
 
    if (seam) PetscCall(ComputeDeclaredSeamMismatch(st, info, seam->declared));
 
    if (st == STATE_B) {
        const char *ufld[] = {"Ucont"};
        const char *cfld[] = {"Ucat"};
        PetscCall(SetDirectUcontField(user, st));
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, ufld));
        PetscCall(VecCopy(user->Ucont, Ubase));
        if (seam) PetscCall(ComputeLocalDuplicateMismatch(user, user->lUcont, seam->ucont_global));
        PetscCall(Contra2Cart(user));
        PetscCall(SynchronizePeriodicCellFields(user, 1, cfld));
        PetscCall(VecCopy(user->Ucat, target));          /* saved recovered Cartesian state */
        PetscCall(Contra2Cart(user));
        PetscCall(SynchronizePeriodicCellFields(user, 1, cfld));
    } else {
        const char *cfld[] = {"Ucat"};
        const char *ufld[] = {"Ucont"};
        PetscCall(SetUcatField(user, st));
        PetscCall(SynchronizePeriodicCellFields(user, 1, cfld));
        PetscCall(VecCopy(user->Ucat, target));          /* saved declared Cartesian state */
        PetscCall(Cart2Contra(user));                    /* global Ucont from lUcat + metrics */
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, ufld));
        PetscCall(VecCopy(user->Ucont, Ubase));
        if (seam) PetscCall(ComputeLocalDuplicateMismatch(user, user->lUcont, seam->ucont_global));
        PetscCall(Contra2Cart(user));                    /* recovered Cartesian from Ucont */
        PetscCall(SynchronizePeriodicCellFields(user, 1, cfld));
    }
 
    if (seam) {
        PetscCall(UpdateLocalGhosts(user, "Ucat"));
        PetscCall(ComputeLocalDuplicateMismatch(user, user->lUcat, seam->ucat_global));
        PetscCall(ComputeLocalOuterGhostMismatch(user, user->lUcat, seam->ucat_ghost));
    }
 
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ur));
    PetscCall(DMDAVecGetArrayRead(user->fda, target, &ut));
    for (PetscInt k = info.zs; k < info.zs+info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys+info.ym; j++)
            for (PetscInt i = info.xs; i < info.xs+info.xm; i++) {
                if (i<1||i>mx-2||j<1||j>my-2||k<1||k>mz-2) continue;   /* interior cells only */
                err = PetscMax(err, PetscAbsReal(ur[k][j][i].x - ut[k][j][i].x));
                err = PetscMax(err, PetscAbsReal(ur[k][j][i].y - ut[k][j][i].y));
                err = PetscMax(err, PetscAbsReal(ur[k][j][i].z - ut[k][j][i].z));
            }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ur));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, target, &ut));
    PetscCall(VecDestroy(&target));
    PetscCall(UpdateLocalGhosts(user, "Ucat"));         /* lUcat now consistent with base Ucont */
 
    PetscCallMPI(MPI_Allreduce(&err, &err_global, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    PetscCall(ComputeMaxDiscreteDivergence(user, maxdiv));
    *repeat_inf = err_global;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeMaxGradientContribution()

static PetscErrorCode ComputeMaxGradientContribution ( UserCtx *  user, PetscReal *  gradmax  )

static

Computes the global maximum Cartesian velocity-gradient row-sum used by Candidate D.

Definition at line 550 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***ucat, ***csi, ***eta, ***zet;
    PetscReal ***aj;
    PetscReal loc = 0.0, glo;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz;
    PetscFunctionBeginUser;
    PetscCall(UpdateLocalGhosts(user, "Ucat"));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lUcat, &ucat));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCsi,  &csi));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lEta,  &eta));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lZet,  &zet));
    PetscCall(DMDAVecGetArrayRead(user->da,  user->lAj,   &aj));
    for (PetscInt k = info.zs; k < info.zs+info.zm; k++)
        for (PetscInt j = info.ys; j < info.ys+info.ym; j++)
            for (PetscInt i = info.xs; i < info.xs+info.xm; i++) {
                if (i<1||i>mx-2||j<1||j>my-2||k<1||k>mz-2) continue;
                const Cmpnts duc = { 0.5*(ucat[k][j][i+1].x-ucat[k][j][i-1].x),
                                     0.5*(ucat[k][j][i+1].y-ucat[k][j][i-1].y),
                                     0.5*(ucat[k][j][i+1].z-ucat[k][j][i-1].z) };
                const Cmpnts due = { 0.5*(ucat[k][j+1][i].x-ucat[k][j-1][i].x),
                                     0.5*(ucat[k][j+1][i].y-ucat[k][j-1][i].y),
                                     0.5*(ucat[k][j+1][i].z-ucat[k][j-1][i].z) };
                const Cmpnts duz = { 0.5*(ucat[k+1][j][i].x-ucat[k-1][j][i].x),
                                     0.5*(ucat[k+1][j][i].y-ucat[k-1][j][i].y),
                                     0.5*(ucat[k+1][j][i].z-ucat[k-1][j][i].z) };
                const Cmpnts C = csi[k][j][i], E = eta[k][j][i], Z = zet[k][j][i];
                const PetscReal Ajc = aj[k][j][i];
#define ROWSUM(cmp) ( \
    PetscAbsReal(Ajc*(C.x*duc.cmp + E.x*due.cmp + Z.x*duz.cmp)) + \
    PetscAbsReal(Ajc*(C.y*duc.cmp + E.y*due.cmp + Z.y*duz.cmp)) + \
    PetscAbsReal(Ajc*(C.z*duc.cmp + E.z*due.cmp + Z.z*duz.cmp)) )
                loc = PetscMax(loc, PetscMax(ROWSUM(x), PetscMax(ROWSUM(y), ROWSUM(z))));
#undef ROWSUM
            }
    PetscCall(DMDAVecRestoreArrayRead(user->da,  user->lAj,   &aj));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lZet,  &zet));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lEta,  &eta));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCsi,  &csi));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lUcat, &ucat));
    PetscCallMPI(MPI_Allreduce(&loc, &glo, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    *gradmax = glo;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DenseSpectralRadius()

static PetscErrorCode DenseSpectralRadius ( const PetscReal *  A, PetscInt  n, PetscReal *  rho, PetscReal *  maxRealPart  )

static

Definition at line 600 of file test_momentum_convective_candidates.c.

{
    PetscBLASInt N, lda, lwork, info;
    PetscReal *Acopy, *wr, *wi, *work, dummy = 0.0;
    PetscFunctionBeginUser;
    PetscCall(PetscBLASIntCast(n, &N)); lda = N; lwork = 8*N;
    PetscCall(PetscMalloc4(n*n, &Acopy, n, &wr, n, &wi, (size_t)lwork, &work));
    for (PetscInt t = 0; t < n*n; t++) Acopy[t] = A[t];
    {
        char nochar = 'N';
        LAPACKgeev_(&nochar, &nochar, &N, Acopy, &lda, wr, wi, &dummy, &lda, &dummy, &lda, work, &lwork, &info);
    }
    PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "LAPACK dgeev failed: info=%d", (int)info);
    *rho = 0.0; *maxRealPart = -PETSC_MAX_REAL;
    for (PetscInt t = 0; t < n; t++) {
        *rho = PetscMax(*rho, PetscSqrtReal(wr[t]*wr[t] + wi[t]*wi[t]));
        *maxRealPart = PetscMax(*maxRealPart, wr[t]);
    }
    PetscCall(PetscFree4(Acopy, wr, wi, work));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DenseMaxRealRightEigenpair()

static PetscErrorCode DenseMaxRealRightEigenpair ( const PetscReal *  A, PetscInt  n, PetscReal *  lamr, PetscReal *  lami, PetscReal *  vr_out, PetscReal *  vi_out  )

static

Extracts the right eigenpair whose eigenvalue has largest real part.

For a complex pair, LAPACK stores real and imaginary parts in adjacent columns of VR.

Definition at line 628 of file test_momentum_convective_candidates.c.

{
    PetscBLASInt N, lda, lwork, info;
    PetscReal *Acopy, *wr, *wi, *vr, *work, dummy = 0.0;
    PetscFunctionBeginUser;
    PetscCall(PetscBLASIntCast(n, &N)); lda = N; lwork = 16*N;
    PetscCall(PetscMalloc5(n*n, &Acopy, n, &wr, n, &wi, n*n, &vr, (size_t)lwork, &work));
    for (PetscInt t = 0; t < n*n; t++) Acopy[t] = A[t];
    {
        char jobvl = 'N', jobvr = 'V';
        LAPACKgeev_(&jobvl, &jobvr, &N, Acopy, &lda, wr, wi, &dummy, &lda, vr, &lda, work, &lwork, &info);
    }
    PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "LAPACK dgeev failed: info=%d", (int)info);
    PetscInt best = 0;
    for (PetscInt t = 1; t < n; t++) if (wr[t] > wr[best]) best = t;
    *lamr = wr[best]; *lami = wi[best];
    if (PetscAbsReal(wi[best]) < 1e-14) {
        for (PetscInt r = 0; r < n; r++) { vr_out[r] = vr[r + best*n]; vi_out[r] = 0.0; }
    } else if (wi[best] > 0.0) {
        for (PetscInt r = 0; r < n; r++) { vr_out[r] = vr[r + best*n]; vi_out[r] = vr[r + (best+1)*n]; }
    } else {
        for (PetscInt r = 0; r < n; r++) { vr_out[r] = vr[r + (best-1)*n]; vi_out[r] = -vr[r + best*n]; }
    }
    PetscCall(PetscFree5(Acopy, wr, wi, vr, work));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DenseRKPolynomialSpectralRadius()

static PetscErrorCode DenseRKPolynomialSpectralRadius ( const PetscReal *  J, PetscInt  n, PetscReal  dtau, PetscReal *  rho  )

static

Computes the spectral radius of the RK polynomial by applying it to eig(J).

Definition at line 660 of file test_momentum_convective_candidates.c.

{
    PetscBLASInt N, lda, lwork, info;
    PetscReal *Jcopy, *wr, *wi, *work, dummy = 0.0;
    PetscFunctionBeginUser;
    PetscCall(PetscBLASIntCast(n, &N)); lda = N; lwork = 8*N;
    PetscCall(PetscMalloc4(n*n, &Jcopy, n, &wr, n, &wi, (size_t)lwork, &work));
    for (PetscInt t = 0; t < n*n; t++) Jcopy[t] = J[t];
    {
        char nochar = 'N';
        LAPACKgeev_(&nochar, &nochar, &N, Jcopy, &lda, wr, wi, &dummy, &lda, &dummy, &lda, work, &lwork, &info);
    }
    PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "LAPACK dgeev failed: info=%d", (int)info);
    *rho = 0.0;
    for (PetscInt t = 0; t < n; t++) {
        const PetscReal zr = dtau*wr[t], zi = dtau*wi[t];
        const PetscReal z2r = zr*zr - zi*zi, z2i = 2.0*zr*zi;
        const PetscReal z3r = z2r*zr - z2i*zi, z3i = z2r*zi + z2i*zr;
        const PetscReal z4r = z3r*zr - z3i*zi, z4i = z3r*zi + z3i*zr;
        const PetscReal pr = 1.0 + zr + 0.5*z2r + z3r/6.0 + z4r/24.0;
        const PetscReal pi =       zi + 0.5*z2i + z3i/6.0 + z4i/24.0;
        *rho = PetscMax(*rho, PetscSqrtReal(pr*pr + pi*pi));
    }
    PetscCall(PetscFree4(Jcopy, wr, wi, work));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DenseSigmaMax()

static PetscErrorCode DenseSigmaMax ( const PetscReal *  A, PetscInt  n, PetscReal *  smax, PetscReal *  v1  )

static

Definition at line 689 of file test_momentum_convective_candidates.c.

{
    PetscBLASInt N, lda, lwork, info;
    PetscReal *Acopy, *S, *VT, *work, ufake = 0.0;
    PetscFunctionBeginUser;
    PetscCall(PetscBLASIntCast(n, &N)); lda = N; lwork = 8*N + 4*N;
    PetscCall(PetscMalloc4(n*n, &Acopy, n, &S, n*n, &VT, (size_t)lwork, &work));
    for (PetscInt t = 0; t < n*n; t++) Acopy[t] = A[t];
    {
        char jobu = 'N', jobvt = v1 ? 'S' : 'N';
        LAPACKgesvd_(&jobu, &jobvt, &N, &N, Acopy, &lda, S, &ufake, &lda, VT, &lda, work, &lwork, &info);
    }
    PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "LAPACK dgesvd failed: info=%d", (int)info);
    *smax = S[0];
    if (v1) { for (PetscInt r = 0; r < n; r++) v1[r] = VT[0 + r*n]; }  /* first row of VT = v1^T */
    PetscCall(PetscFree4(Acopy, S, VT, work));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DenseNonNormality()

static PetscReal DenseNonNormality ( const PetscReal *  A, PetscInt  n  )

static

Definition at line 709 of file test_momentum_convective_candidates.c.

{
    PetscReal fro2 = 0.0, comm = 0.0;
    for (PetscInt t = 0; t < n*n; t++) fro2 += A[t]*A[t];
    for (PetscInt p = 0; p < n; p++)
        for (PetscInt q = 0; q < n; q++) {
            PetscReal ata = 0.0, aat = 0.0;
            for (PetscInt r = 0; r < n; r++) { ata += A[p + r*n]*A[q + r*n]; aat += A[r + p*n]*A[r + q*n]; }
            const PetscReal d = ata - aat; comm += d*d;
        }
    return PetscSqrtReal(comm) / PetscMax(fro2, PETSC_MACHINE_EPSILON);
}

Here is the caller graph for this function:

DenseSkewnessDefect()

static PetscReal DenseSkewnessDefect ( const PetscReal *  A, PetscInt  n  )

static

Computes the normalized Frobenius defect from skew symmetry.

Definition at line 725 of file test_momentum_convective_candidates.c.

{
    PetscReal fro2 = 0.0, sym2 = 0.0;
    for (PetscInt i = 0; i < n*n; i++) fro2 += A[i]*A[i];
    for (PetscInt c = 0; c < n; c++)
        for (PetscInt r = 0; r < n; r++) {
            const PetscReal s = A[r + c*n] + A[c + r*n];
            sym2 += s*s;
        }
    return PetscSqrtReal(sym2) / PetscMax(PetscSqrtReal(fro2), PETSC_MACHINE_EPSILON);
}

Here is the caller graph for this function:

PrintSpectrumSummary()

static PetscErrorCode PrintSpectrumSummary ( const char *  name, const PetscReal *  J, PetscInt  n  )

static

Prints eigenvalue and norm summary for one dense operator.

Definition at line 740 of file test_momentum_convective_candidates.c.

{
    PetscReal rho, maxre, smax;
    PetscFunctionBeginUser;
    PetscCall(DenseSpectralRadius(J, n, &rho, &maxre));
    PetscCall(DenseSigmaMax(J, n, &smax, NULL));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "    %-7s rho=%.6e sigma=%.6e max_real=%.6e skew=%.3e nonnormality=%.3e\n",
        name, (double)rho, (double)smax, (double)maxre,
        (double)DenseSkewnessDefect(J, n), (double)DenseNonNormality(J, n)));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DenseMatVec()

static void DenseMatVec ( const PetscReal *  J, const PetscReal *  x, PetscReal *  y, PetscInt  n  )

static

Applies a dense column-major matrix to an active-space vector.

Definition at line 1440 of file test_momentum_convective_candidates.c.

{
    for (PetscInt r = 0; r < n; r++) {
        PetscReal s = 0.0;
        for (PetscInt c = 0; c < n; c++) s += J[r + c*n]*x[c];
        y[r] = s;
    }
}

Here is the caller graph for this function:

FillAuditDirection()

static void FillAuditDirection ( const DofMap *  map, PetscInt  kind, PetscReal *  x  )

static

Fills one deterministic active-space vector used by dense/matrix-free checks.

Definition at line 1452 of file test_momentum_convective_candidates.c.

{
    PetscReal n2 = 0.0;
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscReal a = 0.31*(PetscReal)(map->ci[m]+1)
                          + 0.47*(PetscReal)(map->cj[m]+1)
                          + 0.59*(PetscReal)(map->ck[m]+1)
                          + 0.23*(PetscReal)(map->comp[m]+1);
        if (kind == 0) x[m] = PetscSinReal(a);
        else if (kind == 1) x[m] = PetscCosReal(1.7*a) + 0.25*PetscSinReal(0.9*(PetscReal)(m+1));
        else x[m] = (map->comp[m] == kind-2) ? PetscSinReal(a) : 0.0;
        n2 += x[m]*x[m];
    }
    n2 = PetscSqrtReal(n2);
    if (n2 > 0.0) for (PetscInt m = 0; m < map->n; m++) x[m] /= n2;
}

Here is the caller graph for this function:

AddActiveVector()

static PetscErrorCode AddActiveVector ( UserCtx *  user, Vec  U, const DofMap *  map, const PetscReal *  x, PetscReal  scale  )

static

Adds a dense active-space vector into a global contravariant vector.

Definition at line 1276 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***a;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArray(user->fda, U, &a));
    for (PetscInt m = 0; m < map->n; m++) {
        PetscReal *p = (PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]];
        p[map->comp[m]] += scale*x[m];
    }
    PetscCall(DMDAVecRestoreArray(user->fda, U, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DenseShiftIdentity()

static PetscErrorCode DenseShiftIdentity ( const PetscReal *  A, PetscInt  n, PetscReal  shift, PetscReal *  B  )

static

Copies a dense matrix and adds a scalar shift to its diagonal.

Definition at line 761 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    PetscCall(PetscArraycpy(B, A, (size_t)n*n));
    for (PetscInt d = 0; d < n; d++) B[d + d*n] += shift;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MatMul()

static void MatMul ( const PetscReal *  A, const PetscReal *  B, PetscReal *  out, PetscInt  n  )

static

Definition at line 770 of file test_momentum_convective_candidates.c.

{
    for (PetscInt c = 0; c < n; c++)
        for (PetscInt r = 0; r < n; r++) {
            PetscReal s = 0.0;
            for (PetscInt t = 0; t < n; t++) s += A[r + t*n]*B[t + c*n];
            out[r + c*n] = s;
        }
}

Here is the caller graph for this function:

RKPolynomial()

static PetscErrorCode RKPolynomial ( const PetscReal *  J, PetscReal  dtau, PetscInt  n, PetscReal *  P  )

static

Definition at line 781 of file test_momentum_convective_candidates.c.

{
    PetscReal *M, *T, *T2;
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc3(n*n, &M, n*n, &T, n*n, &T2));
    for (PetscInt t = 0; t < n*n; t++) M[t] = dtau*J[t];
    /* start S = I + M/4 */
    for (PetscInt t = 0; t < n*n; t++) T[t] = M[t]/4.0;
    for (PetscInt d = 0; d < n; d++) T[d + d*n] += 1.0;
    const PetscReal coef[3] = {3.0, 2.0, 1.0};   /* divide by 3, then 2, then 1 */
    for (int s = 0; s < 3; s++) {
        MatMul(M, T, T2, n);                      /* T2 = M*S */
        for (PetscInt t = 0; t < n*n; t++) T[t] = T2[t]/coef[s];
        for (PetscInt d = 0; d < n; d++) T[d + d*n] += 1.0;   /* S = I + M/coef * S */
    }
    for (PetscInt t = 0; t < n*n; t++) P[t] = T[t];
    PetscCall(PetscFree3(M, T, T2));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

AmplificationMetric()

static PetscErrorCode AmplificationMetric ( const PetscReal *  J, PetscInt  n, PetscReal  dtau, PMetric  which, PetscReal *  metric  )

static

Evaluates either spectral-radius or 2-norm amplification for one pseudo-time step.

Definition at line 806 of file test_momentum_convective_candidates.c.

{
    PetscReal *Pm;
    PetscFunctionBeginUser;
    if (which == METRIC_RHO) {
        PetscCall(DenseRKPolynomialSpectralRadius(J, n, dtau, metric));
        PetscFunctionReturn(0);
    }
    PetscCall(PetscMalloc1(n*n, &Pm));
    PetscCall(RKPolynomial(J, dtau, n, Pm));
    PetscCall(DenseSigmaMax(Pm, n, metric, NULL));
    PetscCall(PetscFree(Pm));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StableCFL()

static PetscErrorCode StableCFL ( const PetscReal *  J, PetscInt  n, PetscReal  lam, PMetric  which, StableCFLResult *  result  )

static

Definition at line 823 of file test_momentum_convective_candidates.c.

{
    const PetscReal tol = 1e-8, probe = 1e-8, scan_step = 0.01, hi = 4.0;
    const PetscReal probe_tol = 1e-12, min_positive_cfl = 1e-6;
    PetscReal met;
    PetscFunctionBeginUser;
    result->status = STABLE_CFL_NONE;
    result->cfl = 0.0;
    if (which == METRIC_RHO) {
        PetscReal rhoJ, maxreJ;
        PetscCall(DenseSpectralRadius(J, n, &rhoJ, &maxreJ));
        if (maxreJ > 1e-8) PetscFunctionReturn(0);
    }
    PetscCall(AmplificationMetric(J, n, probe/lam, which, &met));
    if (met > 1.0 + probe_tol) PetscFunctionReturn(0);
 
    PetscReal stable = probe, cross_b = -1.0;
    for (PetscReal cfl = scan_step; cfl <= hi + 1e-12; cfl += scan_step) {
        PetscCall(AmplificationMetric(J, n, cfl/lam, which, &met));
        if (met > 1.0 + tol) { cross_b = cfl; break; }
        stable = cfl;
    }
    if (cross_b < 0.0) {
        result->status = STABLE_CFL_EXCEEDS_SCAN;
        result->cfl = hi;
        PetscFunctionReturn(0);
    }
    PetscReal cross_a = stable;
    for (int it = 0; it < 40; it++) {
        const PetscReal mid = 0.5*(cross_a + cross_b);
        PetscCall(AmplificationMetric(J, n, mid/lam, which, &met));
        if (met > 1.0 + tol) cross_b = mid; else cross_a = mid;
    }
    if (cross_a < min_positive_cfl) PetscFunctionReturn(0);
    result->status = STABLE_CFL_FINITE;
    result->cfl = cross_a;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StableCFLStatusText()

static const char * StableCFLStatusText ( StableCFLResult  r )

static

Returns human-readable text for a stable-CFL search result.

Definition at line 866 of file test_momentum_convective_candidates.c.

{
    if (r.status == STABLE_CFL_NONE) return "no positive stable interval connected to zero";
    if (r.status == STABLE_CFL_EXCEEDS_SCAN) return "stable through CFL >= 4.0";
    return "stable CFL";
}

Here is the caller graph for this function:

PrintStableCFLLine()

static PetscErrorCode PrintStableCFLLine ( const char *  candidate, StableCFLResult  eig, StableCFLResult  norm  )

static

Prints one candidate's eigenvalue and norm stable-CFL statuses.

Definition at line 876 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "    cand %s: eig %s", candidate, StableCFLStatusText(eig)));
    if (eig.status == STABLE_CFL_FINITE) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " = %.4f", (double)eig.cfl));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, " | norm %s", StableCFLStatusText(norm)));
    if (norm.status == STABLE_CFL_FINITE) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " = %.4f", (double)norm.cfl));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyP()

static void ApplyP ( const PetscReal *  P, const PetscReal *  x, PetscReal *  out, PetscInt  n  )

static

Definition at line 890 of file test_momentum_convective_candidates.c.

{
    for (PetscInt r = 0; r < n; r++) { PetscReal s = 0.0; for (PetscInt c = 0; c < n; c++) s += P[r + c*n]*x[c]; out[r] = s; }
}

Here is the caller graph for this function:

PrintFrozenAmplificationTable()

static PetscErrorCode PrintFrozenAmplificationTable ( const char *  title, const PetscReal *  J, PetscInt  n, const PetscReal  lams[3], const char *  cn[3]  )

static

Prints frozen RK amplification tables for the supplied operator and candidates.

Definition at line 898 of file test_momentum_convective_candidates.c.

{
    const PetscReal cfls[5] = {0.25, 0.50, 1.00, 1.50, 2.00};
    PetscReal *Pm;
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc1(n*n, &Pm));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- %s ---\n"
        "    candidate  CFL    dtau          rho(P)        sigma_max(P)\n", title));
    for (int c = 0; c < 3; c++) {
        if (!(lams[c] > 0.0)) continue;
        for (int q = 0; q < 5; q++) {
            const PetscReal dtau = cfls[q]/lams[c];
            PetscReal rhoP, smaxP;
            PetscCall(RKPolynomial(J, dtau, n, Pm));
            PetscCall(DenseRKPolynomialSpectralRadius(J, n, dtau, &rhoP));
            PetscCall(DenseSigmaMax(Pm, n, &smaxP, NULL));
            PetscCall(PetscPrintf(PETSC_COMM_WORLD,
                "    %-9s %.2f   %.6e  %.6e  %.6e\n",
                cn[c], (double)cfls[q], (double)dtau, (double)rhoP, (double)smaxP));
        }
    }
    PetscCall(PetscFree(Pm));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

VecNorm2Array()

static PetscReal VecNorm2Array ( const PetscReal *  x, PetscInt  n  )

static

Computes the Euclidean norm of a dense vector.

Definition at line 929 of file test_momentum_convective_candidates.c.

{
    PetscReal s = 0.0;
    for (PetscInt i = 0; i < n; i++) s += x[i]*x[i];
    return PetscSqrtReal(s);
}

Here is the caller graph for this function:

VecNormInfArray()

static PetscReal VecNormInfArray ( const PetscReal *  x, PetscInt  n  )

static

Computes the infinity norm of a dense vector.

Definition at line 939 of file test_momentum_convective_candidates.c.

{
    PetscReal s = 0.0;
    for (PetscInt i = 0; i < n; i++) s = PetscMax(s, PetscAbsReal(x[i]));
    return s;
}

Here is the caller graph for this function:

DenseFrobenius()

static PetscReal DenseFrobenius ( const PetscReal *  A, PetscInt  n  )

static

Computes the Frobenius norm of a dense column-major matrix.

Definition at line 949 of file test_momentum_convective_candidates.c.

{
    PetscReal s = 0.0;
    for (PetscInt i = 0; i < n*n; i++) s += A[i]*A[i];
    return PetscSqrtReal(s);
}

Here is the caller graph for this function:

DofTouchesPeriodicRepresentative()

static PetscBool DofTouchesPeriodicRepresentative ( const DofMap *  map, PetscInt  m, DMDALocalInfo  info  )

static

True if an active representative lies on a plane adjacent to periodic duplicates.

Definition at line 959 of file test_momentum_convective_candidates.c.

{
    return (PetscBool)(map->ci[m] == 1 || map->ci[m] == info.mx-2 ||
                       map->cj[m] == 1 || map->cj[m] == info.my-2 ||
                       map->ck[m] == 1 || map->ck[m] == info.mz-2);
}

Here is the caller graph for this function:

PrintPeriodicSpaceAudit()

static PetscErrorCode PrintPeriodicSpaceAudit ( UserCtx *  user, const DofMap *  map  )

static

Prints the component-wise periodic storage count actually used by the active map.

Definition at line 969 of file test_momentum_convective_candidates.c.

{
    PetscInt count[3] = {0,0,0};
    PetscInt min_i[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
    PetscInt min_j[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
    PetscInt min_k[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
    PetscInt max_i[3] = {-PETSC_MAX_INT,-PETSC_MAX_INT,-PETSC_MAX_INT};
    PetscInt max_j[3] = {-PETSC_MAX_INT,-PETSC_MAX_INT,-PETSC_MAX_INT};
    PetscInt max_k[3] = {-PETSC_MAX_INT,-PETSC_MAX_INT,-PETSC_MAX_INT};
    PetscFunctionBeginUser;
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscInt c = map->comp[m];
        count[c]++;
        min_i[c] = PetscMin(min_i[c], map->ci[m]); max_i[c] = PetscMax(max_i[c], map->ci[m]);
        min_j[c] = PetscMin(min_j[c], map->cj[m]); max_j[c] = PetscMax(max_j[c], map->cj[m]);
        min_k[c] = PetscMin(min_k[c], map->ck[m]); max_k[c] = PetscMax(max_k[c], map->ck[m]);
    }
    const PetscInt per_comp = PeriodicRepCount(user->info.mx) *
                              PeriodicRepCount(user->info.my) *
                              PeriodicRepCount(user->info.mz);
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- periodic independent-space audit ---\n"
        "    synchronization convention: duplicate planes 0<-m-2 and m-1<-1 in x/y/z for Ucont.\n"
        "    active rows and columns both use the same representatives i,j,k=1..m-2.\n"
        "    comp  expected  actual   i-range   j-range   k-range\n"));
    for (PetscInt c = 0; c < 3; c++) {
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "      %d   %7d  %6d   [%d,%d]   [%d,%d]   [%d,%d]\n",
            (int)c, (int)per_comp, (int)count[c],
            (int)min_i[c], (int)max_i[c], (int)min_j[c], (int)max_j[c], (int)min_k[c], (int)max_k[c]));
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PrintBlockAndSymmetricLocalization()

static PetscErrorCode PrintBlockAndSymmetricLocalization ( UserCtx *  user, const DofMap *  map, const PetscReal *  J  )

static

Prints 3x3 component block norms and localized symmetric rows.

Definition at line 1006 of file test_momentum_convective_candidates.c.

{
    PetscReal block2[3][3] = {{0.0}}, sym2[3][3] = {{0.0}};
    PetscReal *row2;
    PetscInt top[8];
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc1(map->n, &row2));
    for (PetscInt r = 0; r < map->n; r++) row2[r] = 0.0;
    for (PetscInt c = 0; c < map->n; c++) {
        for (PetscInt r = 0; r < map->n; r++) {
            const PetscInt rb = map->comp[r], cb = map->comp[c];
            const PetscReal a = J[r + c*map->n];
            const PetscReal s = 0.5*(J[r + c*map->n] + J[c + r*map->n]);
            block2[rb][cb] += a*a;
            sym2[rb][cb] += s*s;
            row2[r] += s*s;
        }
    }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- component block norms J_ab = dR_a/dU_b ---\n"
        "    row-comp col-comp   ||J_ab||F    ||0.5(J+J^T)_ab||F\n"));
    for (PetscInt rb = 0; rb < 3; rb++)
        for (PetscInt cb = 0; cb < 3; cb++)
            PetscCall(PetscPrintf(PETSC_COMM_WORLD,
                "       %d        %d      %.6e      %.6e\n",
                (int)rb, (int)cb, (double)PetscSqrtReal(block2[rb][cb]),
                (double)PetscSqrtReal(sym2[rb][cb])));
 
    for (PetscInt q = 0; q < 8; q++) {
        top[q] = -1;
        for (PetscInt r = 0; r < map->n; r++) {
            PetscBool used = PETSC_FALSE;
            for (PetscInt p = 0; p < q; p++) if (top[p] == r) used = PETSC_TRUE;
            if (!used && (top[q] < 0 || row2[r] > row2[top[q]])) top[q] = r;
        }
    }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- largest rows of S=0.5*(J+J^T) ---\n"
        "    row comp (i,j,k) seam-adj  ||S_row||2   largest symmetric columns\n"));
    for (PetscInt q = 0; q < 8; q++) {
        const PetscInt r = top[q];
        PetscInt best[3] = {-1,-1,-1};
        for (PetscInt pass = 0; pass < 3; pass++) {
            for (PetscInt c = 0; c < map->n; c++) {
                PetscBool used = PETSC_FALSE;
                for (PetscInt p = 0; p < pass; p++) if (best[p] == c) used = PETSC_TRUE;
                const PetscReal mag = PetscAbsReal(0.5*(J[r + c*map->n] + J[c + r*map->n]));
                if (!used && (best[pass] < 0 ||
                    mag > PetscAbsReal(0.5*(J[r + best[pass]*map->n] + J[best[pass] + r*map->n])))) best[pass] = c;
            }
        }
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    %3d  %d   (%d,%d,%d)    %s    %.6e",
            (int)r, (int)map->comp[r], (int)map->ci[r], (int)map->cj[r], (int)map->ck[r],
            DofTouchesPeriodicRepresentative(map, r, user->info) ? "yes" : "no ",
            (double)PetscSqrtReal(row2[r])));
        for (PetscInt p = 0; p < 3; p++) {
            const PetscInt c = best[p];
            const PetscReal s = 0.5*(J[r + c*map->n] + J[c + r*map->n]);
            PetscCall(PetscPrintf(PETSC_COMM_WORLD, " | c%d:%d(%d,%d,%d)=%.3e",
                (int)p, (int)map->comp[c], (int)map->ci[c], (int)map->cj[c],
                (int)map->ck[c], (double)s));
        }
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
    }
    PetscCall(PetscFree(row2));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FillStaggeredFourierBasis()

static void FillStaggeredFourierBasis ( const DofMap *  map, PetscInt  wx, PetscInt  wy, PetscInt  wz, PetscReal *  Q  )

static

Fills the six real basis vectors for one staggered same-wavevector subspace.

Definition at line 1079 of file test_momentum_convective_candidates.c.

{
    PetscInt nrep = 0;
    for (PetscInt m = 0; m < map->n; m++) nrep = PetscMax(nrep, map->ci[m]);
    for (PetscInt q = 0; q < 6*map->n; q++) Q[q] = 0.0;
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscInt comp = map->comp[m];
        const PetscReal x = (PetscReal)(map->ci[m]-1) - (comp == 0 ? 0.5 : 0.0);
        const PetscReal y = (PetscReal)(map->cj[m]-1) - (comp == 1 ? 0.5 : 0.0);
        const PetscReal z = (PetscReal)(map->ck[m]-1) - (comp == 2 ? 0.5 : 0.0);
        const PetscReal phase = 2.0*PETSC_PI*((PetscReal)wx*x + (PetscReal)wy*y + (PetscReal)wz*z)/(PetscReal)nrep;
        Q[m + (2*comp+0)*map->n] = PetscCosReal(phase);
        Q[m + (2*comp+1)*map->n] = PetscSinReal(phase);
    }
    for (PetscInt q = 0; q < 6; q++) {
        PetscReal n2 = 0.0;
        for (PetscInt m = 0; m < map->n; m++) n2 += Q[m + q*map->n]*Q[m + q*map->n];
        n2 = PetscSqrtReal(n2);
        if (n2 > 0.0) for (PetscInt m = 0; m < map->n; m++) Q[m + q*map->n] /= n2;
    }
}

Here is the caller graph for this function:

StaggeredFourierSymbol()

static PetscErrorCode StaggeredFourierSymbol ( const DofMap *  map, const PetscReal *  J, PetscInt  wx, PetscInt  wy, PetscInt  wz, PetscReal  A6[36], PetscReal *  leak  )

static

Builds the 6x6 projected real symbol and leakage for one wavevector.

Definition at line 1105 of file test_momentum_convective_candidates.c.

{
    PetscReal *Q, *JQ;
    PetscReal all2 = 0.0, leak2 = 0.0;
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc2((size_t)6*map->n, &Q, (size_t)6*map->n, &JQ));
    FillStaggeredFourierBasis(map, wx, wy, wz, Q);
    for (PetscInt q = 0; q < 6; q++) DenseMatVec(J, &Q[q*map->n], &JQ[q*map->n], map->n);
    for (PetscInt c = 0; c < 6; c++)
        for (PetscInt r = 0; r < 6; r++) {
            PetscReal s = 0.0;
            for (PetscInt m = 0; m < map->n; m++) s += Q[m + r*map->n] * JQ[m + c*map->n];
            A6[r + c*6] = s;
        }
    for (PetscInt c = 0; c < 6; c++) {
        for (PetscInt m = 0; m < map->n; m++) {
            PetscReal proj = 0.0;
            for (PetscInt r = 0; r < 6; r++) proj += Q[m + r*map->n] * A6[r + c*6];
            const PetscReal d = JQ[m + c*map->n] - proj;
            leak2 += d*d;
            all2 += JQ[m + c*map->n]*JQ[m + c*map->n];
        }
    }
    *leak = PetscSqrtReal(leak2) / PetscMax(PETSC_MACHINE_EPSILON, PetscSqrtReal(all2));
    PetscCall(PetscFree2(Q, JQ));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SymbolEigenSummary()

static PetscErrorCode SymbolEigenSummary ( const PetscReal  A6[36], PetscReal *  maxre, PetscReal  wr_out[6], PetscReal  wi_out[6]  )

static

Eigenvalue summary for a 6x6 real symbol.

Definition at line 1138 of file test_momentum_convective_candidates.c.

{
    PetscBLASInt N = 6, lda = 6, lwork = 128, info;
    PetscReal Acopy[36], work[128], dummy = 0.0;
    PetscFunctionBeginUser;
    for (PetscInt t = 0; t < 36; t++) Acopy[t] = A6[t];
    {
        char nochar = 'N';
        LAPACKgeev_(&nochar, &nochar, &N, Acopy, &lda, wr_out, wi_out, &dummy, &lda, &dummy, &lda, work, &lwork, &info);
    }
    PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "LAPACK dgeev failed for 6x6 symbol: info=%d", (int)info);
    *maxre = -PETSC_MAX_REAL;
    for (PetscInt q = 0; q < 6; q++) *maxre = PetscMax(*maxre, wr_out[q]);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PrintFourierChecks()

static PetscErrorCode PrintFourierChecks ( const DofMap *  map, const PetscReal *  J  )

static

Projects production J on full staggered six-dimensional Fourier subspaces.

Definition at line 1158 of file test_momentum_convective_candidates.c.

{
    typedef struct { const char *name; PetscInt wx, wy, wz; } Mode;
    const Mode modes[4] = {{"x",1,0,0}, {"y",0,1,0}, {"xy",1,1,0}, {"x-y",1,2,0}};
    PetscReal best_re = -PETSC_MAX_REAL, best_leak = 0.0;
    PetscInt best_wx = 0, best_wy = 0, best_wz = 0, nrep = 0;
    PetscFunctionBeginUser;
    for (PetscInt m = 0; m < map->n; m++) nrep = PetscMax(nrep, map->ci[m]);
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- staggered 6x6 Fourier-symbol check ---\n"
        "    mode      max_real    leakage    eigenvalues (real,imag)\n"));
    for (PetscInt im = 0; im < 4; im++) {
        PetscReal A6[36], wr[6], wi[6], maxre, leak;
        PetscCall(StaggeredFourierSymbol(map, J, modes[im].wx, modes[im].wy, modes[im].wz, A6, &leak));
        PetscCall(SymbolEigenSummary(A6, &maxre, wr, wi));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    %-5s(%d,%d,%d) %.6e  %.3e",
            modes[im].name, (int)modes[im].wx, (int)modes[im].wy, (int)modes[im].wz,
            (double)maxre, (double)leak));
        for (PetscInt q = 0; q < 6; q++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  %.3e%+.3ei", (double)wr[q], (double)wi[q]));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
    }
    for (PetscInt wz = 0; wz < nrep; wz++)
        for (PetscInt wy = 0; wy < nrep; wy++)
            for (PetscInt wx = 0; wx < nrep; wx++) {
                if (wx == 0 && wy == 0 && wz == 0) continue;
                PetscReal A6[36], wr[6], wi[6], maxre, leak;
                PetscCall(StaggeredFourierSymbol(map, J, wx, wy, wz, A6, &leak));
                PetscCall(SymbolEigenSummary(A6, &maxre, wr, wi));
                if (maxre > best_re) {
                    best_re = maxre; best_leak = leak; best_wx = wx; best_wy = wy; best_wz = wz;
                }
            }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "    best scanned wavevector (%d,%d,%d): max_real=%.6e leakage=%.3e\n",
        (int)best_wx, (int)best_wy, (int)best_wz, (double)best_re, (double)best_leak));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TraceStateAResidualStages()

static PetscErrorCode TraceStateAResidualStages ( UserCtx *  user, Vec  Ubase, Vec  Rhs, const DofMap *  map, const PetscReal *  J  )

static

Trace observable production residual stages for one deterministic perturbation.

Definition at line 1200 of file test_momentum_convective_candidates.c.

{
    Vec Upert, Conv;
    PetscReal *v, *ract;
    PetscReal seam_ucont[3], seam_ucat[3], conv2, convinf, rhs2, rhsinf;
    PetscInt active_rows[4] = {0,0,0,0};
    PetscFunctionBeginUser;
    PetscCall(VecDuplicate(Ubase, &Upert));
    PetscCall(VecDuplicate(user->lUcont, &Conv));
    PetscCall(PetscMalloc2(map->n, &v, map->n, &ract));
    FillAuditDirection(map, 0, v);
    PetscCall(VecCopy(Ubase, Upert));
    PetscCall(AddActiveVector(user, Upert, map, v, 1e-6));
    PetscCall(VecCopy(Upert, user->Ucont));
    {
        const char *fld[] = {"Ucont"};
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fld));
    }
    PetscCall(ComputeLocalDuplicateMismatch(user, user->lUcont, seam_ucont));
    PetscCall(Contra2Cart(user));
    PetscCall(UpdateLocalGhosts(user, "Ucat"));
    PetscCall(ComputeLocalDuplicateMismatch(user, user->lUcat, seam_ucat));
    PetscCall(Convection(user, user->lUcont, user->lUcat, Conv));
    PetscCall(VecNorm(Conv, NORM_2, &conv2));
    PetscCall(VecNorm(Conv, NORM_INFINITY, &convinf));
    PetscCall(EvalConvResidual(user, Upert, Rhs, map, ract));
    rhs2 = VecNorm2Array(ract, map->n);
    rhsinf = VecNormInfArray(ract, map->n);
    {
        PetscReal ***nvert;
        PetscCall(DMDAVecGetArrayRead(user->da, user->lNvert, &nvert));
        for (PetscInt m = 0; m < map->n; m++) {
            const PetscInt rows = MomCellActiveRows(nvert, map->ck[m], map->cj[m], map->ci[m],
                                                    user->info.mx, user->info.my, user->info.mz,
                                                    PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, 0);
            if (rows & (1 << map->comp[m])) active_rows[map->comp[m]]++;
            else active_rows[3]++;
        }
        PetscCall(DMDAVecRestoreArrayRead(user->da, user->lNvert, &nvert));
    }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- observable State A residual trace (mixed perturbation, amp=1e-6) ---\n"
        "    after Ucont periodic sync duplicate mismatch (x,y,z)=%.3e %.3e %.3e\n"
        "    after Contra2Cart Ucat duplicate mismatch       (x,y,z)=%.3e %.3e %.3e\n"
        "    direct Convection() Cartesian norm: ||Conv||2=%.6e ||Conv||inf=%.6e\n"
        "    final active ComputeRHS() norm:     ||R||2=%.6e ||R||inf=%.6e\n"
        "    active-row mask counts by component: u=%d v=%d w=%d discarded=%d\n"
        "    final active-space skewness defect already reported from J: %.3e\n",
        (double)seam_ucont[0], (double)seam_ucont[1], (double)seam_ucont[2],
        (double)seam_ucat[0], (double)seam_ucat[1], (double)seam_ucat[2],
        (double)conv2, (double)convinf, (double)rhs2, (double)rhsinf,
        (int)active_rows[0], (int)active_rows[1], (int)active_rows[2], (int)active_rows[3],
        (double)DenseSkewnessDefect(J, map->n)));
    PetscCall(PetscFree2(v, ract));
    PetscCall(VecDestroy(&Upert));
    PetscCall(VecDestroy(&Conv));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DenseRelativeDiff()

static PetscReal DenseRelativeDiff ( const PetscReal *  A, const PetscReal *  B, PetscInt  n, PetscReal  denom_ref  )

static

Computes a normalized Frobenius difference between two dense matrices.

Definition at line 1263 of file test_momentum_convective_candidates.c.

{
    PetscReal s = 0.0;
    for (PetscInt i = 0; i < n*n; i++) {
        const PetscReal d = A[i] - B[i];
        s += d*d;
    }
    return PetscSqrtReal(s) / PetscMax(1.0, denom_ref);
}

Here is the caller graph for this function:

ExtractActiveVector()

static PetscErrorCode ExtractActiveVector ( UserCtx *  user, Vec  U, const DofMap *  map, PetscReal *  x  )

static

Extracts active-space entries from a global contravariant vector.

Definition at line 1293 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***a;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->fda, U, &a));
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscReal *p = (const PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]];
        x[m] = p[map->comp[m]];
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, U, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DofWeight()

static PetscReal DofWeight ( const DofMap *  map, PetscInt  m  )

static

Returns a deterministic checksum weight for an active DOF.

Definition at line 1309 of file test_momentum_convective_candidates.c.

{
    return 1.0 + 0.013*(PetscReal)(map->comp[m]+1)
               + 0.017*(PetscReal)map->ci[m]
               + 0.019*(PetscReal)map->cj[m]
               + 0.023*(PetscReal)map->ck[m];
}

Here is the caller graph for this function:

ActiveStats()

static PetscErrorCode ActiveStats ( const DofMap *  map, const PetscReal *  x, GlobalVecStats *  stats  )

static

Computes global active-vector norms and checksum.

Definition at line 1320 of file test_momentum_convective_candidates.c.

{
    PetscReal loc2 = 0.0, locinf = 0.0, locsum = 0.0;
    PetscReal glo2, gloinf, glosum;
    PetscFunctionBeginUser;
    for (PetscInt m = 0; m < map->n; m++) {
        loc2 += x[m]*x[m];
        locinf = PetscMax(locinf, PetscAbsReal(x[m]));
        locsum += DofWeight(map, m)*x[m];
    }
    PetscCallMPI(MPI_Allreduce(&loc2, &glo2, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
    PetscCallMPI(MPI_Allreduce(&locinf, &gloinf, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
    PetscCallMPI(MPI_Allreduce(&locsum, &glosum, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
    stats->n2 = PetscSqrtReal(glo2); stats->ninf = gloinf; stats->checksum = glosum;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FillDeterministicDirection()

static PetscErrorCode FillDeterministicDirection ( UserCtx *  user, Vec  V, const DofMap *  map, PetscReal *  x  )

static

Fills a globally normalized deterministic active-space perturbation direction.

Definition at line 1340 of file test_momentum_convective_candidates.c.

{
    PetscReal loc2 = 0.0, glo2;
    PetscFunctionBeginUser;
    PetscCall(VecSet(V, 0.0));
    for (PetscInt m = 0; m < map->n; m++) {
        x[m] = PetscSinReal(0.37*(PetscReal)(map->ci[m]+1)
                          + 0.51*(PetscReal)(map->cj[m]+1)
                          + 0.73*(PetscReal)(map->ck[m]+1)
                          + 0.29*(PetscReal)(map->comp[m]+1));
        loc2 += x[m]*x[m];
    }
    PetscCallMPI(MPI_Allreduce(&loc2, &glo2, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
    const PetscReal invn = 1.0/PetscSqrtReal(glo2);
    for (PetscInt m = 0; m < map->n; m++) x[m] *= invn;
    PetscCall(AddActiveVector(user, V, map, x, 1.0));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FourStage()

static PetscErrorCode FourStage ( UserCtx *  user, Vec  U0full, PetscReal  dtau, Vec  Rhs, const DofMap *  map, PetscReal *  Rscratch, Vec  Uwork, Vec  Uout  )

static

Definition at line 1363 of file test_momentum_convective_candidates.c.

{
    const PetscReal alfa[4] = {0.25, 1.0/3.0, 0.5, 1.0};
    PetscFunctionBeginUser;
    PetscCall(VecCopy(U0full, Uwork));                 /* stage state */
    for (int s = 0; s < 4; s++) {
        PetscCall(EvalConvResidual(user, Uwork, Rhs, map, Rscratch));   /* R(U^{(s)}) into active rows */
        PetscCall(VecCopy(U0full, Uout));               /* U^{(s+1)} = U0 + alfa*dtau*R */
        Cmpnts ***a;
        PetscCall(DMDAVecGetArray(user->fda, Uout, &a));
        for (PetscInt m = 0; m < map->n; m++) {
            PetscReal *p = (PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]];
            p[map->comp[m]] += alfa[s]*dtau*Rscratch[m];
        }
        PetscCall(DMDAVecRestoreArray(user->fda, Uout, &a));
    PetscCall(VecCopy(Uout, Uwork));
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnchoredStage()

static PetscErrorCode SetAnchoredStage ( UserCtx *  user, Vec  U0full, PetscReal  scale, const DofMap *  map, const PetscReal *  Ract, Vec  Ustage  )

static

Forms one anchored RK stage state from the base state and active residual.

Definition at line 1387 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    PetscCall(VecCopy(U0full, Ustage));
    PetscCall(AddActiveVector(user, Ustage, map, Ract, scale));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildStageStates()

static PetscErrorCode BuildStageStates ( UserCtx *  user, Vec  U0full, PetscReal  dtau, Vec  Rhs, const DofMap *  map, PetscReal *  Rscratch, Vec  Y1, Vec  Y2, Vec  Y3  )

static

Builds the first three anchored RK stage states for a base vector.

Definition at line 1399 of file test_momentum_convective_candidates.c.

{
    const PetscReal alfa[3] = {0.25, 1.0/3.0, 0.5};
    PetscFunctionBeginUser;
    PetscCall(EvalConvResidual(user, U0full, Rhs, map, Rscratch));
    PetscCall(SetAnchoredStage(user, U0full, alfa[0]*dtau, map, Rscratch, Y1));
    PetscCall(EvalConvResidual(user, Y1, Rhs, map, Rscratch));
    PetscCall(SetAnchoredStage(user, U0full, alfa[1]*dtau, map, Rscratch, Y2));
    PetscCall(EvalConvResidual(user, Y2, Rhs, map, Rscratch));
    PetscCall(SetAnchoredStage(user, U0full, alfa[2]*dtau, map, Rscratch, Y3));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildFDJacobian()

static PetscErrorCode BuildFDJacobian ( UserCtx *  user, Vec  Ucenter, PetscReal  epsrel, Vec  Rhs, const DofMap *  map, PetscReal *  Rp, PetscReal *  Rm, Vec  Uwork, PetscReal *  J  )

static

Builds a centered finite-difference Jacobian of the production convective residual.

Definition at line 1417 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    for (PetscInt col = 0; col < map->n; col++) {
        PetscReal u0;
        PetscCall(GetDof(user, Ucenter, map, col, &u0));
        const PetscReal eps = epsrel*PetscMax(1.0, PetscAbsReal(u0));
        PetscCall(VecCopy(Ucenter, Uwork));
        PetscCall(PerturbDof(user, Uwork, map, col, +eps));
        PetscCall(EvalConvResidual(user, Uwork, Rhs, map, Rp));
        PetscCall(VecCopy(Ucenter, Uwork));
        PetscCall(PerturbDof(user, Uwork, map, col, -eps));
        PetscCall(EvalConvResidual(user, Uwork, Rhs, map, Rm));
        for (PetscInt row = 0; row < map->n; row++) J[row + col*map->n] = (Rp[row]-Rm[row])/(2.0*eps);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CheckDenseJacobianAction()

static PetscErrorCode CheckDenseJacobianAction ( UserCtx *  user, Vec  Ubase, Vec  Rhs, Vec  Uwork, const DofMap *  map, const PetscReal *  J, PetscReal  epsrel  )

static

Verifies that the assembled dense Jacobian has the same action as production FD Jv.

Definition at line 1472 of file test_momentum_convective_candidates.c.

{
    Vec Up, Um;
    PetscReal *v, *jd, *Rp, *Rm, *jmf;
    PetscFunctionBeginUser;
    PetscCall(VecDuplicate(Ubase, &Up));
    PetscCall(VecDuplicate(Ubase, &Um));
    PetscCall(PetscMalloc5(map->n, &v, map->n, &jd, map->n, &Rp, map->n, &Rm, map->n, &jmf));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- dense Jacobian action check against matrix-free production Jv ---\n"
        "    direction        rel_L2       rel_Linf      ||Jv||2      ||Jv||inf\n"));
    const char *names[5] = {"mixed-sin", "mixed-cos", "u-only", "v-only", "w-only"};
    for (PetscInt kind = 0; kind < 5; kind++) {
        FillAuditDirection(map, kind, v);
        DenseMatVec(J, v, jd, map->n);
        PetscCall(VecCopy(Ubase, Up));
        PetscCall(VecCopy(Ubase, Um));
        PetscCall(AddActiveVector(user, Up, map, v, +epsrel));
        PetscCall(AddActiveVector(user, Um, map, v, -epsrel));
        PetscCall(EvalConvResidual(user, Up, Rhs, map, Rp));
        PetscCall(EvalConvResidual(user, Um, Rhs, map, Rm));
        for (PetscInt m = 0; m < map->n; m++) jmf[m] = (Rp[m]-Rm[m])/(2.0*epsrel);
        PetscReal e2 = 0.0, einf = 0.0;
        for (PetscInt m = 0; m < map->n; m++) {
            const PetscReal d = jd[m] - jmf[m];
            e2 += d*d; einf = PetscMax(einf, PetscAbsReal(d));
        }
        const PetscReal n2 = VecNorm2Array(jmf, map->n);
        const PetscReal ni = VecNormInfArray(jmf, map->n);
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    %-13s %.3e    %.3e    %.6e  %.6e\n",
            names[kind], (double)(PetscSqrtReal(e2)/PetscMax(PETSC_MACHINE_EPSILON,n2)),
            (double)(einf/PetscMax(PETSC_MACHINE_EPSILON,ni)), (double)n2, (double)ni));
    }
    (void)Uwork;
    PetscCall(PetscFree5(v, jd, Rp, Rm, jmf));
    PetscCall(VecDestroy(&Up));
    PetscCall(VecDestroy(&Um));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MatrixFreeJv()

static PetscErrorCode MatrixFreeJv ( UserCtx *  user, Vec  Ubase, Vec  Rhs, const DofMap *  map, const PetscReal *  v, PetscReal  eps, PetscReal *  jv  )

static

Matrix-free production Jacobian-vector product on the active space.

Definition at line 1518 of file test_momentum_convective_candidates.c.

{
    Vec Up, Um;
    PetscReal *Rp, *Rm;
    PetscFunctionBeginUser;
    PetscCall(VecDuplicate(Ubase, &Up));
    PetscCall(VecDuplicate(Ubase, &Um));
    PetscCall(PetscMalloc2(map->n, &Rp, map->n, &Rm));
    PetscCall(VecCopy(Ubase, Up));
    PetscCall(VecCopy(Ubase, Um));
    PetscCall(AddActiveVector(user, Up, map, v, +eps));
    PetscCall(AddActiveVector(user, Um, map, v, -eps));
    PetscCall(EvalConvResidual(user, Up, Rhs, map, Rp));
    PetscCall(EvalConvResidual(user, Um, Rhs, map, Rm));
    for (PetscInt m = 0; m < map->n; m++) jv[m] = (Rp[m]-Rm[m])/(2.0*eps);
    PetscCall(PetscFree2(Rp, Rm));
    PetscCall(VecDestroy(&Up));
    PetscCall(VecDestroy(&Um));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DenseAdditivityError()

static PetscReal DenseAdditivityError ( const PetscReal *  A, const PetscReal *  B, const PetscReal *  C, PetscInt  n  )

static

Normalized Frobenius norm of A - (B+C).

Definition at line 1543 of file test_momentum_convective_candidates.c.

{
    PetscReal num = 0.0, den = 0.0;
    for (PetscInt t = 0; t < n*n; t++) {
        const PetscReal d = A[t] - B[t] - C[t];
        num += d*d; den += A[t]*A[t];
    }
    return PetscSqrtReal(num) / PetscMax(1.0, PetscSqrtReal(den));
}

Here is the caller graph for this function:

DenseCommutatorNorm()

static PetscReal DenseCommutatorNorm ( const PetscReal *  A, const PetscReal *  B, PetscInt  n  )

static

Computes C = A*B - B*A and returns its normalized Frobenius norm.

Definition at line 1557 of file test_momentum_convective_candidates.c.

{
    PetscReal num = 0.0;
    for (PetscInt c = 0; c < n; c++) {
        for (PetscInt r = 0; r < n; r++) {
            PetscReal ab = 0.0, ba = 0.0;
            for (PetscInt q = 0; q < n; q++) {
                ab += A[r + q*n] * B[q + c*n];
                ba += B[r + q*n] * A[q + c*n];
            }
            const PetscReal d = ab - ba;
            num += d*d;
        }
    }
    const PetscReal den = PetscMax(DenseFrobenius(A, n)*DenseFrobenius(B, n), PETSC_MACHINE_EPSILON);
    return PetscSqrtReal(num) / den;
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunStateADirectionalMechanismAudit()

static PetscErrorCode RunStateADirectionalMechanismAudit ( const PetscReal *  Jxy, UserCtx *  user_xy, Vec  Ubase_xy, Vec  Rhs_xy, const DofMap *  map_xy, PetscReal  epsrel  )

static

Builds Jx, Jy, Jxy and prints additivity/commutator/eigenpair diagnostics.

Definition at line 1578 of file test_momentum_convective_candidates.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; DofMap map;
    Vec Ux, Uy, Rhs, Uwork;
    PetscReal *Rp, *Rm, *Jx, *Jy, *Jsum, *vr, *vi, *jdr, *jdi, *jmfr, *jmfi;
    PetscReal repeat_err, maxdiv, lamr, lami;
    SeamDiagnostics seam;
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 4, 4, 4, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    PetscCall(ConfigureCandidateFixture(simCtx, user));
    PetscCall(DofMapBuild(user, &map));
    PetscCheck(map.n == map_xy->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "directional audit map size mismatch");
    PetscCall(VecDuplicate(user->Ucont, &Ux));
    PetscCall(VecDuplicate(user->Ucont, &Uy));
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(VecDuplicate(user->Ucont, &Uwork));
    PetscCall(PetscMalloc5(map.n, &Rp, map.n, &Rm, (size_t)map.n*map.n, &Jx,
                           (size_t)map.n*map.n, &Jy, (size_t)map.n*map.n, &Jsum));
    PetscCall(BuildBaseState(user, STATE_A_X, Ux, &repeat_err, &maxdiv, &seam));
    PetscCall(BuildFDJacobian(user, Ux, epsrel, Rhs, &map, Rp, Rm, Uwork, Jx));
    PetscCall(BuildBaseState(user, STATE_A_Y, Uy, &repeat_err, &maxdiv, &seam));
    PetscCall(BuildFDJacobian(user, Uy, epsrel, Rhs, &map, Rp, Rm, Uwork, Jy));
    for (PetscInt t = 0; t < map.n*map.n; t++) Jsum[t] = Jx[t] + Jy[t];
 
    const PetscReal add_err = DenseAdditivityError(Jxy, Jx, Jy, map.n);
    const PetscReal comm = DenseCommutatorNorm(Jx, Jy, map.n);
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- State A directional additivity / noncommutation audit ---\n"
        "    ||Jxy-(Jx+Jy)||F/max(1,||Jxy||F) = %.3e\n"
        "    ||Jx Jy - Jy Jx||F/(||Jx||F||Jy||F) = %.3e\n"
        "    spectra:\n", (double)add_err, (double)comm));
    PetscCall(PrintSpectrumSummary("Jx", Jx, map.n));
    PetscCall(PrintSpectrumSummary("Jy", Jy, map.n));
    PetscCall(PrintSpectrumSummary("Jx+Jy", Jsum, map.n));
    PetscCall(PrintSpectrumSummary("Jxy", Jxy, map.n));
 
    PetscCall(PetscMalloc6(map.n, &vr, map.n, &vi, map.n, &jdr, map.n, &jdi, map.n, &jmfr, map.n, &jmfi));
    PetscCall(DenseMaxRealRightEigenpair(Jxy, map.n, &lamr, &lami, vr, vi));
    DenseMatVec(Jxy, vr, jdr, map.n);
    DenseMatVec(Jxy, vi, jdi, map.n);
    PetscReal ed2 = 0.0, em2 = 0.0, v2 = 0.0;
    for (PetscInt m = 0; m < map.n; m++) {
        const PetscReal rr = jdr[m] - (lamr*vr[m] - lami*vi[m]);
        const PetscReal ri = jdi[m] - (lami*vr[m] + lamr*vi[m]);
        ed2 += rr*rr + ri*ri;
        v2 += vr[m]*vr[m] + vi[m]*vi[m];
    }
    PetscCall(MatrixFreeJv(user_xy, Ubase_xy, Rhs_xy, map_xy, vr, epsrel, jmfr));
    PetscCall(MatrixFreeJv(user_xy, Ubase_xy, Rhs_xy, map_xy, vi, epsrel, jmfi));
    for (PetscInt m = 0; m < map.n; m++) {
        const PetscReal rr = jmfr[m] - (lamr*vr[m] - lami*vi[m]);
        const PetscReal ri = jmfi[m] - (lami*vr[m] + lamr*vi[m]);
        em2 += rr*rr + ri*ri;
    }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- State A unstable eigenpair verification ---\n"
        "    lambda_max_real = %.12e%+.12ei\n"
        "    ||Jdense v-lambda v||2/||v||2 = %.3e\n"
        "    ||Jmf    v-lambda v||2/||v||2 = %.3e\n",
        (double)lamr, (double)lami, (double)(PetscSqrtReal(ed2)/PetscSqrtReal(v2)),
        (double)(PetscSqrtReal(em2)/PetscSqrtReal(v2))));
 
    PetscCall(PetscFree6(vr, vi, jdr, jdi, jmfr, jmfi));
    PetscCall(PetscFree5(Rp, Rm, Jx, Jy, Jsum));
    PetscCall(VecDestroy(&Ux)); PetscCall(VecDestroy(&Uy)); PetscCall(VecDestroy(&Rhs)); PetscCall(VecDestroy(&Uwork));
    PetscCall(DofMapDestroy(&map));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ExtractLocalActiveVector()

static PetscErrorCode ExtractLocalActiveVector ( UserCtx *  user, Vec  local, const DofMap *  map, PetscReal *  x  )

static

Extract active-space entries from a local vector.

Definition at line 1660 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***a;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->fda, local, &a));
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscReal *p = (const PetscReal*)&a[map->ck[m]][map->cj[m]][map->ci[m]];
        x[m] = p[map->comp[m]];
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, local, &a));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

MapCartesianResidualToRct()

static PetscErrorCode MapCartesianResidualToRct ( UserCtx *  user, Vec  Rc, Vec  Rct  )

static

Mirrors ComputeRHS's Cartesian residual -> contravariant local Rct mapping.

Definition at line 1676 of file test_momentum_convective_candidates.c.

{
    DMDALocalInfo info = user->info;
    const PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
    const PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
    const PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
    const PetscInt lxs = (xs==0) ? xs+1 : xs, lxe = (xe==mx) ? xe-1 : xe;
    const PetscInt lys = (ys==0) ? ys+1 : ys, lye = (ye==my) ? ye-1 : ye;
    const PetscInt lzs = (zs==0) ? zs+1 : zs, lze = (ze==mz) ? ze-1 : ze;
    Cmpnts ***csi, ***eta, ***zet, ***rc, ***rct;
    PetscReal ***aj;
    PetscFunctionBeginUser;
    PetscCall(VecSet(Rct, 0.0));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecGetArrayRead(user->fda, Rc, &rc));
    PetscCall(DMDAVecGetArray(user->fda, Rct, &rct));
    for (PetscInt k = lzs; k < lze; k++)
        for (PetscInt j = lys; j < lye; j++)
            for (PetscInt i = lxs; i < lxe; i++) {
                rct[k][j][i].x = aj[k][j][i] *
                    (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
                     0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
                     0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
                rct[k][j][i].y = aj[k][j][i] *
                    (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
                     0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
                     0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
                rct[k][j][i].z = aj[k][j][i] *
                    (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
                     0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
                     0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
            }
    PetscCall(DMDAVecRestoreArray(user->fda, Rct, &rct));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, Rc, &rc));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(SynchronizePeriodicLocalStaggeredField(user, Rct));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

AverageRctToFinalActive()

static PetscErrorCode AverageRctToFinalActive ( UserCtx *  user, Vec  Rct, const DofMap *  map, PetscReal *  out  )

static

Mirrors ComputeRHS's final Rct face averaging and active cleanup for P=0/no body force.

Definition at line 1724 of file test_momentum_convective_candidates.c.

{
    Cmpnts ***rct;
    PetscFunctionBeginUser;
    PetscCall(DMDAVecGetArrayRead(user->fda, Rct, &rct));
    for (PetscInt m = 0; m < map->n; m++) {
        const PetscInt i = map->ci[m], j = map->cj[m], k = map->ck[m], c = map->comp[m];
        if (c == 0) out[m] = 0.5*(rct[k][j][i].x + rct[k][j][i+1].x);
        else if (c == 1) out[m] = 0.5*(rct[k][j][i].y + rct[k][j+1][i].y);
        else out[m] = 0.5*(rct[k][j][i].z + rct[k+1][j][i].z);
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, Rct, &rct));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

EvalResidualStage()

static PetscErrorCode EvalResidualStage ( UserCtx *  user, Vec  Ucont_in, Vec  Rhs, const DofMap *  map, ResidualStage  stage, PetscReal *  out  )

static

Evaluates one observable/mirrored residual-path stage for a Ucont input.

Definition at line 1742 of file test_momentum_convective_candidates.c.

{
    Vec Conv, Rc, Rct;
    PetscFunctionBeginUser;
    if (stage == STAGE_FINAL) {
        PetscCall(EvalConvResidual(user, Ucont_in, Rhs, map, out));
        PetscFunctionReturn(0);
    }
    PetscCall(VecCopy(Ucont_in, user->Ucont));
    {
        const char *fld[] = {"Ucont"};
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fld));
    }
    PetscCall(Contra2Cart(user));
    PetscCall(UpdateLocalGhosts(user, "Ucat"));
    if (stage == STAGE_UCAT) {
        PetscCall(ExtractActiveVector(user, user->Ucat, map, out));
        PetscFunctionReturn(0);
    }
    PetscCall(VecDuplicate(user->lUcont, &Conv));
    PetscCall(VecDuplicate(user->lUcont, &Rc));
    PetscCall(VecDuplicate(user->lUcont, &Rct));
    PetscCall(Convection(user, user->lUcont, user->lUcat, Conv));
    PetscCall(VecSet(Rc, 0.0));
    PetscCall(VecAXPY(Rc, -1.0, Conv));
    if (stage == STAGE_CART_RESID) PetscCall(ExtractLocalActiveVector(user, Rc, map, out));
    else {
        PetscCall(MapCartesianResidualToRct(user, Rc, Rct));
        if (stage == STAGE_RCT) PetscCall(ExtractLocalActiveVector(user, Rct, map, out));
        else PetscCall(AverageRctToFinalActive(user, Rct, map, out));
    }
    PetscCall(VecDestroy(&Conv));
    PetscCall(VecDestroy(&Rc));
    PetscCall(VecDestroy(&Rct));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildStageJacobian()

static PetscErrorCode BuildStageJacobian ( UserCtx *  user, Vec  Ucenter, PetscReal  epsrel, Vec  Rhs, const DofMap *  map, ResidualStage  stage, PetscReal *  Rp, PetscReal *  Rm, Vec  Uwork, PetscReal *  J  )

static

Builds a finite-difference Jacobian for one residual-path stage.

Definition at line 1783 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    for (PetscInt col = 0; col < map->n; col++) {
        PetscReal u0;
        PetscCall(GetDof(user, Ucenter, map, col, &u0));
        const PetscReal eps = epsrel*PetscMax(1.0, PetscAbsReal(u0));
        PetscCall(VecCopy(Ucenter, Uwork));
        PetscCall(PerturbDof(user, Uwork, map, col, +eps));
        PetscCall(EvalResidualStage(user, Uwork, Rhs, map, stage, Rp));
        PetscCall(VecCopy(Ucenter, Uwork));
        PetscCall(PerturbDof(user, Uwork, map, col, -eps));
        PetscCall(EvalResidualStage(user, Uwork, Rhs, map, stage, Rm));
        for (PetscInt row = 0; row < map->n; row++) J[row + col*map->n] = (Rp[row]-Rm[row])/(2.0*eps);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunStateAResidualPathIsolation()

static PetscErrorCode RunStateAResidualPathIsolation ( UserCtx *  user, Vec  Ubase, Vec  Rhs, Vec  Uwork, const DofMap *  map, PetscReal  epsrel  )

static

Prints the stage where positive-real spectrum first appears.

Definition at line 1806 of file test_momentum_convective_candidates.c.

{
    const char *names[4] = {"Contra2Cart Ucat", "Cartesian Rc=-Conv", "mapped local Rct", "final active RHS"};
    PetscReal *Rp, *Rm, *J;
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc3(map->n, &Rp, map->n, &Rm, (size_t)map->n*map->n, &J));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- State A residual-path stage Jacobians ---\n"
        "    stage                 rho          sigma        max_real     skew       nonnormal\n"));
    for (PetscInt s = 0; s < 4; s++) {
        PetscReal rho, maxre, smax;
        PetscCall(BuildStageJacobian(user, Ubase, epsrel, Rhs, map, (ResidualStage)s, Rp, Rm, Uwork, J));
        PetscCall(DenseSpectralRadius(J, map->n, &rho, &maxre));
        PetscCall(DenseSigmaMax(J, map->n, &smax, NULL));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    %-20s %.6e  %.6e  %.6e  %.3e  %.3e\n",
            names[s], (double)rho, (double)smax, (double)maxre,
            (double)DenseSkewnessDefect(J, map->n), (double)DenseNonNormality(J, map->n)));
    }
    PetscCall(PetscFree3(Rp, Rm, J));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildStageTangent()

static PetscErrorCode BuildStageTangent ( const PetscReal *  J0, const PetscReal *  J1, const PetscReal *  J2, const PetscReal *  J3, PetscReal  dtau, PetscInt  n, PetscReal *  T4  )

static

Builds the exact four-stage tangent from stage-dependent Jacobians.

Definition at line 1833 of file test_momentum_convective_candidates.c.

{
    const PetscReal alfa[4] = {0.25, 1.0/3.0, 0.5, 1.0};
    PetscReal *T1, *T2, *T3, *Tmp;
    PetscFunctionBeginUser;
    PetscCall(PetscMalloc4(n*n, &T1, n*n, &T2, n*n, &T3, n*n, &Tmp));
 
    for (PetscInt i = 0; i < n*n; i++) T1[i] = alfa[0]*dtau*J0[i];
    for (PetscInt d = 0; d < n; d++) T1[d + d*n] += 1.0;
 
    MatMul(J1, T1, Tmp, n);
    for (PetscInt i = 0; i < n*n; i++) T2[i] = alfa[1]*dtau*Tmp[i];
    for (PetscInt d = 0; d < n; d++) T2[d + d*n] += 1.0;
 
    MatMul(J2, T2, Tmp, n);
    for (PetscInt i = 0; i < n*n; i++) T3[i] = alfa[2]*dtau*Tmp[i];
    for (PetscInt d = 0; d < n; d++) T3[d + d*n] += 1.0;
 
    MatMul(J3, T3, Tmp, n);
    for (PetscInt i = 0; i < n*n; i++) T4[i] = alfa[3]*dtau*Tmp[i];
    for (PetscInt d = 0; d < n; d++) T4[d + d*n] += 1.0;
 
    PetscCall(PetscFree4(T1, T2, T3, Tmp));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BuildPhiJacobian()

static PetscErrorCode BuildPhiJacobian ( UserCtx *  user, Vec  Ucenter, PetscReal  dtau, Vec  Rhs, const DofMap *  map, PetscReal *  Rscratch, Vec  Upert, Vec  Ustage, Vec  PhiP, Vec  PhiM, PetscReal *  xp, PetscReal *  xm, PetscReal  epsrel, PetscReal *  JPhi  )

static

Builds a finite-difference Jacobian of the complete nonlinear four-stage map.

Definition at line 1864 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    for (PetscInt col = 0; col < map->n; col++) {
        PetscReal u0;
        PetscCall(GetDof(user, Ucenter, map, col, &u0));
        const PetscReal eps = epsrel*PetscMax(1.0, PetscAbsReal(u0));
        PetscCall(VecCopy(Ucenter, Upert));
        PetscCall(PerturbDof(user, Upert, map, col, +eps));
        PetscCall(FourStage(user, Upert, dtau, Rhs, map, Rscratch, Ustage, PhiP));
        PetscCall(VecCopy(Ucenter, Upert));
        PetscCall(PerturbDof(user, Upert, map, col, -eps));
        PetscCall(FourStage(user, Upert, dtau, Rhs, map, Rscratch, Ustage, PhiM));
        PetscCall(ExtractActiveVector(user, PhiP, map, xp));
        PetscCall(ExtractActiveVector(user, PhiM, map, xm));
        for (PetscInt row = 0; row < map->n; row++) JPhi[row + col*map->n] = (xp[row]-xm[row])/(2.0*eps);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunRKTangentDiagnostics()

static PetscErrorCode RunRKTangentDiagnostics ( UserCtx *  user, CandState  st, Vec  Ubase, Vec  Rhs, const DofMap *  map, PetscReal  epsrel, const PetscReal  lams[3], const char *  cn[3], const PetscReal *  J0  )

static

Runs stage-dependent RK tangent and direct nonlinear perturbation diagnostics.

Definition at line 1891 of file test_momentum_convective_candidates.c.

{
    const PetscReal cflsB[4] = {0.1, 0.25, 0.5, 1.0};
    const PetscReal cflsOther[1] = {0.5};
    const PetscReal *cfls = (st == STATE_B) ? cflsB : cflsOther;
    const PetscInt ncfl = (st == STATE_B) ? 4 : 1;
    const PetscReal amps[3] = {1e-4, 1e-5, 1e-6};
    Vec Y1, Y2, Y3, Upert, Ustage, PhiP, PhiM, Phi0;
    PetscReal *Rtmp, *J1, *J2, *J3, *T4, *Pm, *JPhi, *v1, *xrand;
    PetscReal *meas, *phi0, *phip, *xscaled, *predP, *predT;
    PetscFunctionBeginUser;
 
    PetscCall(VecDuplicate(Ubase, &Y1));
    PetscCall(VecDuplicate(Ubase, &Y2));
    PetscCall(VecDuplicate(Ubase, &Y3));
    PetscCall(VecDuplicate(Ubase, &Upert));
    PetscCall(VecDuplicate(Ubase, &Ustage));
    PetscCall(VecDuplicate(Ubase, &PhiP));
    PetscCall(VecDuplicate(Ubase, &PhiM));
    PetscCall(VecDuplicate(Ubase, &Phi0));
    PetscCall(PetscMalloc5(map->n, &Rtmp, (size_t)map->n*map->n, &J1,
                           (size_t)map->n*map->n, &J2, (size_t)map->n*map->n, &J3,
                           (size_t)map->n*map->n, &T4));
    PetscCall(PetscMalloc5((size_t)map->n*map->n, &Pm, (size_t)map->n*map->n, &JPhi,
                           map->n, &v1, map->n, &xrand, map->n, &meas));
    PetscCall(PetscMalloc5(map->n, &phi0, map->n, &phip, map->n, &xscaled,
                           map->n, &predP, map->n, &predT));
 
    PetscReal nrm = 0.0;
    for (PetscInt m = 0; m < map->n; m++) {
        xrand[m] = PetscSinReal((PetscReal)(m+1)*1.2345);
        nrm += xrand[m]*xrand[m];
    }
    nrm = PetscSqrtReal(nrm);
    for (PetscInt m = 0; m < map->n; m++) xrand[m] /= nrm;
 
    for (int cand = 0; cand < 3; cand++) {
      if (!(lams[cand] > 0.0)) continue;
      for (PetscInt icfl = 0; icfl < ncfl; icfl++) {
        const PetscReal cfl = cfls[icfl], dtau = cfl/lams[cand];
        PetscCall(BuildStageStates(user, Ubase, dtau, Rhs, map, Rtmp, Y1, Y2, Y3));
        PetscCall(BuildFDJacobian(user, Y1, epsrel, Rhs, map, phi0, phip, Upert, J1));
        PetscCall(BuildFDJacobian(user, Y2, epsrel, Rhs, map, phi0, phip, Upert, J2));
        PetscCall(BuildFDJacobian(user, Y3, epsrel, Rhs, map, phi0, phip, Upert, J3));
        PetscCall(BuildStageTangent(J0, J1, J2, J3, dtau, map->n, T4));
        PetscCall(RKPolynomial(J0, dtau, map->n, Pm));
        PetscCall(BuildPhiJacobian(user, Ubase, dtau, Rhs, map, Rtmp,
                                   Upert, Ustage, PhiP, PhiM, phip, phi0, epsrel, JPhi));
 
        const PetscReal froT = DenseFrobenius(T4, map->n);
        const PetscReal froPhi = DenseFrobenius(JPhi, map->n);
        const PetscReal relTP = DenseRelativeDiff(T4, Pm, map->n, froT);
        const PetscReal relPhiT = DenseRelativeDiff(JPhi, T4, map->n, froPhi);
        const PetscReal relPhiP = DenseRelativeDiff(JPhi, Pm, map->n, froPhi);
        PetscReal smaxT, rhoT, dummyT;
        PetscCall(DenseSpectralRadius(T4, map->n, &rhoT, &dummyT));
        PetscCall(DenseSigmaMax(T4, map->n, &smaxT, v1));
 
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "  --- convection-only stage-dependent RK tangent (cand %s, CFL=%.2f, dtau=%.6e) ---\n"
            "    rho(T4)=%.6e  sigma_max(T4)=%.6e\n"
            "    ||T4-P(hJ0)||F/max(1,||T4||F)      = %.3e\n"
            "    ||J_Phi-T4||F/max(1,||J_Phi||F)    = %.3e\n"
            "    ||J_Phi-P(hJ0)||F/max(1,||J_Phi||F)= %.3e\n",
            cn[cand], (double)cfl, (double)dtau, (double)rhoT, (double)smaxT,
            (double)relTP, (double)relPhiT, (double)relPhiP));
 
        PetscCall(FourStage(user, Ubase, dtau, Rhs, map, Rtmp, Ustage, Phi0));
        PetscCall(ExtractActiveVector(user, Phi0, map, phi0));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    direction  amp    amp(meas)  amp(frozen)  amp(stage)  err(frozen)  err(stage)\n"));
        for (int dir = 0; dir < 2; dir++) {
            const PetscReal *xv = (dir == 0) ? xrand : v1;
            const char *dname = (dir == 0) ? "random" : "v1(T4)";
            for (int a = 0; a < 3; a++) {
                const PetscReal amp = amps[a];
                PetscCall(VecCopy(Ubase, Upert));
                PetscCall(AddActiveVector(user, Upert, map, xv, amp));
                PetscCall(FourStage(user, Upert, dtau, Rhs, map, Rtmp, Ustage, PhiP));
                PetscCall(ExtractActiveVector(user, PhiP, map, phip));
                for (PetscInt m = 0; m < map->n; m++) meas[m] = phip[m] - phi0[m];
                for (PetscInt m = 0; m < map->n; m++) xscaled[m] = amp*xv[m];
                ApplyP(Pm, xscaled, predP, map->n);
                ApplyP(T4, xscaled, predT, map->n);
 
                PetscReal eP = 0.0, eT = 0.0;
                for (PetscInt m = 0; m < map->n; m++) {
                    const PetscReal dP = meas[m] - predP[m];
                    const PetscReal dT = meas[m] - predT[m];
                    eP += dP*dP; eT += dT*dT;
                }
                const PetscReal nm = VecNorm2Array(meas, map->n);
                const PetscReal nP = VecNorm2Array(predP, map->n);
                const PetscReal nT = VecNorm2Array(predT, map->n);
                PetscCall(PetscPrintf(PETSC_COMM_WORLD,
                    "    %-8s %.0e  %.6e  %.6e  %.6e  %.3e  %.3e\n",
                    dname, (double)amp, (double)(nm/amp), (double)(nP/amp), (double)(nT/amp),
                    (double)(PetscSqrtReal(eP)/PetscMax(PETSC_MACHINE_EPSILON, nP)),
                    (double)(PetscSqrtReal(eT)/PetscMax(PETSC_MACHINE_EPSILON, nT))));
            }
        }
 
        PetscCall(PicurvAssertBool((PetscBool)(relPhiT < 1e-6),
                                   "stage tangent matches direct finite-difference RK map"));
        if (st == STATE_B) {
            PetscCall(PicurvAssertBool((PetscBool)(relTP > 1e-5),
                                       "B: non-steady base makes frozen-Jacobian RK map measurably different"));
            PetscCall(PicurvAssertBool((PetscBool)(relPhiP > 1e-5),
                                       "B: direct RK map confirms frozen-Jacobian error"));
        } else {
            PetscCall(PicurvAssertBool((PetscBool)(relTP < 1e-8),
                                       "steady base reduces stage tangent to frozen RK polynomial"));
        }
      }
    }
 
    PetscCall(PetscFree5(Rtmp, J1, J2, J3, T4));
    PetscCall(PetscFree5(Pm, JPhi, v1, xrand, meas));
    PetscCall(PetscFree5(phi0, phip, xscaled, predP, predT));
    PetscCall(VecDestroy(&Y1)); PetscCall(VecDestroy(&Y2)); PetscCall(VecDestroy(&Y3));
    PetscCall(VecDestroy(&Upert)); PetscCall(VecDestroy(&Ustage));
    PetscCall(VecDestroy(&PhiP)); PetscCall(VecDestroy(&PhiM)); PetscCall(VecDestroy(&Phi0));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunState()

static PetscErrorCode RunState ( CandState  st, const char *  name  )

static

Definition at line 2022 of file test_momentum_convective_candidates.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; DofMap map;
    Vec Ubase, Rhs, Uwork, Uout;
    PetscReal *Rref, *Rrep, *Jbest;
    PetscReal repeat_err, maxdiv, det_err;
    SeamDiagnostics seam;
    MomStabilityReport rep;
    const PetscInt N = (st == STATE_C) ? 5 : 4;   /* C uses one extra point for canonical shear. */
    PetscFunctionBeginUser;
 
    /* periodic Cartesian fixture, inviscid + centered + P=0, no LES/RANS/Clark/IB/body force. */
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, N, N, N, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    PetscCall(ConfigureCandidateFixture(simCtx, user));
 
    PetscCall(DofMapBuild(user, &map));
    PetscCall(VecDuplicate(user->Ucont, &Ubase));
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(VecDuplicate(user->Ucont, &Uwork));
    PetscCall(VecDuplicate(user->Ucont, &Uout));
    PetscCall(PetscMalloc3(map.n, &Rref, map.n, &Rrep, (size_t)map.n*map.n, &Jbest));
 
    PetscCall(BuildBaseState(user, st, Ubase, &repeat_err, &maxdiv, &seam));
 
    /* residual determinism: evaluate the base residual twice. */
    PetscCall(EvalConvResidual(user, Ubase, Rhs, &map, Rref));
    PetscCall(EvalConvResidual(user, Ubase, Rhs, &map, Rrep));
    det_err = 0.0; for (PetscInt m = 0; m < map.n; m++) det_err = PetscMax(det_err, PetscAbsReal(Rref[m]-Rrep[m]));
    const PetscReal R0_2 = VecNorm2Array(Rref, map.n);
    const PetscReal R0_inf = VecNormInfArray(Rref, map.n);
 
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "\n================ STATE %s ================\n"
        "  grid: DMDA %d^3 (periodic) | independent face DOFs: actual=%d expected=%d\n"
        "  declared endpoint mismatch       (x,y,z) = %.3e %.3e %.3e\n"
        "  actual Ucat duplicate mismatch   (x,y,z) = %.3e %.3e %.3e\n"
        "  local lUcat ghost mismatch       (x,y,z) = %.3e %.3e %.3e\n"
        "  actual Ucont duplicate mismatch  (x,y,z) = %.3e %.3e %.3e\n"
        "  ||Ucat_repeat - Ucat_reference||inf = %.3e\n"
        "  max|div_h Ucont|                    = %.3e\n"
        "  ||R(U0)||2 = %.6e  ||R(U0)||inf = %.6e\n"
        "  residual determinism ||R_rep-R_ref||inf = %.3e\n",
        name, (int)(N+1), (int)map.n,
        (int)map.expected_n,
        (double)seam.declared[0], (double)seam.declared[1], (double)seam.declared[2],
        (double)seam.ucat_global[0], (double)seam.ucat_global[1], (double)seam.ucat_global[2],
        (double)seam.ucat_ghost[0], (double)seam.ucat_ghost[1], (double)seam.ucat_ghost[2],
        (double)seam.ucont_global[0], (double)seam.ucont_global[1], (double)seam.ucont_global[2],
        (double)repeat_err, (double)maxdiv,
        (double)R0_2, (double)R0_inf, (double)det_err));
    if (st == STATE_A) PetscCall(PrintPeriodicSpaceAudit(user, &map));
 
    /* ---- epsilon-convergence study: build J at several eps, compare to next-finer ---- */
    const PetscReal epsrel[5] = {1e-4, 1e-5, 1e-6, 1e-7, 1e-8};
    PetscReal *Jprev, *Jcur;
    PetscCall(PetscMalloc2((size_t)map.n*map.n, &Jprev, (size_t)map.n*map.n, &Jcur));
    PetscReal best_rel = PETSC_MAX_REAL; PetscInt best_e = 2;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  epsilon-convergence (||J_e - J_e/10||_F / ||J||_F):\n"));
    for (int e = 0; e < 5; e++) {
        /* build column-major J at epsrel[e] */
        for (PetscInt col = 0; col < map.n; col++) {
            PetscReal u0; PetscCall(GetDof(user, Ubase, &map, col, &u0));
            const PetscReal eps = epsrel[e]*PetscMax(1.0, PetscAbsReal(u0));
            PetscReal *Rp = Rref, *Rm = Rrep;   /* reuse scratch */
            PetscCall(VecCopy(Ubase, Uwork));
            PetscCall(PerturbDof(user, Uwork, &map, col, +eps));
            PetscCall(EvalConvResidual(user, Uwork, Rhs, &map, Rp));
            PetscCall(VecCopy(Ubase, Uwork));
            PetscCall(PerturbDof(user, Uwork, &map, col, -eps));
            PetscCall(EvalConvResidual(user, Uwork, Rhs, &map, Rm));
            for (PetscInt row = 0; row < map.n; row++) Jcur[row + col*map.n] = (Rp[row]-Rm[row])/(2.0*eps);
        }
        if (e > 0) {
            PetscReal num = 0.0, den = 0.0;
            for (PetscInt t = 0; t < map.n*map.n; t++) { const PetscReal d = Jprev[t]-Jcur[t]; num += d*d; den += Jcur[t]*Jcur[t]; }
            const PetscReal rel = PetscSqrtReal(num)/PetscMax(1.0, PetscSqrtReal(den));
            PetscCall(PetscPrintf(PETSC_COMM_WORLD, "    eps=%.0e -> %.0e : rel=%.3e\n",
                                  (double)epsrel[e-1], (double)epsrel[e], (double)rel));
            if (rel < best_rel) { best_rel = rel; best_e = e; }
        }
        if (st == STATE_A) {
            PetscReal erho, emaxre, esmax;
            PetscCall(DenseSpectralRadius(Jcur, map.n, &erho, &emaxre));
            PetscCall(DenseSigmaMax(Jcur, map.n, &esmax, NULL));
            PetscCall(PetscPrintf(PETSC_COMM_WORLD,
                "    eps=%.0e metrics: rho=%.6e sigma=%.6e max_real=%.6e fro=%.6e skew=%.3e\n",
                (double)epsrel[e], (double)erho, (double)esmax, (double)emaxre,
                (double)DenseFrobenius(Jcur, map.n), (double)DenseSkewnessDefect(Jcur, map.n)));
        }
        PetscCall(PetscArraycpy(Jprev, Jcur, (size_t)map.n*map.n));   /* prev <- current (no pointer swap) */
    }
    /* rebuild J at the plateau epsilon into Jbest */
    {
        const PetscReal er = epsrel[best_e];
        for (PetscInt col = 0; col < map.n; col++) {
            PetscReal u0; PetscCall(GetDof(user, Ubase, &map, col, &u0));
            const PetscReal eps = er*PetscMax(1.0, PetscAbsReal(u0));
            PetscCall(VecCopy(Ubase, Uwork)); PetscCall(PerturbDof(user, Uwork, &map, col, +eps));
            PetscCall(EvalConvResidual(user, Uwork, Rhs, &map, Rref));
            PetscCall(VecCopy(Ubase, Uwork)); PetscCall(PerturbDof(user, Uwork, &map, col, -eps));
            PetscCall(EvalConvResidual(user, Uwork, Rhs, &map, Rrep));
            for (PetscInt row = 0; row < map.n; row++) Jbest[row + col*map.n] = (Rref[row]-Rrep[row])/(2.0*eps);
        }
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  selected plateau eps = %.0e\n", (double)er));
    }
    PetscCall(PetscFree2(Jprev, Jcur));
 
    /* The FD probes leave the production vectors at the final perturbed evaluation;
       restore the base state before every downstream diagnostic. */
    PetscCall(VecCopy(Ubase, user->Ucont));
    {
        const char *fld[] = {"Ucont"};
        PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fld));
    }
 
    /* base state must be restored after the Jacobian build. */
    PetscBool eqbase; { Vec chk; PetscCall(VecDuplicate(Ubase,&chk)); PetscCall(VecCopy(user->Ucont,chk));
        PetscCall(VecAXPY(chk, -1.0, Ubase)); PetscReal nb; PetscCall(VecNorm(chk, NORM_INFINITY, &nb));
        eqbase = (PetscBool)(nb < 1e-12); PetscCall(VecDestroy(&chk));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  base Ucont restored after Jacobian: %s\n", eqbase?"yes":"NO")); }
 
    /* The FD operator acts on Ucont; J is on the convective residual (sign per ComputeRHS). */
    PetscReal rho, maxre, smax;
    PetscCall(DenseSpectralRadius(Jbest, map.n, &rho, &maxre));
    PetscCall(DenseSigmaMax(Jbest, map.n, &smax, NULL));
    const PetscReal eta = DenseNonNormality(Jbest, map.n);
 
    /* ---- candidate estimates (convection-only: lambda_cX = lambda_X - lambda_t) ---- */
    PetscCall(EvalConvResidual(user, Ubase, Rhs, &map, Rref));   /* restore lUcat consistent w/ base */
    PetscCall(UpdateLocalGhosts(user, "Ucont"));
    PetscCall(Contra2Cart(user)); PetscCall(UpdateLocalGhosts(user, "Ucat"));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
    const PetscReal lcB = rep.lambda_B - rep.lambda_t;
    const PetscReal lcC = rep.lambda_C - rep.lambda_t;
    const PetscReal lcD = rep.lambda_D - rep.lambda_t;
    PetscReal gradmax = 0.0;
    PetscCall(ComputeMaxGradientContribution(user, &gradmax));
 
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  --- convective Jacobian spectrum ---\n"
        "    rho(J)=%.6e  sigma_max(J)=%.6e  nonnormality=%.3e  max_real_eig=%.3e\n"
        "  --- candidate convective estimates ---\n"
        "    lambda_cB=%.6e  lambda_cC=%.6e  lambda_cD=%.6e\n"
        "    max local |grad u| row-sum contribution = %.6e\n"
        "    rB/rho=%.3f rC/rho=%.3f rD/rho=%.3f | rB/smax=%.3f rC/smax=%.3f rD/smax=%.3f\n",
        (double)rho, (double)smax, (double)eta, (double)maxre,
        (double)lcB, (double)lcC, (double)lcD,
        (double)gradmax,
        (double)(lcB/rho), (double)(lcC/rho), (double)(lcD/rho),
        (double)(lcB/smax), (double)(lcC/smax), (double)(lcD/smax)));
    if (st == STATE_A) {
        PetscCall(CheckDenseJacobianAction(user, Ubase, Rhs, Uwork, &map, Jbest, epsrel[best_e]));
        PetscCall(RunStateADirectionalMechanismAudit(Jbest, user, Ubase, Rhs, &map, epsrel[best_e]));
        PetscCall(PrintBlockAndSymmetricLocalization(user, &map, Jbest));
        PetscCall(PrintFourierChecks(&map, Jbest));
        PetscCall(TraceStateAResidualStages(user, Ubase, Rhs, &map, Jbest));
        PetscCall(RunStateAResidualPathIsolation(user, Ubase, Rhs, Uwork, &map, epsrel[best_e]));
    }
 
    /* ---- RK pseudo-CFL stability per candidate (convective Jacobian) ---- */
    const PetscReal lams[3] = {lcB, lcC, lcD}; const char *cn[3] = {"B","C","D"};
    PetscReal *Jpseudo;
    PetscCall(PetscMalloc1((size_t)map.n*map.n, &Jpseudo));
    PetscCall(DenseShiftIdentity(Jbest, map.n, -rep.lambda_t, Jpseudo));
    const PetscReal lams_full[3] = {rep.lambda_B, rep.lambda_C, rep.lambda_D};
    PetscCall(PrintFrozenAmplificationTable("CONVECTION-ONLY FROZEN JACOBIAN", Jbest, map.n, lams, cn));
    PetscCall(PrintFrozenAmplificationTable("COMPLETE PSEUDO-TIME FROZEN OPERATOR", Jpseudo, map.n, lams_full, cn));
    PetscCall(PetscFree(Jpseudo));
 
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  --- RK stable pseudo-CFL (convective J) ---\n"));
    for (int c = 0; c < 3; c++) {
        if (!(lams[c] > 0.0)) continue;
        StableCFLResult cfl_rho, cfl_nrm;
        PetscCall(StableCFL(Jbest, map.n, lams[c], METRIC_RHO,  &cfl_rho));
        PetscCall(StableCFL(Jbest, map.n, lams[c], METRIC_NORM, &cfl_nrm));
        PetscCall(PrintStableCFLLine(cn[c], cfl_rho, cfl_nrm));
    }
 
    /* ---- exact stage-dependent tangent map and direct nonlinear RK cross-check ---- */
    if (lcC > 0.0) PetscCall(RunRKTangentDiagnostics(user, st, Ubase, Rhs, &map,
                                                     epsrel[best_e], lams, cn, Jbest));
 
    /* ---- robust automated assertions per state ---- */
    PetscCall(PicurvAssertBool(eqbase, "base Ucont restored after Jacobian build"));
    PetscCall(PicurvAssertBool((PetscBool)(det_err < 1e-10), "residual evaluation deterministic"));
    PetscCall(PicurvAssertBool((PetscBool)(rho > 0.0 && PetscIsNormalReal(rho)), "finite rho(J)"));
    PetscCall(PicurvAssertBool((PetscBool)(smax > 0.0 && PetscIsNormalReal(smax)), "finite sigma_max(J)"));
    if (st == STATE_A) {
        for (int d = 0; d < 3; d++) {
            PetscCall(PicurvAssertRealNear(seam.declared[d], 0.0, 1e-12, "A: declared periodic seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_global[d], 0.0, 1e-12, "A: Ucat duplicate seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_ghost[d], 0.0, 1e-12, "A: lUcat ghost seam"));
            PetscCall(PicurvAssertRealNear(seam.ucont_global[d], 0.0, 1e-12, "A: Ucont duplicate seam"));
        }
        PetscCall(PicurvAssertRealNear(repeat_err, 0.0, 1e-9, "A: recovered Ucat repeat near roundoff"));
        PetscCall(PicurvAssertRealNear(R0_2, 0.0, 1e-12, "A: base residual 2-norm near roundoff"));
        PetscCall(PicurvAssertRealNear(R0_inf, 0.0, 1e-12, "A: base residual inf-norm near roundoff"));
        PetscCall(PicurvAssertRealNear(lcB, lcC, 1e-9, "A: B == C (no divergence)"));
        PetscCall(PicurvAssertRealNear(lcC, lcD, 1e-9, "A: C == D (no shear)"));
    } else if (st == STATE_B) {
        for (int d = 0; d < 3; d++) {
            PetscCall(PicurvAssertRealNear(seam.declared[d], 0.0, 1e-12, "B: declared periodic seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_global[d], 0.0, 1e-12, "B: Ucat duplicate seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_ghost[d], 0.0, 1e-12, "B: lUcat ghost seam"));
            PetscCall(PicurvAssertRealNear(seam.ucont_global[d], 0.0, 1e-12, "B: Ucont duplicate seam"));
        }
        PetscCall(PicurvAssertRealNear(repeat_err, 0.0, 1e-9, "B: recovered Ucat repeat near roundoff"));
        PetscCall(PicurvAssertBool((PetscBool)(maxdiv > 1e-3), "B: nonzero discrete divergence"));
        PetscCall(PicurvAssertBool((PetscBool)(R0_2 > 1e-3), "B: materially nonzero base residual 2-norm"));
        PetscCall(PicurvAssertBool((PetscBool)(R0_inf > 1e-3), "B: materially nonzero base residual inf-norm"));
        PetscCall(PicurvAssertBool((PetscBool)(lcC > lcB + 1e-9), "B: C > B"));
    } else {
        for (int d = 0; d < 3; d++) {
            PetscCall(PicurvAssertRealNear(seam.declared[d], 0.0, 1e-12, "C: declared periodic seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_global[d], 0.0, 1e-12, "C: Ucat duplicate seam"));
            PetscCall(PicurvAssertRealNear(seam.ucat_ghost[d], 0.0, 1e-12, "C: lUcat ghost seam"));
            PetscCall(PicurvAssertRealNear(seam.ucont_global[d], 0.0, 1e-12, "C: Ucont duplicate seam"));
        }
        PetscCall(PicurvAssertRealNear(repeat_err, 0.0, 1e-9, "C: recovered Ucat repeat near roundoff"));
        PetscCall(PicurvAssertBool((PetscBool)(maxdiv < 1e-6), "C: divergence near zero"));
        PetscCall(PicurvAssertRealNear(R0_2, 0.0, 1e-12, "C: base residual 2-norm near roundoff"));
        PetscCall(PicurvAssertRealNear(R0_inf, 0.0, 1e-12, "C: base residual inf-norm near roundoff"));
        PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(lcC-lcB) < 1e-6), "C: C ~= B"));
        PetscCall(PicurvAssertBool((PetscBool)(gradmax > 1e-6), "C: canonical shear has nonzero local gradient contribution"));
    }
 
    PetscCall(PetscFree3(Rref, Rrep, Jbest));
    PetscCall(VecDestroy(&Ubase)); PetscCall(VecDestroy(&Rhs));
    PetscCall(VecDestroy(&Uwork)); PetscCall(VecDestroy(&Uout));
    PetscCall(DofMapDestroy(&map));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestStateA()

static PetscErrorCode TestStateA ( void  )

static

Runs the State A candidate harness.

Definition at line 2259 of file test_momentum_convective_candidates.c.

{ return RunState(STATE_A, "A (uniform div-free)"); }

Here is the call graph for this function:

Here is the caller graph for this function:

TestStateB()

static PetscErrorCode TestStateB ( void  )

static

Runs the State B candidate harness.

Definition at line 2264 of file test_momentum_convective_candidates.c.

{ return RunState(STATE_B, "B (nonzero divergence)"); }

Here is the call graph for this function:

Here is the caller graph for this function:

TestStateC()

static PetscErrorCode TestStateC ( void  )

static

Runs the State C candidate harness.

Definition at line 2269 of file test_momentum_convective_candidates.c.

{ return RunState(STATE_C, "C (div-free shear)"); }

Here is the call graph for this function:

Here is the caller graph for this function:

RunStateAGridAuditOne()

static PetscErrorCode RunStateAGridAuditOne ( PetscInt  N )

static

Runs one State A grid-size audit case.

Definition at line 2274 of file test_momentum_convective_candidates.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; DofMap map;
    Vec Ubase, Rhs, Uwork;
    PetscReal *Rp, *Rm, *J;
    PetscReal repeat_err, maxdiv;
    SeamDiagnostics seam;
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, N, N, N, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    PetscCall(ConfigureCandidateFixture(simCtx, user));
    PetscCall(DofMapBuild(user, &map));
    PetscCall(VecDuplicate(user->Ucont, &Ubase));
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(VecDuplicate(user->Ucont, &Uwork));
    PetscCall(PetscMalloc3(map.n, &Rp, map.n, &Rm, (size_t)map.n*map.n, &J));
    PetscCall(BuildBaseState(user, STATE_A, Ubase, &repeat_err, &maxdiv, &seam));
    PetscCall(BuildFDJacobian(user, Ubase, 1e-5, Rhs, &map, Rp, Rm, Uwork, J));
    PetscReal rho, maxre, smax;
    PetscCall(DenseSpectralRadius(J, map.n, &rho, &maxre));
    PetscCall(DenseSigmaMax(J, map.n, &smax, NULL));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  State A grid audit: DMDA %d^3 independent=%d expected=%d max_real=%.6e rho=%.6e sigma=%.6e skew=%.3e nonnormality=%.3e repeat=%.3e\n",
        (int)(N+1), (int)map.n, (int)map.expected_n, (double)maxre, (double)rho, (double)smax,
        (double)DenseSkewnessDefect(J, map.n), (double)DenseNonNormality(J, map.n), (double)repeat_err));
    if (N >= 5) {
        PetscReal *vr, *vi;
        PetscReal lamr, lami, v2 = 0.0, seam2 = 0.0, maxabs = 0.0;
        PetscCall(PetscMalloc2(map.n, &vr, map.n, &vi));
        PetscCall(DenseMaxRealRightEigenpair(J, map.n, &lamr, &lami, vr, vi));
        for (PetscInt m = 0; m < map.n; m++) {
            const PetscReal a2 = vr[m]*vr[m] + vi[m]*vi[m];
            v2 += a2;
            if (DofTouchesPeriodicRepresentative(&map, m, user->info)) seam2 += a2;
            maxabs = PetscMax(maxabs, PetscSqrtReal(a2));
        }
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "    larger-grid eigenvector: lambda=%.6e%+.6ei seam_energy=%.3f max_entry/||v||=%.3e\n",
            (double)lamr, (double)lami, (double)(seam2/PetscMax(v2, PETSC_MACHINE_EPSILON)),
            (double)(maxabs/PetscSqrtReal(PetscMax(v2, PETSC_MACHINE_EPSILON)))));
        PetscCall(PrintBlockAndSymmetricLocalization(user, &map, J));
        PetscCall(PetscFree2(vr, vi));
    }
    PetscCall(PetscFree3(Rp, Rm, J));
    PetscCall(VecDestroy(&Ubase)); PetscCall(VecDestroy(&Rhs)); PetscCall(VecDestroy(&Uwork));
    PetscCall(DofMapDestroy(&map));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunStateASplitOne()

static PetscErrorCode RunStateASplitOne ( CandState  st, const char *  label  )

static

Runs one State A transport-direction split.

Definition at line 2326 of file test_momentum_convective_candidates.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL; DofMap map;
    Vec Ubase, Rhs, Uwork;
    PetscReal *Rp, *Rm, *J;
    PetscReal repeat_err, maxdiv;
    SeamDiagnostics seam;
    PetscFunctionBeginUser;
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, 4, 4, 4, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    PetscCall(ConfigureCandidateFixture(simCtx, user));
    PetscCall(DofMapBuild(user, &map));
    PetscCall(VecDuplicate(user->Ucont, &Ubase));
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(VecDuplicate(user->Ucont, &Uwork));
    PetscCall(PetscMalloc3(map.n, &Rp, map.n, &Rm, (size_t)map.n*map.n, &J));
    PetscCall(BuildBaseState(user, st, Ubase, &repeat_err, &maxdiv, &seam));
    PetscCall(BuildFDJacobian(user, Ubase, 1e-5, Rhs, &map, Rp, Rm, Uwork, J));
    PetscReal rho, maxre, smax;
    PetscCall(DenseSpectralRadius(J, map.n, &rho, &maxre));
    PetscCall(DenseSigmaMax(J, map.n, &smax, NULL));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  %-4s rho=%.6e sigma=%.6e max_real=%.6e skew=%.3e nonnormality=%.3e div=%.3e repeat=%.3e\n",
        label, (double)rho, (double)smax, (double)maxre,
        (double)DenseSkewnessDefect(J, map.n), (double)DenseNonNormality(J, map.n),
        (double)maxdiv, (double)repeat_err));
    PetscCall(PetscFree3(Rp, Rm, J));
    PetscCall(VecDestroy(&Ubase)); PetscCall(VecDestroy(&Rhs)); PetscCall(VecDestroy(&Uwork));
    PetscCall(DofMapDestroy(&map));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestStateAGridAudit()

static PetscErrorCode TestStateAGridAudit ( void  )

static

Runs the State A grid-dependence and active-space audit.

Definition at line 2361 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "\n================ STATE A GRID / ACTIVE-SPACE AUDIT ================\n"
        "  residual sign convention: ComputeRHS returns the production convection residual used by pseudo-time updates.\n"
        "  component-staggered periodic duplicate planes: 0<-m-2 and m-1<-1 in each direction.\n"
        "  independent map: all three Ucont components use representatives i,j,k=1..m-2; count = 3*(m-2)^3.\n"));
    PetscCall(RunStateAGridAuditOne(4));
    PetscCall(RunStateAGridAuditOne(5));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  State A transport-direction split on DMDA 5^3:\n"));
    PetscCall(RunStateASplitOne(STATE_A_X, "A-x"));
    PetscCall(RunStateASplitOne(STATE_A_Y, "A-y"));
    PetscCall(RunStateASplitOne(STATE_A,   "A-xy"));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteStateADecompReference()

static PetscErrorCode WriteStateADecompReference ( GlobalVecStats  r0, GlobalVecStats  jv, GlobalVecStats  phi, const MomStabilityReport *  rep  )

static

Writes the one-rank State A matrix-free decomposition reference.

Definition at line 2382 of file test_momentum_convective_candidates.c.

{
    PetscMPIInt rank;
    PetscFunctionBeginUser;
    if (!g_ref_path_set) PetscFunctionReturn(0);
    PetscCheck(g_ref_token_set, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "-candidate_ref_token is required when -candidate_ref_path is set");
    PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    if (rank == 0) {
        FILE *fp = fopen(g_ref_path, "w");
        PetscCheck(fp != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "could not write State A MPI reference");
        fprintf(fp, "PICURV_CANDIDATE_STATEA_REF_V2 %s\n", g_ref_token);
        fprintf(fp, "%.17e %.17e %.17e\n", (double)r0.n2, (double)r0.ninf, (double)r0.checksum);
        fprintf(fp, "%.17e %.17e %.17e\n", (double)jv.n2, (double)jv.ninf, (double)jv.checksum);
        fprintf(fp, "%.17e %.17e %.17e\n", (double)phi.n2, (double)phi.ninf, (double)phi.checksum);
        fprintf(fp, "%.17e %.17e %.17e %d %d %d %d\n",
                (double)(rep->lambda_B - rep->lambda_t),
                (double)(rep->lambda_C - rep->lambda_t),
                (double)(rep->lambda_D - rep->lambda_t),
                (int)rep->active_cells, (int)rep->cblock, (int)rep->ci, (int)rep->cj);
        fclose(fp);
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ReadAndCompareStateADecompReference()

static PetscErrorCode ReadAndCompareStateADecompReference ( GlobalVecStats  r0, GlobalVecStats  jv, GlobalVecStats  phi, const MomStabilityReport *  rep  )

static

Compares a distributed State A matrix-free check against the one-rank reference.

Definition at line 2411 of file test_momentum_convective_candidates.c.

{
    PetscMPIInt rank;
    PetscReal vals[16] = {0.0};
    PetscFunctionBeginUser;
    PetscCheck(g_ref_path_set && g_ref_token_set, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "two-rank State A decomp check requires matching -candidate_ref_path and -candidate_ref_token from a preceding one-rank run");
    PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    if (rank == 0) {
        char magic[64], token[128];
        FILE *fp = fopen(g_ref_path, "r");
        PetscCheck(fp != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
                   "State A one-rank decomp reference missing; run the Makefile target so -n 1 precedes -n 2");
        PetscCheck(fscanf(fp, "%63s %127s", magic, token) == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "bad reference header");
        PetscCheck(strcmp(magic, "PICURV_CANDIDATE_STATEA_REF_V2") == 0, PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                   "bad reference magic");
        PetscCheck(strcmp(token, g_ref_token) == 0, PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                   "State A one-rank decomp reference token mismatch");
        int active, cblock, ci, cj;
        PetscCheck(fscanf(fp, "%le %le %le", &vals[0], &vals[1], &vals[2]) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "bad reference R0");
        PetscCheck(fscanf(fp, "%le %le %le", &vals[3], &vals[4], &vals[5]) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "bad reference Jv");
        PetscCheck(fscanf(fp, "%le %le %le", &vals[6], &vals[7], &vals[8]) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "bad reference Phi");
        PetscCheck(fscanf(fp, "%le %le %le %d %d %d %d", &vals[9], &vals[10], &vals[11],
                          &active, &cblock, &ci, &cj) == 7, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "bad reference estimator");
        vals[12] = (PetscReal)active; vals[13] = (PetscReal)cblock; vals[14] = (PetscReal)ci; vals[15] = (PetscReal)cj;
        fclose(fp);
        PetscCheck(remove(g_ref_path) == 0, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
                   "could not remove consumed State A MPI reference '%s'", g_ref_path);
    }
    PetscCallMPI(MPI_Bcast(vals, 16, MPIU_REAL, 0, PETSC_COMM_WORLD));
    PetscCall(PicurvAssertRealNear(r0.n2, vals[0], 1e-10, "MPI State A R0 L2"));
    PetscCall(PicurvAssertRealNear(r0.ninf, vals[1], 1e-10, "MPI State A R0 Linf"));
    PetscCall(PicurvAssertRealNear(r0.checksum, vals[2], 1e-10, "MPI State A R0 checksum"));
    PetscCall(PicurvAssertRealNear(jv.n2, vals[3], 1e-8, "MPI State A Jv L2"));
    PetscCall(PicurvAssertRealNear(jv.ninf, vals[4], 1e-8, "MPI State A Jv Linf"));
    PetscCall(PicurvAssertRealNear(jv.checksum, vals[5], 1e-8, "MPI State A Jv checksum"));
    PetscCall(PicurvAssertRealNear(phi.n2, vals[6], 1e-8, "MPI State A DPhi v L2"));
    PetscCall(PicurvAssertRealNear(phi.ninf, vals[7], 1e-8, "MPI State A DPhi v Linf"));
    PetscCall(PicurvAssertRealNear(phi.checksum, vals[8], 1e-8, "MPI State A DPhi v checksum"));
    PetscCall(PicurvAssertRealNear(rep->lambda_B - rep->lambda_t, vals[9], 1e-10, "MPI State A lambda_cB"));
    PetscCall(PicurvAssertRealNear(rep->lambda_C - rep->lambda_t, vals[10], 1e-10, "MPI State A lambda_cC"));
    PetscCall(PicurvAssertRealNear(rep->lambda_D - rep->lambda_t, vals[11], 1e-10, "MPI State A lambda_cD"));
    PetscCall(PicurvAssertBool((PetscBool)(rep->active_cells == (PetscInt)vals[12]), "MPI State A active cell count"));
    PetscCall(PicurvAssertBool((PetscBool)(rep->cblock == (PetscInt)vals[13]), "MPI State A controlling block"));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StateADecompDirectionalCheck()

static PetscErrorCode StateADecompDirectionalCheck ( UserCtx *  user, Vec  Ubase, Vec  Rhs, PetscReal  lcC, const MomStabilityReport *  rep  )

static

Runs State A matrix-free residual, Jv, and four-stage MPI decomposition checks.

Definition at line 2462 of file test_momentum_convective_candidates.c.

{
    DofMap map;
    Vec V, Up, Um, PhiP, PhiM, Ustage;
    PetscReal *v, *R0, *Rp, *Rm, *ActP, *ActM, *Out;
    GlobalVecStats sR0, sJv, sPhi;
    const PetscReal eps = 1e-6, dtau = 0.5/lcC;
    PetscMPIInt size;
    PetscFunctionBeginUser;
    PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    PetscCall(DofMapBuildOwned(user, &map));
    PetscCall(VecDuplicate(Ubase, &V)); PetscCall(VecDuplicate(Ubase, &Up));
    PetscCall(VecDuplicate(Ubase, &Um)); PetscCall(VecDuplicate(Ubase, &PhiP));
    PetscCall(VecDuplicate(Ubase, &PhiM)); PetscCall(VecDuplicate(Ubase, &Ustage));
    PetscCall(PetscMalloc6(map.n, &v, map.n, &R0, map.n, &Rp, map.n, &Rm, map.n, &ActP, map.n, &ActM));
    PetscCall(PetscMalloc1(map.n, &Out));
    PetscCall(FillDeterministicDirection(user, V, &map, v));
 
    PetscCall(EvalConvResidual(user, Ubase, Rhs, &map, R0));
    PetscCall(ActiveStats(&map, R0, &sR0));
    PetscCall(VecCopy(Ubase, Up)); PetscCall(VecAXPY(Up, eps, V));
    PetscCall(VecCopy(Ubase, Um)); PetscCall(VecAXPY(Um, -eps, V));
    PetscCall(EvalConvResidual(user, Up, Rhs, &map, Rp));
    PetscCall(EvalConvResidual(user, Um, Rhs, &map, Rm));
    for (PetscInt m = 0; m < map.n; m++) Out[m] = (Rp[m]-Rm[m])/(2.0*eps);
    PetscCall(ActiveStats(&map, Out, &sJv));
 
    PetscCall(FourStage(user, Up, dtau, Rhs, &map, Rp, Ustage, PhiP));
    PetscCall(FourStage(user, Um, dtau, Rhs, &map, Rm, Ustage, PhiM));
    PetscCall(ExtractActiveVector(user, PhiP, &map, ActP));
    PetscCall(ExtractActiveVector(user, PhiM, &map, ActM));
    for (PetscInt m = 0; m < map.n; m++) Out[m] = (ActP[m]-ActM[m])/(2.0*eps);
    PetscCall(ActiveStats(&map, Out, &sPhi));
 
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "  State A matrix-free MPI check: R0 L2=%.6e Linf=%.6e checksum=%.6e\n"
        "                                Jv L2=%.6e Linf=%.6e checksum=%.6e\n"
        "                           DPhi v L2=%.6e Linf=%.6e checksum=%.6e\n",
        (double)sR0.n2, (double)sR0.ninf, (double)sR0.checksum,
        (double)sJv.n2, (double)sJv.ninf, (double)sJv.checksum,
        (double)sPhi.n2, (double)sPhi.ninf, (double)sPhi.checksum));
    if (size == 1) PetscCall(WriteStateADecompReference(sR0, sJv, sPhi, rep));
    else           PetscCall(ReadAndCompareStateADecompReference(sR0, sJv, sPhi, rep));
 
    PetscCall(PetscFree6(v, R0, Rp, Rm, ActP, ActM)); PetscCall(PetscFree(Out));
    PetscCall(VecDestroy(&V)); PetscCall(VecDestroy(&Up)); PetscCall(VecDestroy(&Um));
    PetscCall(VecDestroy(&PhiP)); PetscCall(VecDestroy(&PhiM)); PetscCall(VecDestroy(&Ustage));
    PetscCall(DofMapDestroy(&map));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunDecompBaseline()

static PetscErrorCode RunDecompBaseline ( CandState  st, const char *  name, PetscReal  expB, PetscReal  expC, PetscReal  expD  )

static

Runs one decomp baseline and compares scalar estimates with regenerated references.

Definition at line 2517 of file test_momentum_convective_candidates.c.

{
    SimCtx *simCtx = NULL; UserCtx *user = NULL;
    Vec Ubase, Rhs;
    PetscReal repeat_err, maxdiv;
    SeamDiagnostics seam;
    MomStabilityReport rep;
    const PetscInt N = 8;
    PetscFunctionBeginUser;
 
    PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx, &user, N, N, N, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE));
    PetscCall(ConfigureCandidateFixture(simCtx, user));
    PetscCall(VecDuplicate(user->Ucont, &Ubase));
    PetscCall(VecDuplicate(user->Ucont, &Rhs));
    PetscCall(BuildBaseState(user, st, Ubase, &repeat_err, &maxdiv, &seam));
    PetscCall(ComputeMomentumStabilityEstimate(user, 1, simCtx->dt, MOM_STAB_CAND_C, &rep));
 
    const PetscReal lcB = rep.lambda_B - rep.lambda_t;
    const PetscReal lcC = rep.lambda_C - rep.lambda_t;
    const PetscReal lcD = rep.lambda_D - rep.lambda_t;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "\n================ DECOMP BASELINE %s ================\n"
        "  active cells: %d | repeat=%.3e | max|div_h Ucont|=%.3e\n"
        "  lambda_cB=%.6e  lambda_cC=%.6e  lambda_cD=%.6e\n",
        name, (int)rep.active_cells, (double)repeat_err, (double)maxdiv,
        (double)lcB, (double)lcC, (double)lcD));
 
    PetscCall(PicurvAssertBool((PetscBool)(rep.active_cells == 343), "decomp: active-cell count invariant"));
    PetscCall(PicurvAssertRealNear(expB, lcB, 5e-7, "decomp: lambda_cB matches one-rank baseline"));
    PetscCall(PicurvAssertRealNear(expC, lcC, 5e-7, "decomp: lambda_cC matches one-rank baseline"));
    PetscCall(PicurvAssertRealNear(expD, lcD, 5e-7, "decomp: lambda_cD matches one-rank baseline"));
    if (st == STATE_A) {
        PetscCall(StateADecompDirectionalCheck(user, Ubase, Rhs, lcC, &rep));
        PetscCall(PicurvAssertRealNear(maxdiv, 0.0, 1e-12, "decomp A: divergence-free"));
        PetscCall(PicurvAssertRealNear(lcB, lcC, 1e-12, "decomp A: B == C"));
        PetscCall(PicurvAssertRealNear(lcC, lcD, 1e-12, "decomp A: C == D"));
    } else if (st == STATE_B) {
        PetscCall(PicurvAssertBool((PetscBool)(maxdiv > 1e-3), "decomp B: nonzero divergence"));
        PetscCall(PicurvAssertBool((PetscBool)(lcC > lcB), "decomp B: C > B"));
    } else {
        PetscCall(PicurvAssertRealNear(maxdiv, 0.0, 1e-12, "decomp C: divergence-free"));
        PetscCall(PicurvAssertRealNear(lcB, lcC, 1e-12, "decomp C: B == C"));
    }
 
    PetscCall(VecDestroy(&Ubase));
    PetscCall(VecDestroy(&Rhs));
    PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TestDecompA()

static PetscErrorCode TestDecompA ( void  )

static

Runs the State A decomposition baseline.

Definition at line 2571 of file test_momentum_convective_candidates.c.

{ return RunDecompBaseline(STATE_A, "A (uniform div-free)", 1.1000000, 1.1000000, 1.1000000); }

Here is the call graph for this function:

Here is the caller graph for this function:

TestDecompB()

static PetscErrorCode TestDecompB ( void  )

static

Runs the State B decomposition baseline.

Definition at line 2576 of file test_momentum_convective_candidates.c.

{ return RunDecompBaseline(STATE_B, "B (nonzero divergence)", 0.878379697325, 0.974927912182, 1.572173303885); }

Here is the call graph for this function:

Here is the caller graph for this function:

TestDecompC()

static PetscErrorCode TestDecompC ( void  )

static

Runs the State C decomposition baseline.

Definition at line 2581 of file test_momentum_convective_candidates.c.

{ return RunDecompBaseline(STATE_C, "C (div-free shear)", 1.5000000, 1.5000000, 1.780330085890); }

Here is the call graph for this function:

Here is the caller graph for this function:

ConfigureMPIReferenceOptions()

static PetscErrorCode ConfigureMPIReferenceOptions ( void  )

static

Reads optional paired-run MPI reference path and token.

Definition at line 2586 of file test_momentum_convective_candidates.c.

{
    PetscFunctionBeginUser;
    PetscCall(PetscOptionsGetString(NULL, NULL, "-candidate_ref_path",
                                    g_ref_path, sizeof(g_ref_path), &g_ref_path_set));
    PetscCall(PetscOptionsGetString(NULL, NULL, "-candidate_ref_token",
                                    g_ref_token, sizeof(g_ref_token), &g_ref_token_set));
    if (g_ref_path_set || g_ref_token_set) {
        PetscCheck(g_ref_path_set && g_ref_token_set, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                   "-candidate_ref_path and -candidate_ref_token must be supplied together");
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

PETSc entry point for the focused convective-candidate harness.

Definition at line 2603 of file test_momentum_convective_candidates.c.

{
    PetscErrorCode ierr;
    PetscMPIInt size;
    const PicurvTestCase cases[] = {
        {"candidate-state-A-uniform", TestStateA},
        {"candidate-state-B-divergence", TestStateB},
        {"candidate-state-C-shear", TestStateC},
        {"candidate-state-A-grid-audit", TestStateAGridAudit},
    };
    const PicurvTestCase decomp_cases[] = {
        {"candidate-decomp-A-uniform", TestDecompA},
        {"candidate-decomp-B-divergence", TestDecompB},
        {"candidate-decomp-C-shear", TestDecompC},
    };
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv A4a convective-candidate study");
    if (ierr) return (int)ierr;
    ierr = ConfigureMPIReferenceOptions(); if (ierr) { PetscFinalize(); return (int)ierr; }
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); if (ierr) { PetscFinalize(); return (int)ierr; }
    if (size == 1) {
        ierr = PicurvRunTests("unit-momentum-candidates", cases, sizeof(cases)/sizeof(cases[0]));
        if (!ierr) ierr = PicurvRunTests("unit-momentum-candidates-decomp", decomp_cases, sizeof(decomp_cases)/sizeof(decomp_cases[0]));
    } else {
        ierr = PicurvRunTests("unit-momentum-candidates-decomp", decomp_cases, sizeof(decomp_cases)/sizeof(decomp_cases[0]));
    }
    if (ierr) { PetscFinalize(); return (int)ierr; }
    return (int)PetscFinalize();
}

Here is the call graph for this function:

Variable Documentation

g_ref_path

char g_ref_path[PETSC_MAX_PATH_LEN] = ""

static

Definition at line 36 of file test_momentum_convective_candidates.c.

g_ref_token

char g_ref_token[128] = ""

static

Definition at line 37 of file test_momentum_convective_candidates.c.

g_ref_path_set

PetscBool g_ref_path_set = PETSC_FALSE

static

Definition at line 38 of file test_momentum_convective_candidates.c.

g_ref_token_set

PetscBool g_ref_token_set = PETSC_FALSE

static

Definition at line 39 of file test_momentum_convective_candidates.c.