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 } 48 for (i=0;i<n;i++) { 49 ierr = VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 50 } 51 PetscFunctionReturn(0); 52 } 53 54 static PetscErrorCode SNESDestroy_ASPIN(SNES snes) 55 { 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = SNESDestroy(&snes->npc);CHKERRQ(ierr); 60 /* reset NEWTONLS and free the data */ 61 ierr = SNESReset(snes);CHKERRQ(ierr); 62 ierr = PetscFree(snes->data);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 /* -------------------------------------------------------------------------- */ 67 /*MC 68 SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton 69 70 Options Database: 71 + -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM) 72 . -npc_sub_snes_ - options prefix of the subdomain nonlinear solves 73 . -npc_sub_ksp_ - options prefix of the subdomain Krylov solver 74 - -npc_sub_pc_ - options prefix of the subdomain preconditioner 75 76 Notes: 77 This routine sets up an instance of NETWONLS with nonlinear left preconditioning. It differs from other 78 similar functionality in SNES as it creates a linear shell matrix that corresponds to the product 79 80 \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b)) 81 82 which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of 83 nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner 84 factorizations are reused on each application of J_b^{-1}. 85 86 Level: 87 intermediate 88 89 References: 90 + 1. - X. C. Cai and D. E. Keyes, "Nonlinearly preconditioned inexact Newton algorithms", SIAM J. Sci. Comput., 24, 2002. 91 - 2. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 92 SIAM Review, 57(4), 2015 93 94 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNASM, SNESGetNPC(), SNESGetNPCSide() 95 96 M*/ 97 PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes) 98 { 99 PetscErrorCode ierr; 100 SNES npc; 101 KSP ksp; 102 PC pc; 103 Mat aspinmat; 104 Vec F; 105 PetscInt n; 106 SNESLineSearch linesearch; 107 108 PetscFunctionBegin; 109 /* set up the solver */ 110 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 111 ierr = SNESSetNPCSide(snes,PC_LEFT);CHKERRQ(ierr); 112 ierr = SNESSetFunctionType(snes,SNES_FUNCTION_PRECONDITIONED);CHKERRQ(ierr); 113 ierr = SNESGetNPC(snes,&npc);CHKERRQ(ierr); 114 ierr = SNESSetType(npc,SNESNASM);CHKERRQ(ierr); 115 ierr = SNESNASMSetType(npc,PC_ASM_BASIC);CHKERRQ(ierr); 116 ierr = SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);CHKERRQ(ierr); 117 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 118 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 119 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 120 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 121 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);CHKERRQ(ierr); 122 123 /* set up the shell matrix */ 124 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 125 ierr = VecGetLocalSize(F,&n);CHKERRQ(ierr); 126 ierr = MatCreateShell(PetscObjectComm((PetscObject)snes),n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);CHKERRQ(ierr); 127 ierr = MatSetType(aspinmat,MATSHELL);CHKERRQ(ierr); 128 ierr = MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);CHKERRQ(ierr); 129 ierr = SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);CHKERRQ(ierr); 130 ierr = MatDestroy(&aspinmat);CHKERRQ(ierr); 131 132 snes->ops->destroy = SNESDestroy_ASPIN; 133 134 PetscFunctionReturn(0); 135 } 136