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 void *ctx; 6 SNES snes; 7 PetscInt n, i; 8 VecScatter *oscatter; 9 SNES *subsnes; 10 PetscBool match; 11 MPI_Comm comm; 12 KSP ksp; 13 Vec *x, *b; 14 Vec W; 15 SNES npc; 16 Mat subJ, subpJ; 17 18 PetscFunctionBegin; 19 PetscCall(MatShellGetContext(m, &ctx)); 20 snes = (SNES)ctx; 21 PetscCall(SNESGetNPC(snes, &npc)); 22 PetscCall(SNESGetFunction(npc, &W, NULL, NULL)); 23 PetscCall(PetscObjectTypeCompare((PetscObject)npc, SNESNASM, &match)); 24 if (!match) { 25 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 26 SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "MatMultASPIN requires that the nonlinear preconditioner be Nonlinear additive Schwarz"); 27 } 28 PetscCall(SNESNASMGetSubdomains(npc, &n, &subsnes, NULL, &oscatter, NULL)); 29 PetscCall(SNESNASMGetSubdomainVecs(npc, &n, &x, &b, NULL, NULL)); 30 31 PetscCall(VecSet(Y, 0)); 32 PetscCall(MatMult(npc->jacobian_pre, X, W)); 33 34 for (i = 0; i < n; i++) PetscCall(VecScatterBegin(oscatter[i], W, b[i], INSERT_VALUES, SCATTER_FORWARD)); 35 for (i = 0; i < n; i++) { 36 PetscCall(VecScatterEnd(oscatter[i], W, b[i], INSERT_VALUES, SCATTER_FORWARD)); 37 PetscCall(VecSet(x[i], 0.)); 38 PetscCall(SNESGetJacobian(subsnes[i], &subJ, &subpJ, NULL, NULL)); 39 PetscCall(SNESGetKSP(subsnes[i], &ksp)); 40 PetscCall(KSPSetOperators(ksp, subJ, subpJ)); 41 PetscCall(KSPSolve(ksp, b[i], x[i])); 42 PetscCall(VecScatterBegin(oscatter[i], x[i], Y, ADD_VALUES, SCATTER_REVERSE)); 43 PetscCall(VecScatterEnd(oscatter[i], x[i], Y, ADD_VALUES, SCATTER_REVERSE)); 44 } 45 PetscFunctionReturn(0); 46 } 47 48 static PetscErrorCode SNESDestroy_ASPIN(SNES snes) { 49 PetscFunctionBegin; 50 PetscCall(SNESDestroy(&snes->npc)); 51 /* reset NEWTONLS and free the data */ 52 PetscCall(SNESReset(snes)); 53 PetscCall(PetscFree(snes->data)); 54 PetscFunctionReturn(0); 55 } 56 57 /* -------------------------------------------------------------------------- */ 58 /*MC 59 SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton 60 61 Options Database: 62 + -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM) 63 . -npc_sub_snes_ - options prefix of the subdomain nonlinear solves 64 . -npc_sub_ksp_ - options prefix of the subdomain Krylov solver 65 - -npc_sub_pc_ - options prefix of the subdomain preconditioner 66 67 Notes: 68 This routine sets up an instance of NETWONLS with nonlinear left preconditioning. It differs from other 69 similar functionality in SNES as it creates a linear shell matrix that corresponds to the product 70 71 \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b)) 72 73 which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of 74 nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner 75 factorizations are reused on each application of J_b^{-1}. 76 77 The Krylov method used in this nonlinear solver is run with NO preconditioner, because the preconditioning is done 78 at the nonlinear level, but the Jacobian for the original function must be provided (or calculated via coloring and 79 finite differences automatically) in the Pmat location of SNESSetJacobian() because the action of the original Jacobian 80 is needed by the shell matrix used to apply the Jacobian of the nonlinear preconditioned problem (see above). 81 Note that since the Pmat is not used to construct a preconditioner it could be provided in a matrix-free form. 82 The code for this implementation is a bit confusing because the Amat of SNESSetJacobian() applies the Jacobian of the 83 nonlinearly preconditioned function Jacobian while the Pmat provides the Jacobian of the original user provided function. 84 Note that the original SNES and nonlinear preconditioner preconditioner (see SNESGetNPC()), in this case NASM, share 85 the same Jacobian matrices. SNESNASM computes the needed Jacobian in SNESNASMComputeFinalJacobian_Private(). 86 87 Level: intermediate 88 89 References: 90 + * - X. C. Cai and D. E. Keyes, "Nonlinearly preconditioned inexact Newton algorithms", SIAM J. Sci. Comput., 24, 2002. 91 - * - 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 SNES npc; 99 KSP ksp; 100 PC pc; 101 Mat aspinmat; 102 Vec F; 103 PetscInt n; 104 SNESLineSearch linesearch; 105 106 PetscFunctionBegin; 107 /* set up the solver */ 108 PetscCall(SNESSetType(snes, SNESNEWTONLS)); 109 PetscCall(SNESSetNPCSide(snes, PC_LEFT)); 110 PetscCall(SNESSetFunctionType(snes, SNES_FUNCTION_PRECONDITIONED)); 111 PetscCall(SNESGetNPC(snes, &npc)); 112 PetscCall(SNESSetType(npc, SNESNASM)); 113 PetscCall(SNESNASMSetType(npc, PC_ASM_BASIC)); 114 PetscCall(SNESNASMSetComputeFinalJacobian(npc, PETSC_TRUE)); 115 PetscCall(SNESGetKSP(snes, &ksp)); 116 PetscCall(KSPGetPC(ksp, &pc)); 117 PetscCall(PCSetType(pc, PCNONE)); 118 PetscCall(SNESGetLineSearch(snes, &linesearch)); 119 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 120 121 /* set up the shell matrix */ 122 PetscCall(SNESGetFunction(snes, &F, NULL, NULL)); 123 PetscCall(VecGetLocalSize(F, &n)); 124 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)snes), n, n, PETSC_DECIDE, PETSC_DECIDE, snes, &aspinmat)); 125 PetscCall(MatSetType(aspinmat, MATSHELL)); 126 PetscCall(MatShellSetOperation(aspinmat, MATOP_MULT, (void (*)(void))MatMultASPIN)); 127 PetscCall(SNESSetJacobian(snes, aspinmat, NULL, NULL, NULL)); 128 PetscCall(MatDestroy(&aspinmat)); 129 130 snes->ops->destroy = SNESDestroy_ASPIN; 131 132 PetscFunctionReturn(0); 133 } 134