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++) {
35: VecScatterBegin(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);
36: }
37: for (i=0;i<n;i++) {
38: VecScatterEnd(oscatter[i],W,b[i],INSERT_VALUES,SCATTER_FORWARD);
39: VecSet(x[i],0.);
40: SNESGetJacobian(subsnes[i],&subJ,&subpJ,NULL,NULL);
41: SNESGetKSP(subsnes[i],&ksp);
42: KSPSetOperators(ksp,subJ,subpJ);
43: KSPSolve(ksp,b[i],x[i]);
44: VecScatterBegin(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);
45: VecScatterEnd(oscatter[i],x[i],Y,ADD_VALUES,SCATTER_REVERSE);
46: }
47: return 0;
48: }
50: static PetscErrorCode SNESDestroy_ASPIN(SNES snes)
51: {
52: SNESDestroy(&snes->npc);
53: /* reset NEWTONLS and free the data */
54: SNESReset(snes);
55: PetscFree(snes->data);
56: return 0;
57: }
59: /* -------------------------------------------------------------------------- */
60: /*MC
61: SNESASPIN - Helper SNES type for Additive-Schwarz Preconditioned Inexact Newton
63: Options Database:
64: + -npc_snes_ - options prefix of the nonlinear subdomain solver (must be of type NASM)
65: . -npc_sub_snes_ - options prefix of the subdomain nonlinear solves
66: . -npc_sub_ksp_ - options prefix of the subdomain Krylov solver
67: - -npc_sub_pc_ - options prefix of the subdomain preconditioner
69: Notes:
70: This routine sets up an instance of NETWONLS 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 NASM, 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) {
122: SNESLineSearchSetType(linesearch,SNESLINESEARCHBT);
123: }
125: /* set up the shell matrix */
126: SNESGetFunction(snes,&F,NULL,NULL);
127: VecGetLocalSize(F,&n);
128: MatCreateShell(PetscObjectComm((PetscObject)snes),n,n,PETSC_DECIDE,PETSC_DECIDE,snes,&aspinmat);
129: MatSetType(aspinmat,MATSHELL);
130: MatShellSetOperation(aspinmat,MATOP_MULT,(void(*)(void))MatMultASPIN);
131: SNESSetJacobian(snes,aspinmat,NULL,NULL,NULL);
132: MatDestroy(&aspinmat);
134: snes->ops->destroy = SNESDestroy_ASPIN;
136: return 0;
137: }