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 21 PetscFunctionBegin; 22 ierr = MatShellGetContext(m,&ctx);CHKERRQ(ierr); 23 snes = (SNES)ctx; 24 ierr = SNESGetPC(snes,&npc);CHKERRQ(ierr); 25 ierr = SNESGetFunction(npc,&W,NULL,NULL);CHKERRQ(ierr); 26 ierr = PetscObjectTypeCompare((PetscObject)npc,SNESNASM,&match);CHKERRQ(ierr); 27 if (!match) { 28 ierr = PetscObjectGetComm((PetscObject)snes,&comm); 29 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz"); 30 } 31 ierr = SNESNASMGetSubdomains(npc,&n,&subsnes,NULL,&oscatter,NULL);CHKERRQ(ierr); 32 ierr = SNESNASMGetSubdomainVecs(npc,&n,&x,&b,NULL,NULL);CHKERRQ(ierr); 33 34 ierr = VecSet(Y,0);CHKERRQ(ierr); 35 ierr = MatMult(npc->jacobian_pre,X,W);CHKERRQ(ierr); 36 37 for (i=0;i<n;i++) { 38 ierr = VecScatterBegin(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 39 } 40 for (i=0;i<n;i++) { 41 ierr = VecScatterEnd(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 42 ierr = VecSet(x[i],0.);CHKERRQ(ierr); 43 ierr = SNESGetKSP(subsnes[i],&ksp);CHKERRQ(ierr); 44 ierr = KSPSolve(ksp,b[i],x[i]);CHKERRQ(ierr); 45 ierr = VecScatterBegin(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 46 } 47 for (i=0;i<n;i++) { 48 ierr = VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 49 } 50 PetscFunctionReturn(0); 51 } 52 53 EXTERN_C_BEGIN 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESCreate_ASPIN" 56 /* -------------------------------------------------------------------------- */ 57 /*MC 58 SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton 59 60 Options Database: 61 + -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM) 62 . -npc_sub_snes_ - options prefix of the subdomain nonlinear solves 63 . -npc_sub_ksp_ - options prefix of the subdomain Krylov solver 64 - -npc_sub_pc_ - options prefix of the subdomain preconditioner 65 66 Notes: This routine sets up an instance of NETWONLS with nonlinear left preconditioning. It differs from other 67 similar functionality in SNES as it creates a linear shell matrix that corresponds to the product: 68 69 \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b)) 70 71 which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of 72 nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner 73 factorizations are reused on each application of J_b^{-1}. 74 75 Level: intermediate 76 77 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNASM, SNESGetPC(), SNESGetPCSide() 78 79 M*/ 80 PetscErrorCode SNESCreate_ASPIN(SNES snes) 81 { 82 PetscErrorCode ierr; 83 SNES npc; 84 KSP ksp; 85 PC pc; 86 Mat aspinmat; 87 MPI_Comm comm; 88 Vec F; 89 PetscInt n; 90 SNESLineSearch linesearch; 91 92 PetscFunctionBegin; 93 /* set up the solver */ 94 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 95 ierr = SNESSetPCSide(snes,PC_LEFT);CHKERRQ(ierr); 96 ierr = SNESGetPC(snes,&npc);CHKERRQ(ierr); 97 ierr = SNESSetType(npc,SNESNASM);CHKERRQ(ierr); 98 ierr = SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);CHKERRQ(ierr); 99 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 100 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 101 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 102 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 103 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 104 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);CHKERRQ(ierr); 105 106 /* set up the shell matrix */ 107 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 108 ierr = VecGetLocalSize(F,&n);CHKERRQ(ierr); 109 ierr = MatCreateShell(comm,n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);CHKERRQ(ierr); 110 ierr = MatSetType(aspinmat,MATSHELL);CHKERRQ(ierr); 111 ierr = MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);CHKERRQ(ierr); 112 113 ierr = SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);CHKERRQ(ierr); 114 ierr = MatDestroy(&aspinmat);CHKERRQ(ierr); 115 116 PetscFunctionReturn(0); 117 } 118 EXTERN_C_END 119