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