test_momentum_newton_boundary_fixedpoint.c File Reference

Opt-in integration regression for the Newton–Krylov deterministic Cartesian-seed correction on the production-sized straight duct. More...

#include "test_support.h"
#include "initialcondition.h"
#include "runloop.h"
#include "solvers.h"
#include "../../src/momentum_newton_krylov.c"

Include dependency graph for test_momentum_newton_boundary_fixedpoint.c:

Go to the source code of this file.

Data Structures
struct	FpSnapshot
struct	SolveMonitorCtx
struct	SolveResult

Macros
#define	MomentumSolver_NewtonKrylov   MomentumSolver_NewtonKrylov_BoundaryFixedpointPrivateCopy
#define	FP_FIELD_COUNT   16
#define	FP_SCALAR_COUNT   10

Functions
static void	FpGetVectors (UserCtx *user, Vec field[16])
	Collects the persistent vectors observed by the residual/boundary path.
static void	FpGetScalars (UserCtx *user, PetscReal scalar[10])
	Collects the persistent scalar diagnostics maintained by boundary processing.
static PetscErrorCode	FpCapture (UserCtx *user, FpSnapshot *snap)
	Deep-copies the audited persistent state into a test-owned snapshot.
static PetscErrorCode	FpRestore (UserCtx *user, const FpSnapshot *snap)
	Restores a previously captured persistent-state snapshot.
static PetscErrorCode	FpDestroy (FpSnapshot *snap)
	Releases all vectors owned by one persistent-state snapshot.
static PetscErrorCode	CollectiveDifference (Vec a, Vec b, PetscReal *norm_inf, PetscReal *norm_2)
	Computes collective global difference norms for global or rank-local vectors.
static PetscErrorCode	DifferenceMaxRank (Vec a, Vec b, PetscMPIInt *max_rank)
	Finds the MPI rank containing the largest local entry difference.
static PetscErrorCode	SeedCartesianFromX (UserCtx *user, Vec X)
	Reconstructs Ucat/lUcat deterministically from X, mirroring the seed now performed by the production callback.
static PetscErrorCode	ResidualNCalls (MomentumNewtonKrylovContext *ctx, Vec X, PetscInt boundary_calls, Vec F)
	Deterministic residual at X using an explicit number of boundary-condition calls (each three internal passes), for the boundary-map compatibility test.
static PetscErrorCode	ConfigureNewtonOptions (void)
	Installs the production Newton, line-search, and GMRES options.
static PetscErrorCode	SolveMonitor (SNES snes, PetscInt it, PetscReal norm, void *vctx)
	Records per-Newton nonlinear norm, KSP iterations/reason, and accepted lambda.
static PetscErrorCode	CreateProductionSNES (UserCtx *user, MomentumNewtonKrylovContext *pctx, PetscBool use_mgs, SNES *snes, Mat *J, KSP *ksp)
	Creates a test-owned matrix-free SNES driven by the production residual callback.
static PetscErrorCode	RunProductionFullSolve (UserCtx *user, const FpSnapshot *entry, Vec X0, PetscBool use_mgs, const char *label, SolveResult *res)
	Runs one complete production-callback step-2 solve from the restored entry state.
static PetscErrorCode	StopAfterOne (SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal fnorm, SNESConvergedReason *reason, void *vctx)
	Stops SNES immediately after its first accepted Newton correction.
static PetscErrorCode	TestNewtonKrylovProductionRegression (void)
	Builds the production-sized duct, advances the true physical step 1, and exercises purity, full-solve convergence, boundary-map compatibility, and projection using the production Newton–Krylov callback and solver.
int	main (int argc, char **argv)
	Runs the opt-in Newton–Krylov production regression executable.

Detailed Description

Opt-in integration regression for the Newton–Krylov deterministic Cartesian-seed correction on the production-sized straight duct.

The diagnostic phase established that the production Newton–Krylov residual must reconstruct Ucat/lUcat from the current SNES trial vector X before the first boundary sweep (see src/momentum_newton_krylov.c). That correction is now in production. This file locks the resulting guarantees in place using the actual production residual callback (included via a privately renamed copy of the production translation unit) – it deliberately does NOT keep a duplicate test-local seeded callback, so it exercises exactly what ships:

1. residual purity at the reproduced physical step-2 state (immediate and after real PETSc MFFD products);
2. a complete production step-2 solve from the true projected step-1 state using PETSc's default classical Gram–Schmidt;
3. boundary-map compatibility: the converged three-pass solution also zeros the 24-pass outlet residual;
4. clean pressure projection of the converged momentum state.

The heavy fixture keeps this target out of make check; the cheap seed-removal detector lives in unit-newton-krylov (default suite). Central checks run on one and four MPI ranks.

Definition in file test_momentum_newton_boundary_fixedpoint.c.

---

Data Structure Documentation

FpSnapshot

struct FpSnapshot

Definition at line 43 of file test_momentum_newton_boundary_fixedpoint.c.

Data Fields
Vec	value[16]
PetscReal	scalar[10]

SolveMonitorCtx

struct SolveMonitorCtx

Definition at line 234 of file test_momentum_newton_boundary_fixedpoint.c.

Data Fields
const char *	label
PetscReal *	initial_norm

SolveResult

struct SolveResult

Definition at line 260 of file test_momentum_newton_boundary_fixedpoint.c.

Data Fields
SNESConvergedReason	reason
PetscInt	its
PetscInt	fevals
PetscInt	ksp_total
PetscReal	initial_norm
PetscReal	final_norm
PetscBool	converged
Vec	Xfinal

Macro Definition Documentation

MomentumSolver_NewtonKrylov

#define MomentumSolver_NewtonKrylov   MomentumSolver_NewtonKrylov_BoundaryFixedpointPrivateCopy

Definition at line 32 of file test_momentum_newton_boundary_fixedpoint.c.

FP_FIELD_COUNT

#define FP_FIELD_COUNT   16

Definition at line 40 of file test_momentum_newton_boundary_fixedpoint.c.

FP_SCALAR_COUNT

#define FP_SCALAR_COUNT   10

Definition at line 41 of file test_momentum_newton_boundary_fixedpoint.c.

Function Documentation

FpGetVectors()

static void FpGetVectors ( UserCtx *  user, Vec  field[16]  )

static

Collects the persistent vectors observed by the residual/boundary path.

Definition at line 49 of file test_momentum_newton_boundary_fixedpoint.c.

{
    field[0] = user->Ucont; field[1] = user->lUcont;
    field[2] = user->Ucat; field[3] = user->lUcat;
    field[4] = user->P; field[5] = user->lP;
    field[6] = user->Bcs.Ubcs; field[7] = user->Bcs.Uch;
    field[8] = user->Rhs;
    field[9] = user->Ucont_o; field[10] = user->lUcont_o;
    field[11] = user->Ucont_rm1; field[12] = user->lUcont_rm1;
    field[13] = user->Ucat_o; field[14] = user->P_o;
    field[15] = user->CellFieldAtCorner;
}

Here is the caller graph for this function:

FpGetScalars()

static void FpGetScalars ( UserCtx *  user, PetscReal  scalar[10]  )

static

Collects the persistent scalar diagnostics maintained by boundary processing.

Definition at line 63 of file test_momentum_newton_boundary_fixedpoint.c.

{
    SimCtx *s = user->simCtx;
    scalar[0] = s->FluxInSum; scalar[1] = s->FluxOutSum;
    scalar[2] = s->FarFluxInSum; scalar[3] = s->FarFluxOutSum;
    scalar[4] = s->Fluxsum; scalar[5] = s->AreaInSum;
    scalar[6] = s->AreaOutSum; scalar[7] = s->bulkVelocityCorrection;
    scalar[8] = s->boundaryVelocityCorrection; scalar[9] = user->FluxIntpSum;
}

Here is the caller graph for this function:

FpCapture()

static PetscErrorCode FpCapture ( UserCtx *  user, FpSnapshot *  snap  )

static

Deep-copies the audited persistent state into a test-owned snapshot.

Definition at line 74 of file test_momentum_newton_boundary_fixedpoint.c.

{
    Vec field[FP_FIELD_COUNT];
    PetscFunctionBeginUser;
    memset(snap, 0, sizeof(*snap));
    FpGetVectors(user, field);
    for (PetscInt i = 0; i < FP_FIELD_COUNT; ++i) {
        if (field[i]) {
            PetscCall(VecDuplicate(field[i], &snap->value[i]));
            PetscCall(VecCopy(field[i], snap->value[i]));
        }
    }
    FpGetScalars(user, snap->scalar);
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FpRestore()

static PetscErrorCode FpRestore ( UserCtx *  user, const FpSnapshot *  snap  )

static

Restores a previously captured persistent-state snapshot.

Definition at line 91 of file test_momentum_newton_boundary_fixedpoint.c.

{
    Vec field[FP_FIELD_COUNT];
    SimCtx *s = user->simCtx;
    PetscFunctionBeginUser;
    FpGetVectors(user, field);
    for (PetscInt i = 0; i < FP_FIELD_COUNT; ++i)
        if (field[i] && snap->value[i]) PetscCall(VecCopy(snap->value[i], field[i]));
    s->FluxInSum = snap->scalar[0]; s->FluxOutSum = snap->scalar[1];
    s->FarFluxInSum = snap->scalar[2]; s->FarFluxOutSum = snap->scalar[3];
    s->Fluxsum = snap->scalar[4]; s->AreaInSum = snap->scalar[5];
    s->AreaOutSum = snap->scalar[6]; s->bulkVelocityCorrection = snap->scalar[7];
    s->boundaryVelocityCorrection = snap->scalar[8]; user->FluxIntpSum = snap->scalar[9];
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FpDestroy()

static PetscErrorCode FpDestroy ( FpSnapshot *  snap )

static

Releases all vectors owned by one persistent-state snapshot.

Definition at line 108 of file test_momentum_newton_boundary_fixedpoint.c.

{
    PetscFunctionBeginUser;
    for (PetscInt i = 0; i < FP_FIELD_COUNT; ++i) PetscCall(VecDestroy(&snap->value[i]));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

CollectiveDifference()

static PetscErrorCode CollectiveDifference ( Vec  a, Vec  b, PetscReal *  norm_inf, PetscReal *  norm_2  )

static

Computes collective global difference norms for global or rank-local vectors.

Definition at line 120 of file test_momentum_newton_boundary_fixedpoint.c.

{
    Vec delta = NULL;
    MPI_Comm comm;
    PetscMPIInt comm_size, world_size;
    PetscReal local_inf, local_2;
    PetscFunctionBeginUser;
    PetscCall(VecDuplicate(a, &delta));
    PetscCall(VecWAXPY(delta, -1.0, a, b));
    PetscCall(VecNorm(delta, NORM_INFINITY, &local_inf));
    PetscCall(VecNorm(delta, NORM_2, &local_2));
    comm = PetscObjectComm((PetscObject)a);
    PetscCallMPI(MPI_Comm_size(comm, &comm_size));
    PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
    if (comm_size == 1 && world_size > 1) {
        PetscReal square = local_2 * local_2;
        PetscCallMPI(MPI_Allreduce(&local_inf, norm_inf, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
        PetscCallMPI(MPI_Allreduce(&square, norm_2, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
        *norm_2 = PetscSqrtReal(*norm_2);
    } else {
        *norm_inf = local_inf;
        *norm_2 = local_2;
    }
    PetscCall(VecDestroy(&delta));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

DifferenceMaxRank()

static PetscErrorCode DifferenceMaxRank ( Vec  a, Vec  b, PetscMPIInt *  max_rank  )

static

Finds the MPI rank containing the largest local entry difference.

Definition at line 148 of file test_momentum_newton_boundary_fixedpoint.c.

{
    const PetscScalar *aa = NULL, *bb = NULL;
    PetscInt nlocal;
    PetscReal local_max = 0.0, *all_max = NULL;
    PetscMPIInt size;
    PetscFunctionBeginUser;
    PetscCall(VecGetLocalSize(a, &nlocal));
    PetscCall(VecGetArrayRead(a, &aa)); PetscCall(VecGetArrayRead(b, &bb));
    for (PetscInt i = 0; i < nlocal; ++i) local_max = PetscMax(local_max, PetscAbsScalar(bb[i] - aa[i]));
    PetscCall(VecRestoreArrayRead(b, &bb)); PetscCall(VecRestoreArrayRead(a, &aa));
    PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    PetscCall(PetscMalloc1(size, &all_max));
    PetscCallMPI(MPI_Allgather(&local_max, 1, MPIU_REAL, all_max, 1, MPIU_REAL, PETSC_COMM_WORLD));
    *max_rank = 0;
    for (PetscMPIInt r = 1; r < size; ++r) if (all_max[r] > all_max[*max_rank]) *max_rank = r;
    PetscCall(PetscFree(all_max));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

SeedCartesianFromX()

static PetscErrorCode SeedCartesianFromX ( UserCtx *  user, Vec  X  )

static

Reconstructs Ucat/lUcat deterministically from X, mirroring the seed now performed by the production callback.

Used only by the N-pass compatibility residual (Experiment 3), which needs a boundary-pass count other than the production default of one call (three passes). For a single call it reproduces the production callback exactly.

Definition at line 180 of file test_momentum_newton_boundary_fixedpoint.c.

{
    const char *staggered[] = {"Ucont"};
    const char *cell[] = {"Ucat"};
    PetscFunctionBeginUser;
    PetscCall(VecCopy(X, user->Ucont));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, staggered));
    PetscCall(Contra2Cart(user));
    PetscCall(SynchronizePeriodicCellFields(user, 1, cell));
    PetscCall(UpdateLocalGhosts(user, "Ucat"));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ResidualNCalls()

static PetscErrorCode ResidualNCalls ( MomentumNewtonKrylovContext *  ctx, Vec  X, PetscInt  boundary_calls, Vec  F  )

static

Deterministic residual at X using an explicit number of boundary-condition calls (each three internal passes), for the boundary-map compatibility test.

Definition at line 197 of file test_momentum_newton_boundary_fixedpoint.c.

{
    UserCtx *user = ctx->user;
    PetscFunctionBeginUser;
    PetscCall(SeedCartesianFromX(user, X));
    for (PetscInt i = 0; i < boundary_calls; ++i) PetscCall(ApplyBoundaryConditions(user));
    PetscCall(ComputeTotalResidual(user));
    PetscCall(VecCopy(user->Rhs, F));
    PetscCall(VecScale(F, -1.0));
    PetscCall(MomentumNewtonKrylov_ApplyConstraints(ctx, X, F));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ConfigureNewtonOptions()

static PetscErrorCode ConfigureNewtonOptions ( void  )

static

Installs the production Newton, line-search, and GMRES options.

Definition at line 216 of file test_momentum_newton_boundary_fixedpoint.c.

{
    PetscFunctionBeginUser;
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_type", "newtonls"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_atol", "1e-10"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_rtol", "1e-8"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_stol", "1e-12"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_max_it", "25"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_snes_linesearch_type", "bt"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_type", "gmres"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_atol", "1e-10"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_rtol", "1e-6"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_max_it", "400"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_ksp_gmres_restart", "80"));
    PetscCall(PetscOptionsSetValue(NULL, "-mom_nk_pc_type", "none"));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

SolveMonitor()

static PetscErrorCode SolveMonitor ( SNES  snes, PetscInt  it, PetscReal  norm, void *  vctx  )

static

Records per-Newton nonlinear norm, KSP iterations/reason, and accepted lambda.

Definition at line 237 of file test_momentum_newton_boundary_fixedpoint.c.

{
    SolveMonitorCtx *m = (SolveMonitorCtx *)vctx;
    SNESLineSearch ls = NULL;
    KSP ksp = NULL;
    KSPConvergedReason kr = KSP_CONVERGED_ITERATING;
    PetscReal lambda = -1.0;
    PetscInt cum_ksp = 0, step_ksp = 0;
    PetscFunctionBeginUser;
    if (it == 0 && m->initial_norm) *m->initial_norm = norm;
    PetscCall(SNESGetLinearSolveIterations(snes, &cum_ksp));
    PetscCall(SNESGetKSP(snes, &ksp));
    PetscCall(KSPGetConvergedReason(ksp, &kr));
    PetscCall(KSPGetIterationNumber(ksp, &step_ksp));
    PetscCall(SNESGetLineSearch(snes, &ls));
    if (it > 0) PetscCall(SNESLineSearchGetLambda(ls, &lambda));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "FULLSOLVE %s newton=%d nonlinear_norm=% .16e last_ksp_its=%d cumulative_ksp_its=%d "
        "last_ksp_reason=%d accepted_lambda=% .6e\n",
        m->label, (int)it, (double)norm, (int)step_ksp, (int)cum_ksp, (int)kr, (double)lambda));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

CreateProductionSNES()

static PetscErrorCode CreateProductionSNES ( UserCtx *  user, MomentumNewtonKrylovContext *  pctx, PetscBool  use_mgs, SNES *  snes, Mat *  J, KSP *  ksp  )

static

Creates a test-owned matrix-free SNES driven by the production residual callback.

Parameters

use_mgs

When false (production default), leaves GMRES orthogonalization at PETSc's default classical Gram–Schmidt – the path the production solver uses.

Definition at line 273 of file test_momentum_newton_boundary_fixedpoint.c.

{
    PC pc = NULL;
    PetscFunctionBeginUser;
    PetscCall(SNESCreate(PETSC_COMM_WORLD, snes));
    PetscCall(SNESSetOptionsPrefix(*snes, "mom_nk_"));
    PetscCall(SNESSetType(*snes, SNESNEWTONLS));
    PetscCall(SNESSetDM(*snes, user->fda));
    PetscCall(SNESSetFunction(*snes, NULL, MomentumNewtonKrylov_FormResidual, pctx));
    PetscCall(MatCreateSNESMF(*snes, J));
    PetscCall(SNESSetJacobian(*snes, *J, *J, MatMFFDComputeJacobian, NULL));
    PetscCall(SNESGetKSP(*snes, ksp));
    PetscCall(KSPSetType(*ksp, KSPGMRES));
    if (use_mgs) PetscCall(KSPGMRESSetOrthogonalization(*ksp, KSPGMRESModifiedGramSchmidtOrthogonalization));
    PetscCall(KSPGMRESSetRestart(*ksp, 80));
    PetscCall(KSPGetPC(*ksp, &pc));
    PetscCall(PCSetType(pc, PCNONE));
    PetscCall(SNESSetFromOptions(*snes));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

RunProductionFullSolve()

static PetscErrorCode RunProductionFullSolve ( UserCtx *  user, const FpSnapshot *  entry, Vec  X0, PetscBool  use_mgs, const char *  label, SolveResult *  res  )

static

Runs one complete production-callback step-2 solve from the restored entry state.

Definition at line 296 of file test_momentum_newton_boundary_fixedpoint.c.

{
    MomentumNewtonKrylovContext pctx = {user, NULL, PETSC_FALSE, 0.0};
    SNES snes = NULL; Mat J = NULL; KSP ksp = NULL; Vec solution = NULL;
    SolveMonitorCtx m;
    PetscFunctionBeginUser;
    memset(res, 0, sizeof(*res));
    res->initial_norm = -1.0;
    PetscCall(FpRestore(user, entry));
    PetscCall(VecDuplicate(X0, &solution)); PetscCall(VecCopy(X0, solution));
    PetscCall(CreateProductionSNES(user, &pctx, use_mgs, &snes, &J, &ksp));
    m.label = label; m.initial_norm = &res->initial_norm;
    PetscCall(SNESMonitorSet(snes, SolveMonitor, &m, NULL));
    PetscCall(SNESSolve(snes, NULL, solution));
    PetscCall(SNESGetConvergedReason(snes, &res->reason));
    PetscCall(SNESGetIterationNumber(snes, &res->its));
    PetscCall(SNESGetNumberFunctionEvals(snes, &res->fevals));
    PetscCall(SNESGetLinearSolveIterations(snes, &res->ksp_total));
    PetscCall(SNESGetFunctionNorm(snes, &res->final_norm));
    res->converged = (PetscBool)(res->reason > 0);
    PetscCall(VecDuplicate(solution, &res->Xfinal)); PetscCall(VecCopy(solution, res->Xfinal));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "FULLSOLVE_SUMMARY %s reason=%s (%d) newton=%d fevals=%d ksp_total=%d orthogonalization=%s "
        "initial_norm=% .16e final_norm=% .16e -> %s\n",
        label, SNESConvergedReasons[res->reason], (int)res->reason, (int)res->its, (int)res->fevals,
        (int)res->ksp_total, use_mgs ? "modified-GS" : "classical-GS(default)",
        (double)res->initial_norm, (double)res->final_norm, res->converged ? "CONVERGED" : "FAILED"));
    PetscCall(VecDestroy(&solution)); PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StopAfterOne()

static PetscErrorCode StopAfterOne ( SNES  snes, PetscInt  it, PetscReal  xnorm, PetscReal  gnorm, PetscReal  fnorm, SNESConvergedReason *  reason, void *  vctx  )

static

Stops SNES immediately after its first accepted Newton correction.

Definition at line 329 of file test_momentum_newton_boundary_fixedpoint.c.

{
    PetscReal *norms = (PetscReal *)vctx;
    PetscFunctionBeginUser;
    (void)snes; (void)xnorm; (void)gnorm;
    if (it == 0) norms[0] = fnorm;
    if (it == 1) norms[1] = fnorm;
    *reason = it >= 1 ? SNES_CONVERGED_ITS : SNES_CONVERGED_ITERATING;
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

TestNewtonKrylovProductionRegression()

static PetscErrorCode TestNewtonKrylovProductionRegression ( void  )

static

Builds the production-sized duct, advances the true physical step 1, and exercises purity, full-solve convergence, boundary-map compatibility, and projection using the production Newton–Krylov callback and solver.

Definition at line 350 of file test_momentum_newton_boundary_fixedpoint.c.

{
    SimCtx *simCtx = NULL;
    UserCtx *user = NULL;
    char tmpdir[PETSC_MAX_PATH_LEN] = "";
    FpSnapshot entry = {0};
    MomentumNewtonKrylovContext pctx;
    Vec X0 = NULL, X1 = NULL, Xfinal = NULL;
    Vec F0 = NULL, F1 = NULL, Fafter = NULL, v = NULL, jv = NULL;
    Vec F3 = NULL, F24 = NULL;
    SolveResult res_cgs = {0};
    PetscReal norms[2] = {0.0, 0.0};
    PetscReal imm_inf = 0.0, imm_2 = 0.0, mffd_inf = 0.0, mffd_2 = 0.0;
    PetscReal f3_2 = 0.0, f24_2 = 0.0, f24_inf = 0.0, f24m3_inf = 0.0, f24m3_2 = 0.0;
    PetscReal max_div = -1.0;
    const char *fields[] = {"Ucont"};
    const char *stag[] = {"Ucont"}; const char *cell[] = {"Ucat"};
    const PetscReal purity_inf_tol = 1e-13, purity_2_tol = 1e-12;
    const PetscReal compat_tol = 1e-8, div_tol = 1e-6;
    PetscMPIInt world, mffd_rank = 0;
    PetscErrorCode prod_ierr;
 
    PetscFunctionBeginUser;
    PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "\n########## NEWTON-KRYLOV PRODUCTION REGRESSION (ranks=%d) ##########\n", (int)world));
    PetscCall(PicurvBuildMomentumPurityRuntimeContext(NULL, &simCtx, &user, tmpdir, sizeof(tmpdir)));
    PetscCall(InitializeEulerianState(simCtx));
    PetscCall(ConfigureNewtonOptions());
 
    /* Complete physical step 1 (momentum + Poisson + projection + finalization)
     * through the real run path, advance histories, then enter step 2. */
    simCtx->step = 1; simCtx->ti = simCtx->dt;
    PetscCall(FlowSolver(simCtx));
    PetscCall(UpdateSolverHistoryVectors(user));
    simCtx->step = 2; simCtx->ti = 2.0 * simCtx->dt;
 
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    pctx.user = user; pctx.history_file = NULL; pctx.have_initial_norm = PETSC_FALSE; pctx.initial_norm = 0.0;
 
    /* Canonicalize the entry exactly as production MomentumSolver_NewtonKrylov does. */
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, fields));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(FpCapture(user, &entry));
    PetscCall(VecDuplicate(user->Ucont, &X0)); PetscCall(VecCopy(user->Ucont, X0));
 
    /* -------- 4.1 residual purity at the step-2 / Newton-1 state -------- */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n===== 4.1 production-callback residual purity =====\n"));
    {   /* Reach the accepted Newton-1 state X1 with the production callback + default CGS. */
        SNES snes = NULL; Mat J = NULL; KSP ksp = NULL;
        PetscCall(VecDuplicate(X0, &X1)); PetscCall(VecCopy(X0, X1));
        PetscCall(CreateProductionSNES(user, &pctx, PETSC_FALSE, &snes, &J, &ksp));
        PetscCall(SNESSetConvergenceTest(snes, StopAfterOne, norms, NULL));
        PetscCall(SNESSolve(snes, NULL, X1));
        PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
    }
    PetscCall(VecDuplicate(X1, &F0)); PetscCall(VecDuplicate(X1, &F1));
    PetscCall(VecDuplicate(X1, &Fafter)); PetscCall(VecDuplicate(X1, &v)); PetscCall(VecDuplicate(X1, &jv));
 
    /* Immediate repeatability with hidden state poisoned between the two calls. */
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, X1, F0, &pctx));
    simCtx->FluxInSum = 1234.0; simCtx->FluxOutSum = -4321.0;
    simCtx->FarFluxInSum = 77.0; simCtx->FarFluxOutSum = -88.0;
    PetscCall(VecSet(user->Ucat, 7.0)); PetscCall(VecSet(user->lUcat, 7.0));
    PetscCall(MomentumNewtonKrylov_FormResidual(NULL, X1, F1, &pctx));
    PetscCall(CollectiveDifference(F0, F1, &imm_inf, &imm_2));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "PURITY immediate inf=% .3e l2=% .3e (tol inf=%.0e l2=%.0e)\n",
        (double)imm_inf, (double)imm_2, (double)purity_inf_tol, (double)purity_2_tol));
 
    /* Post-MFFD repeatability: real PETSc MFFD products, then reevaluate. */
    {
        SNES snes = NULL; Mat J = NULL; KSP ksp = NULL;
        PetscInt lo, hi; PetscScalar *a = NULL; PetscReal vn;
        PetscCall(VecGetOwnershipRange(v, &lo, &hi));
        PetscCall(VecGetArray(v, &a));
        for (PetscInt i = lo; i < hi; ++i) a[i - lo] = PetscSinReal(0.17 * (PetscReal)(i + 1));
        PetscCall(VecRestoreArray(v, &a));
        PetscCall(VecNorm(v, NORM_2, &vn)); PetscCall(VecScale(v, 1.0 / vn));
        PetscCall(CreateProductionSNES(user, &pctx, PETSC_FALSE, &snes, &J, &ksp));
        PetscCall(MatMFFDSetBase(J, X1, F0));
        for (PetscInt i = 0; i < 40; ++i) PetscCall(MatMult(J, v, jv));
        PetscCall(MomentumNewtonKrylov_FormResidual(NULL, X1, Fafter, &pctx));
        PetscCall(CollectiveDifference(F0, Fafter, &mffd_inf, &mffd_2));
        PetscCall(DifferenceMaxRank(F0, Fafter, &mffd_rank));
        PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
    }
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "PURITY post_mffd inf=% .3e l2=% .3e max_rank=%d\n",
        (double)mffd_inf, (double)mffd_2, (int)mffd_rank));
 
    /* -------- 4.2 full step-2 solve with default classical Gram-Schmidt -------- */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n===== 4.2 full production step-2 solve (default classical GS) =====\n"));
    PetscCall(RunProductionFullSolve(user, &entry, X0, PETSC_FALSE, "CGS-default", &res_cgs));
 
    /* Also drive the ACTUAL production solver wrapper (default GS) and take X_final
     * from it. The wrapper requires user->Rhs unallocated on entry. */
    PetscCall(FpRestore(user, &entry));
    PetscCall(VecDestroy(&user->Rhs));
    prod_ierr = MomentumSolver_NewtonKrylov_BoundaryFixedpointPrivateCopy(user, NULL, NULL);
    PetscCall(VecDuplicate(user->Ucont, &user->Rhs));
    PetscCall(VecDuplicate(user->Ucont, &Xfinal)); PetscCall(VecCopy(user->Ucont, Xfinal));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "PROD_SOLVER ierr=%d mom_last_converged=%d\n",
        (int)prod_ierr, (int)simCtx->mom_last_converged));
    {
        PetscReal dinf, d2; PetscMPIInt r;
        PetscCall(CollectiveDifference(res_cgs.Xfinal, Xfinal, &dinf, &d2));
        PetscCall(DifferenceMaxRank(res_cgs.Xfinal, Xfinal, &r));
        PetscCall(PetscPrintf(PETSC_COMM_WORLD,
            "PROD_SOLVER test_snes_vs_wrapper inf=% .3e l2=% .3e max_rank=%d\n",
            (double)dinf, (double)d2, (int)r));
    }
 
    /* -------- 4.3 boundary-map compatibility at the converged solution -------- */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n===== 4.3 boundary-map compatibility (3 vs 24 passes) =====\n"));
    PetscCall(VecDuplicate(Xfinal, &F3)); PetscCall(VecDuplicate(Xfinal, &F24));
    PetscCall(FpRestore(user, &entry));
    PetscCall(ResidualNCalls(&pctx, Xfinal, 1, F3));    /* 3 passes == production callback */
    PetscCall(VecNorm(F3, NORM_2, &f3_2));
    PetscCall(FpRestore(user, &entry));
    PetscCall(ResidualNCalls(&pctx, Xfinal, 8, F24));   /* 24 passes */
    PetscCall(VecNorm(F24, NORM_2, &f24_2));
    PetscCall(VecNorm(F24, NORM_INFINITY, &f24_inf));
    PetscCall(CollectiveDifference(F3, F24, &f24m3_inf, &f24m3_2));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "COMPAT F3_l2=% .6e F24_l2=% .6e F24_inf=% .6e F24_minus_F3_l2=% .6e (tol %.0e)\n",
        (double)f3_2, (double)f24_2, (double)f24_inf, (double)f24m3_2, (double)compat_tol));
 
    /* -------- 4.4 projection compatibility -------- */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n===== 4.4 projection compatibility =====\n"));
    PetscCall(FpRestore(user, &entry));
    PetscCall(VecCopy(Xfinal, user->Ucont));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, stag));
    PetscCall(Contra2Cart(user));
    PetscCall(SynchronizePeriodicCellFields(user, 1, cell));
    PetscCall(UpdateLocalGhosts(user, "Ucat"));
    PetscCall(ApplyBoundaryConditions(user));
    PetscCall(PoissonSolver_MG(&simCtx->usermg));
    PetscCall(UpdatePressure(user));
    PetscCall(Projection(user));
    PetscCall(UpdateLocalGhosts(user, "P"));
    PetscCall(ComputeDivergence(user));
    max_div = simCtx->MaxDiv;
    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
        "PROJECTION max_divergence=% .6e net_divergence_sum=% .6e (tol %.0e)\n",
        (double)max_div, (double)simCtx->summationRHS, (double)div_tol));
 
    /* -------- consolidated verdicts (printed before assertions) -------- */
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nREGRESSION_VERDICTS (ranks=%d)\n", (int)world));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PRODUCTION_PURITY: %s\n",
        (imm_inf <= purity_inf_tol && imm_2 <= purity_2_tol &&
         mffd_inf <= purity_inf_tol && mffd_2 <= purity_2_tol) ? "PASS" : "FAIL"));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FULL_STEP2_DEFAULT_GS: %s\n",
        res_cgs.converged ? "PASS" : "FAIL"));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PRODUCTION_SOLVER_WRAPPER: %s\n",
        (prod_ierr == 0 && simCtx->mom_last_converged) ? "PASS" : "FAIL"));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  THREE_PASS_24_PASS_COMPATIBLE: %s\n",
        (f24_2 <= compat_tol) ? "PASS" : "FAIL"));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PROJECTION_DIVERGENCE_FREE: %s\n",
        (max_div <= div_tol) ? "PASS" : "FAIL"));
 
    /* -------- assertions -------- */
    PetscCheck(imm_inf <= purity_inf_tol && imm_2 <= purity_2_tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "production residual not immediately repeatable (inf=%g l2=%g)", (double)imm_inf, (double)imm_2);
    PetscCheck(mffd_inf <= purity_inf_tol && mffd_2 <= purity_2_tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "production residual changed after MFFD products (inf=%g l2=%g)", (double)mffd_inf, (double)mffd_2);
    PetscCheck(res_cgs.converged, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "full step-2 solve (default GS) did not converge (reason=%d)", (int)res_cgs.reason);
    PetscCheck(res_cgs.reason != SNES_DIVERGED_LINE_SEARCH && res_cgs.reason != SNES_DIVERGED_LINEAR_SOLVE,
        PETSC_COMM_WORLD, PETSC_ERR_PLIB, "full step-2 solve diverged via line search or linear solve (reason=%d)",
        (int)res_cgs.reason);
    PetscCheck(prod_ierr == 0 && simCtx->mom_last_converged, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "production Newton-Krylov solver wrapper did not converge (ierr=%d)", (int)prod_ierr);
    PetscCheck(f24_2 <= compat_tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "converged three-pass solution does not satisfy the 24-pass residual (||F24||2=%g)", (double)f24_2);
    PetscCheck(max_div <= div_tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
        "projection of the converged momentum state is not divergence free (max_div=%g)", (double)max_div);
 
    PetscCall(VecDestroy(&F24)); PetscCall(VecDestroy(&F3));
    PetscCall(VecDestroy(&Xfinal)); PetscCall(VecDestroy(&res_cgs.Xfinal));
    PetscCall(VecDestroy(&jv)); PetscCall(VecDestroy(&v));
    PetscCall(VecDestroy(&Fafter)); PetscCall(VecDestroy(&F1)); PetscCall(VecDestroy(&F0));
    PetscCall(VecDestroy(&X1)); PetscCall(VecDestroy(&X0));
    PetscCall(FpDestroy(&entry));
    PetscCall(VecDestroy(&user->Rhs));
    PetscCall(PicurvDestroyRuntimeContext(&simCtx));
    PetscCall(PicurvRemoveTempDir(tmpdir));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

Runs the opt-in Newton–Krylov production regression executable.

Definition at line 542 of file test_momentum_newton_boundary_fixedpoint.c.

{
    PetscErrorCode ierr;
    const PicurvTestCase cases[] = {
        {"newton-krylov-production-regression", TestNewtonKrylovProductionRegression}
    };
 
    ierr = PetscInitialize(&argc, &argv, NULL, "PICurv Newton-Krylov production regression");
    if (ierr) return (int)ierr;
    ierr = PicurvRunTests("unit-momentum-newton-boundary-fixedpoint", cases, 1);
    if (PetscFinalize()) return 1;
    return (int)ierr;
}

Here is the call graph for this function: