1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdm.h> 3 4 PetscErrorCode MatMultASPIN(Mat m,Vec X,Vec Y) 5 { 6 PetscErrorCode ierr; 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; 19 20 PetscFunctionBegin; 21 ierr = MatShellGetContext(m,&ctx);CHKERRQ(ierr); 22 snes = (SNES)ctx; 23 ierr = SNESGetNPC(snes,&npc);CHKERRQ(ierr); 24 ierr = SNESGetFunction(npc,&W,NULL,NULL);CHKERRQ(ierr); 25 ierr = PetscObjectTypeCompare((PetscObject)npc,SNESNASM,&match);CHKERRQ(ierr); 26 if (!match) { 27 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 28 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz"); 29 } 30 ierr = SNESNASMGetSubdomains(npc,&n,&subsnes,NULL,&oscatter,NULL);CHKERRQ(ierr); 31 ierr = SNESNASMGetSubdomainVecs(npc,&n,&x,&b,NULL,NULL);CHKERRQ(ierr); 32 33 ierr = VecSet(Y,0);CHKERRQ(ierr); 34 ierr = MatMult(npc->jacobian_pre,X,W);CHKERRQ(ierr); 35 36 for (i=0;i<n;i++) { 37 ierr = VecScatterBegin(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 38 } 39 for (i=0;i<n;i++) { 40 ierr = VecScatterEnd(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 41 ierr = VecSet(x[i],0.);CHKERRQ(ierr); 42 ierr = SNESGetJacobian(subsnes[i],&subJ,&subpJ,NULL,NULL);CHKERRQ(ierr); 43 ierr = SNESGetKSP(subsnes[i],&ksp);CHKERRQ(ierr); 44 ierr = KSPSetOperators(ksp,subJ,subpJ);CHKERRQ(ierr); 45 ierr = KSPSolve(ksp,b[i],x[i]);CHKERRQ(ierr); 46 ierr = VecScatterBegin(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 47 ierr = VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 48 } 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode SNESDestroy_ASPIN(SNES snes) 53 { 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = SNESDestroy(&snes->npc);CHKERRQ(ierr); 58 /* reset NEWTONLS and free the data */ 59 ierr = SNESReset(snes);CHKERRQ(ierr); 60 ierr = PetscFree(snes->data);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 /* -------------------------------------------------------------------------- */ 65 /*MC 66 SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton 67 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 73 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 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 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}. 83 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(). 93 94 Level: intermediate 95 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 100 101 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNASM, SNESGetNPC(), SNESGetNPCSide() 102 103 M*/ 104 PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes) 105 { 106 PetscErrorCode ierr; 107 SNES npc; 108 KSP ksp; 109 PC pc; 110 Mat aspinmat; 111 Vec F; 112 PetscInt n; 113 SNESLineSearch linesearch; 114 115 PetscFunctionBegin; 116 /* set up the solver */ 117 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 118 ierr = SNESSetNPCSide(snes,PC_LEFT);CHKERRQ(ierr); 119 ierr = SNESSetFunctionType(snes,SNES_FUNCTION_PRECONDITIONED);CHKERRQ(ierr); 120 ierr = SNESGetNPC(snes,&npc);CHKERRQ(ierr); 121 ierr = SNESSetType(npc,SNESNASM);CHKERRQ(ierr); 122 ierr = SNESNASMSetType(npc,PC_ASM_BASIC);CHKERRQ(ierr); 123 ierr = SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);CHKERRQ(ierr); 124 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 125 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 126 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 127 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 128 if (!((PetscObject)linesearch)->type_name) { 129 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);CHKERRQ(ierr); 130 } 131 132 /* set up the shell matrix */ 133 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 134 ierr = VecGetLocalSize(F,&n);CHKERRQ(ierr); 135 ierr = MatCreateShell(PetscObjectComm((PetscObject)snes),n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);CHKERRQ(ierr); 136 ierr = MatSetType(aspinmat,MATSHELL);CHKERRQ(ierr); 137 ierr = MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);CHKERRQ(ierr); 138 ierr = SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);CHKERRQ(ierr); 139 ierr = MatDestroy(&aspinmat);CHKERRQ(ierr); 140 141 snes->ops->destroy = SNESDestroy_ASPIN; 142 143 PetscFunctionReturn(0); 144 } 145