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: {
      7:   void           *ctx;
      8:   SNES           snes;
      9:   PetscInt       n,i;
     10:   VecScatter     *oscatter;
     11:   SNES           *subsnes;
     12:   PetscBool      match;
     13:   MPI_Comm       comm;
     14:   KSP            ksp;
     15:   Vec            *x,*b;
     16:   Vec            W;
     17:   SNES           npc;
     18:   Mat            subJ,subpJ;

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

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

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

     52: static PetscErrorCode SNESDestroy_ASPIN(SNES snes)
     53: {

     57:   SNESDestroy(&snes->npc);
     58:   /* reset NEWTONLS and free the data */
     59:   SNESReset(snes);
     60:   PetscFree(snes->data);
     61:   return(0);
     62: }

     64: /* -------------------------------------------------------------------------- */
     65: /*MC
     66:       SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton

     68:    Options Database:
     69: +  -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM)
     70: .  -npc_sub_snes_ - options prefix of the subdomain nonlinear solves
     71: .  -npc_sub_ksp_ - options prefix of the subdomain Krylov solver
     72: -  -npc_sub_pc_ - options prefix of the subdomain preconditioner

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

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

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

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

     94:    Level: intermediate

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

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

    103: M*/
    104: PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes)
    105: {
    107:   SNES           npc;
    108:   KSP            ksp;
    109:   PC             pc;
    110:   Mat            aspinmat;
    111:   Vec            F;
    112:   PetscInt       n;
    113:   SNESLineSearch linesearch;

    116:   /* set up the solver */
    117:   SNESSetType(snes,SNESNEWTONLS);
    118:   SNESSetNPCSide(snes,PC_LEFT);
    119:   SNESSetFunctionType(snes,SNES_FUNCTION_PRECONDITIONED);
    120:   SNESGetNPC(snes,&npc);
    121:   SNESSetType(npc,SNESNASM);
    122:   SNESNASMSetType(npc,PC_ASM_BASIC);
    123:   SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);
    124:   SNESGetKSP(snes,&ksp);
    125:   KSPGetPC(ksp,&pc);
    126:   PCSetType(pc,PCNONE);
    127:   SNESGetLineSearch(snes,&linesearch);
    128:   if (!((PetscObject)linesearch)->type_name) {
    129:     SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);
    130:   }

    132:   /* set up the shell matrix */
    133:   SNESGetFunction(snes,&F,NULL,NULL);
    134:   VecGetLocalSize(F,&n);
    135:   MatCreateShell(PetscObjectComm((PetscObject)snes),n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);
    136:   MatSetType(aspinmat,MATSHELL);
    137:   MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);
    138:   SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);
    139:   MatDestroy(&aspinmat);

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

    143:   return(0);
    144: }