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	MomentumNewtonKrylovContext

Macros
#define	__FUNCT__   "MomentumNewtonKrylov_Validate"
#define	__FUNCT__   "MomentumNewtonKrylov_ClassifyRow"
#define	__FUNCT__   "MomentumNewtonKrylov_ApplyConstraints"
#define	__FUNCT__   "MomentumNewtonKrylov_FormResidual"
#define	__FUNCT__   "MomentumSolver_NewtonKrylov"

Enumerations
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_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.
PetscErrorCode	MomentumSolver_NewtonKrylov (UserCtx *user, IBMNodes *ibm, FSInfo *fsi)
	Runs one per-call matrix-free Newton–Krylov momentum solve.

---

Data Structure Documentation

MomentumNewtonKrylovContext

struct MomentumNewtonKrylovContext

Definition at line 3 of file momentum_newton_krylov.c.

Collaboration diagram for MomentumNewtonKrylovContext:

Data Fields
UserCtx *	user

Macro Definition Documentation

__FUNCT__ [1/5]

#define __FUNCT__   "MomentumNewtonKrylov_Validate"

Definition at line 23 of file momentum_newton_krylov.c.

__FUNCT__ [2/5]

#define __FUNCT__   "MomentumNewtonKrylov_ClassifyRow"

Definition at line 23 of file momentum_newton_krylov.c.

__FUNCT__ [3/5]

#define __FUNCT__   "MomentumNewtonKrylov_ApplyConstraints"

Definition at line 23 of file momentum_newton_krylov.c.

__FUNCT__ [4/5]

#define __FUNCT__   "MomentumNewtonKrylov_FormResidual"

Definition at line 23 of file momentum_newton_krylov.c.

__FUNCT__ [5/5]

#define __FUNCT__   "MomentumSolver_NewtonKrylov"

Definition at line 23 of file momentum_newton_krylov.c.

Enumeration Type Documentation

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 7 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 29 of file momentum_newton_krylov.c.

{
    SimCtx   *simCtx;
    PetscReal mask_max = 0.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.");
    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.");
    PetscCheck(simCtx->StartStep == 0, PETSC_COMM_WORLD, PETSC_ERR_SUP,
               "Newton Krylov version one supports fresh starts only (StartStep must be zero)." );
 
    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.

Supported handlers may overwrite flux totals and other diagnostics on every call, but repeatability tests require those values not to affect a later call at the same X. No histories, pressure, viscosity, or controller state advance is performed here.

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 259 of file momentum_newton_krylov.c.

{
    MomentumNewtonKrylovContext *ctx = (MomentumNewtonKrylovContext *)vctx;
    UserCtx *user = ctx->user;
    const char *staggered_fields[] = {"Ucont"};
 
    PetscFunctionBeginUser;
    (void)snes;
    PetscCall(VecCopy(X, user->Ucont));
    PetscCall(SynchronizePeriodicStaggeredFields(user, 1, staggered_fields));
    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_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 195 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 129 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:

MomentumSolver_NewtonKrylov()

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

Runs one per-call matrix-free Newton–Krylov momentum solve.

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

Parameters

user	Single-block momentum context.
ibm	Must be NULL in version one.
fsi	Must be NULL in version one.

Returns

0 on convergence or an error after restoring the canonical entry state.

Definition at line 286 of file momentum_newton_krylov.c.

{
    PetscErrorCode ierr = PETSC_SUCCESS, cleanup_ierr;
    SimCtx        *simCtx;
    SNES           snes = NULL;
    Mat            J = NULL;
    Vec            solution = NULL, entry_backup = NULL;
    KSP            ksp = NULL;
    PC             pc = NULL;
    const char    *pc_type = NULL;
    PetscBool      pc_is_none = PETSC_FALSE;
    PetscBool      restore_entry = PETSC_FALSE;
    PetscBool      rhs_created = 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;
    const char *staggered_fields[] = {"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 = 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 = MatCreateSNESMF(snes, &J); if (ierr) goto cleanup;
    ierr = SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, NULL); 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 = PCSetType(pc, PCNONE); if (ierr) goto cleanup;
    ierr = SNESSetFromOptions(snes); if (ierr) goto cleanup;
    ierr = PCGetType(pc, &pc_type); if (ierr) goto cleanup;
    ierr = PetscStrcmp(pc_type, PCNONE, &pc_is_none); if (ierr) goto cleanup;
    if (!pc_is_none) {
        LOG(GLOBAL, LOG_ERROR,
            "Newton Krylov version one requires PCNONE; option processing selected PC type '%s'.\n",
            pc_type ? pc_type : "(unset)");
        ierr = PETSC_ERR_SUP;
        goto cleanup;
    }
 
    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;
 
    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:
    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;
    }
    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 = MatDestroy(&J); if (!ierr) ierr = cleanup_ierr;
    cleanup_ierr = SNESDestroy(&snes); if (!ierr) ierr = cleanup_ierr;
    PetscFunctionReturn(ierr);
}

Here is the call graph for this function: