1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdm.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatMultASPIN" 6 PetscErrorCode MatMultASPIN(Mat m,Vec X,Vec Y) 7 { 8 PetscErrorCode ierr; 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; 21 22 PetscFunctionBegin; 23 ierr = MatShellGetContext(m,&ctx);CHKERRQ(ierr); 24 snes = (SNES)ctx; 25 ierr = SNESGetNPC(snes,&npc);CHKERRQ(ierr); 26 ierr = SNESGetFunction(npc,&W,NULL,NULL);CHKERRQ(ierr); 27 ierr = PetscObjectTypeCompare((PetscObject)npc,SNESNASM,&match);CHKERRQ(ierr); 28 if (!match) { 29 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 30 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz"); 31 } 32 ierr = SNESNASMGetSubdomains(npc,&n,&subsnes,NULL,&oscatter,NULL);CHKERRQ(ierr); 33 ierr = SNESNASMGetSubdomainVecs(npc,&n,&x,&b,NULL,NULL);CHKERRQ(ierr); 34 35 ierr = VecSet(Y,0);CHKERRQ(ierr); 36 ierr = MatMult(npc->jacobian_pre,X,W);CHKERRQ(ierr); 37 38 for (i=0;i<n;i++) { 39 ierr = VecScatterBegin(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 40 } 41 for (i=0;i<n;i++) { 42 ierr = VecScatterEnd(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 43 ierr = VecSet(x[i],0.);CHKERRQ(ierr); 44 ierr = SNESGetJacobian(subsnes[i],&subJ,&subpJ,NULL,NULL);CHKERRQ(ierr); 45 ierr = SNESGetKSP(subsnes[i],&ksp);CHKERRQ(ierr); 46 ierr = KSPSetOperators(ksp,subJ,subpJ);CHKERRQ(ierr); 47 ierr = KSPSolve(ksp,b[i],x[i]);CHKERRQ(ierr); 48 ierr = VecScatterBegin(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 49 } 50 for (i=0;i<n;i++) { 51 ierr = VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 52 } 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "SNESDestroy_ASPIN" 58 PetscErrorCode SNESDestroy_ASPIN(SNES snes) 59 { 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 64 /* reset NEWTONLS and free the data */ 65 ierr = SNESReset(snes);CHKERRQ(ierr); 66 ierr = PetscFree(snes->data);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "SNESCreate_ASPIN" 72 /* -------------------------------------------------------------------------- */ 73 /*MC 74 SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton 75 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 81 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: 84 85 \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b)) 86 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}. 90 91 Level: intermediate 92 93 References: 94 95 X.-C. Cai and D. E. Keyes, "Nonlinearly preconditioned inexact {Newton} algorithms", SIAM J. Sci. Comput., 24, 2002. 96 97 Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 98 SIAM Review, 57(4), 2015 99 100 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNASM, SNESGetNPC(), SNESGetNPCSide() 101 102 M*/ 103 PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes) 104 { 105 PetscErrorCode ierr; 106 SNES npc; 107 KSP ksp; 108 PC pc; 109 Mat aspinmat; 110 Vec F; 111 PetscInt n; 112 SNESLineSearch linesearch; 113 114 PetscFunctionBegin; 115 /* set up the solver */ 116 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 117 ierr = SNESSetNPCSide(snes,PC_LEFT);CHKERRQ(ierr); 118 ierr = SNESSetFunctionType(snes,SNES_FUNCTION_PRECONDITIONED);CHKERRQ(ierr); 119 ierr = SNESGetNPC(snes,&npc);CHKERRQ(ierr); 120 ierr = SNESSetType(npc,SNESNASM);CHKERRQ(ierr); 121 ierr = SNESNASMSetType(npc,PC_ASM_BASIC);CHKERRQ(ierr); 122 ierr = SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);CHKERRQ(ierr); 123 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 124 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 125 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 126 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 127 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);CHKERRQ(ierr); 128 129 /* set up the shell matrix */ 130 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 131 ierr = VecGetLocalSize(F,&n);CHKERRQ(ierr); 132 ierr = MatCreateShell(PetscObjectComm((PetscObject)snes),n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);CHKERRQ(ierr); 133 ierr = MatSetType(aspinmat,MATSHELL);CHKERRQ(ierr); 134 ierr = MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);CHKERRQ(ierr); 135 ierr = SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);CHKERRQ(ierr); 136 ierr = MatDestroy(&aspinmat);CHKERRQ(ierr); 137 138 snes->ops->destroy = SNESDestroy_ASPIN; 139 140 PetscFunctionReturn(0); 141 } 142