1: #include <petsc/private/snesimpl.h>             /*I   "petscsnes.h"   I*/
  2: #include <petscdm.h>

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

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

 35:   VecSet(Y,0);
 36:   MatMult(npc->jacobian_pre,X,W);

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

 58: PetscErrorCode SNESDestroy_ASPIN(SNES snes)
 59: {

 63:   SNESDestroy(&snes->pc);
 64:   /* reset NEWTONLS and free the data */
 65:   SNESReset(snes);
 66:   PetscFree(snes->data);
 67:   return(0);
 68: }

 72: /* -------------------------------------------------------------------------- */
 73: /*MC
 74:       SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton

 76:    Options Database:
 77: +  -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM)
 78: .  -npc_sub_snes_ - options prefix of the subdomain nonlinear solves
 79: .  -npc_sub_ksp_ - options prefix of the subdomain Krylov solver
 80: -  -npc_sub_pc_ - options prefix of the subdomain preconditioner

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

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

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

 91:    Level: intermediate

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

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

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

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

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

133:   return(0);
134: }