momentum_newton_krylov.c File Reference

#include "momentumsolvers.h"

Include dependency graph for momentum_newton_krylov.c:

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

Go to the source code of this file.

Data Structures
struct	MomentumNewtonJacobian
struct	MomentumPreconditionerDescription
struct	MomentumPreconditionerModelOps
struct	MomentumPreconditionerEngine
struct	MomentumNewtonKrylovContext

Macros
#define	__FUNCT__   "MomentumNewtonKrylov_Validate"
#define	__FUNCT__   "MomentumNewtonKrylov_ReadLinearizationConfig"
#define	__FUNCT__   "FrozenMomentumJacobian_DescribePointBlock"
#define	__FUNCT__   "FrozenMomentumJacobian_AssemblePointBlocks"
#define	__FUNCT__   "MomentumPreconditionerEngine_ApplyConstraintRows"
#define	__FUNCT__   "MomentumNewtonJacobian_Create"
#define	__FUNCT__   "MomentumPreconditionerEngine_Create"
#define	__FUNCT__   "MomentumPreconditionerEngine_Assemble"
#define	__FUNCT__   "MomentumNewtonKrylov_FormJacobian"
#define	__FUNCT__   "MomentumNewtonKrylov_ClassifyRow"
#define	__FUNCT__   "MomentumPreconditionerEngine_CreateExactPointBlockMatrix"
#define	__FUNCT__   "MomentumNewtonKrylov_ApplyConstraints"
#define	__FUNCT__   "MomentumNewtonKrylov_FormResidual"
#define	__FUNCT__   "MomentumSolver_NewtonKrylov"

Enumerations
enum	MomentumNewtonJacobianType { MOM_NK_JACOBIAN_FINITE_DIFFERENCE }
enum	MomentumNewtonFiniteDifferenceMode { MOM_NK_FD_MODE_MATRIX_FREE }
enum	MomentumPreconditionerModel { MOM_NK_PC_MODEL_NONE , MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN }
enum	MomentumPreconditionerStructure { MOM_NK_PC_STRUCTURE_NONE , MOM_NK_PC_STRUCTURE_POINT_BLOCK }
enum	MomentumNewtonKrylovRowType { MOM_NK_ROW_PHYSICAL = 0 , MOM_NK_ROW_FIXED_CONDITIONED , MOM_NK_ROW_FIXED_HOMOGENEOUS , MOM_NK_ROW_PERIODIC_DUPLICATE }

Functions
static PetscErrorCode	MomentumNewtonKrylov_Validate (UserCtx *user)
	Rejects configurations outside the audited version-one feature set.
static PetscErrorCode	MomentumNewtonKrylov_FormResidual (SNES snes, Vec X, Vec F, void *vctx)
	Adapts a PETSc trial vector to the existing momentum residual path.
static PetscErrorCode	MomentumNewtonKrylov_Monitor (SNES snes, PetscInt iteration, PetscReal norm, void *vctx)
	Captures SNES iteration norms and optionally writes PICurv history rows.
static void	MomentumNewtonKrylov_OpenHistory (MomentumNewtonKrylovContext *ctx)
	Opens the optional rank-zero Newton iteration-history file.
static void	MomentumNewtonKrylov_WriteSummary (const MomentumNewtonKrylovContext *ctx, SNESConvergedReason reason, PetscInt nonlinear_its, PetscInt function_evals, PetscInt linear_its, PetscReal final_norm, PetscBool committed)
	Appends one rank-zero structured Newton result for a physical step.
static PetscErrorCode	MomentumNewtonKrylov_ApplyConstraints (MomentumNewtonKrylovContext *ctx, Vec X, Vec F)
	Replaces every non-independent residual row with an explicit equation.
static MomentumNewtonKrylovRowType	MomentumNewtonKrylov_ClassifyRow (UserCtx *user, PetscInt i, PetscInt j, PetscInt k, PetscInt component, PetscInt *ri, PetscInt *rj, PetscInt *rk)
	Classifies one stored staggered component row and its periodic representative.
static PetscErrorCode	MomentumNewtonKrylov_ReadLinearizationConfig (MomentumNewtonJacobian *jacobian, MomentumPreconditionerDescription *description)
	Reads application-owned Jacobian and preconditioner mathematics.
static PetscErrorCode	MomentumNewtonKrylov_FormJacobian (SNES snes, Vec current_solution, Mat jacobian_operator, Mat preconditioning_matrix, void *vctx)
	Updates the Jacobian and then assembles any separate preconditioning matrix.
static PetscErrorCode	FrozenMomentumJacobian_AssemblePointBlocks (UserCtx *user, Vec current_solution, Mat preconditioning_matrix)
	Inserts only the audited interior frozen-momentum point blocks.
static PetscErrorCode	MomentumPreconditionerEngine_CreateExactPointBlockMatrix (UserCtx *user, Mat *preconditioning_matrix)
	Creates the frozen point-block P matrix with its exact scalar pattern.
static PetscReal	FrozenMomentumJacobian_MetricNormSquared (Cmpnts metric)
	Returns the squared Euclidean norm of one metric vector.
static void	FrozenMomentumJacobian_PointBlock (const SimCtx *simCtx, const Cmpnts ***ucont, const Cmpnts ***csi, const Cmpnts ***eta, const Cmpnts ***zet, const PetscReal ***aj, PetscInt i, PetscInt j, PetscInt k, PetscScalar block[9])
	Returns the audited frozen-momentum point block in modern residual sign.
static PetscErrorCode	FrozenMomentumJacobian_DescribePointBlock (UserCtx *user, MomentumPreconditionerDescription *description)
	Describes the audited frozen-coefficient point-block model.
static PetscErrorCode	MomentumPreconditionerEngine_ApplyConstraintRows (UserCtx *user, Mat preconditioning_matrix)
	Inserts all common fixed, homogeneous, and periodic-duplicate rows.
static PetscErrorCode	MomentumNewtonJacobian_Create (SNES snes, MomentumNewtonJacobian *jacobian)
	Creates the selected Jacobian operator; currently PETSc MFFD only.
static PetscErrorCode	MomentumNewtonJacobian_Update (SNES snes, Vec current_solution, MomentumNewtonJacobian *jacobian)
	Updates the matrix-free finite-difference operator base.
static PetscErrorCode	MomentumNewtonJacobian_Register (SNES snes, MomentumNewtonJacobian *jacobian, MomentumPreconditionerEngine *engine, MomentumNewtonKrylovContext *ctx)
	Registers the application orchestration callback and both SNES matrices.
static PetscErrorCode	MomentumNewtonJacobian_Destroy (MomentumNewtonJacobian *jacobian)
	Destroys a partially or fully created Jacobian operator.
static PetscErrorCode	MomentumPreconditionerEngine_Create (UserCtx *user, Mat jacobian_operator, const MomentumPreconditionerDescription *requested, MomentumPreconditionerEngine *engine)
	Validates a model/structure and creates or aliases its matrix.
static PetscErrorCode	MomentumPreconditionerEngine_Assemble (MomentumPreconditionerEngine *engine, UserCtx *user, Vec current_solution)
	Runs model insertion, common row handling, and final assembly.
static PetscErrorCode	MomentumPreconditionerEngine_ConfigurePetscPC (MomentumPreconditionerEngine *engine, PC pc)
	Applies the validated model/structure-to-PETSc-PC mapping.
static PetscErrorCode	MomentumPreconditionerEngine_ValidatePetscPC (MomentumPreconditionerEngine *engine, PC pc)
	Rejects raw options that select an unvalidated PETSc PC backend.
static PetscErrorCode	MomentumPreconditionerEngine_Destroy (MomentumPreconditionerEngine *engine)
	Destroys only a separately owned preconditioning matrix.
PetscErrorCode	MomentumSolver_NewtonKrylov (UserCtx *user, IBMNodes *ibm, FSInfo *fsi)
	Solves one physical momentum step with matrix-free Newton–Krylov.

Variables
static const MomentumPreconditionerModelOps	frozen_momentum_point_block_ops

---

Data Structure Documentation

MomentumNewtonJacobian

struct MomentumNewtonJacobian

Definition at line 11 of file momentum_newton_krylov.c.

Data Fields
MomentumNewtonJacobianType	type
MomentumNewtonFiniteDifferenceMode	finite_difference_mode
Mat	jacobian_operator

MomentumPreconditionerDescription

struct MomentumPreconditionerDescription

Definition at line 27 of file momentum_newton_krylov.c.

Data Fields
MomentumPreconditionerModel	model
MomentumPreconditionerStructure	structure
PetscInt	block_size
PetscInt	stencil_width

MomentumPreconditionerEngine

struct MomentumPreconditionerEngine

Definition at line 39 of file momentum_newton_krylov.c.

Collaboration diagram for MomentumPreconditionerEngine:

Data Fields
MomentumPreconditionerDescription	description
const MomentumPreconditionerModelOps *	model_ops
Mat	preconditioning_matrix
PetscBool	aliases_jacobian_operator
PetscBool	owns_preconditioning_matrix
const char *	petsc_pc_type

MomentumNewtonKrylovContext

struct MomentumNewtonKrylovContext

Definition at line 48 of file momentum_newton_krylov.c.

Collaboration diagram for MomentumNewtonKrylovContext:

Data Fields
UserCtx *	user
FILE *	history_file
PetscBool	have_initial_norm
PetscReal	initial_norm
MomentumNewtonJacobian	jacobian
MomentumPreconditionerEngine	preconditioning_engine

Macro Definition Documentation

__FUNCT__ [1/14]

#define __FUNCT__   "MomentumNewtonKrylov_Validate"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [2/14]

#define __FUNCT__   "MomentumNewtonKrylov_ReadLinearizationConfig"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [3/14]

#define __FUNCT__   "FrozenMomentumJacobian_DescribePointBlock"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [4/14]

#define __FUNCT__   "FrozenMomentumJacobian_AssemblePointBlocks"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [5/14]

#define __FUNCT__   "MomentumPreconditionerEngine_ApplyConstraintRows"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [6/14]

#define __FUNCT__   "MomentumNewtonJacobian_Create"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [7/14]

#define __FUNCT__   "MomentumPreconditionerEngine_Create"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [8/14]

#define __FUNCT__   "MomentumPreconditionerEngine_Assemble"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [9/14]

#define __FUNCT__   "MomentumNewtonKrylov_FormJacobian"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [10/14]

#define __FUNCT__   "MomentumNewtonKrylov_ClassifyRow"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [11/14]

#define __FUNCT__   "MomentumPreconditionerEngine_CreateExactPointBlockMatrix"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [12/14]

#define __FUNCT__   "MomentumNewtonKrylov_ApplyConstraints"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [13/14]

#define __FUNCT__   "MomentumNewtonKrylov_FormResidual"

Definition at line 203 of file momentum_newton_krylov.c.

__FUNCT__ [14/14]

#define __FUNCT__   "MomentumSolver_NewtonKrylov"

Definition at line 203 of file momentum_newton_krylov.c.

Enumeration Type Documentation

MomentumNewtonJacobianType

enum MomentumNewtonJacobianType

Enumerator
MOM_NK_JACOBIAN_FINITE_DIFFERENCE

Definition at line 3 of file momentum_newton_krylov.c.

             {
    MOM_NK_JACOBIAN_FINITE_DIFFERENCE
} MomentumNewtonJacobianType;

MomentumNewtonFiniteDifferenceMode

enum MomentumNewtonFiniteDifferenceMode

Enumerator
MOM_NK_FD_MODE_MATRIX_FREE

Definition at line 7 of file momentum_newton_krylov.c.

             {
    MOM_NK_FD_MODE_MATRIX_FREE
} MomentumNewtonFiniteDifferenceMode;

MomentumPreconditionerModel

enum MomentumPreconditionerModel

Enumerator
MOM_NK_PC_MODEL_NONE
MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN

Definition at line 17 of file momentum_newton_krylov.c.

             {
    MOM_NK_PC_MODEL_NONE,
    MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN
} MomentumPreconditionerModel;

MomentumPreconditionerStructure

enum MomentumPreconditionerStructure

Enumerator
MOM_NK_PC_STRUCTURE_NONE
MOM_NK_PC_STRUCTURE_POINT_BLOCK

Definition at line 22 of file momentum_newton_krylov.c.

             {
    MOM_NK_PC_STRUCTURE_NONE,
    MOM_NK_PC_STRUCTURE_POINT_BLOCK
} MomentumPreconditionerStructure;

MomentumNewtonKrylovRowType

enum MomentumNewtonKrylovRowType

Enumerator
MOM_NK_ROW_PHYSICAL
MOM_NK_ROW_FIXED_CONDITIONED
MOM_NK_ROW_FIXED_HOMOGENEOUS
MOM_NK_ROW_PERIODIC_DUPLICATE

Definition at line 59 of file momentum_newton_krylov.c.

             {
    MOM_NK_ROW_PHYSICAL = 0,
    MOM_NK_ROW_FIXED_CONDITIONED,
    MOM_NK_ROW_FIXED_HOMOGENEOUS,
    MOM_NK_ROW_PERIODIC_DUPLICATE
} MomentumNewtonKrylovRowType;

Function Documentation

MomentumNewtonKrylov_Validate()

static PetscErrorCode MomentumNewtonKrylov_Validate ( UserCtx *  user )

static

Rejects configurations outside the audited version-one feature set.

Parameters

user	Single-block momentum context to validate.

Returns

PetscErrorCode 0 when the configuration is supported.

Definition at line 209 of file momentum_newton_krylov.c.

{
    SimCtx   *simCtx;
    PetscReal mask_max = 0.0;
    PetscInt  velocity_dof = 0;
 
    PetscFunctionBeginUser;
    PetscCheck(user != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Newton Krylov requires a non-NULL UserCtx.");
    simCtx = user->simCtx;
    PetscCheck(simCtx != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Newton Krylov requires UserCtx::simCtx.");
    PetscCall(DMDAGetInfo(user->fda, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                           &velocity_dof, NULL, NULL, NULL, NULL, NULL));
    PetscCheck(velocity_dof == 3, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "Newton Krylov requires a three-component velocity DMDA (got dof=%d).",
               velocity_dof);
    PetscCheck(simCtx->block_number == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one supports exactly one block (got %d).",
               simCtx->block_number);
    PetscCheck(!simCtx->immersed, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support immersed boundaries.");
    PetscCheck(!simCtx->movefsi && !simCtx->rotatefsi, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support moving or rotating bodies/FSI.");
    PetscCheck(!simCtx->moveframe && !simCtx->rotateframe, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support moving or rotating reference frames.");
    PetscCheck(!simCtx->rans, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support RANS.");
    PetscCheck(!simCtx->clark, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support the Clark model.");
    PetscCheck(!simCtx->TwoD, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support TwoD component masking.");
    PetscCheck(!simCtx->wallfunction, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not support wall functions.");
    for (PetscInt face = 0; face < 6; ++face) {
        const BoundaryFaceConfig *cfg = &user->boundary_faces[face];
        PetscBool supported = PETSC_FALSE;
 
        switch (cfg->handler_type) {
            case BC_HANDLER_WALL_NOSLIP:
                supported = (PetscBool)(cfg->mathematical_type == WALL);
                break;
            case BC_HANDLER_INLET_CONSTANT_VELOCITY:
            case BC_HANDLER_INLET_PARABOLIC:
            case BC_HANDLER_INLET_PROFILE_FROM_FILE:
                supported = (PetscBool)(cfg->mathematical_type == INLET);
                break;
            case BC_HANDLER_OUTLET_CONSERVATION:
                supported = (PetscBool)(cfg->mathematical_type == OUTLET);
                break;
            case BC_HANDLER_PERIODIC_GEOMETRIC:
                supported = (PetscBool)(cfg->mathematical_type == PERIODIC);
                break;
            default:
                supported = PETSC_FALSE;
                break;
        }
        PetscCheck(supported, PETSC_COMM_WORLD, PETSC_ERR_SUP,
                   "Newton Krylov version one does not support boundary face %d with mathematical type %d and handler %d.",
                   face, (PetscInt)cfg->mathematical_type, (PetscInt)cfg->handler_type);
    }
 
    PetscCheck((user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC) ==
               (user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC),
               PETSC_COMM_WORLD, PETSC_ERR_SUP, "Newton Krylov requires paired x-periodic faces.");
    PetscCheck((user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC) ==
               (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC),
               PETSC_COMM_WORLD, PETSC_ERR_SUP, "Newton Krylov requires paired y-periodic faces.");
    PetscCheck((user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC) ==
               (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC),
               PETSC_COMM_WORLD, PETSC_ERR_SUP, "Newton Krylov requires paired z-periodic faces.");
 
    PetscCheck(user->Nvert != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Newton Krylov requires the cell-mask vector Nvert.");
    PetscCall(VecMax(user->Nvert, NULL, &mask_max));
    PetscCheck(mask_max <= 0.1, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not define equations for masked solid cells (max Nvert=%g).",
               (double)mask_max);
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumNewtonKrylov_FormResidual()

static PetscErrorCode MomentumNewtonKrylov_FormResidual ( SNES  snes, Vec  X, Vec  F, void *  vctx  )

static

Adapts a PETSc trial vector to the existing momentum residual path.

A matrix-free SNES residual must be a deterministic function of the trial vector X alone: F(X) may not depend on any state left by a previous residual or MFFD evaluation, or finite-difference Jacobian actions become inconsistent. To honor that contract this callback fully derives the Cartesian velocity state (Ucat/lUcat) from X before the first boundary sweep – see the inline comment below for why ApplyBoundaryConditions()'s own internal reconstruction is not sufficient for the first outlet pass.

State invariants:

• On entry, X is the only input that determines the result; user->Ucont, user->Ucat and their local ghosts are treated as scratch and are fully overwritten from X.
• Supported handlers may overwrite flux totals and other diagnostics on every call, but those values must not affect a later call at the same X; the deterministic seed guarantees this.
• No histories, pressure, viscosity, or controller state advance here.

Side effects: overwrites user->Ucont/lUcont, user->Ucat/lUcat, user->Rhs, the boundary Ubcs targets, and boundary flux/area diagnostics; writes F.

Parameters

snes	Calling nonlinear solver.
X	Trial solution (read-only).
F	Residual output.
vctx	Pointer to MomentumNewtonKrylovContext.

Returns

PetscErrorCode 0 on success.

Definition at line 1066 of file momentum_newton_krylov.c.

{
    MomentumNewtonKrylovContext *ctx = (MomentumNewtonKrylovContext *)vctx;
    UserCtx *user = ctx->user;
    const FieldId staggered_fields[] = {FIELD_ID_UCONT};
    const FieldId cell_fields[] = {FIELD_ID_UCAT};
 
    PetscFunctionBeginUser;
    (void)snes;
    PetscCall(VecCopy(X, user->Ucont));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, staggered_fields));
 
    /* Deterministic pre-boundary Cartesian seed. Establish the full Ucat/lUcat
     * state from the current X before any boundary handler runs:
     *
     *     X -> Ucont/lUcont -> Ucat -> periodic Ucat -> lUcat -> boundaries
     *
     * Why this is required, and why it is NOT redundant with the reconstruction
     * already performed inside ApplyBoundaryConditions():
     *
     *  1. A matrix-free SNES residual must be a deterministic function of X. If
     *     the Cartesian state is left over from a previous residual/MFFD call,
     *     F(X) becomes history dependent and the finite-difference Jacobian
     *     action Jv = (F(X+hv)-F(X))/h is invalidated.
     *  2. The conservation-outlet handler reads lUcat during the FIRST boundary
     *     sweep (it measures the uncorrected outflow and builds the outlet
     *     profile from the Cartesian field). Without this seed it would read the
     *     stale lUcat from the preceding evaluation.
     *  3. ApplyBoundaryConditions() does reconstruct Ucat/lUcat, but only AFTER
     *     each handler sweep (Contra2Cart runs after BoundarySystem_ExecuteStep
     *     within every pass). Those internal updates therefore prepare passes 2
     *     and 3 -- they cannot prepare the very first outlet read of pass 1.
     *  4. Contra2Cart() rebuilds the global Ucat interior from the current
     *     lUcont, but it does not by itself refresh lUcat (nor lUcont; that was
     *     done by the SynchronizePeriodicStaggeredFields call above).
     *  5. SynchronizePeriodicCellFields(FIELD_ID_UCAT) must run before the ghost
     *     scatter so periodic duplicate planes are finalized consistently (it is
     *     a no-op when no direction is periodic, as on the straight duct).
     *  6. UpdateLocalGhosts(FIELD_ID_UCAT) is required because the outlet handler reads
     *     lUcat -- the local ghosted vector -- not merely the global Ucat.
     *
     * Do NOT "simplify" this to a bare Contra2Cart(user), and do NOT delete it
     * as apparently redundant with ApplyBoundaryConditions(): the three internal
     * boundary passes remain necessary (they refresh the Cartesian state after
     * each boundary correction), but only this sequence makes pass 1's input a
     * deterministic function of X. */
    PetscCall(Contra2Cart(user));
    PetscCall(SynchronizePeriodicCellFields(user, 1, cell_fields));
    PetscCall(UpdateLocalGhosts(user, FIELD_ID_UCAT));
 
    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:

MomentumNewtonKrylov_Monitor()

static PetscErrorCode MomentumNewtonKrylov_Monitor ( SNES  snes, PetscInt  iteration, PetscReal  norm, void *  vctx  )

static

Captures SNES iteration norms and optionally writes PICurv history rows.

SNES supplies the already-computed norm, so this monitor never causes an additional nonlinear residual evaluation. PETSc monitors selected through -mom_nk_snes_monitor remain independent and may run alongside this callback.

Definition at line 98 of file momentum_newton_krylov.c.

{
    MomentumNewtonKrylovContext *ctx = (MomentumNewtonKrylovContext *)vctx;
    SimCtx *simCtx = ctx->user->simCtx;
 
    (void)snes;
    PetscFunctionBeginUser;
    if (iteration == 0 && !ctx->have_initial_norm) {
        ctx->initial_norm = norm;
        ctx->have_initial_norm = PETSC_TRUE;
    }
    if (ctx->history_file) {
        (void)fprintf(ctx->history_file,
                      "step: %d | block: %d | newton: %d | nonlinear_norm: %.16e\n",
                      (int)simCtx->step, (int)ctx->user->_this, (int)iteration,
                      (double)norm);
        (void)fflush(ctx->history_file);
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumNewtonKrylov_OpenHistory()

static void MomentumNewtonKrylov_OpenHistory ( MomentumNewtonKrylovContext *  ctx )

static

Opens the optional rank-zero Newton iteration-history file.

Definition at line 121 of file momentum_newton_krylov.c.

{
    SimCtx *simCtx = ctx->user->simCtx;
    char path[PETSC_MAX_PATH_LEN + 128];
    const char *mode;
 
    if (!simCtx->mom_nk_monitor_history || simCtx->rank != 0) return;
    if (PetscSNPrintf(path, sizeof(path),
                     "%s/Momentum_Solver_Newton_Krylov_History_Block_%d.log",
                     simCtx->log_dir, (int)ctx->user->_this)) return;
    mode = (simCtx->step == simCtx->StartStep + 1 && !simCtx->continueMode) ? "w" : "a";
    ctx->history_file = fopen(path, mode);
    if (!ctx->history_file) {
        LOG(GLOBAL, LOG_WARNING, "Could not open Newton iteration-history log '%s'.\n", path);
        return;
    }
    if (mode[0] == 'w') {
        (void)fprintf(ctx->history_file,
                      "# step | block | Newton iteration | nonlinear residual norm\n");
    } else if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
        (void)fprintf(ctx->history_file, "# Continuation from step %d\n", (int)simCtx->StartStep);
    }
}

Here is the caller graph for this function:

MomentumNewtonKrylov_WriteSummary()

static void MomentumNewtonKrylov_WriteSummary ( const MomentumNewtonKrylovContext *  ctx, SNESConvergedReason  reason, PetscInt  nonlinear_its, PetscInt  function_evals, PetscInt  linear_its, PetscReal  final_norm, PetscBool  committed  )

static

Appends one rank-zero structured Newton result for a physical step.

File failures are deliberately diagnostic-only: rollback and PETSc cleanup must retain their original error behavior.

Definition at line 150 of file momentum_newton_krylov.c.

{
    SimCtx *simCtx = ctx->user->simCtx;
    char path[PETSC_MAX_PATH_LEN + 128];
    const char *mode;
    const char *reason_name;
    FILE *file;
 
    if (simCtx->rank != 0) return;
    if (PetscSNPrintf(path, sizeof(path),
                     "%s/Momentum_Solver_Newton_Krylov_Summary_Block_%d.log",
                     simCtx->log_dir, (int)ctx->user->_this)) return;
    mode = (simCtx->step == simCtx->StartStep + 1 && !simCtx->continueMode) ? "w" : "a";
    file = fopen(path, mode);
    if (!file) {
        LOG(GLOBAL, LOG_WARNING, "Could not open Newton summary log '%s'.\n", path);
        return;
    }
    if (mode[0] == 'w') {
        (void)fprintf(file,
                      "# step | block | solver | Jacobian | preconditioner | SNES reason | "
                      "reason code | Newton iterations | "
                      "residual evaluations | Krylov iterations | initial nonlinear norm | "
                      "final nonlinear norm | state\n");
    } else if (simCtx->continueMode && simCtx->step == simCtx->StartStep + 1) {
        (void)fprintf(file, "# Continuation from step %d\n", (int)simCtx->StartStep);
    }
    reason_name = reason == SNES_CONVERGED_ITERATING
                    ? "SNES_CONVERGED_ITERATING" : SNESConvergedReasons[reason];
    (void)fprintf(file,
                  "step: %d | block: %d | solver: Newton Krylov | "
                  "Jacobian: finite_difference / matrix_free | Preconditioner: %s | "
                  "reason: %s | reason_code: %d | "
                  "newton: %d | evals: %d | krylov: %d | initial: ",
                  (int)simCtx->step, (int)ctx->user->_this,
                  ctx->preconditioning_engine.description.model == MOM_NK_PC_MODEL_NONE ?
                      "none" : "frozen_momentum_jacobian / point_block",
                  reason_name, (int)reason,
                  (int)nonlinear_its, (int)function_evals, (int)linear_its);
    if (ctx->have_initial_norm) (void)fprintf(file, "%.16e", (double)ctx->initial_norm);
    else                        (void)fprintf(file, "unavailable");
    (void)fprintf(file, " | final: %.16e | state: %s\n", (double)final_norm,
                  committed ? "committed" : "rolled_back");
    (void)fclose(file);
}

Here is the caller graph for this function:

MomentumNewtonKrylov_ApplyConstraints()

static PetscErrorCode MomentumNewtonKrylov_ApplyConstraints ( MomentumNewtonKrylovContext *  ctx, Vec  X, Vec  F  )

static

Replaces every non-independent residual row with an explicit equation.

Conditioned face-normal rows use F=X-Uconditioned, unconditioned legacy dummy/tangential rows use F=X, and periodic duplicates use Fdup=Xdup-Xrep. These equations prevent the zero Jacobian rows produced by simply retaining EnforceRHSBoundaryConditions() zeros in a matrix-free Newton operator. Immersed, masked, TwoD, and interface rows are rejected before this callback is installed.

Parameters

ctx	Active solve context.
X	Unconditioned PETSc trial state.
F	Residual vector to update in place.

Returns

PetscErrorCode 0 on success.

Definition at line 986 of file momentum_newton_krylov.c.

{
    UserCtx       *user = ctx->user;
    DMDALocalInfo  info = user->info;
    Vec            local_x = NULL;
    Cmpnts       ***x = NULL, ***conditioned = NULL, ***f = NULL, ***lx = NULL;
    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;
 
    PetscFunctionBeginUser;
    PetscCall(DMGetLocalVector(user->fda, &local_x));
    PetscCall(DMGlobalToLocalBegin(user->fda, X, INSERT_VALUES, local_x));
    PetscCall(DMGlobalToLocalEnd(user->fda, X, INSERT_VALUES, local_x));
    PetscCall(DMDAVecGetArrayRead(user->fda, X, &x));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &conditioned));
    PetscCall(DMDAVecGetArray(user->fda, F, &f));
    PetscCall(DMDAVecGetArrayRead(user->fda, local_x, &lx));
 
    for (PetscInt k = zs; k < ze; ++k) {
        for (PetscInt j = ys; j < ye; ++j) {
            for (PetscInt i = xs; i < xe; ++i) {
                PetscScalar *fv = &f[k][j][i].x;
                const PetscScalar *xv = &x[k][j][i].x;
                const PetscScalar *cv = &conditioned[k][j][i].x;
 
                for (PetscInt component = 0; component < 3; ++component) {
                    PetscInt ri, rj, rk;
                    MomentumNewtonKrylovRowType row = MomentumNewtonKrylov_ClassifyRow(
                        user, i, j, k, component, &ri, &rj, &rk);
                    const PetscScalar *rv = &lx[rk][rj][ri].x;
 
                    if (row == MOM_NK_ROW_FIXED_CONDITIONED) fv[component] = xv[component] - cv[component];
                    else if (row == MOM_NK_ROW_FIXED_HOMOGENEOUS) fv[component] = xv[component];
                    else if (row == MOM_NK_ROW_PERIODIC_DUPLICATE) fv[component] = xv[component] - rv[component];
                }
            }
        }
    }
 
    PetscCall(DMDAVecRestoreArrayRead(user->fda, local_x, &lx));
    PetscCall(DMDAVecRestoreArray(user->fda, F, &f));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &conditioned));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, X, &x));
    PetscCall(DMRestoreLocalVector(user->fda, &local_x));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumNewtonKrylov_ClassifyRow()

static MomentumNewtonKrylovRowType MomentumNewtonKrylov_ClassifyRow ( UserCtx *  user, PetscInt  i, PetscInt  j, PetscInt  k, PetscInt  component, PetscInt *  ri, PetscInt *  rj, PetscInt *  rk  )

static

Classifies one stored staggered component row and its periodic representative.

Negative nonperiodic planes and positive dummy planes are zeroed by the legacy residual treatment. Only a face-normal row actually written by its boundary handler is conditioned; the remaining zeroed rows use F=X. Periodic endpoint rows use the same wrapped representatives as SynchronizePeriodicStaggeredFields(). At periodic/nonperiodic intersections the periodic equation owns rows that a shrunken boundary-handler loop does not condition.

Parameters

user	Block context containing boundary metadata and global dimensions.
i	Global i index.
j	Global j index.
k	Global k index.
component	Stored component, 0=x, 1=y, 2=z.
ri	Returned periodic representative i index.
rj	Returned periodic representative j index.
rk	Returned periodic representative k index.

Returns

Exact version-one row category.

Definition at line 769 of file momentum_newton_krylov.c.

{
    const PetscInt mx = user->info.mx, my = user->info.my, mz = user->info.mz;
    const PetscInt coord[3] = {i, j, k};
    const PetscInt size[3] = {mx, my, mz};
    const BCFace neg_face[3] = {BC_FACE_NEG_X, BC_FACE_NEG_Y, BC_FACE_NEG_Z};
    PetscBool periodic[3], periodic_duplicate = PETSC_FALSE;
    PetscBool residual_zeroed = PETSC_FALSE, conditioned = PETSC_FALSE;
 
    *ri = i; *rj = j; *rk = k;
    for (PetscInt axis = 0; axis < 3; ++axis) {
        periodic[axis] = (PetscBool)(
            user->boundary_faces[neg_face[axis]].mathematical_type == PERIODIC);
        if (periodic[axis] && coord[axis] == 0) {
            periodic_duplicate = PETSC_TRUE;
            if (axis == 0) *ri = -2;
            else if (axis == 1) *rj = -2;
            else *rk = -2;
        }
        if (periodic[axis] && coord[axis] == size[axis] - 1) {
            periodic_duplicate = PETSC_TRUE;
            if (axis == 0) *ri = mx + 1;
            else if (axis == 1) *rj = my + 1;
            else *rk = mz + 1;
        }
 
        if (!periodic[axis] && coord[axis] == 0) residual_zeroed = PETSC_TRUE;
        if (coord[axis] == size[axis] - 1) residual_zeroed = PETSC_TRUE;
        if (!periodic[axis] && coord[axis] == size[axis] - 2 && component == axis)
            residual_zeroed = PETSC_TRUE;
    }
 
    if (!periodic[component] &&
        (coord[component] == 0 || coord[component] == size[component] - 2)) {
        PetscBool tangential_interior = PETSC_TRUE;
        for (PetscInt axis = 0; axis < 3; ++axis) {
            if (axis == component) continue;
            if (coord[axis] < 1 || coord[axis] > size[axis] - 2)
                tangential_interior = PETSC_FALSE;
        }
        conditioned = tangential_interior;
    }
 
    if (conditioned) return MOM_NK_ROW_FIXED_CONDITIONED;
    if (periodic_duplicate) return MOM_NK_ROW_PERIODIC_DUPLICATE;
    if (residual_zeroed) return MOM_NK_ROW_FIXED_HOMOGENEOUS;
    return MOM_NK_ROW_PHYSICAL;
}

Here is the caller graph for this function:

MomentumNewtonKrylov_ReadLinearizationConfig()

static PetscErrorCode MomentumNewtonKrylov_ReadLinearizationConfig ( MomentumNewtonJacobian *  jacobian, MomentumPreconditionerDescription *  description  )

static

Reads application-owned Jacobian and preconditioner mathematics.

Definition at line 295 of file momentum_newton_krylov.c.

{
    char type[48] = "finite_difference";
    char finite_difference_mode[32] = "matrix_free";
    char model[48] = "none";
    char structure[32] = "none";
    PetscBool set = PETSC_FALSE, match = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCheck(jacobian != NULL && description != NULL, PETSC_COMM_SELF,
               PETSC_ERR_ARG_NULL, "Newton Krylov linearization configuration is NULL.");
    PetscCall(PetscOptionsGetString(NULL, NULL, "-mom_nk_jacobian_type",
                                    type, sizeof(type), &set));
    PetscCall(PetscStrcasecmp(type, "finite_difference", &match));
    PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "-mom_nk_jacobian_type must be 'finite_difference' (got '%s').", type);
    jacobian->type = MOM_NK_JACOBIAN_FINITE_DIFFERENCE;
    PetscCall(PetscOptionsGetString(NULL, NULL, "-mom_nk_jacobian_fd_mode",
                                    finite_difference_mode,
                                    sizeof(finite_difference_mode), &set));
    PetscCall(PetscStrcasecmp(finite_difference_mode, "matrix_free", &match));
    PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "-mom_nk_jacobian_fd_mode must be 'matrix_free' (got '%s').",
               finite_difference_mode);
    jacobian->finite_difference_mode = MOM_NK_FD_MODE_MATRIX_FREE;
 
    PetscCall(PetscOptionsGetString(NULL, NULL, "-mom_nk_preconditioner_model",
                                    model, sizeof(model), &set));
    PetscCall(PetscStrcasecmp(model, "none", &match));
    if (match) description->model = MOM_NK_PC_MODEL_NONE;
    if (!match) {
        PetscCall(PetscStrcasecmp(model, "frozen_momentum_jacobian", &match));
        if (match) description->model = MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN;
    }
    PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "-mom_nk_preconditioner_model must be 'none' or "
               "'frozen_momentum_jacobian' (got '%s').", model);
    PetscCall(PetscOptionsGetString(NULL, NULL, "-mom_nk_preconditioner_structure",
                                    structure, sizeof(structure), &set));
    PetscCall(PetscStrcasecmp(structure, "none", &match));
    if (match) description->structure = MOM_NK_PC_STRUCTURE_NONE;
    if (!match) {
        PetscCall(PetscStrcasecmp(structure, "point_block", &match));
        if (match) description->structure = MOM_NK_PC_STRUCTURE_POINT_BLOCK;
    }
    PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
               "-mom_nk_preconditioner_structure must be 'none' or 'point_block' "
               "(got '%s').", structure);
    PetscCheck((description->model == MOM_NK_PC_MODEL_NONE &&
                description->structure == MOM_NK_PC_STRUCTURE_NONE) ||
               (description->model == MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN &&
                description->structure == MOM_NK_PC_STRUCTURE_POINT_BLOCK),
               PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Unsupported Newton Krylov preconditioner model/structure combination: "
               "model='%s', structure='%s'.", model, structure);
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumNewtonKrylov_FormJacobian()

static PetscErrorCode MomentumNewtonKrylov_FormJacobian ( SNES  snes, Vec  current_solution, Mat  jacobian_operator, Mat  preconditioning_matrix, void *  vctx  )

static

Updates the Jacobian and then assembles any separate preconditioning matrix.

Definition at line 734 of file momentum_newton_krylov.c.

{
    MomentumNewtonKrylovContext *ctx = (MomentumNewtonKrylovContext *)vctx;
    PetscFunctionBeginUser;
    (void)jacobian_operator;
    (void)preconditioning_matrix;
    PetscCall(MomentumNewtonJacobian_Update(snes, current_solution, &ctx->jacobian));
    PetscCall(MomentumPreconditionerEngine_Assemble(
        &ctx->preconditioning_engine, ctx->user, current_solution));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FrozenMomentumJacobian_AssemblePointBlocks()

static PetscErrorCode FrozenMomentumJacobian_AssemblePointBlocks ( UserCtx *  user, Vec  current_solution, Mat  preconditioning_matrix  )

static

Inserts only the audited interior frozen-momentum point blocks.

Definition at line 485 of file momentum_newton_krylov.c.

{
    DMDALocalInfo info = user->info;
    Cmpnts ***ucont = NULL, ***csi = NULL, ***eta = NULL, ***zet = NULL;
    PetscReal ***aj = NULL;
    PetscErrorCode ierr = PETSC_SUCCESS, cleanup_ierr;
 
    PetscFunctionBeginUser;
    /*
     * current_solution is the current SNES trial Ucont Vec: it is layout-compatible
     * with user->fda/user->Ucont, but is not necessarily the canonical user->Ucont
     * selected by UpdateLocalGhosts(FIELD_ID_UCONT). Scatter that trial Vec directly with
     * user->fda so PETSc applies its MPI ownership, periodic topology, component
     * ordering, and ghost mapping without canonical-field synchronization or mutation
     * of current_solution. UpdateLocalGhosts additionally repairs the component-normal
     * staggered buffers Uxi(i=-1,mx), Ueta(j=-1,my), and Uzeta(k=-1,mz); the audited
     * point-block velocity stencil reads none of those planes, so the repair cannot
     * change a coefficient. Revisit this choice and extend the periodic-localization
     * tests if the model stencil is expanded to read any repaired normal buffer.
     */
    ierr = DMGlobalToLocalBegin(user->fda, current_solution, INSERT_VALUES, user->lUcont);
    if (ierr) goto cleanup;
    ierr = DMGlobalToLocalEnd(user->fda, current_solution, INSERT_VALUES, user->lUcont);
    if (ierr) goto cleanup;
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcont, &ucont); if (ierr) goto cleanup;
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi); if (ierr) goto cleanup;
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta); if (ierr) goto cleanup;
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet); if (ierr) goto cleanup;
    ierr = DMDAVecGetArrayRead(user->da, user->lAj, &aj); if (ierr) goto cleanup;
    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) {
                for (PetscInt component = 0; component < 3; ++component) {
                    MatStencil row = {.i = i, .j = j, .k = k, .c = component};
                    PetscInt ri, rj, rk;
                    MomentumNewtonKrylovRowType type = MomentumNewtonKrylov_ClassifyRow(
                        user, i, j, k, component, &ri, &rj, &rk);
                    if (type == MOM_NK_ROW_PHYSICAL) {
                        PetscScalar block[9];
                        MatStencil cols[3] = {
                            {.i = i, .j = j, .k = k, .c = 0},
                            {.i = i, .j = j, .k = k, .c = 1},
                            {.i = i, .j = j, .k = k, .c = 2}
                        };
                        FrozenMomentumJacobian_PointBlock(user->simCtx, (const Cmpnts ***)ucont,
                            (const Cmpnts ***)csi, (const Cmpnts ***)eta, (const Cmpnts ***)zet,
                            (const PetscReal ***)aj, i, j, k, block);
                        ierr = MatSetValuesStencil(preconditioning_matrix, 1, &row, 3, cols,
                                                  &block[3 * component], INSERT_VALUES);
                        if (ierr) goto cleanup;
                    }
                }
            }
        }
    }
cleanup:
    if (aj) { cleanup_ierr = DMDAVecRestoreArrayRead(user->da, user->lAj, &aj); if (!ierr) ierr = cleanup_ierr; }
    if (zet) { cleanup_ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet); if (!ierr) ierr = cleanup_ierr; }
    if (eta) { cleanup_ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta); if (!ierr) ierr = cleanup_ierr; }
    if (csi) { cleanup_ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi); if (!ierr) ierr = cleanup_ierr; }
    if (ucont) { cleanup_ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcont, &ucont); if (!ierr) ierr = cleanup_ierr; }
    PetscFunctionReturn(ierr);
}

Here is the call graph for this function:

MomentumPreconditionerEngine_CreateExactPointBlockMatrix()

static PetscErrorCode MomentumPreconditionerEngine_CreateExactPointBlockMatrix ( UserCtx *  user, Mat *  preconditioning_matrix  )

static

Creates the frozen point-block P matrix with its exact scalar pattern.

Row and column ownership follows the velocity DMDA global vector. Physical rows reserve the three same-point components, fixed rows reserve only their diagonal, and periodic duplicate rows additionally reserve their wrapped representative. DMDA AO, local mapping, and stencil metadata are retained so the existing insertion paths keep their exact ordering.

Definition at line 830 of file momentum_newton_krylov.c.

{
    DMDALocalInfo info;
    ISLocalToGlobalMapping local_to_global = NULL;
    Mat matrix = NULL;
    MPI_Comm comm;
    PetscInt local_size, global_size, ownership_start, ownership_end;
    PetscInt *diagonal_nnz = NULL, *offdiagonal_nnz = NULL;
    PetscInt ghost_starts[4] = {0, 0, 0, 0}, ghost_sizes[3] = {0, 0, 0};
    PetscErrorCode ierr = PETSC_SUCCESS, cleanup_ierr;
 
    PetscFunctionBeginUser;
    PetscCheck(preconditioning_matrix != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Point-block matrix output is NULL.");
    *preconditioning_matrix = NULL;
    comm = PetscObjectComm((PetscObject)user->fda);
    PetscCall(DMDAGetLocalInfo(user->fda, &info));
    PetscCall(VecGetLocalSize(user->Ucont, &local_size));
    PetscCall(VecGetSize(user->Ucont, &global_size));
    PetscCall(VecGetOwnershipRange(user->Ucont, &ownership_start, &ownership_end));
    PetscCheck(local_size == ownership_end - ownership_start, comm, PETSC_ERR_PLIB,
               "Velocity ownership range does not match its local size.");
    ierr = PetscCalloc2(local_size, &diagonal_nnz,
                        local_size, &offdiagonal_nnz); if (ierr) goto cleanup;
    ierr = DMGetLocalToGlobalMapping(user->fda, &local_to_global); if (ierr) goto cleanup;
    ierr = DMDAGetGhostCorners(user->fda,
                               &ghost_starts[0], &ghost_starts[1], &ghost_starts[2],
                               &ghost_sizes[0], &ghost_sizes[1], &ghost_sizes[2]);
    if (ierr) goto cleanup;
 
    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) {
                for (PetscInt component = 0; component < 3; ++component) {
                    PetscInt ri, rj, rk, column_count;
                    MatStencil row_stencil = {.i = i, .j = j, .k = k, .c = component};
                    MatStencil column_stencils[3];
                    PetscInt row_local, row, column_locals[3], columns[3];
                    MomentumNewtonKrylovRowType type = MomentumNewtonKrylov_ClassifyRow(
                        user, i, j, k, component, &ri, &rj, &rk);
 
                    if (type == MOM_NK_ROW_PHYSICAL) {
                        column_count = 3;
                        for (PetscInt column_component = 0; column_component < 3;
                             ++column_component) {
                            column_stencils[column_component] = (MatStencil){
                                .i = i, .j = j, .k = k, .c = column_component
                            };
                        }
                    } else {
                        column_count = 1;
                        column_stencils[0] = row_stencil;
                        if (type == MOM_NK_ROW_PERIODIC_DUPLICATE) {
                            column_stencils[column_count++] = (MatStencil){
                                .i = ri, .j = rj, .k = rk, .c = component
                            };
                        }
                    }
                    row_local = component + 3 * (
                        (i - ghost_starts[0]) + ghost_sizes[0] * (
                        (j - ghost_starts[1]) + ghost_sizes[1] *
                        (k - ghost_starts[2])));
                    for (PetscInt column_index = 0; column_index < column_count;
                         ++column_index) {
                        const MatStencil column = column_stencils[column_index];
                        if (!(column.i >= ghost_starts[0] &&
                              column.i < ghost_starts[0] + ghost_sizes[0] &&
                              column.j >= ghost_starts[1] &&
                              column.j < ghost_starts[1] + ghost_sizes[1] &&
                              column.k >= ghost_starts[2] &&
                              column.k < ghost_starts[2] + ghost_sizes[2])) {
                            ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__,
                                              PETSC_ERR_ARG_OUTOFRANGE,
                                              PETSC_ERROR_INITIAL,
                                              "Point-block column lies outside the DMDA ghost stencil.");
                            goto cleanup;
                        }
                        column_locals[column_index] = column.c + 3 * (
                            (column.i - ghost_starts[0]) + ghost_sizes[0] * (
                            (column.j - ghost_starts[1]) + ghost_sizes[1] *
                            (column.k - ghost_starts[2])));
                    }
                    ierr = ISLocalToGlobalMappingApply(local_to_global, 1,
                                                        &row_local, &row);
                    if (ierr) goto cleanup;
                    ierr = ISLocalToGlobalMappingApply(local_to_global, column_count,
                                                        column_locals, columns);
                    if (ierr) goto cleanup;
                    if (!(row >= ownership_start && row < ownership_end)) {
                        ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__,
                                          PETSC_ERR_PLIB, PETSC_ERROR_INITIAL,
                                          "DMDA-mapped point-block row is not locally owned.");
                        goto cleanup;
                    }
                    for (PetscInt column_index = 0; column_index < column_count;
                         ++column_index) {
                        PetscBool duplicate = PETSC_FALSE;
                        for (PetscInt previous = 0; previous < column_index; ++previous)
                            if (columns[previous] == columns[column_index]) duplicate = PETSC_TRUE;
                        if (columns[column_index] < 0) {
                            ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__,
                                              PETSC_ERR_PLIB, PETSC_ERROR_INITIAL,
                                              "DMDA-mapped point-block column is invalid.");
                            goto cleanup;
                        }
                        if (duplicate) continue;
                        if (columns[column_index] >= ownership_start &&
                            columns[column_index] < ownership_end)
                            ++diagonal_nnz[row - ownership_start];
                        else
                            ++offdiagonal_nnz[row - ownership_start];
                    }
                }
            }
        }
    }
 
    ierr = MatCreateAIJ(comm, local_size, local_size, global_size, global_size,
                        0, diagonal_nnz, 0, offdiagonal_nnz, &matrix);
    if (ierr) goto cleanup;
    ierr = MatSetBlockSize(matrix, 3); if (ierr) goto cleanup;
    ierr = MatSetLocalToGlobalMapping(matrix, local_to_global, local_to_global);
    if (ierr) goto cleanup;
    ierr = MatSetStencil(matrix, 3, ghost_sizes, ghost_starts, 3);
    if (ierr) goto cleanup;
    ierr = MatSetDM(matrix, user->fda); if (ierr) goto cleanup;
    ierr = MatSetOption(matrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
    if (ierr) goto cleanup;
    ierr = PetscFree2(diagonal_nnz, offdiagonal_nnz); if (ierr) goto cleanup;
 
    *preconditioning_matrix = matrix;
    matrix = NULL;
 
cleanup:
    cleanup_ierr = MatDestroy(&matrix); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = PetscFree2(diagonal_nnz, offdiagonal_nnz);
    if (!ierr) ierr = cleanup_ierr;
    PetscFunctionReturn(ierr);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FrozenMomentumJacobian_MetricNormSquared()

static PetscReal FrozenMomentumJacobian_MetricNormSquared ( Cmpnts  metric )

static

Returns the squared Euclidean norm of one metric vector.

Definition at line 355 of file momentum_newton_krylov.c.

{
    return metric.x * metric.x + metric.y * metric.y + metric.z * metric.z;
}

Here is the caller graph for this function:

FrozenMomentumJacobian_PointBlock()

static void FrozenMomentumJacobian_PointBlock ( const SimCtx *  simCtx, const Cmpnts ***  ucont, const Cmpnts ***  csi, const Cmpnts ***  eta, const Cmpnts ***  zet, const PetscReal ***  aj, PetscInt  i, PetscInt  j, PetscInt  k, PetscScalar  block[9]  )

static

Returns the audited frozen-momentum point block in modern residual sign.

Definition at line 361 of file momentum_newton_krylov.c.

{
    const PetscReal dtc = MomentumBDFCoefficient((SimCtx *)simCtx) / simCtx->dt;
    const PetscReal inverse_reynolds = simCtx->ren > 0.0 ? 1.0 / simCtx->ren : 0.0;
    const PetscReal AJip = 0.5 * (aj[k][j][i] + aj[k][j][i + 1]);
    const PetscReal AJjp = 0.5 * (aj[k][j][i] + aj[k][j + 1][i]);
    const PetscReal AJkp = 0.5 * (aj[k][j][i] + aj[k + 1][j][i]);
    const PetscReal g11ip = csi[k][j][i].x * csi[k][j][i].x +
                             csi[k][j][i].y * csi[k][j][i].y +
                             csi[k][j][i].z * csi[k][j][i].z;
    const PetscReal g22ip = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j][i + 1]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j - 1][i]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j - 1][i + 1]));
    const PetscReal g33ip = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(zet[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k][j][i + 1]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k - 1][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k - 1][j][i + 1]));
    const PetscReal g11jp = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j + 1][i]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j][i - 1]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j + 1][i - 1]));
    const PetscReal g22jp = eta[k][j][i].x * eta[k][j][i].x +
                             eta[k][j][i].y * eta[k][j][i].y +
                             eta[k][j][i].z * eta[k][j][i].z;
    const PetscReal g33jp = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(zet[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k][j + 1][i]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k - 1][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(zet[k - 1][j + 1][i]));
    const PetscReal g11kp = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k + 1][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k][j][i - 1]) +
        FrozenMomentumJacobian_MetricNormSquared(csi[k + 1][j][i - 1]));
    const PetscReal g22kp = 0.25 * (
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k + 1][j][i]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k][j - 1][i]) +
        FrozenMomentumJacobian_MetricNormSquared(eta[k + 1][j - 1][i]));
    const PetscReal g33kp = zet[k][j][i].x * zet[k][j][i].x +
                             zet[k][j][i].y * zet[k][j][i].y +
                             zet[k][j][i].z * zet[k][j][i].z;
    const PetscReal U0jp = 0.25 * (ucont[k][j][i].x + ucont[k][j][i - 1].x +
                                   ucont[k][j + 1][i].x + ucont[k][j + 1][i - 1].x);
    const PetscReal U0kp = 0.25 * (ucont[k][j][i].x + ucont[k][j][i - 1].x +
                                   ucont[k + 1][j][i].x + ucont[k + 1][j][i - 1].x);
    const PetscReal U1ip = 0.25 * (ucont[k][j][i].y + ucont[k][j - 1][i].y +
                                   ucont[k][j][i + 1].y + ucont[k][j - 1][i + 1].y);
    const PetscReal U1kp = 0.25 * (ucont[k][j][i].y + ucont[k][j - 1][i].y +
                                   ucont[k + 1][j][i].y + ucont[k + 1][j - 1][i].y);
    const PetscReal U2ip = 0.25 * (ucont[k][j][i].z + ucont[k - 1][j][i].z +
                                   ucont[k][j][i + 1].z + ucont[k - 1][j][i + 1].z);
    const PetscReal U2jp = 0.25 * (ucont[k][j][i].z + ucont[k - 1][j][i].z +
                                   ucont[k][j + 1][i].z + ucont[k - 1][j + 1][i].z);
    PetscReal A[6][4] = {{0.0}};
    PetscReal Su, Sv, Sw, nui, nuj, nuk;
 
    A[0][0] =  0.125 * aj[k][j][i]         * ucont[k][j][i].y;
    A[0][1] = -0.125 * aj[k][j - 1][i]     * ucont[k][j - 1][i].y;
    A[0][2] =  0.125 * aj[k][j][i + 1]     * ucont[k][j][i + 1].y;
    A[0][3] = -0.125 * aj[k][j - 1][i + 1] * ucont[k][j - 1][i + 1].y;
    A[1][0] =  0.125 * aj[k][j][i]         * ucont[k][j][i].z;
    A[1][1] = -0.125 * aj[k - 1][j][i]     * ucont[k - 1][j][i].z;
    A[1][2] =  0.125 * aj[k][j][i + 1]     * ucont[k][j][i + 1].z;
    A[1][3] = -0.125 * aj[k - 1][j][i + 1] * ucont[k - 1][j][i + 1].z;
    A[2][0] = -0.125 * aj[k][j + 1][i - 1] * ucont[k][j + 1][i - 1].x;
    A[2][1] = -0.125 * aj[k][j][i - 1]     * ucont[k][j][i - 1].x;
    A[2][2] =  0.125 * aj[k][j + 1][i]     * ucont[k][j + 1][i].x;
    A[2][3] =  0.125 * aj[k][j][i]         * ucont[k][j][i].x;
    A[3][0] =  0.125 * aj[k][j][i]         * ucont[k][j][i].z;
    A[3][1] = -0.125 * aj[k - 1][j][i]     * ucont[k - 1][j][i].z;
    A[3][2] =  0.125 * aj[k][j + 1][i]     * ucont[k][j + 1][i].z;
    A[3][3] = -0.125 * aj[k - 1][j + 1][i] * ucont[k - 1][j + 1][i].z;
    A[4][0] = -0.125 * aj[k + 1][j][i - 1] * ucont[k + 1][j][i - 1].x;
    A[4][1] = -0.125 * aj[k][j][i - 1]     * ucont[k][j][i - 1].x;
    A[4][2] =  0.125 * aj[k + 1][j][i]     * ucont[k + 1][j][i].x;
    A[4][3] =  0.125 * aj[k][j][i]         * ucont[k][j][i].x;
    A[5][0] = -0.125 * aj[k + 1][j - 1][i] * ucont[k + 1][j - 1][i].y;
    A[5][1] = -0.125 * aj[k][j - 1][i]     * ucont[k][j - 1][i].y;
    A[5][2] =  0.125 * aj[k + 1][j][i]     * ucont[k + 1][j][i].y;
    A[5][3] =  0.125 * aj[k][j][i]         * ucont[k][j][i].y;
    Su = A[0][0] + A[0][1] + A[0][2] + A[0][3] +
         A[1][0] + A[1][1] + A[1][2] + A[1][3];
    Sv = A[2][0] + A[2][1] + A[2][2] + A[2][3] +
         A[3][0] + A[3][1] + A[3][2] + A[3][3];
    Sw = A[4][0] + A[4][1] + A[4][2] + A[4][3] +
         A[5][0] + A[5][1] + A[5][2] + A[5][3];
    nui = AJip * AJip * (g11ip + g22ip + g33ip) * inverse_reynolds;
    nuj = AJjp * AJjp * (g11jp + g22jp + g33jp) * inverse_reynolds;
    nuk = AJkp * AJkp * (g11kp + g22kp + g33kp) * inverse_reynolds;
 
    /* The modern residual is the negative of the legacy residual. */
    block[0] = dtc + nui + Su; block[1] = 0.5 * AJip * U1ip; block[2] = 0.5 * AJip * U2ip;
    block[3] = 0.5 * AJjp * U0jp; block[4] = dtc + nuj + Sv; block[5] = 0.5 * AJjp * U2jp;
    block[6] = 0.5 * AJkp * U0kp; block[7] = 0.5 * AJkp * U1kp; block[8] = dtc + nuk + Sw;
}

Here is the call graph for this function:

Here is the caller graph for this function:

FrozenMomentumJacobian_DescribePointBlock()

static PetscErrorCode FrozenMomentumJacobian_DescribePointBlock ( UserCtx *  user, MomentumPreconditionerDescription *  description  )

static

Describes the audited frozen-coefficient point-block model.

Definition at line 468 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    (void)user;
    PetscCheck(description != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Preconditioner description is NULL.");
    description->model = MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN;
    description->structure = MOM_NK_PC_STRUCTURE_POINT_BLOCK;
    description->block_size = 3;
    description->stencil_width = 0;
    PetscFunctionReturn(PETSC_SUCCESS);
}

MomentumPreconditionerEngine_ApplyConstraintRows()

static PetscErrorCode MomentumPreconditionerEngine_ApplyConstraintRows ( UserCtx *  user, Mat  preconditioning_matrix  )

static

Inserts all common fixed, homogeneous, and periodic-duplicate rows.

Definition at line 553 of file momentum_newton_krylov.c.

{
    DMDALocalInfo info = user->info;
 
    PetscFunctionBeginUser;
    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) {
                for (PetscInt component = 0; component < 3; ++component) {
                    PetscInt ri, rj, rk;
                    MomentumNewtonKrylovRowType type = MomentumNewtonKrylov_ClassifyRow(
                        user, i, j, k, component, &ri, &rj, &rk);
                    if (type != MOM_NK_ROW_PHYSICAL) {
                        MatStencil row = {.i = i, .j = j, .k = k, .c = component};
                        MatStencil columns[2] = {row, row};
                        PetscScalar values[2] = {1.0, -1.0};
                        PetscInt column_count = 1;
                        if (type == MOM_NK_ROW_PERIODIC_DUPLICATE) {
                            columns[column_count++] = (MatStencil){
                                .i = ri, .j = rj, .k = rk, .c = component
                            };
                        }
                        PetscCall(MatSetValuesStencil(preconditioning_matrix, 1, &row,
                                                      column_count, columns, values,
                                                      INSERT_VALUES));
                    }
                }
            }
        }
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumNewtonJacobian_Create()

static PetscErrorCode MomentumNewtonJacobian_Create ( SNES  snes, MomentumNewtonJacobian *  jacobian  )

static

Creates the selected Jacobian operator; currently PETSc MFFD only.

Definition at line 595 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    PetscCheck(jacobian->type == MOM_NK_JACOBIAN_FINITE_DIFFERENCE &&
               jacobian->finite_difference_mode == MOM_NK_FD_MODE_MATRIX_FREE,
               PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Unsupported Newton Krylov Jacobian type/finite-difference-mode combination.");
    PetscCall(MatCreateSNESMF(snes, &jacobian->jacobian_operator));
    PetscCall(MatSetOptionsPrefix(jacobian->jacobian_operator, "mom_nk_"));
    PetscCall(PetscObjectSetName((PetscObject)jacobian->jacobian_operator,
                                 "momentum_jacobian_finite_difference_matrix_free"));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumNewtonJacobian_Update()

static PetscErrorCode MomentumNewtonJacobian_Update ( SNES  snes, Vec  current_solution, MomentumNewtonJacobian *  jacobian  )

static

Updates the matrix-free finite-difference operator base.

Definition at line 611 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    PetscCall(MatMFFDComputeJacobian(snes, current_solution, jacobian->jacobian_operator,
                                     jacobian->jacobian_operator, NULL));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumNewtonJacobian_Register()

static PetscErrorCode MomentumNewtonJacobian_Register ( SNES  snes, MomentumNewtonJacobian *  jacobian, MomentumPreconditionerEngine *  engine, MomentumNewtonKrylovContext *  ctx  )

static

Registers the application orchestration callback and both SNES matrices.

Definition at line 621 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    PetscCall(SNESSetJacobian(snes, jacobian->jacobian_operator,
                              engine->preconditioning_matrix,
                              MomentumNewtonKrylov_FormJacobian, ctx));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumNewtonJacobian_Destroy()

static PetscErrorCode MomentumNewtonJacobian_Destroy ( MomentumNewtonJacobian *  jacobian )

static

Destroys a partially or fully created Jacobian operator.

Definition at line 633 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    PetscCall(MatDestroy(&jacobian->jacobian_operator));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumPreconditionerEngine_Create()

static PetscErrorCode MomentumPreconditionerEngine_Create ( UserCtx *  user, Mat  jacobian_operator, const MomentumPreconditionerDescription *  requested, MomentumPreconditionerEngine *  engine  )

static

Validates a model/structure and creates or aliases its matrix.

Definition at line 643 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    engine->description = *requested;
    if (requested->model == MOM_NK_PC_MODEL_NONE &&
        requested->structure == MOM_NK_PC_STRUCTURE_NONE) {
        engine->description.block_size = 0;
        engine->description.stencil_width = 0;
        engine->preconditioning_matrix = jacobian_operator;
        engine->aliases_jacobian_operator = PETSC_TRUE;
        engine->owns_preconditioning_matrix = PETSC_FALSE;
        engine->petsc_pc_type = PCNONE;
    } else if (requested->model == MOM_NK_PC_MODEL_FROZEN_MOMENTUM_JACOBIAN &&
               requested->structure == MOM_NK_PC_STRUCTURE_POINT_BLOCK) {
        engine->model_ops = &frozen_momentum_point_block_ops;
        PetscCall(engine->model_ops->Describe(user, &engine->description));
        PetscCall(MomentumPreconditionerEngine_CreateExactPointBlockMatrix(
            user, &engine->preconditioning_matrix));
        engine->owns_preconditioning_matrix = PETSC_TRUE;
        PetscCall(PetscObjectSetName((PetscObject)engine->preconditioning_matrix,
                                     "momentum_preconditioner_frozen_point_block"));
        engine->aliases_jacobian_operator = PETSC_FALSE;
        engine->petsc_pc_type = PCPBJACOBI;
    } else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,
                "Unsupported Newton Krylov preconditioner model/structure combination.");
    }
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumPreconditionerEngine_Assemble()

static PetscErrorCode MomentumPreconditionerEngine_Assemble ( MomentumPreconditionerEngine *  engine, UserCtx *  user, Vec  current_solution  )

static

Runs model insertion, common row handling, and final assembly.

Definition at line 678 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    if (engine->aliases_jacobian_operator) PetscFunctionReturn(PETSC_SUCCESS);
    PetscCall(MatZeroEntries(engine->preconditioning_matrix));
    PetscCall(engine->model_ops->AssembleInterior(user, current_solution,
                                                   engine->preconditioning_matrix));
    PetscCall(MomentumPreconditionerEngine_ApplyConstraintRows(
        user, engine->preconditioning_matrix));
    PetscCall(MatAssemblyBegin(engine->preconditioning_matrix, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(engine->preconditioning_matrix, MAT_FINAL_ASSEMBLY));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the call graph for this function:

Here is the caller graph for this function:

MomentumPreconditionerEngine_ConfigurePetscPC()

static PetscErrorCode MomentumPreconditionerEngine_ConfigurePetscPC ( MomentumPreconditionerEngine *  engine, PC  pc  )

static

Applies the validated model/structure-to-PETSc-PC mapping.

Definition at line 694 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    PetscCall(PCSetType(pc, engine->petsc_pc_type));
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumPreconditionerEngine_ValidatePetscPC()

static PetscErrorCode MomentumPreconditionerEngine_ValidatePetscPC ( MomentumPreconditionerEngine *  engine, PC  pc  )

static

Rejects raw options that select an unvalidated PETSc PC backend.

Definition at line 703 of file momentum_newton_krylov.c.

{
    const char *actual_type = NULL;
    PetscBool matches = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PetscCall(PCGetType(pc, &actual_type));
    PetscCall(PetscStrcmp(actual_type, engine->petsc_pc_type, &matches));
    PetscCheck(matches, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov preconditioner model/structure requires internal PETSc PC "
               "type '%s', but raw option processing selected '%s'.",
               engine->petsc_pc_type, actual_type ? actual_type : "(unset)");
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumPreconditionerEngine_Destroy()

static PetscErrorCode MomentumPreconditionerEngine_Destroy ( MomentumPreconditionerEngine *  engine )

static

Destroys only a separately owned preconditioning matrix.

Definition at line 720 of file momentum_newton_krylov.c.

{
    PetscFunctionBeginUser;
    if (!engine->owns_preconditioning_matrix) engine->preconditioning_matrix = NULL;
    PetscCall(MatDestroy(&engine->preconditioning_matrix));
    engine->aliases_jacobian_operator = PETSC_FALSE;
    engine->owns_preconditioning_matrix = PETSC_FALSE;
    PetscFunctionReturn(PETSC_SUCCESS);
}

Here is the caller graph for this function:

MomentumSolver_NewtonKrylov()

PetscErrorCode MomentumSolver_NewtonKrylov	(	UserCtx *	user,
		IBMNodes *	ibm,
		FSInfo *	fsi
	)

Solves one physical momentum step with matrix-free Newton–Krylov.

The Jacobian operator is the finite-difference, matrix-free operator of the complete deterministic residual. With no mathematical preconditioner its matrix is also passed to PETSc's preconditioning slot and PCNONE is derived. The optional provisional frozen-momentum, point-block model approximates same-cell frozen coupling, owns a separate matrix, and derives a block-Jacobi PETSc backend. It does not relax the version-one physics or boundary-condition restrictions. All PETSc solver objects are local to this call. Rows removed by legacy boundary residual enforcement are made explicit: conditioned normal rows use X-Uconditioned, untouched dummy/tangential rows use X, and periodic duplicates use Xdup-Xrepresentative. Unsupported masked, interface, and component-disabled rows are rejected before setup.

Parameters

user	Single-block momentum context.
ibm	Must be NULL; immersed boundaries are not supported in version one.
fsi	Must be NULL; moving-body coupling is not supported in version one.

Returns

0 on convergence, PETSC_ERR_CONV_FAILED after rollback on nonconvergence.

Definition at line 1131 of file momentum_newton_krylov.c.

{
    PetscErrorCode ierr = PETSC_SUCCESS, cleanup_ierr;
    SimCtx        *simCtx;
    SNES           snes = NULL;
    Vec            solution = NULL, entry_backup = NULL;
    KSP            ksp = NULL;
    PC             pc = NULL;
    MomentumPreconditionerDescription preconditioner_description = {0};
    PetscBool      restore_entry = PETSC_FALSE;
    PetscBool      rhs_created = PETSC_FALSE;
    PetscBool      solve_started = PETSC_FALSE;
    PetscBool      committed = PETSC_FALSE;
    SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
    PetscInt       nonlinear_its = 0, function_evals = 0, linear_its = 0;
    PetscReal      final_norm = PETSC_MAX_REAL;
    MomentumNewtonKrylovContext ctx = {0};
    const FieldId staggered_fields[] = {FIELD_ID_UCONT};
 
    PetscFunctionBeginUser;
    PetscCall(MomentumNewtonKrylov_Validate(user));
    PetscCheck(ibm == NULL && fsi == NULL, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one does not accept IBM or FSI objects.");
    PetscCheck(user->Rhs == NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "Newton Krylov requires UserCtx::Rhs to be unallocated on entry.");
    simCtx = user->simCtx;
 
    ierr = VecDuplicate(user->Ucont, &solution); if (ierr) goto cleanup;
    ierr = VecDuplicate(user->Ucont, &entry_backup); if (ierr) goto cleanup;
    ierr = VecDuplicate(user->Ucont, &user->Rhs); if (ierr) goto cleanup;
    rhs_created = PETSC_TRUE;
    ierr = SNESCreate(PetscObjectComm((PetscObject)user->Ucont), &snes); if (ierr) goto cleanup;
 
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); if (ierr) goto cleanup;
    ierr = ApplyBoundaryConditions(user); if (ierr) goto cleanup;
    ierr = VecCopy(user->Ucont, entry_backup); if (ierr) goto cleanup;
    restore_entry = PETSC_TRUE;
    ierr = VecCopy(user->Ucont, solution); if (ierr) goto cleanup;
 
    ctx.user = user;
    ierr = MomentumNewtonKrylov_ReadLinearizationConfig(
        &ctx.jacobian, &preconditioner_description); if (ierr) goto cleanup;
    ierr = SNESSetOptionsPrefix(snes, "mom_nk_"); if (ierr) goto cleanup;
    ierr = SNESSetType(snes, SNESNEWTONLS); if (ierr) goto cleanup;
    ierr = SNESSetDM(snes, user->fda); if (ierr) goto cleanup;
    ierr = SNESSetFunction(snes, NULL, MomentumNewtonKrylov_FormResidual, &ctx); if (ierr) goto cleanup;
    ierr = MomentumNewtonJacobian_Create(snes, &ctx.jacobian); if (ierr) goto cleanup;
    ierr = MomentumPreconditionerEngine_Create(user, ctx.jacobian.jacobian_operator,
                                                &preconditioner_description,
                                                &ctx.preconditioning_engine); if (ierr) goto cleanup;
    ierr = MomentumNewtonJacobian_Register(snes, &ctx.jacobian,
                                            &ctx.preconditioning_engine, &ctx); if (ierr) goto cleanup;
    ierr = SNESGetKSP(snes, &ksp); if (ierr) goto cleanup;
    ierr = KSPSetType(ksp, KSPGMRES); if (ierr) goto cleanup;
    ierr = KSPGetPC(ksp, &pc); if (ierr) goto cleanup;
    ierr = MomentumPreconditionerEngine_ConfigurePetscPC(&ctx.preconditioning_engine, pc); if (ierr) goto cleanup;
    ierr = SNESSetFromOptions(snes); if (ierr) goto cleanup;
    ierr = SNESMonitorSet(snes, MomentumNewtonKrylov_Monitor, &ctx, NULL); if (ierr) goto cleanup;
    ierr = MomentumPreconditionerEngine_ValidatePetscPC(&ctx.preconditioning_engine, pc); if (ierr) goto cleanup;
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Newton Krylov Jacobian: finite_difference / matrix_free; "
              "Preconditioner: %s; PETSc Jacobian matrix type: MATMFFD; "
              "PETSc PC type: %s.\n",
              ctx.preconditioning_engine.description.model == MOM_NK_PC_MODEL_NONE ?
                  "none" : "frozen_momentum_jacobian / point_block",
              ctx.preconditioning_engine.petsc_pc_type);
 
    MomentumNewtonKrylov_OpenHistory(&ctx);
    solve_started = PETSC_TRUE;
    ierr = SNESSolve(snes, NULL, solution);
    if (ierr) goto cleanup;
    ierr = SNESGetConvergedReason(snes, &reason); if (ierr) goto cleanup;
    ierr = SNESGetIterationNumber(snes, &nonlinear_its); if (ierr) goto cleanup;
    ierr = SNESGetNumberFunctionEvals(snes, &function_evals); if (ierr) goto cleanup;
    ierr = SNESGetLinearSolveIterations(snes, &linear_its); if (ierr) goto cleanup;
    ierr = SNESGetFunctionNorm(snes, &final_norm); if (ierr) goto cleanup;
 
    if (reason > 0) {
        ierr = VecCopy(solution, user->Ucont); if (ierr) goto cleanup;
        simCtx->mom_last_converged = PETSC_TRUE;
    } else {
        ierr = VecCopy(entry_backup, user->Ucont); if (ierr) goto cleanup;
        simCtx->mom_last_converged = PETSC_FALSE;
    }
    ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); if (ierr) goto cleanup;
    ierr = ApplyBoundaryConditions(user); if (ierr) goto cleanup;
    restore_entry = PETSC_FALSE;
    committed = (PetscBool)(reason > 0);
 
    LOG_ALLOW(GLOBAL, LOG_INFO,
              "Newton Krylov momentum solve: reason=%s (%d), Newton iterations=%d, residual evaluations=%d, Krylov iterations=%d, final norm=%.6e, state=%s.\n",
              SNESConvergedReasons[reason], (PetscInt)reason, nonlinear_its, function_evals,
              linear_its, (double)final_norm, reason > 0 ? "committed" : "rolled back");
    if (reason <= 0) ierr = PETSC_ERR_CONV_FAILED;
 
cleanup:
    /* A PETSc solve error can bypass the normal statistics path. Query whatever
       SNES retained without replacing the primary error so a failed attempt is
       still represented in the structured log. */
    if (solve_started && snes) {
        (void)SNESGetConvergedReason(snes, &reason);
        (void)SNESGetIterationNumber(snes, &nonlinear_its);
        (void)SNESGetNumberFunctionEvals(snes, &function_evals);
        (void)SNESGetLinearSolveIterations(snes, &linear_its);
        (void)SNESGetFunctionNorm(snes, &final_norm);
    }
    if (restore_entry && entry_backup) {
        cleanup_ierr = VecCopy(entry_backup, user->Ucont);
        if (!ierr) ierr = cleanup_ierr;
        cleanup_ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields);
        if (!ierr) ierr = cleanup_ierr;
        cleanup_ierr = ApplyBoundaryConditions(user);
        if (!ierr) ierr = cleanup_ierr;
        simCtx->mom_last_converged = PETSC_FALSE;
        committed = PETSC_FALSE;
    }
    if (ctx.history_file) {
        (void)fclose(ctx.history_file);
        ctx.history_file = NULL;
    }
    if (solve_started) {
        MomentumNewtonKrylov_WriteSummary(&ctx, reason, nonlinear_its, function_evals,
                                          linear_its, final_norm, committed);
    }
    if (rhs_created) {
        cleanup_ierr = VecDestroy(&user->Rhs);
        if (!ierr) ierr = cleanup_ierr;
    }
    cleanup_ierr = VecDestroy(&entry_backup); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = VecDestroy(&solution); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = MomentumPreconditionerEngine_Destroy(&ctx.preconditioning_engine); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = MomentumNewtonJacobian_Destroy(&ctx.jacobian); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = SNESDestroy(&snes); if (!ierr) ierr = cleanup_ierr;
    PetscFunctionReturn(ierr);
}

Here is the call graph for this function:

Variable Documentation

frozen_momentum_point_block_ops

const MomentumPreconditionerModelOps frozen_momentum_point_block_ops

static

Initial value:

= {
    FrozenMomentumJacobian_DescribePointBlock,
    FrozenMomentumJacobian_AssemblePointBlocks
}

Definition at line 587 of file momentum_newton_krylov.c.

                                                                              {
    FrozenMomentumJacobian_DescribePointBlock,
    FrozenMomentumJacobian_AssemblePointBlocks
};