Actual source code: aspin.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
4: 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: Notes:
69: This solver transform the given nonlinear problem to a new form and then runs matrix-free Newton-Krylov with no
70: preconditioner on that transformed problem.
72: This routine sets up an instance of `SNESNETWONLS` with nonlinear left preconditioning. It differs from other
73: similar functionality in `SNES` as it creates a linear shell matrix that corresponds to the product
75: \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b))
77: which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of
78: nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner
79: factorizations are reused on each application of J_b^{-1}.
81: The Krylov method used in this nonlinear solver is run with NO preconditioner, because the preconditioning is done
82: at the nonlinear level, but the Jacobian for the original function must be provided (or calculated via coloring and
83: finite differences automatically) in the Pmat location of `SNESSetJacobian()` because the action of the original Jacobian
84: is needed by the shell matrix used to apply the Jacobian of the nonlinear preconditioned problem (see above).
85: Note that since the Pmat is not used to construct a preconditioner it could be provided in a matrix-free form.
86: The code for this implementation is a bit confusing because the Amat of `SNESSetJacobian()` applies the Jacobian of the
87: nonlinearly preconditioned function Jacobian while the Pmat provides the Jacobian of the original user provided function.
88: Note that the original `SNES` and nonlinear preconditioner preconditioner (see `SNESGetNPC()`), in this case `SNESNASM`, share
89: the same Jacobian matrices. `SNESNASM` computes the needed Jacobian in `SNESNASMComputeFinalJacobian_Private()`.
91: Level: intermediate
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: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNASM`, `SNESGetNPC()`, `SNESGetNPCSide()`
100: M*/
101: PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes)
102: {
103: SNES npc;
104: KSP ksp;
105: PC pc;
106: Mat aspinmat;
107: Vec F;
108: PetscInt n;
109: SNESLineSearch linesearch;
111: PetscFunctionBegin;
112: /* set up the solver */
113: PetscCall(SNESSetType(snes, SNESNEWTONLS));
114: PetscCall(SNESSetNPCSide(snes, PC_LEFT));
115: PetscCall(SNESSetFunctionType(snes, SNES_FUNCTION_PRECONDITIONED));
116: PetscCall(SNESGetNPC(snes, &npc));
117: PetscCall(SNESSetType(npc, SNESNASM));
118: PetscCall(SNESNASMSetType(npc, PC_ASM_BASIC));
119: PetscCall(SNESNASMSetComputeFinalJacobian(npc, PETSC_TRUE));
120: PetscCall(SNESGetKSP(snes, &ksp));
121: PetscCall(KSPGetPC(ksp, &pc));
122: PetscCall(PCSetType(pc, PCNONE));
123: PetscCall(SNESGetLineSearch(snes, &linesearch));
124: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
126: /* set up the shell matrix */
127: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
128: PetscCall(VecGetLocalSize(F, &n));
129: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)snes), n, n, PETSC_DECIDE, PETSC_DECIDE, snes, &aspinmat));
130: PetscCall(MatSetType(aspinmat, MATSHELL));
131: PetscCall(MatShellSetOperation(aspinmat, MATOP_MULT, (void (*)(void))MatMultASPIN));
132: PetscCall(SNESSetJacobian(snes, aspinmat, NULL, NULL, NULL));
133: PetscCall(MatDestroy(&aspinmat));
135: snes->ops->destroy = SNESDestroy_ASPIN;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }