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 SNES npc; 19 20 PetscFunctionBegin; 21 ierr = MatShellGetContext(m,&ctx);CHKERRQ(ierr); 22 snes = (SNES)ctx; 23 ierr = SNESGetPC(snes,&npc);CHKERRQ(ierr); 24 ierr = PetscObjectTypeCompare((PetscObject)npc,SNESNASM,&match);CHKERRQ(ierr); 25 if (!match) { 26 ierr = PetscObjectGetComm((PetscObject)snes,&comm); 27 SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz"); 28 } 29 ierr = SNESNASMGetSubdomains(npc,&n,&subsnes,NULL,&oscatter,NULL);CHKERRQ(ierr); 30 ierr = SNESNASMGetSubdomainVecs(npc,&n,&x,&b,NULL,NULL);CHKERRQ(ierr); 31 ierr = MatMult(npc->jacobian_pre,X,Y);CHKERRQ(ierr); 32 for (i=0;i<n;i++) { 33 ierr = VecScatterBegin(oscatter[i],Y,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 34 ierr = VecScatterEnd(oscatter[i],Y,b[i],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 35 } 36 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 37 for (i=0;i<n;i++) { 38 ierr = VecSet(x[i],0);CHKERRQ(ierr); 39 ierr = SNESGetKSP(subsnes[i],&ksp);CHKERRQ(ierr); 40 ierr = KSPSolve(ksp,b[i],x[i]);CHKERRQ(ierr); 41 } 42 ierr = VecSet(Y,0);CHKERRQ(ierr); 43 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 44 for (i=0;i<n;i++) { 45 ierr = VecScatterBegin(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 46 ierr = VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 47 } 48 49 ierr = MPI_Barrier(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 50 51 PetscFunctionReturn(0); 52 } 53 54 EXTERN_C_BEGIN 55 #undef __FUNCT__ 56 #define __FUNCT__ "SNESCreate_ASPIN" 57 PetscErrorCode SNESCreate_ASPIN(SNES snes) 58 { 59 PetscErrorCode ierr; 60 SNES npc; 61 KSP ksp; 62 PC pc; 63 Mat aspinmat; 64 MPI_Comm comm; 65 Vec F; 66 PetscInt n; 67 68 PetscFunctionBegin; 69 /* set up the solver */ 70 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 71 ierr = SNESSetPCSide(snes,PC_LEFT);CHKERRQ(ierr); 72 ierr = SNESGetPC(snes,&npc);CHKERRQ(ierr); 73 ierr = SNESSetType(npc,SNESNASM);CHKERRQ(ierr); 74 ierr = SNESNASMSetComputeFinalJacobian(npc,PETSC_TRUE);CHKERRQ(ierr); 75 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 76 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 77 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 78 ierr = PetscObjectGetComm((PetscObject)snes,&comm); 79 80 /* set up the shell matrix */ 81 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr); 82 ierr = VecGetLocalSize(F,&n);CHKERRQ(ierr); 83 ierr = MatCreateShell(comm,n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);CHKERRQ(ierr); 84 ierr = MatSetType(aspinmat,MATSHELL);CHKERRQ(ierr); 85 ierr = MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);CHKERRQ(ierr); 86 87 ierr = SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);CHKERRQ(ierr); 88 ierr = MatDestroy(&aspinmat);CHKERRQ(ierr); 89 90 PetscFunctionReturn(0); 91 } 92 EXTERN_C_END 93