Actual source code: aspin.c

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

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

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

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

 48: static PetscErrorCode SNESDestroy_ASPIN(SNES snes)
 49: {
 50:   SNESDestroy(&snes->npc);
 51:   /* reset NEWTONLS and free the data */
 52:   SNESReset(snes);
 53:   PetscFree(snes->data);
 54:   return 0;
 55: }

 57: /*MC
 58:       SNESASPIN - Helper `SNES` type for Additive-Schwarz Preconditioned Inexact Newton

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

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

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

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

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

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

 89:    Level: intermediate

 91:    References:
 92: +  * - X. C. Cai and D. E. Keyes, "Nonlinearly preconditioned inexact Newton algorithms",  SIAM J. Sci. Comput., 24, 2002.
 93: -  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
 94:    SIAM Review, 57(4), 2015

 96: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNASM`, `SNESGetNPC()`, `SNESGetNPCSide()`

 98: M*/
 99: PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes)
100: {
101:   SNES           npc;
102:   KSP            ksp;
103:   PC             pc;
104:   Mat            aspinmat;
105:   Vec            F;
106:   PetscInt       n;
107:   SNESLineSearch linesearch;

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

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

132:   snes->ops->destroy = SNESDestroy_ASPIN;

134:   return 0;
135: }