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