Actual source code: aspin.c

  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdm.h>

  4: static PetscErrorCode MatMultASPIN(Mat m, Vec X, Vec Y)
  5: {
  6:   void       *ctx;
  7:   SNES        snes;
  8:   PetscInt    n, i;
  9:   VecScatter *oscatter;
 10:   SNES       *subsnes;
 11:   PetscBool   match;
 12:   MPI_Comm    comm;
 13:   KSP         ksp;
 14:   Vec        *x, *b;
 15:   Vec         W;
 16:   SNES        npc;
 17:   Mat         subJ, subpJ;

 19:   PetscFunctionBegin;
 20:   PetscCall(MatShellGetContext(m, &ctx));
 21:   snes = (SNES)ctx;
 22:   PetscCall(SNESGetNPC(snes, &npc));
 23:   PetscCall(SNESGetFunction(npc, &W, NULL, NULL));
 24:   PetscCall(PetscObjectTypeCompare((PetscObject)npc, SNESNASM, &match));
 25:   if (!match) {
 26:     PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
 27:     SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz");
 28:   }
 29:   PetscCall(SNESNASMGetSubdomains(npc, &n, &subsnes, NULL, &oscatter, NULL));
 30:   PetscCall(SNESNASMGetSubdomainVecs(npc, &n, &x, &b, NULL, NULL));

 32:   PetscCall(VecSet(Y, 0));
 33:   PetscCall(MatMult(npc->jacobian_pre, X, W));

 35:   for (i = 0; i < n; i++) PetscCall(VecScatterBegin(oscatter[i], W, b[i], INSERT_VALUES, SCATTER_FORWARD));
 36:   for (i = 0; i < n; i++) {
 37:     PetscCall(VecScatterEnd(oscatter[i], W, b[i], INSERT_VALUES, SCATTER_FORWARD));
 38:     PetscCall(VecSet(x[i], 0.));
 39:     PetscCall(SNESGetJacobian(subsnes[i], &subJ, &subpJ, NULL, NULL));
 40:     PetscCall(SNESGetKSP(subsnes[i], &ksp));
 41:     PetscCall(KSPSetOperators(ksp, subJ, subpJ));
 42:     PetscCall(KSPSolve(ksp, b[i], x[i]));
 43:     PetscCall(VecScatterBegin(oscatter[i], x[i], Y, ADD_VALUES, SCATTER_REVERSE));
 44:     PetscCall(VecScatterEnd(oscatter[i], x[i], Y, ADD_VALUES, SCATTER_REVERSE));
 45:   }
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: static PetscErrorCode SNESDestroy_ASPIN(SNES snes)
 50: {
 51:   PetscFunctionBegin;
 52:   PetscCall(SNESDestroy(&snes->npc));
 53:   /* reset NEWTONLS and free the data */
 54:   PetscCall(SNESReset(snes));
 55:   PetscCall(PetscFree(snes->data));
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }

 59: /*MC
 60:    SNESASPIN - Helper `SNES` type for Additive-Schwarz Preconditioned Inexact Newton {cite}`ck02`, {cite}`bruneknepleysmithtu15`

 62:    Options Database Keys:
 63: +  -npc_snes_     - options prefix of the nonlinear subdomain solver (must be of type `NASM`)
 64: .  -npc_sub_snes_ - options prefix of the subdomain nonlinear solves
 65: .  -npc_sub_ksp_  - options prefix of the subdomain Krylov solver
 66: -  -npc_sub_pc_   - options prefix of the subdomain preconditioner

 68:      Level: intermediate

 70:     Notes:
 71:     This solver transform the given nonlinear problem to a new form and then runs matrix-free Newton-Krylov with no
 72:     preconditioner on that transformed problem.

 74:     This routine sets up an instance of `SNESNETWONLS` with nonlinear left preconditioning.  It differs from other
 75:     similar functionality in `SNES` as it creates a linear shell matrix that corresponds to the product

 77:     $$
 78:     \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b))
 79:     $$

 81:     which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of
 82:     nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner
 83:     factorizations are reused on each application of $J_b^{-1}$.

 85:     The Krylov method used in this nonlinear solver is run with NO preconditioner, because the preconditioning is done
 86:     at the nonlinear level, but the Jacobian for the original function must be provided (or calculated via coloring and
 87:     finite differences automatically) in the Pmat location of `SNESSetJacobian()` because the action of the original Jacobian
 88:     is needed by the shell matrix used to apply the Jacobian of the nonlinear preconditioned problem (see above).
 89:     Note that since the Pmat is not used to construct a preconditioner it could be provided in a matrix-free form.
 90:     The code for this implementation is a bit confusing because the Amat of `SNESSetJacobian()` applies the Jacobian of the
 91:     nonlinearly preconditioned function Jacobian while the Pmat provides the Jacobian of the original user provided function.
 92:     Note that the original `SNES` and nonlinear preconditioner (see `SNESGetNPC()`), in this case `SNESNASM`, share
 93:     the same Jacobian matrices. `SNESNASM` computes the needed Jacobian in `SNESNASMComputeFinalJacobian_Private()`.

 95: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNASM`, `SNESGetNPC()`, `SNESGetNPCSide()`
 96: M*/
 97: PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes)
 98: {
 99:   SNES           npc;
100:   KSP            ksp;
101:   PC             pc;
102:   Mat            aspinmat;
103:   Vec            F;
104:   PetscInt       n;
105:   SNESLineSearch linesearch;

107:   PetscFunctionBegin;
108:   /* set up the solver */
109:   PetscCall(SNESSetType(snes, SNESNEWTONLS));
110:   PetscCall(SNESSetNPCSide(snes, PC_LEFT));
111:   PetscCall(SNESSetFunctionType(snes, SNES_FUNCTION_PRECONDITIONED));
112:   PetscCall(SNESGetNPC(snes, &npc));
113:   PetscCall(SNESSetType(npc, SNESNASM));
114:   PetscCall(SNESNASMSetType(npc, PC_ASM_BASIC));
115:   PetscCall(SNESNASMSetComputeFinalJacobian(npc, PETSC_TRUE));
116:   PetscCall(SNESGetKSP(snes, &ksp));
117:   PetscCall(KSPGetPC(ksp, &pc));
118:   PetscCall(PCSetType(pc, PCNONE));
119:   PetscCall(SNESGetLineSearch(snes, &linesearch));
120:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));

122:   /* set up the shell matrix */
123:   PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
124:   PetscCall(VecGetLocalSize(F, &n));
125:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)snes), n, n, PETSC_DECIDE, PETSC_DECIDE, snes, &aspinmat));
126:   PetscCall(MatSetType(aspinmat, MATSHELL));
127:   PetscCall(MatShellSetOperation(aspinmat, MATOP_MULT, (void (*)(void))MatMultASPIN));
128:   PetscCall(SNESSetJacobian(snes, aspinmat, NULL, NULL, NULL));
129:   PetscCall(MatDestroy(&aspinmat));

131:   snes->ops->destroy = SNESDestroy_ASPIN;
132:   PetscCall(PetscObjectChangeTypeName((PetscObject)snes, SNESASPIN));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }