Actual source code: nasm.c

  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdm.h>

  4: typedef struct {
  5:   PetscInt    n;             /* local subdomains */
  6:   SNES       *subsnes;       /* nonlinear solvers for each subdomain */
  7:   Vec        *x;             /* solution vectors */
  8:   Vec        *xl;            /* solution local vectors */
  9:   Vec        *y;             /* step vectors */
 10:   Vec        *b;             /* rhs vectors */
 11:   Vec         weight;        /* weighting for adding updates on overlaps, in global space */
 12:   VecScatter *oscatter;      /* scatter from global space to the subdomain global space */
 13:   VecScatter *oscatter_copy; /* copy of the above */
 14:   VecScatter *iscatter;      /* scatter from global space to the nonoverlapping subdomain space */
 15:   VecScatter *gscatter;      /* scatter from global space to the subdomain local space */
 16:   PCASMType   type;          /* ASM type */
 17:   PetscBool   usesdm;        /* use the DM for setting up the subproblems */
 18:   PetscBool   finaljacobian; /* compute the jacobian of the converged solution */
 19:   PetscReal   damping;       /* damping parameter for updates from the blocks */
 20:   PetscBool   weight_set;    /* use a weight in the overlap updates */

 22:   /* logging events */
 23:   PetscLogEvent eventrestrictinterp;
 24:   PetscLogEvent eventsubsolve;

 26:   PetscInt fjtype; /* type of computed jacobian */
 27:   Vec      xinit;  /* initial solution in case the final jacobian type is computed as first */
 28: } SNES_NASM;

 30: const char *const SNESNASMTypes[]   = {"NONE", "RESTRICT", "INTERPOLATE", "BASIC", "PCASMType", "PC_ASM_", NULL};
 31: const char *const SNESNASMFJTypes[] = {"FINALOUTER", "FINALINNER", "INITIAL"};

 33: static PetscErrorCode SNESReset_NASM(SNES snes)
 34: {
 35:   SNES_NASM *nasm = (SNES_NASM *)snes->data;
 36:   PetscInt   i;

 38:   for (i = 0; i < nasm->n; i++) {
 39:     if (nasm->xl) VecDestroy(&nasm->xl[i]);
 40:     if (nasm->x) VecDestroy(&nasm->x[i]);
 41:     if (nasm->y) VecDestroy(&nasm->y[i]);
 42:     if (nasm->b) VecDestroy(&nasm->b[i]);

 44:     if (nasm->subsnes) SNESDestroy(&nasm->subsnes[i]);
 45:     if (nasm->oscatter) VecScatterDestroy(&nasm->oscatter[i]);
 46:     if (nasm->oscatter_copy) VecScatterDestroy(&nasm->oscatter_copy[i]);
 47:     if (nasm->iscatter) VecScatterDestroy(&nasm->iscatter[i]);
 48:     if (nasm->gscatter) VecScatterDestroy(&nasm->gscatter[i]);
 49:   }

 51:   PetscFree(nasm->x);
 52:   PetscFree(nasm->xl);
 53:   PetscFree(nasm->y);
 54:   PetscFree(nasm->b);

 56:   if (nasm->xinit) VecDestroy(&nasm->xinit);

 58:   PetscFree(nasm->subsnes);
 59:   PetscFree(nasm->oscatter);
 60:   PetscFree(nasm->oscatter_copy);
 61:   PetscFree(nasm->iscatter);
 62:   PetscFree(nasm->gscatter);

 64:   if (nasm->weight_set) VecDestroy(&nasm->weight);

 66:   nasm->eventrestrictinterp = 0;
 67:   nasm->eventsubsolve       = 0;
 68:   return 0;
 69: }

 71: static PetscErrorCode SNESDestroy_NASM(SNES snes)
 72: {
 73:   SNESReset_NASM(snes);
 74:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetType_C", NULL);
 75:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetType_C", NULL);
 76:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", NULL);
 77:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", NULL);
 78:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetDamping_C", NULL);
 79:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetDamping_C", NULL);
 80:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", NULL);
 81:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", NULL);
 82:   PetscFree(snes->data);
 83:   return 0;
 84: }

 86: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
 87: {
 88:   Vec bcs = (Vec)ctx;

 90:   VecCopy(bcs, l);
 91:   return 0;
 92: }

 94: static PetscErrorCode SNESSetUp_NASM(SNES snes)
 95: {
 96:   SNES_NASM  *nasm = (SNES_NASM *)snes->data;
 97:   DM          dm, subdm;
 98:   DM         *subdms;
 99:   PetscInt    i;
100:   const char *optionsprefix;
101:   Vec         F;
102:   PetscMPIInt size;
103:   KSP         ksp;
104:   PC          pc;

106:   if (!nasm->subsnes) {
107:     SNESGetDM(snes, &dm);
108:     if (dm) {
109:       nasm->usesdm = PETSC_TRUE;
110:       DMCreateDomainDecomposition(dm, &nasm->n, NULL, NULL, NULL, &subdms);
112:       DMCreateDomainDecompositionScatters(dm, nasm->n, subdms, &nasm->iscatter, &nasm->oscatter, &nasm->gscatter);
113:       PetscMalloc1(nasm->n, &nasm->oscatter_copy);
114:       for (i = 0; i < nasm->n; i++) VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);

116:       SNESGetOptionsPrefix(snes, &optionsprefix);
117:       PetscMalloc1(nasm->n, &nasm->subsnes);
118:       for (i = 0; i < nasm->n; i++) {
119:         SNESCreate(PETSC_COMM_SELF, &nasm->subsnes[i]);
120:         PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
121:         SNESAppendOptionsPrefix(nasm->subsnes[i], optionsprefix);
122:         SNESAppendOptionsPrefix(nasm->subsnes[i], "sub_");
123:         SNESSetDM(nasm->subsnes[i], subdms[i]);
124:         MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]), &size);
125:         if (size == 1) {
126:           SNESGetKSP(nasm->subsnes[i], &ksp);
127:           KSPGetPC(ksp, &pc);
128:           KSPSetType(ksp, KSPPREONLY);
129:           PCSetType(pc, PCLU);
130:         }
131:         SNESSetFromOptions(nasm->subsnes[i]);
132:         DMDestroy(&subdms[i]);
133:       }
134:       PetscFree(subdms);
135:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Cannot construct local problems automatically without a DM!");
136:   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must set subproblems manually if there is no DM!");
137:   /* allocate the global vectors */
138:   if (!nasm->x) PetscCalloc1(nasm->n, &nasm->x);
139:   if (!nasm->xl) PetscCalloc1(nasm->n, &nasm->xl);
140:   if (!nasm->y) PetscCalloc1(nasm->n, &nasm->y);
141:   if (!nasm->b) PetscCalloc1(nasm->n, &nasm->b);

143:   for (i = 0; i < nasm->n; i++) {
144:     SNESGetFunction(nasm->subsnes[i], &F, NULL, NULL);
145:     if (!nasm->x[i]) VecDuplicate(F, &nasm->x[i]);
146:     if (!nasm->y[i]) VecDuplicate(F, &nasm->y[i]);
147:     if (!nasm->b[i]) VecDuplicate(F, &nasm->b[i]);
148:     if (!nasm->xl[i]) {
149:       SNESGetDM(nasm->subsnes[i], &subdm);
150:       DMCreateLocalVector(subdm, &nasm->xl[i]);
151:       DMGlobalToLocalHookAdd(subdm, DMGlobalToLocalSubDomainDirichletHook_Private, NULL, nasm->xl[i]);
152:     }
153:   }
154:   if (nasm->finaljacobian) {
155:     SNESSetUpMatrices(snes);
156:     if (nasm->fjtype == 2) VecDuplicate(snes->vec_sol, &nasm->xinit);
157:     for (i = 0; i < nasm->n; i++) SNESSetUpMatrices(nasm->subsnes[i]);
158:   }
159:   return 0;
160: }

162: static PetscErrorCode SNESSetFromOptions_NASM(SNES snes, PetscOptionItems *PetscOptionsObject)
163: {
164:   PCASMType  asmtype;
165:   PetscBool  flg, monflg;
166:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

168:   PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Additive Schwarz options");
169:   PetscOptionsEnum("-snes_nasm_type", "Type of restriction/extension", "", SNESNASMTypes, (PetscEnum)nasm->type, (PetscEnum *)&asmtype, &flg);
170:   if (flg) SNESNASMSetType(snes, asmtype);
171:   flg    = PETSC_FALSE;
172:   monflg = PETSC_TRUE;
173:   PetscOptionsReal("-snes_nasm_damping", "The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)", "SNESNASMSetDamping", nasm->damping, &nasm->damping, &flg);
174:   if (flg) SNESNASMSetDamping(snes, nasm->damping);
175:   PetscOptionsDeprecated("-snes_nasm_sub_view", NULL, "3.15", "Use -snes_view ::ascii_info_detail");
176:   PetscOptionsBool("-snes_nasm_finaljacobian", "Compute the global jacobian of the final iterate (for ASPIN)", "", nasm->finaljacobian, &nasm->finaljacobian, NULL);
177:   PetscOptionsEList("-snes_nasm_finaljacobian_type", "The type of the final jacobian computed.", "", SNESNASMFJTypes, 3, SNESNASMFJTypes[0], &nasm->fjtype, NULL);
178:   PetscOptionsBool("-snes_nasm_log", "Log times for subSNES solves and restriction", "", monflg, &monflg, &flg);
179:   if (flg) {
180:     PetscLogEventRegister("SNESNASMSubSolve", ((PetscObject)snes)->classid, &nasm->eventsubsolve);
181:     PetscLogEventRegister("SNESNASMRestrict", ((PetscObject)snes)->classid, &nasm->eventrestrictinterp);
182:   }
183:   PetscOptionsHeadEnd();
184:   return 0;
185: }

187: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
188: {
189:   SNES_NASM        *nasm = (SNES_NASM *)snes->data;
190:   PetscMPIInt       rank, size;
191:   PetscInt          i, N, bsz;
192:   PetscBool         iascii, isstring;
193:   PetscViewer       sviewer;
194:   MPI_Comm          comm;
195:   PetscViewerFormat format;
196:   const char       *prefix;

198:   PetscObjectGetComm((PetscObject)snes, &comm);
199:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
200:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
201:   MPI_Comm_rank(comm, &rank);
202:   MPI_Comm_size(comm, &size);
203:   MPIU_Allreduce(&nasm->n, &N, 1, MPIU_INT, MPI_SUM, comm);
204:   if (iascii) {
205:     PetscViewerASCIIPrintf(viewer, "  total subdomain blocks = %" PetscInt_FMT "\n", N);
206:     PetscViewerGetFormat(viewer, &format);
207:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
208:       if (nasm->subsnes) {
209:         PetscViewerASCIIPrintf(viewer, "  Local solver information for first block on rank 0:\n");
210:         SNESGetOptionsPrefix(snes, &prefix);
211:         PetscViewerASCIIPrintf(viewer, "  Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : "");
212:         PetscViewerASCIIPushTab(viewer);
213:         PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
214:         if (rank == 0) {
215:           PetscViewerASCIIPushTab(viewer);
216:           SNESView(nasm->subsnes[0], sviewer);
217:           PetscViewerASCIIPopTab(viewer);
218:         }
219:         PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
220:         PetscViewerASCIIPopTab(viewer);
221:       }
222:     } else {
223:       /* print the solver on each block */
224:       PetscViewerASCIIPushSynchronized(viewer);
225:       PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] number of local blocks = %" PetscInt_FMT "\n", (int)rank, nasm->n);
226:       PetscViewerFlush(viewer);
227:       PetscViewerASCIIPopSynchronized(viewer);
228:       PetscViewerASCIIPrintf(viewer, "  Local solver information for each block is in the following SNES objects:\n");
229:       PetscViewerASCIIPushTab(viewer);
230:       PetscViewerASCIIPrintf(viewer, "- - - - - - - - - - - - - - - - - -\n");
231:       PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
232:       for (i = 0; i < nasm->n; i++) {
233:         VecGetLocalSize(nasm->x[i], &bsz);
234:         PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT ", size = %" PetscInt_FMT "\n", (int)rank, i, bsz);
235:         SNESView(nasm->subsnes[i], sviewer);
236:         PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n");
237:       }
238:       PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
239:       PetscViewerFlush(viewer);
240:       PetscViewerASCIIPopTab(viewer);
241:     }
242:   } else if (isstring) {
243:     PetscViewerStringSPrintf(viewer, " blocks=%" PetscInt_FMT ",type=%s", N, SNESNASMTypes[nasm->type]);
244:     PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
245:     if (nasm->subsnes && rank == 0) SNESView(nasm->subsnes[0], sviewer);
246:     PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
247:   }
248:   return 0;
249: }

251: /*@
252:    SNESNASMSetType - Set the type of subdomain update used for the nonlinear additive Schwarz solver

254:    Logically Collective

256:    Input Parameters:
257: +  snes - the `SNES` context
258: -  type - the type of update, `PC_ASM_BASIC` or `PC_ASM_RESTRICT`

260:    Level: intermediate

262: .seealso: `SNESNASM`, `SNESNASMGetType()`, `PCASMSetType()`, `PC_ASM_BASIC`, `PC_ASM_RESTRICT`, `PCASMType`
263: @*/
264: PetscErrorCode SNESNASMSetType(SNES snes, PCASMType type)
265: {
266:   PetscErrorCode (*f)(SNES, PCASMType);

268:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetType_C", &f);
269:   if (f) (f)(snes, type);
270:   return 0;
271: }

273: static PetscErrorCode SNESNASMSetType_NASM(SNES snes, PCASMType type)
274: {
275:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

278:   nasm->type = type;
279:   return 0;
280: }

282: /*@
283:    SNESNASMGetType - Get the type of subdomain update used for the nonlinear additive Schwarz solver

285:    Logically Collective

287:    Input Parameters:
288: .  snes - the `SNES` context

290:    Output Parameters:
291: .  type - the type of update

293:    Level: intermediate

295: .seealso: `SNESNASM`, `SNESNASMSetType()`, `PCASMGetType()`, `PC_ASM_BASIC`, `PC_ASM_RESTRICT`, `PCASMType`
296: @*/
297: PetscErrorCode SNESNASMGetType(SNES snes, PCASMType *type)
298: {
299:   PetscUseMethod(snes, "SNESNASMGetType_C", (SNES, PCASMType *), (snes, type));
300:   return 0;
301: }

303: static PetscErrorCode SNESNASMGetType_NASM(SNES snes, PCASMType *type)
304: {
305:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

307:   *type = nasm->type;
308:   return 0;
309: }

311: /*@
312:    SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems in the nonlinear additive Schwarz solver

314:    Logically Collective

316:    Input Parameters:
317: +  snes - the SNES context
318: .  n - the number of local subdomains
319: .  subsnes - solvers defined on the local subdomains
320: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
321: .  oscatter - scatters into the overlapping portions of the local subdomains
322: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

324:    Level: intermediate

326: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
327: @*/
328: PetscErrorCode SNESNASMSetSubdomains(SNES snes, PetscInt n, SNES subsnes[], VecScatter iscatter[], VecScatter oscatter[], VecScatter gscatter[])
329: {
330:   PetscErrorCode (*f)(SNES, PetscInt, SNES *, VecScatter *, VecScatter *, VecScatter *);

332:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", &f);
333:   if (f) (f)(snes, n, subsnes, iscatter, oscatter, gscatter);
334:   return 0;
335: }

337: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes, PetscInt n, SNES subsnes[], VecScatter iscatter[], VecScatter oscatter[], VecScatter gscatter[])
338: {
339:   PetscInt   i;
340:   SNES_NASM *nasm = (SNES_NASM *)snes->data;


344:   /* tear down the previously set things */
345:   SNESReset(snes);

347:   nasm->n = n;
348:   if (oscatter) {
349:     for (i = 0; i < n; i++) PetscObjectReference((PetscObject)oscatter[i]);
350:   }
351:   if (iscatter) {
352:     for (i = 0; i < n; i++) PetscObjectReference((PetscObject)iscatter[i]);
353:   }
354:   if (gscatter) {
355:     for (i = 0; i < n; i++) PetscObjectReference((PetscObject)gscatter[i]);
356:   }
357:   if (oscatter) {
358:     PetscMalloc1(n, &nasm->oscatter);
359:     PetscMalloc1(n, &nasm->oscatter_copy);
360:     for (i = 0; i < n; i++) {
361:       nasm->oscatter[i] = oscatter[i];
362:       VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
363:     }
364:   }
365:   if (iscatter) {
366:     PetscMalloc1(n, &nasm->iscatter);
367:     for (i = 0; i < n; i++) nasm->iscatter[i] = iscatter[i];
368:   }
369:   if (gscatter) {
370:     PetscMalloc1(n, &nasm->gscatter);
371:     for (i = 0; i < n; i++) nasm->gscatter[i] = gscatter[i];
372:   }

374:   if (subsnes) {
375:     PetscMalloc1(n, &nasm->subsnes);
376:     for (i = 0; i < n; i++) nasm->subsnes[i] = subsnes[i];
377:   }
378:   return 0;
379: }

381: /*@
382:    SNESNASMGetSubdomains - Get the local subdomain contexts for the nonlinear additive Schwarz solver

384:    Not Collective but some of the objects returned will be parallel

386:    Input Parameter:
387: .  snes - the `SNES` context

389:    Output Parameters:
390: +  n - the number of local subdomains
391: .  subsnes - solvers defined on the local subdomains
392: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
393: .  oscatter - scatters into the overlapping portions of the local subdomains
394: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

396:    Level: intermediate

398: .seealso: `SNESNASM`, `SNESNASMSetSubdomains()`
399: @*/
400: PetscErrorCode SNESNASMGetSubdomains(SNES snes, PetscInt *n, SNES *subsnes[], VecScatter *iscatter[], VecScatter *oscatter[], VecScatter *gscatter[])
401: {
402:   PetscErrorCode (*f)(SNES, PetscInt *, SNES **, VecScatter **, VecScatter **, VecScatter **);

404:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", &f);
405:   if (f) (f)(snes, n, subsnes, iscatter, oscatter, gscatter);
406:   return 0;
407: }

409: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes, PetscInt *n, SNES *subsnes[], VecScatter *iscatter[], VecScatter *oscatter[], VecScatter *gscatter[])
410: {
411:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

413:   if (n) *n = nasm->n;
414:   if (oscatter) *oscatter = nasm->oscatter;
415:   if (iscatter) *iscatter = nasm->iscatter;
416:   if (gscatter) *gscatter = nasm->gscatter;
417:   if (subsnes) *subsnes = nasm->subsnes;
418:   return 0;
419: }

421: /*@
422:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors for the nonlinear additive Schwarz solver

424:    Not Collective

426:    Input Parameter:
427: .  snes - the `SNES` context

429:    Output Parameters:
430: +  n - the number of local subdomains
431: .  x - The subdomain solution vector
432: .  y - The subdomain step vector
433: .  b - The subdomain RHS vector
434: -  xl - The subdomain local vectors (ghosted)

436:    Level: developer

438: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
439: @*/
440: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes, PetscInt *n, Vec **x, Vec **y, Vec **b, Vec **xl)
441: {
442:   PetscErrorCode (*f)(SNES, PetscInt *, Vec **, Vec **, Vec **, Vec **);

444:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", &f);
445:   if (f) (f)(snes, n, x, y, b, xl);
446:   return 0;
447: }

449: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes, PetscInt *n, Vec **x, Vec **y, Vec **b, Vec **xl)
450: {
451:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

453:   if (n) *n = nasm->n;
454:   if (x) *x = nasm->x;
455:   if (y) *y = nasm->y;
456:   if (b) *b = nasm->b;
457:   if (xl) *xl = nasm->xl;
458:   return 0;
459: }

461: /*@
462:    SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence for the
463:    nonlinear additive Schwarz solver

465:    Collective

467:    Input Parameters:
468: +  snes - the SNES context
469: -  flg - `PETSC_TRUE` to compute the Jacobians

471:    Level: developer

473:    Notes:
474:    This is used almost exclusively in the implementation of `SNESASPIN`, where the converged subdomain and global Jacobian
475:    is needed at each linear iteration.

477: .seealso: `SNESNASM`, `SNESNASMGetSubdomains()`
478: @*/
479: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes, PetscBool flg)
480: {
481:   PetscErrorCode (*f)(SNES, PetscBool);

483:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", &f);
484:   if (f) (f)(snes, flg);
485:   return 0;
486: }

488: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes, PetscBool flg)
489: {
490:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

492:   nasm->finaljacobian = flg;
493:   return 0;
494: }

496: /*@
497:    SNESNASMSetDamping - Sets the update damping for `SNESNASM` the nonlinear additive Schwarz solver

499:    Logically collective

501:    Input Parameters:
502: +  snes - the `SNES` context
503: -  dmp - damping

505:    Level: intermediate

507:    Notes:
508:     The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)

510: .seealso: `SNESNASM`, `SNESNASMGetDamping()`
511: @*/
512: PetscErrorCode SNESNASMSetDamping(SNES snes, PetscReal dmp)
513: {
514:   PetscErrorCode (*f)(SNES, PetscReal);

516:   PetscObjectQueryFunction((PetscObject)snes, "SNESNASMSetDamping_C", (void (**)(void)) & f);
517:   if (f) (f)(snes, dmp);
518:   return 0;
519: }

521: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes, PetscReal dmp)
522: {
523:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

525:   nasm->damping = dmp;
526:   return 0;
527: }

529: /*@
530:    SNESNASMGetDamping - Gets the update damping for `SNESNASM` the nonlinear additive Schwarz solver

532:    Not Collective

534:    Input Parameter:
535: .  snes - the SNES context

537:    Output Parameter:
538: .  dmp - damping

540:    Level: intermediate

542: .seealso: `SNESNASM`, `SNESNASMSetDamping()`
543: @*/
544: PetscErrorCode SNESNASMGetDamping(SNES snes, PetscReal *dmp)
545: {
546:   PetscUseMethod(snes, "SNESNASMGetDamping_C", (SNES, PetscReal *), (snes, dmp));
547:   return 0;
548: }

550: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes, PetscReal *dmp)
551: {
552:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

554:   *dmp = nasm->damping;
555:   return 0;
556: }

558: /*
559:   Input Parameters:
560: + snes - The solver
561: . B - The RHS vector
562: - X - The initial guess

564:   Output Parameters:
565: . Y - The solution update

567:   TODO: All scatters should be packed into one
568: */
569: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes, Vec B, Vec Y, Vec X)
570: {
571:   SNES_NASM *nasm = (SNES_NASM *)snes->data;
572:   SNES       subsnes;
573:   PetscInt   i;
574:   PetscReal  dmp;
575:   Vec        Xl, Bl, Yl, Xlloc;
576:   VecScatter iscat, oscat, gscat, oscat_copy;
577:   DM         dm, subdm;
578:   PCASMType  type;

580:   SNESNASMGetType(snes, &type);
581:   SNESGetDM(snes, &dm);
582:   VecSet(Y, 0);
583:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
584:   for (i = 0; i < nasm->n; i++) {
585:     /* scatter the solution to the global solution and the local solution */
586:     Xl         = nasm->x[i];
587:     Xlloc      = nasm->xl[i];
588:     oscat      = nasm->oscatter[i];
589:     oscat_copy = nasm->oscatter_copy[i];
590:     gscat      = nasm->gscatter[i];
591:     VecScatterBegin(oscat, X, Xl, INSERT_VALUES, SCATTER_FORWARD);
592:     VecScatterBegin(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
593:     if (B) {
594:       /* scatter the RHS to the local RHS */
595:       Bl = nasm->b[i];
596:       VecScatterBegin(oscat_copy, B, Bl, INSERT_VALUES, SCATTER_FORWARD);
597:     }
598:   }
599:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);

601:   if (nasm->eventsubsolve) PetscLogEventBegin(nasm->eventsubsolve, snes, 0, 0, 0);
602:   for (i = 0; i < nasm->n; i++) {
603:     Xl      = nasm->x[i];
604:     Xlloc   = nasm->xl[i];
605:     Yl      = nasm->y[i];
606:     subsnes = nasm->subsnes[i];
607:     SNESGetDM(subsnes, &subdm);
608:     iscat      = nasm->iscatter[i];
609:     oscat      = nasm->oscatter[i];
610:     oscat_copy = nasm->oscatter_copy[i];
611:     gscat      = nasm->gscatter[i];
612:     VecScatterEnd(oscat, X, Xl, INSERT_VALUES, SCATTER_FORWARD);
613:     VecScatterEnd(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
614:     if (B) {
615:       Bl = nasm->b[i];
616:       VecScatterEnd(oscat_copy, B, Bl, INSERT_VALUES, SCATTER_FORWARD);
617:     } else Bl = NULL;

619:     DMSubDomainRestrict(dm, oscat, gscat, subdm);
620:     VecCopy(Xl, Yl);
621:     SNESSolve(subsnes, Bl, Xl);
622:     VecAYPX(Yl, -1.0, Xl);
623:     VecScale(Yl, nasm->damping);
624:     if (type == PC_ASM_BASIC) {
625:       VecScatterBegin(oscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
626:       VecScatterEnd(oscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
627:     } else if (type == PC_ASM_RESTRICT) {
628:       VecScatterBegin(iscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
629:       VecScatterEnd(iscat, Yl, Y, ADD_VALUES, SCATTER_REVERSE);
630:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Only basic and restrict types are supported for SNESNASM");
631:   }
632:   if (nasm->eventsubsolve) PetscLogEventEnd(nasm->eventsubsolve, snes, 0, 0, 0);
633:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
634:   if (nasm->weight_set) VecPointwiseMult(Y, Y, nasm->weight);
635:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);
636:   SNESNASMGetDamping(snes, &dmp);
637:   VecAXPY(X, dmp, Y);
638:   return 0;
639: }

641: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
642: {
643:   Vec        X    = Xfinal;
644:   SNES_NASM *nasm = (SNES_NASM *)snes->data;
645:   SNES       subsnes;
646:   PetscInt   i, lag = 1;
647:   Vec        Xlloc, Xl, Fl, F;
648:   VecScatter oscat, gscat;
649:   DM         dm, subdm;

651:   if (nasm->fjtype == 2) X = nasm->xinit;
652:   F = snes->vec_func;
653:   if (snes->normschedule == SNES_NORM_NONE) SNESComputeFunction(snes, X, F);
654:   SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
655:   SNESGetDM(snes, &dm);
656:   if (nasm->eventrestrictinterp) PetscLogEventBegin(nasm->eventrestrictinterp, snes, 0, 0, 0);
657:   if (nasm->fjtype != 1) {
658:     for (i = 0; i < nasm->n; i++) {
659:       Xlloc = nasm->xl[i];
660:       gscat = nasm->gscatter[i];
661:       VecScatterBegin(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
662:     }
663:   }
664:   if (nasm->eventrestrictinterp) PetscLogEventEnd(nasm->eventrestrictinterp, snes, 0, 0, 0);
665:   for (i = 0; i < nasm->n; i++) {
666:     Fl      = nasm->subsnes[i]->vec_func;
667:     Xl      = nasm->x[i];
668:     Xlloc   = nasm->xl[i];
669:     subsnes = nasm->subsnes[i];
670:     oscat   = nasm->oscatter[i];
671:     gscat   = nasm->gscatter[i];
672:     if (nasm->fjtype != 1) VecScatterEnd(gscat, X, Xlloc, INSERT_VALUES, SCATTER_FORWARD);
673:     SNESGetDM(subsnes, &subdm);
674:     DMSubDomainRestrict(dm, oscat, gscat, subdm);
675:     if (nasm->fjtype != 1) {
676:       DMLocalToGlobalBegin(subdm, Xlloc, INSERT_VALUES, Xl);
677:       DMLocalToGlobalEnd(subdm, Xlloc, INSERT_VALUES, Xl);
678:     }
679:     if (subsnes->lagjacobian == -1) subsnes->lagjacobian = -2;
680:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
681:     SNESComputeFunction(subsnes, Xl, Fl);
682:     SNESComputeJacobian(subsnes, Xl, subsnes->jacobian, subsnes->jacobian_pre);
683:     if (lag > 1) subsnes->lagjacobian = lag;
684:   }
685:   return 0;
686: }

688: static PetscErrorCode SNESSolve_NASM(SNES snes)
689: {
690:   Vec              F;
691:   Vec              X;
692:   Vec              B;
693:   Vec              Y;
694:   PetscInt         i;
695:   PetscReal        fnorm = 0.0;
696:   SNESNormSchedule normschedule;
697:   SNES_NASM       *nasm = (SNES_NASM *)snes->data;



702:   PetscCitationsRegister(SNESCitation, &SNEScite);
703:   X = snes->vec_sol;
704:   Y = snes->vec_sol_update;
705:   F = snes->vec_func;
706:   B = snes->vec_rhs;

708:   PetscObjectSAWsTakeAccess((PetscObject)snes);
709:   snes->iter = 0;
710:   snes->norm = 0.;
711:   PetscObjectSAWsGrantAccess((PetscObject)snes);
712:   snes->reason = SNES_CONVERGED_ITERATING;
713:   SNESGetNormSchedule(snes, &normschedule);
714:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
715:     /* compute the initial function and preconditioned update delX */
716:     if (!snes->vec_func_init_set) {
717:       SNESComputeFunction(snes, X, F);
718:     } else snes->vec_func_init_set = PETSC_FALSE;

720:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
721:     SNESCheckFunctionNorm(snes, fnorm);
722:     PetscObjectSAWsTakeAccess((PetscObject)snes);
723:     snes->iter = 0;
724:     snes->norm = fnorm;
725:     PetscObjectSAWsGrantAccess((PetscObject)snes);
726:     SNESLogConvergenceHistory(snes, snes->norm, 0);
727:     SNESMonitor(snes, 0, snes->norm);

729:     /* test convergence */
730:     PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
731:     if (snes->reason) return 0;
732:   } else {
733:     PetscObjectSAWsGrantAccess((PetscObject)snes);
734:     SNESLogConvergenceHistory(snes, snes->norm, 0);
735:     SNESMonitor(snes, 0, snes->norm);
736:   }

738:   /* Call general purpose update function */
739:   PetscTryTypeMethod(snes, update, snes->iter);
740:   /* copy the initial solution over for later */
741:   if (nasm->fjtype == 2) VecCopy(X, nasm->xinit);

743:   for (i = 0; i < snes->max_its; i++) {
744:     SNESNASMSolveLocal_Private(snes, B, Y, X);
745:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
746:       SNESComputeFunction(snes, X, F);
747:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
748:       SNESCheckFunctionNorm(snes, fnorm);
749:     }
750:     /* Monitor convergence */
751:     PetscObjectSAWsTakeAccess((PetscObject)snes);
752:     snes->iter = i + 1;
753:     snes->norm = fnorm;
754:     PetscObjectSAWsGrantAccess((PetscObject)snes);
755:     SNESLogConvergenceHistory(snes, snes->norm, 0);
756:     SNESMonitor(snes, snes->iter, snes->norm);
757:     /* Test for convergence */
758:     if (normschedule == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
759:     if (snes->reason) break;
760:     /* Call general purpose update function */
761:     PetscTryTypeMethod(snes, update, snes->iter);
762:   }
763:   if (nasm->finaljacobian) {
764:     SNESNASMComputeFinalJacobian_Private(snes, X);
765:     SNESCheckJacobianDomainerror(snes);
766:   }
767:   if (normschedule == SNES_NORM_ALWAYS) {
768:     if (i == snes->max_its) {
769:       PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
770:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
771:     }
772:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a nonlinear preconditioner */
773:   return 0;
774: }

776: /*MC
777:   SNESNASM - Nonlinear Additive Schwarz solver

779:    Options Database Keys:
780: +  -snes_nasm_log - enable logging events for the communication and solve stages
781: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
782: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
783: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
784: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
785: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
786: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
787: -  -sub_pc_ - options prefix of the subdomain preconditioner

789:    Level: advanced

791:    Note:
792:    This is not often used directly as a solver, it converges too slowly. However it works well as a nonlinear preconditioner for
793:    the `SNESASPIN` solver

795:    Developer Note:
796:    This is a non-Newton based nonlinear solver that does not directly require a Jacobian; hence the flag snes->usesksp is set to
797:        false and `SNESView()` and -snes_view do not display a `KSP` object. However, if the flag nasm->finaljacobian is set (for example, if
798:        `SNESNASM` is used as a nonlinear preconditioner for  `SNESASPIN`) then `SNESSetUpMatrices()` is called to generate the
799:        Jacobian (needed by `SNESASPIN`)
800:        and this utilizes the inner `KSP` object for storing the matrices, but the `KSP` is never used for solving a linear system. When `SNESNASM` is
801:        used by `SNESASPIN` they share the same Jacobian matrices because `SNESSetUp()` (called on the outer `SNESASPIN`) causes the inner `SNES`
802:        object (in this case `SNESNASM`) to inherit the outer Jacobian matrices.

804:    References:
805: .  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
806:    SIAM Review, 57(4), 2015

808: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`, `SNESNASMSetType()`, `SNESNASMGetType()`, `SNESNASMSetSubdomains()`, `SNESNASMGetSubdomains()`, `SNESNASMGetSubdomainVecs()`, `SNESNASMSetComputeFinalJacobian()`, `SNESNASMSetDamping()`, `SNESNASMGetDamping()`
809: M*/

811: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
812: {
813:   SNES_NASM *nasm;

815:   PetscNew(&nasm);
816:   snes->data = (void *)nasm;

818:   nasm->n             = PETSC_DECIDE;
819:   nasm->subsnes       = NULL;
820:   nasm->x             = NULL;
821:   nasm->xl            = NULL;
822:   nasm->y             = NULL;
823:   nasm->b             = NULL;
824:   nasm->oscatter      = NULL;
825:   nasm->oscatter_copy = NULL;
826:   nasm->iscatter      = NULL;
827:   nasm->gscatter      = NULL;
828:   nasm->damping       = 1.;

830:   nasm->type          = PC_ASM_BASIC;
831:   nasm->finaljacobian = PETSC_FALSE;
832:   nasm->weight_set    = PETSC_FALSE;

834:   snes->ops->destroy        = SNESDestroy_NASM;
835:   snes->ops->setup          = SNESSetUp_NASM;
836:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
837:   snes->ops->view           = SNESView_NASM;
838:   snes->ops->solve          = SNESSolve_NASM;
839:   snes->ops->reset          = SNESReset_NASM;

841:   snes->usesksp = PETSC_FALSE;
842:   snes->usesnpc = PETSC_FALSE;

844:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

846:   nasm->fjtype              = 0;
847:   nasm->xinit               = NULL;
848:   nasm->eventrestrictinterp = 0;
849:   nasm->eventsubsolve       = 0;

851:   if (!snes->tolerancesset) {
852:     snes->max_its   = 10000;
853:     snes->max_funcs = 10000;
854:   }

856:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetType_C", SNESNASMSetType_NASM);
857:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetType_C", SNESNASMGetType_NASM);
858:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetSubdomains_C", SNESNASMSetSubdomains_NASM);
859:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomains_C", SNESNASMGetSubdomains_NASM);
860:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetDamping_C", SNESNASMSetDamping_NASM);
861:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetDamping_C", SNESNASMGetDamping_NASM);
862:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMGetSubdomainVecs_C", SNESNASMGetSubdomainVecs_NASM);
863:   PetscObjectComposeFunction((PetscObject)snes, "SNESNASMSetComputeFinalJacobian_C", SNESNASMSetComputeFinalJacobian_NASM);
864:   return 0;
865: }

867: /*@
868:    SNESNASMGetSNES - Gets a subsolver

870:    Not collective

872:    Input Parameters:
873: +  snes - the SNES context
874: -  i - the number of the subsnes to get

876:    Output Parameters:
877: .  subsnes - the subsolver context

879:    Level: intermediate

881: .seealso: `SNESNASM`, `SNESNASMGetNumber()`
882: @*/
883: PetscErrorCode SNESNASMGetSNES(SNES snes, PetscInt i, SNES *subsnes)
884: {
885:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

888:   *subsnes = nasm->subsnes[i];
889:   return 0;
890: }

892: /*@
893:    SNESNASMGetNumber - Gets number of subsolvers

895:    Not collective

897:    Input Parameters:
898: .  snes - the SNES context

900:    Output Parameters:
901: .  n - the number of subsolvers

903:    Level: intermediate

905: .seealso: `SNESNASM`, `SNESNASMGetSNES()`
906: @*/
907: PetscErrorCode SNESNASMGetNumber(SNES snes, PetscInt *n)
908: {
909:   SNES_NASM *nasm = (SNES_NASM *)snes->data;

911:   *n = nasm->n;
912:   return 0;
913: }

915: /*@
916:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

918:    Collective

920:    Input Parameters:
921: +  snes - the SNES context
922: -  weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)

924:    Level: intermediate

926: .seealso: `SNESNASM`
927: @*/
928: PetscErrorCode SNESNASMSetWeight(SNES snes, Vec weight)
929: {
930:   SNES_NASM *nasm = (SNES_NASM *)snes->data;


933:   VecDestroy(&nasm->weight);
934:   nasm->weight_set = PETSC_TRUE;
935:   nasm->weight     = weight;
936:   PetscObjectReference((PetscObject)nasm->weight);

938:   return 0;
939: }