Actual source code: aspin.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
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;
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));
32: PetscCall(VecSet(Y, 0));
33: PetscCall(MatMult(npc->jacobian_pre, X, W));
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: }
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: }
59: /*MC
60: SNESASPIN - Helper `SNES` type for Additive-Schwarz Preconditioned Inexact Newton
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
68: Level: intermediate
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.
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
77: \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b))
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}.
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()`.
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
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;
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));
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));
134: snes->ops->destroy = SNESDestroy_ASPIN;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }