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, (PetscErrorCodeFn *)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: }