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 {cite}`ck02`, {cite}`bruneknepleysmithtu15`
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: $$
78: \sum_{i=0}^{N_b}J_b({X^b_{converged}})^{-1}J(X + \sum_{i=0}^{N_b}(X^b_{converged} - X^b))
79: $$
81: which is the ASPIN preconditioned matrix. Similar solvers may be constructed by having matrix-free differencing of
82: nonlinear solves per linear iteration, but this is far more efficient when subdomain sparse-direct preconditioner
83: factorizations are reused on each application of $J_b^{-1}$.
85: The Krylov method used in this nonlinear solver is run with NO preconditioner, because the preconditioning is done
86: at the nonlinear level, but the Jacobian for the original function must be provided (or calculated via coloring and
87: finite differences automatically) in the Pmat location of `SNESSetJacobian()` because the action of the original Jacobian
88: is needed by the shell matrix used to apply the Jacobian of the nonlinear preconditioned problem (see above).
89: Note that since the Pmat is not used to construct a preconditioner it could be provided in a matrix-free form.
90: The code for this implementation is a bit confusing because the Amat of `SNESSetJacobian()` applies the Jacobian of the
91: nonlinearly preconditioned function Jacobian while the Pmat provides the Jacobian of the original user provided function.
92: Note that the original `SNES` and nonlinear preconditioner (see `SNESGetNPC()`), in this case `SNESNASM`, share
93: the same Jacobian matrices. `SNESNASM` computes the needed Jacobian in `SNESNASMComputeFinalJacobian_Private()`.
95: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNASM`, `SNESGetNPC()`, `SNESGetNPCSide()`
96: M*/
97: PETSC_EXTERN PetscErrorCode SNESCreate_ASPIN(SNES snes)
98: {
99: SNES npc;
100: KSP ksp;
101: PC pc;
102: Mat aspinmat;
103: Vec F;
104: PetscInt n;
105: SNESLineSearch linesearch;
107: PetscFunctionBegin;
108: /* set up the solver */
109: PetscCall(SNESSetType(snes, SNESNEWTONLS));
110: PetscCall(SNESSetNPCSide(snes, PC_LEFT));
111: PetscCall(SNESSetFunctionType(snes, SNES_FUNCTION_PRECONDITIONED));
112: PetscCall(SNESGetNPC(snes, &npc));
113: PetscCall(SNESSetType(npc, SNESNASM));
114: PetscCall(SNESNASMSetType(npc, PC_ASM_BASIC));
115: PetscCall(SNESNASMSetComputeFinalJacobian(npc, PETSC_TRUE));
116: PetscCall(SNESGetKSP(snes, &ksp));
117: PetscCall(KSPGetPC(ksp, &pc));
118: PetscCall(PCSetType(pc, PCNONE));
119: PetscCall(SNESGetLineSearch(snes, &linesearch));
120: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
122: /* set up the shell matrix */
123: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
124: PetscCall(VecGetLocalSize(F, &n));
125: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)snes), n, n, PETSC_DECIDE, PETSC_DECIDE, snes, &aspinmat));
126: PetscCall(MatSetType(aspinmat, MATSHELL));
127: PetscCall(MatShellSetOperation(aspinmat, MATOP_MULT, (void (*)(void))MatMultASPIN));
128: PetscCall(SNESSetJacobian(snes, aspinmat, NULL, NULL, NULL));
129: PetscCall(MatDestroy(&aspinmat));
131: snes->ops->destroy = SNESDestroy_ASPIN;
132: PetscCall(PetscObjectChangeTypeName((PetscObject)snes, SNESASPIN));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }