PETSc

https://petsc.org/release/
Admin
PETSc (Portable, Extensible Toolkit for Scientific Computation) is a software library for the...
Tokens:10,888,198
Snippets:76,388
Trust Score:8.8
Update:2 weeks ago
Show doc for...
Tokens:
Raw
### Setup: help, includes, and material parameters

Source: https://petsc.org/release/src/snes/tutorials/ex56.c.html

Initial setup for the PETSc example: help string, required includes, and material parameters for the elasticity problem.

```c
static char help[] = "3D tensor hexahedra & 3D Laplacian displacement finite element formulation\n\
of linear elasticity.  E=1.0, nu=1/3.\n\
Unit cube domain with Dirichlet boundary\n\n";
```

```c
#include <petscdmplex.h>
#include <petscsnes.h>
#include <petscds.h>
#include <petscdmforest.h>
```

```c
static PetscReal s_soft_alpha = 0.01;
static PetscReal s_mu         = 0.4;
static PetscReal s_lambda     = 0.4;
```

--------------------------------

### Stokes MMS setup: ProcessOptions, CreateMesh, SetupParameters

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Setup routines for the Stokes MMS example. ProcessOptions registers the -sol and -bc command-line options using PetscOptionsEList, with defaults SOL_QUADRATIC and BC_ESSENTIAL. CreateMesh creates a DMPLEX DM and optionally views it with -dm_view. SetupParameters creates a PETSc parameter bag registering mu (default 1.0) and eta (default 100.0) as scalars, then applies command-line overrides.

```C
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  PetscInt sol, bc;

  PetscFunctionBeginUser;
  options->sol = SOL_QUADRATIC;
  options->bc  = BC_ESSENTIAL;
  PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
  sol = options->sol;
  PetscCall(PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, PETSC_STATIC_ARRAY_LENGTH(SolTypes) - 3, SolTypes[options->sol], &sol, NULL));
  options->sol = (SolType)sol;
  bc           = options->bc;
  PetscCall(PetscOptionsEList("-bc", "The boundary condition type", "ex62.c", BCTypes, PETSC_STATIC_ARRAY_LENGTH(BCTypes) - 3, BCTypes[options->bc], &bc, NULL));
  options->bc = (BCType)bc;
  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}
```

```C
static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
{
  PetscFunctionBeginUser;
  PetscCall(DMCreate(comm, dm));
  PetscCall(DMSetType(*dm, DMPLEX));
  PetscCall(DMSetFromOptions(*dm));
  PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
  PetscFunctionReturn(PETSC_SUCCESS);
}
```

```C
static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
{
  Parameter *p;

  PetscFunctionBeginUser;
  /* setup PETSc parameter bag */
  PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &ctx->bag));
  PetscCall(PetscBagGetData(ctx->bag, &p));
  PetscCall(PetscBagSetName(ctx->bag, "par", "Stokes Parameters"));
  PetscCall(PetscBagRegisterScalar(ctx->bag, &p->mu, 1.0, "mu", "Dynamic Shear Viscosity, Pa s"));
  PetscCall(PetscBagRegisterScalar(ctx->bag, &p->eta, 100.0, "eta", "Nitsche penalty parameter (dimensionless)"));
  PetscCall(PetscBagSetFromOptions(ctx->bag));
  {
    PetscViewer       viewer;
    PetscViewerFormat format;
    PetscBool         flg;

```

--------------------------------

### Run ex69 with -snes_view to inspect solver setup

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Runs the P2-P1 Stokes example with -snes_view to inspect the solver configuration. The output shows a Newton line search (type bt) with GMRES KSP and a fieldsplit PC using a Schur preconditioner with FULL factorization, with the Schur complement formed from A11.

```shell
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_view"
SNES Object: 1 MPI process
  type: newtonls
  maximum iterations=50, maximum function evaluations=10000
  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
  total number of linear solver iterations=1
  total number of function evaluations=2
  norm schedule ALWAYS
  SNESLineSearch Object: 1 MPI process
    type: bt
      interpolation: cubic
      alpha=1.000000e-04
    maxstep=1.000000e+08, minlambda=1.000000e-12
    tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
    maximum iterations=40
  KSP Object: 1 MPI process
    type: gmres
      restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
      happy breakdown tolerance=1e-30
    maximum iterations=10000, initial guess is zero
    tolerances: relative=1e-09, absolute=1e-10, divergence=10000.
    left preconditioning
    using PRECONDITIONED norm type for convergence test
  PC Object: 1 MPI process
    type: fieldsplit
      FieldSplit with Schur preconditioner, factorization FULL
      Preconditioner for the Schur complement formed from A11
      Split info:
      Split number 0 Defined by IS
      Split number 1 Defined by IS
      KSP solver for A00 block
        KSP Object: (fieldsplit_velocity_) 1 MPI process
          type: gmres
            restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
            happy breakdown tolerance=1e-30
          maximum iterations=10000, initial guess is zero
          tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
          left preconditioning
```

--------------------------------

### Monitor solver convergence with -snes_monitor and -ksp_monitor_true_residual

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Runs the example with monitoring options to check convergence. The output shows both the linear and nonlinear solvers converge in one iteration, with the inner Schur complement solver converging in seven iterations each time.

```shell
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_monitor -ksp_monitor_true_residual -ksp_converged_reason -fieldsplit_pressure_ksp_converged_reason"
0 SNES Function norm 1.170604545948e-01
  Linear fieldsplit_pressure_ solve converged due to CONVERGED_RTOL iterations 7
  0 KSP preconditioned resid norm 4.965098891419e-01 true resid norm 1.170604545948e-01 ||r(i)||/||b|| 1.000000000000e+00
  Linear fieldsplit_pressure_ solve converged due to CONVERGED_RTOL iterations 7
  1 KSP preconditioned resid norm 9.236813926190e-11 true resid norm 1.460072673561e-12 ||r(i)||/||b|| 1.247280884579e-11
Linear solve converged due to CONVERGED_ATOL iterations 1
1 SNES Function norm 1.460070661322e-12

```

--------------------------------

### Main function declaration

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Starts the main function for the MMS example, declaring the SNES, DM, Vec, and AppCtx variables and a flag for testing patch facet residuals.

```c
int main(int argc, char **argv)
{
  SNES      snes;
  DM        dm;
  Vec       u;
  AppCtx    user;
  PetscBool testPatchFacetResidual = PETSC_FALSE;

```

--------------------------------

### ProcessOptions function start

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Starts the ProcessOptions function that reads command-line options into the AppCtx structure. The body is incomplete in the source; it declares local variables sol and bc before the source cuts off.

```C
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  PetscInt sol, bc;

```

--------------------------------

### View KSP and PC configuration with -snes_view

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Run the PETSc example with -snes_view to display the KSP and PC configuration. The output shows GMRES with restart=30 and a field split preconditioner using a Schur complement with FULL factorization.

```bash
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_view"
  KSP Object: 1 MPI process
    type: gmres
      restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
      happy breakdown tolerance=1e-30
    maximum iterations=10000, initial guess is zero
    tolerances: relative=1e-09, absolute=1e-10, divergence=10000.
    left preconditioning
    using PRECONDITIONED norm type for convergence test
  PC Object: 1 MPI process
    type: fieldsplit
      FieldSplit with Schur preconditioner, factorization FULL
      Preconditioner for the Schur complement formed from A11
      Split info:
      Split number 0 Defined by IS
      Split number 1 Defined by IS

```

--------------------------------

### Main driver for MMS Stokes example

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Main driver for the MMS Stokes example. Initializes PETSc, optionally runs a patch facet residual test, otherwise sets up the mesh, SNES, and user context, creates a pressure null space, solves, and cleans up.

```C
PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
  if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
  else {
    PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
    PetscCall(SNESSetDM(snes, dm));
    PetscCall(DMSetApplicationContext(dm, &user));

    PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    PetscCall(SetupProblem(dm, SetupEqn, &user));
    PetscCall(DMPlexCreateClosureIndex(dm, NULL));

    PetscCall(DMCreateGlobalVector(dm, &u));
    PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    PetscCall(SNESSetFromOptions(snes));
    PetscCall(DMSNESCheckFromOptions(snes, u));
    PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
    {
      Mat          J;
      MatNullSpace sp;

      PetscCall(SNESSetUp(snes));
      PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
      PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
      PetscCall(MatSetNullSpace(J, sp));
      PetscCall(MatNullSpaceDestroy(&sp));
      PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
      PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
    }
    PetscCall(SNESSolve(snes, NULL, u));

    PetscCall(VecDestroy(&u));
    PetscCall(SNESDestroy(&snes));
    PetscCall(DMDestroy(&dm));
    PetscCall(PetscBagDestroy(&user.bag));
  }
  PetscCall(PetscFinalize());
  return 0;
```

--------------------------------

### View system and preconditioner matrices with -snes_view

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Runs the PETSc example snes_tutorials-ex69_p2p1 with -snes_view to display the linear system matrix and the preconditioner matrix. The output shows both matrices have the same nonzero pattern (653 nonzeros) because zeros are inserted to define the structure.

```shell
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_view"
    linear system matrix followed by matrix used to construct the preconditioner:
    Mat Object: 1 MPI process
      type: seqaij
      rows=39, cols=39
      total: nonzeros=653, allocated nonzeros=653
      total number of mallocs used during MatSetValues calls=0
        has attached null space
        using I-node routines: found 24 nodes, limit used is 5
    Mat Object: (prec_) 1 MPI process
      type: seqaij
      rows=39, cols=39
      total: nonzeros=653, allocated nonzeros=653
      total number of mallocs used during MatSetValues calls=0
        using I-node routines: found 24 nodes, limit used is 5

```

--------------------------------

### Main driver for Stokes MMS example

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Main driver for the Stokes MMS example. It initializes PETSc, optionally runs a patch facet residual test, otherwise sets up the mesh, SNES, and user context, configures the FEM problem, solves, and cleans up. The code sets a pressure null space on the Jacobian and names the solution and Jacobian for viewing.

```c
  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
  if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
  else {
    PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
    PetscCall(SNESSetDM(snes, dm));
    PetscCall(DMSetApplicationContext(dm, &user));

    PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    PetscCall(SetupProblem(dm, SetupEqn, &user));
    PetscCall(DMPlexCreateClosureIndex(dm, NULL));

    PetscCall(DMCreateGlobalVector(dm, &u));
    PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    PetscCall(SNESSetFromOptions(snes));
    PetscCall(DMSNESCheckFromOptions(snes, u));
    PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
    {
      Mat          J;
      MatNullSpace sp;

      PetscCall(SNESSetUp(snes));
      PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
      PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
      PetscCall(MatSetNullSpace(J, sp));
      PetscCall(MatNullSpaceDestroy(&sp));
      PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
      PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
    }
    PetscCall(SNESSolve(snes, NULL, u));

    PetscCall(VecDestroy(&u));
    PetscCall(SNESDestroy(&snes));
    PetscCall(DMDestroy(&dm));
    PetscCall(PetscBagDestroy(&user.bag));
  }
  PetscCall(PetscFinalize());
  return 0;
```

--------------------------------

### ProcessOptions: register -sol and -bc options

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Processes command-line options for the Stokes example. Sets defaults SOL_QUADRATIC and BC_ESSENTIAL, then registers the -sol and -bc options via PetscOptionsEList with the SolTypes and BCTypes arrays.

```C
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  PetscInt sol, bc;

  PetscFunctionBeginUser;
  options->sol = SOL_QUADRATIC;
  options->bc  = BC_ESSENTIAL;
  PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
  sol = options->sol;
  PetscCall(PetscOptionsEList("-sol", "The MMS solution", "ex62.c", SolTypes, PETSC_STATIC_ARRAY_LENGTH(SolTypes) - 3, SolTypes[options->sol], &sol, NULL));
  options->sol = (SolType)sol;
  bc           = options->bc;
  PetscCall(PetscOptionsEList("-bc", "The boundary condition type", "ex62.c", BCTypes, PETSC_STATIC_ARRAY_LENGTH(BCTypes) - 3, BCTypes[options->bc], &bc, NULL));
  options->bc = (BCType)bc;
  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}
```

--------------------------------

### Main driver for MMS Stokes example

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Main driver for the MMS Stokes example. It initializes PETSc, optionally runs a patch facet residual test, otherwise sets up the mesh, SNES, and FEM problem, solves, and cleans up. The code sets a pressure null space on the Jacobian and names the solution and Jacobian for viewing.

```c
  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
  if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
  else {
    PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
    PetscCall(SNESSetDM(snes, dm));
    PetscCall(DMSetApplicationContext(dm, &user));

    PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    PetscCall(SetupProblem(dm, SetupEqn, &user));
    PetscCall(DMPlexCreateClosureIndex(dm, NULL));

    PetscCall(DMCreateGlobalVector(dm, &u));
    PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    PetscCall(SNESSetFromOptions(snes));
    PetscCall(DMSNESCheckFromOptions(snes, u));
    PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
    {
      Mat          J;
      MatNullSpace sp;

      PetscCall(SNESSetUp(snes));
      PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
      PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
      PetscCall(MatSetNullSpace(J, sp));
      PetscCall(MatNullSpaceDestroy(&sp));
      PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
      PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
    }
    PetscCall(SNESSolve(snes, NULL, u));

    PetscCall(VecDestroy(&u));
    PetscCall(SNESDestroy(&snes));
    PetscCall(DMDestroy(&dm));
    PetscCall(PetscBagDestroy(&user.bag));
  }
  PetscCall(PetscFinalize());
  return 0;
}
```

--------------------------------

### Run PETSc ex69 with -snes_view to inspect fieldsplit solver

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Runs the PETSc example snes_tutorials-ex69_p2p1 with -snes_view to inspect the fieldsplit preconditioner. The output shows the Schur complement S = A11 - A10 inv(A00) A01, with the pressure block using GMRES/LU and the velocity block using GMRES/LU.

```shell
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_view"
      KSP solver for S = A11 - A10 inv(A00) A01
        KSP Object: (fieldsplit_pressure_) 1 MPI process
          type: gmres
            restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
            happy breakdown tolerance=1e-30
          maximum iterations=10000, initial guess is zero
          tolerances: relative=1e-09, absolute=1e-50, divergence=10000.
          left preconditioning
          using PRECONDITIONED norm type for convergence test
        PC Object: (fieldsplit_pressure_) 1 MPI process
          type: lu
            out-of-place factorization
            tolerance for zero pivot 2.22045e-14
            matrix ordering: nd
            factor fill ratio given 5., needed 1.2439
              Factored matrix:
                Mat Object: 1 MPI process
                  type: seqaij
                  rows=9, cols=9
                  package used to perform factorization: petsc
                  total: nonzeros=51, allocated nonzeros=51
                    not using I-node routines
          linear system matrix followed by matrix used to construct the preconditioner:
          Mat Object: (fieldsplit_pressure_) 1 MPI process
            type: schurcomplement
            rows=9, cols=9
              has attached null space
              Schur complement A11 - A10 inv(A00) A01
              A11
                Mat Object: 1 MPI process
                  type: seqaij
                  rows=9, cols=9
                  total: nonzeros=41, allocated nonzeros=41
                  total number of mallocs used during MatSetValues calls=0
                    has attached null space
                    not using I-node routines
              A10
                Mat Object: 1 MPI process
                  type: seqaij
                  rows=9, cols=30
                  total: nonzeros=122, allocated nonzeros=122
                  total number of mallocs used during MatSetValues calls=0
                    not using I-node routines
              KSP of A00
                KSP Object: (fieldsplit_velocity_) 1 MPI process
                  type: gmres
                    restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
                    happy breakdown tolerance=1e-30
                  maximum iterations=10000, initial guess is zero
                  tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
                  left preconditioning
                  using PRECONDITIONED norm type for convergence test
                PC Object: (fieldsplit_velocity_) 1 MPI process
                  type: lu
                    out-of-place factorization
                    tolerance for zero pivot 2.22045e-14
                    matrix ordering: nd
                    factor fill ratio given 5., needed 1.15761
                      Factored matrix:
                        Mat Object: 1 MPI process
                          type: seqaij
                          rows=30, cols=30
                          package used to perform factorization: petsc
                          total: nonzeros=426, allocated nonzeros=426
                            using I-node routines: found 17 nodes, limit used is 5
                  linear system matrix followed by matrix used to construct the preconditioner:
                  Mat Object: 1 MPI process
                    type: seqaij
                    rows=30, cols=30
                    total: nonzeros=368, allocated nonzeros=368
                    total number of mallocs used during MatSetValues calls=0
                      using I-node routines: found 20 nodes, limit used is 5
                  Mat Object: (fieldsplit_velocity_) 1 MPI process
                    type: seqaij
                    rows=30, cols=30
                    total: nonzeros=368, allocated nonzeros=368
                    total number of mallocs used during MatSetValues calls=0
                      using I-node routines: found 20 nodes, limit used is 5
              A01
                Mat Object: 1 MPI process
                  type: seqaij
                  rows=30, cols=9
                  total: nonzeros=122, allocated nonzeros=122
                  total number of mallocs used during MatSetValues calls=0
                    using I-node routines: found 20 nodes, limit used is 5
          Mat Object: (fieldsplit_pressure_) 1 MPI process
            type: seqaij
            rows=9, cols=9
            total: nonzeros=41, allocated nonzeros=41
            total number of mallocs used during MatSetValues calls=0
              not using I-node routines

```

--------------------------------

### Run PETSc ex69 with -snes_view to inspect velocity block solver

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Run the PETSc example with -snes_view to inspect the velocity block solver. The output shows the GMRES/LU solver for the fieldsplit_velocity_ block, including the KSP and PC objects and the two identical matrices (original and preconditioner).

```text
$ make -f ./gmakefile test search="snes_tutorials-ex69_p2p1" EXTRA_OPTIONS="-snes_view"
      KSP solver for A00 block
        KSP Object: (fieldsplit_velocity_) 1 MPI process
          type: gmres
            restart=30, using classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement
            happy breakdown tolerance=1e-30
          maximum iterations=10000, initial guess is zero
          tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
          left preconditioning
          using PRECONDITIONED norm type for convergence test
        PC Object: (fieldsplit_velocity_) 1 MPI process
          type: lu
            out-of-place factorization
            tolerance for zero pivot 2.22045e-14
            matrix ordering: nd
            factor fill ratio given 5., needed 1.15761
              Factored matrix:
                Mat Object: 1 MPI process
                  type: seqaij
                  rows=30, cols=30
                  package used to perform factorization: petsc
                  total: nonzeros=426, allocated nonzeros=426
                    using I-node routines: found 17 nodes, limit used is 5
          linear system matrix followed by matrix used to construct the preconditioner:
          Mat Object: 1 MPI process
            type: seqaij
            rows=30, cols=30
            total: nonzeros=368, allocated nonzeros=368
            total number of mallocs used during MatSetValues calls=0
              using I-node routines: found 20 nodes, limit used is 5
          Mat Object: (fieldsplit_velocity_) 1 MPI process
            type: seqaij
            rows=30, cols=30
            total: nonzeros=368, allocated nonzeros=368
            total number of mallocs used during MatSetValues calls=0
              using I-node routines: found 20 nodes, limit used is 5

```

--------------------------------

### Configure p2p1_gmg test with -dm_refine_hierarchy and -fieldsplit_velocity_pc_type mg

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Test configuration for the p2p1_gmg example. Uses -dm_refine_hierarchy 2 to create a geometric multigrid hierarchy and -fieldsplit_velocity_pc_type mg for the velocity block. Marked as broken (requires subDMs hooks) and requires the triangle package.

```yaml
    suffix: p2p1_gmg
    TODO: broken (requires subDMs hooks)
    requires: triangle
    args: -dm_plex_separate_marker -dm_refine_hierarchy 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
      -snes_error_if_not_converged -dmsnes_check .001 \
      -ksp_type fgmres -ksp_rtol 1.e-9 -ksp_error_if_not_converged -pc_use_amat \
      -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition a11 \
        -fieldsplit_velocity_pc_type mg \
        -fieldsplit_pressure_ksp_rtol 1.e-9 -fieldsplit_pressure_pc_type lu
```

--------------------------------

### Setup PETSc DMPlex Stokes problem with pressure null space

Source: https://petsc.org/release/tutorials/physics/guide_to_stokes

Sets up a PETSc DMPlex finite element problem with a pressure null space, creates the SNES solver, and solves the Stokes system. The code attaches a null space to the pressure field and to the Jacobian, and supports a -test_patch_facet_residual option.

```c
PetscCall(DMGetField(dm, field, NULL, &pressure));
PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
PetscCall(MatNullSpaceDestroy(&nullspacePres));
}
PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SetupProblem(DM dm, PetscErrorCode (*setupEqn)(DM, AppCtx *), AppCtx *user)
{
  DM              cdm = dm;
  PetscQuadrature q   = NULL;
  PetscBool       simplex;
  PetscInt        dim, Nf = 2, f, Nc[2];
  const char     *name[2]   = {"velocity", "pressure"};
  const char     *prefix[2] = {"vel_", "pres_"};

  PetscFunctionBegin;
  PetscCall(DMGetDimension(dm, &dim));
  PetscCall(DMPlexIsSimplex(dm, &simplex));
  Nc[0] = dim;
  Nc[1] = 1;
  for (f = 0; f < Nf; ++f) {
    PetscFE fe;

    PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc[f], simplex, prefix[f], -1, &fe));
    PetscCall(PetscObjectSetName((PetscObject)fe, name[f]));
    if (!q) PetscCall(PetscFEGetQuadrature(fe, &q));
    PetscCall(PetscFESetQuadrature(fe, q));
    PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
    PetscCall(PetscFEDestroy(&fe));
  }
  PetscCall(DMCreateDS(dm));
  PetscCall((*setupEqn)(dm, user));
  while (cdm) {
    PetscCall(DMCopyDisc(dm, cdm));
    PetscCall(DMSetNullSpaceConstructor(cdm, 1, CreatePressureNullSpace));
    PetscCall(DMGetCoarseDM(cdm, &cdm));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  SNES      snes;
  DM        dm;
  Vec       u;
  AppCtx    user;
  PetscBool testPatchFacetResidual = PETSC_FALSE;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_patch_facet_residual", &testPatchFacetResidual, NULL));
  if (testPatchFacetResidual) PetscCall(TestPatchFacetResidual());
  else {
    PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    PetscCall(SNESCreate(PetscObjectComm((PetscObject)dm), &snes));
    PetscCall(SNESSetDM(snes, dm));
    PetscCall(DMSetApplicationContext(dm, &user));

    PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    PetscCall(SetupProblem(dm, SetupEqn, &user));
    PetscCall(DMPlexCreateClosureIndex(dm, NULL));

    PetscCall(DMCreateGlobalVector(dm, &u));
    PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    PetscCall(SNESSetFromOptions(snes));
    PetscCall(DMSNESCheckFromOptions(snes, u));
    PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
    {
      Mat          J;
      MatNullSpace sp;

      PetscCall(SNESSetUp(snes));
      PetscCall(CreatePressureNullSpace(dm, 1, 1, &sp));
      PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
      PetscCall(MatSetNullSpace(J, sp));
      PetscCall(MatNullSpaceDestroy(&sp));
      PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
      PetscCall(MatViewFromOptions(J, NULL, "-J_view"));
    }
    PetscCall(SNESSolve(snes, NULL, u));

    PetscCall(VecDestroy(&u));
    PetscCall(SNESDestroy(&snes));
    PetscCall(DMDestroy(&dm));
    PetscCall(PetscBagDestroy(&user.bag));
  }
  PetscCall(PetscFinalize());
  return 0;
}
```

--------------------------------

### PETSc ex16.c: Large-deformation elasticity buckling example setup

Source: https://petsc.org/release/src/snes/tutorials/ex16.c.html

This is the beginning of the PETSc example ex16.c, a 3D large deformation elasticity buckling problem. It includes the help string, problem description, includes, constants, type definitions, function prototypes, and the start of main() that initializes PETSc, creates a 3D DMDA, and sets up the SNES object.

```c
static char help[] = "Large-deformation Elasticity Buckling Example";
```

```c
/*F-----------------------------------------------------------------------

    This example solves the 3D large deformation elasticity problem

\begin{equation}
 \int_{\Omega}F \cdot S : \nabla v d\Omega + \int_{\Omega} (loading)\mathbf{e}_y \cdot v d\Omega = 0
\end{equation}

    F is the deformation gradient, and S is the second Piola-Kirchhoff tensor from the Saint Venant-Kirchhoff model of
    hyperelasticity.  \Omega is a (arc) angle subsection of a cylindrical shell of thickness (height), inner radius
    (rad) and width (width).  The problem is discretized using Q1 finite elements on a logically structured grid.
    Homogeneous Dirichlet boundary conditions are applied at the centers of the ends of the sphere.

    This example is tunable with the following options:
    -rad : the radius of the circle
    -arc : set the angle of the arch represented
    -loading : set the bulk loading (the mass)
    -ploading : set the point loading at the top
    -height : set the height of the arch
    -width : set the width of the arch
    -view_line : print initial and final offsets of the centerline of the
                 beam along the x direction

    The material properties may be modified using either:
    -mu      : the first lame' parameter
    -lambda  : the second lame' parameter

    Or:
    -young   : the Young's modulus
    -poisson : the poisson ratio

    This example is meant to show the strain placed upon the nonlinear solvers when trying to "snap through" the arch
    using the loading.  Under certain parameter regimes, the arch will invert under the load, and the number of Newton
    steps will jump considerably.  Composed nonlinear solvers may be used to mitigate this difficulty.  This example
    also demonstrates the use of the arc length continuation method NEWTONAL, which avoids the numerical difficulties
    of the snap-through via tracing the equilibrium path through load increments.

    The initial setup follows the example in pg. 268 of "Nonlinear Finite Element Methods" by Peter Wriggers, but is a
    3D extension.

------------------------------------------------------------------------F*/
```

```c
#include <petscsnes.h>
#include <petscdm.h>
#include <petscdmda.h>
```

```c
#define QP0 0.2113248654051871
#define QP1 0.7886751345948129
#define NQ  2
#define NB  2
#define NEB 8
#define NEQ 8
#define NPB 24
```

```c
#define NVALS NEB *NEQ
const PetscReal pts[NQ] = {QP0, QP1};
const PetscReal wts[NQ] = {0.5, 0.5};
```

```c
PetscScalar vals[NVALS];
PetscScalar grad[3 * NVALS];
```

```c
typedef PetscScalar Field[3];
typedef PetscScalar CoordField[3];
```

```c
typedef PetscScalar JacField[9];
```

```c
PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field ***, Field ***, void *);
PetscErrorCode FormJacobianLocal(DMDALocalInfo *, Field ***, Mat, Mat, void *);
PetscErrorCode DisplayLine(SNES, Vec);
PetscErrorCode FormElements(void);
```

```c
typedef struct {
  PetscReal loading;
  PetscReal mu;
  PetscReal lambda;
  PetscReal rad;
  PetscReal height;
  PetscReal width;
  PetscReal arc;
  PetscReal ploading;
  PetscReal load_factor;
} AppCtx;
```

```c
PetscErrorCode InitialGuess(DM, AppCtx *, Vec);
PetscErrorCode ArcLengthScaling(DM, AppCtx *, Vec);
PetscErrorCode FormRHS(DM, AppCtx *, Vec);
PetscErrorCode FormCoordinates(DM, AppCtx *);
PetscErrorCode TangentLoad(SNES, Vec, Vec, void *);
```

```c
int main(int argc, char **argv)
{
  AppCtx    user; /* user-defined work context */
  PetscInt  mx, my, its;
  MPI_Comm  comm;
  SNES      snes;
  DM        da;
  Vec       x, X, b;
  PetscBool youngflg, poissonflg, muflg, lambdaflg, alflg, view = PETSC_FALSE, viewline = PETSC_FALSE, al_rescale = PETSC_FALSE;
  PetscReal poisson = 0.2, young = 4e4;
  char      filename[PETSC_MAX_PATH_LEN]     = "ex16.vts";
  char      filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  PetscCall(FormElements());
  comm = PETSC_COMM_WORLD;
  PetscCall(SNESCreate(comm, &snes));
  PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 11, 2, 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, NULL, &da));
  PetscCall(DMSetFromOptions(da));
  PetscCall(DMSetUp(da));
  PetscCall(SNESSetDM(snes, (DM)da));
```
Tokens:
Raw