Actual source code: snescomposite.c


  2: /*
  3:       Defines a SNES that can consist of a collection of SNESes
  4: */
  5: #include <petsc/private/snesimpl.h>
  6: #include <petscblaslapack.h>

  8: const char *const SNESCompositeTypes[] = {"ADDITIVE", "MULTIPLICATIVE", "ADDITIVEOPTIMAL", "SNESCompositeType", "SNES_COMPOSITE", NULL};

 10: typedef struct _SNES_CompositeLink *SNES_CompositeLink;
 11: struct _SNES_CompositeLink {
 12:   SNES               snes;
 13:   PetscReal          dmp;
 14:   Vec                X;
 15:   SNES_CompositeLink next;
 16:   SNES_CompositeLink previous;
 17: };

 19: typedef struct {
 20:   SNES_CompositeLink head;
 21:   PetscInt           nsnes;
 22:   SNESCompositeType  type;
 23:   Vec                Xorig;
 24:   PetscInt           innerFailures; /* the number of inner failures we've seen */

 26:   /* context for ADDITIVEOPTIMAL */
 27:   Vec         *Xes, *Fes; /* solution and residual vectors for the subsolvers */
 28:   PetscReal   *fnorms;    /* norms of the solutions */
 29:   PetscScalar *h;         /* the matrix formed as q_ij = (rdot_i, rdot_j) */
 30:   PetscScalar *g;         /* the dotproducts of the previous function with the candidate functions */
 31:   PetscBLASInt n;         /* matrix dimension -- nsnes */
 32:   PetscBLASInt nrhs;      /* the number of right hand sides */
 33:   PetscBLASInt lda;       /* the padded matrix dimension */
 34:   PetscBLASInt ldb;       /* the padded vector dimension */
 35:   PetscReal   *s;         /* the singular values */
 36:   PetscScalar *beta;      /* the RHS and combination */
 37:   PetscReal    rcond;     /* the exit condition */
 38:   PetscBLASInt rank;      /* the effective rank */
 39:   PetscScalar *work;      /* the work vector */
 40:   PetscReal   *rwork;     /* the real work vector used for complex */
 41:   PetscBLASInt lwork;     /* the size of the work vector */
 42:   PetscBLASInt info;      /* the output condition */

 44:   PetscReal rtol; /* restart tolerance for accepting the combination */
 45:   PetscReal stol; /* restart tolerance for the combination */
 46: } SNES_Composite;

 48: static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
 49: {
 50:   SNES_Composite     *jac  = (SNES_Composite *)snes->data;
 51:   SNES_CompositeLink  next = jac->head;
 52:   Vec                 FSub;
 53:   SNESConvergedReason reason;

 56:   if (snes->normschedule == SNES_NORM_ALWAYS) SNESSetInitialFunction(next->snes, F);
 57:   SNESSolve(next->snes, B, X);
 58:   SNESGetConvergedReason(next->snes, &reason);
 59:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 60:     jac->innerFailures++;
 61:     if (jac->innerFailures >= snes->maxFailures) {
 62:       snes->reason = SNES_DIVERGED_INNER;
 63:       return 0;
 64:     }
 65:   }

 67:   while (next->next) {
 68:     /* only copy the function over in the case where the functions correspond */
 69:     if (next->snes->npcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
 70:       SNESGetFunction(next->snes, &FSub, NULL, NULL);
 71:       next = next->next;
 72:       SNESSetInitialFunction(next->snes, FSub);
 73:     } else {
 74:       next = next->next;
 75:     }
 76:     SNESSolve(next->snes, B, X);
 77:     SNESGetConvergedReason(next->snes, &reason);
 78:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 79:       jac->innerFailures++;
 80:       if (jac->innerFailures >= snes->maxFailures) {
 81:         snes->reason = SNES_DIVERGED_INNER;
 82:         return 0;
 83:       }
 84:     }
 85:   }
 86:   if (next->snes->npcside == PC_RIGHT) {
 87:     SNESGetFunction(next->snes, &FSub, NULL, NULL);
 88:     VecCopy(FSub, F);
 89:     if (fnorm) {
 90:       if (snes->xl && snes->xu) {
 91:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
 92:       } else {
 93:         VecNorm(F, NORM_2, fnorm);
 94:       }
 95:       SNESCheckFunctionNorm(snes, *fnorm);
 96:     }
 97:   } else if (snes->normschedule == SNES_NORM_ALWAYS) {
 98:     SNESComputeFunction(snes, X, F);
 99:     if (fnorm) {
100:       if (snes->xl && snes->xu) {
101:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
102:       } else {
103:         VecNorm(F, NORM_2, fnorm);
104:       }
105:       SNESCheckFunctionNorm(snes, *fnorm);
106:     }
107:   }
108:   return 0;
109: }

111: static PetscErrorCode SNESCompositeApply_Additive(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
112: {
113:   SNES_Composite     *jac  = (SNES_Composite *)snes->data;
114:   SNES_CompositeLink  next = jac->head;
115:   Vec                 Y, Xorig;
116:   SNESConvergedReason reason;

118:   Y = snes->vec_sol_update;
119:   if (!jac->Xorig) VecDuplicate(X, &jac->Xorig);
120:   Xorig = jac->Xorig;
121:   VecCopy(X, Xorig);
123:   if (snes->normschedule == SNES_NORM_ALWAYS) {
124:     SNESSetInitialFunction(next->snes, F);
125:     while (next->next) {
126:       next = next->next;
127:       SNESSetInitialFunction(next->snes, F);
128:     }
129:   }
130:   next = jac->head;
131:   VecCopy(Xorig, Y);
132:   SNESSolve(next->snes, B, Y);
133:   SNESGetConvergedReason(next->snes, &reason);
134:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
135:     jac->innerFailures++;
136:     if (jac->innerFailures >= snes->maxFailures) {
137:       snes->reason = SNES_DIVERGED_INNER;
138:       return 0;
139:     }
140:   }
141:   VecAXPY(Y, -1.0, Xorig);
142:   VecAXPY(X, next->dmp, Y);
143:   while (next->next) {
144:     next = next->next;
145:     VecCopy(Xorig, Y);
146:     SNESSolve(next->snes, B, Y);
147:     SNESGetConvergedReason(next->snes, &reason);
148:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
149:       jac->innerFailures++;
150:       if (jac->innerFailures >= snes->maxFailures) {
151:         snes->reason = SNES_DIVERGED_INNER;
152:         return 0;
153:       }
154:     }
155:     VecAXPY(Y, -1.0, Xorig);
156:     VecAXPY(X, next->dmp, Y);
157:   }
158:   if (snes->normschedule == SNES_NORM_ALWAYS) {
159:     SNESComputeFunction(snes, X, F);
160:     if (fnorm) {
161:       if (snes->xl && snes->xu) {
162:         SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
163:       } else {
164:         VecNorm(F, NORM_2, fnorm);
165:       }
166:       SNESCheckFunctionNorm(snes, *fnorm);
167:     }
168:   }
169:   return 0;
170: }

172: /* approximately solve the overdetermined system:

174:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
175:  \alpha_i                      = 1

177:  Which minimizes the L2 norm of the linearization of:
178:  ||F(\sum_i \alpha_i*x_i)||^2

180:  With the constraint that \sum_i\alpha_i = 1
181:  Where x_i is the solution from the ith subsolver.
182:  */
183: static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
184: {
185:   SNES_Composite     *jac  = (SNES_Composite *)snes->data;
186:   SNES_CompositeLink  next = jac->head;
187:   Vec                *Xes = jac->Xes, *Fes = jac->Fes;
188:   PetscInt            i, j;
189:   PetscScalar         tot, total, ftf;
190:   PetscReal           min_fnorm;
191:   PetscInt            min_i;
192:   SNESConvergedReason reason;


196:   if (snes->normschedule == SNES_NORM_ALWAYS) {
197:     next = jac->head;
198:     SNESSetInitialFunction(next->snes, F);
199:     while (next->next) {
200:       next = next->next;
201:       SNESSetInitialFunction(next->snes, F);
202:     }
203:   }

205:   next = jac->head;
206:   i    = 0;
207:   VecCopy(X, Xes[i]);
208:   SNESSolve(next->snes, B, Xes[i]);
209:   SNESGetConvergedReason(next->snes, &reason);
210:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
211:     jac->innerFailures++;
212:     if (jac->innerFailures >= snes->maxFailures) {
213:       snes->reason = SNES_DIVERGED_INNER;
214:       return 0;
215:     }
216:   }
217:   while (next->next) {
218:     i++;
219:     next = next->next;
220:     VecCopy(X, Xes[i]);
221:     SNESSolve(next->snes, B, Xes[i]);
222:     SNESGetConvergedReason(next->snes, &reason);
223:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
224:       jac->innerFailures++;
225:       if (jac->innerFailures >= snes->maxFailures) {
226:         snes->reason = SNES_DIVERGED_INNER;
227:         return 0;
228:       }
229:     }
230:   }

232:   /* all the solutions are collected; combine optimally */
233:   for (i = 0; i < jac->n; i++) {
234:     for (j = 0; j < i + 1; j++) VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]);
235:     VecDotBegin(Fes[i], F, &jac->g[i]);
236:   }

238:   for (i = 0; i < jac->n; i++) {
239:     for (j = 0; j < i + 1; j++) {
240:       VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]);
241:       if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
242:     }
243:     VecDotEnd(Fes[i], F, &jac->g[i]);
244:   }

246:   ftf = (*fnorm) * (*fnorm);

248:   for (i = 0; i < jac->n; i++) {
249:     for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
250:   }

252:   for (i = 0; i < jac->n; i++) {
253:     for (j = 0; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[i + j * jac->n] - jac->g[j] - jac->g[i] + ftf;
254:     jac->beta[i] = ftf - jac->g[i];
255:   }

257:   jac->info  = 0;
258:   jac->rcond = -1.;
259:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
260: #if defined(PETSC_USE_COMPLEX)
261:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, jac->rwork, &jac->info));
262: #else
263:   PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, &jac->info));
264: #endif
265:   PetscFPTrapPop();
268:   tot   = 0.;
269:   total = 0.;
270:   for (i = 0; i < jac->n; i++) {
272:     PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i]));
273:     tot += jac->beta[i];
274:     total += PetscAbsScalar(jac->beta[i]);
275:   }
276:   VecScale(X, (1. - tot));
277:   VecMAXPY(X, jac->n, jac->beta, Xes);
278:   SNESComputeFunction(snes, X, F);

280:   if (snes->xl && snes->xu) {
281:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
282:   } else {
283:     VecNorm(F, NORM_2, fnorm);
284:   }

286:   /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
287:   min_fnorm = jac->fnorms[0];
288:   min_i     = 0;
289:   for (i = 0; i < jac->n; i++) {
290:     if (jac->fnorms[i] < min_fnorm) {
291:       min_fnorm = jac->fnorms[i];
292:       min_i     = i;
293:     }
294:   }

296:   /* stagnation or divergence restart to the solution of the solver that failed the least */
297:   if (PetscRealPart(total) < jac->stol || min_fnorm * jac->rtol < *fnorm) {
298:     VecCopy(jac->Xes[min_i], X);
299:     VecCopy(jac->Fes[min_i], F);
300:     *fnorm = min_fnorm;
301:   }
302:   return 0;
303: }

305: static PetscErrorCode SNESSetUp_Composite(SNES snes)
306: {
307:   DM                 dm;
308:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
309:   SNES_CompositeLink next = jac->head;
310:   PetscInt           n    = 0, i;
311:   Vec                F;

313:   SNESGetDM(snes, &dm);

315:   if (snes->ops->computevariablebounds) {
316:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
317:     if (!snes->xl) VecDuplicate(snes->vec_sol, &snes->xl);
318:     if (!snes->xu) VecDuplicate(snes->vec_sol, &snes->xu);
319:     PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
320:   }

322:   while (next) {
323:     n++;
324:     SNESSetDM(next->snes, dm);
325:     SNESSetApplicationContext(next->snes, snes->user);
326:     if (snes->xl && snes->xu) {
327:       if (snes->ops->computevariablebounds) {
328:         SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
329:       } else {
330:         SNESVISetVariableBounds(next->snes, snes->xl, snes->xu);
331:       }
332:     }

334:     next = next->next;
335:   }
336:   jac->nsnes = n;
337:   SNESGetFunction(snes, &F, NULL, NULL);
338:   if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
339:     VecDuplicateVecs(F, jac->nsnes, &jac->Xes);
340:     PetscMalloc1(n, &jac->Fes);
341:     PetscMalloc1(n, &jac->fnorms);
342:     next = jac->head;
343:     i    = 0;
344:     while (next) {
345:       SNESGetFunction(next->snes, &F, NULL, NULL);
346:       jac->Fes[i] = F;
347:       PetscObjectReference((PetscObject)F);
348:       next = next->next;
349:       i++;
350:     }
351:     /* allocate the subspace direct solve area */
352:     jac->nrhs = 1;
353:     jac->lda  = jac->nsnes;
354:     jac->ldb  = jac->nsnes;
355:     jac->n    = jac->nsnes;

357:     PetscMalloc1(jac->n * jac->n, &jac->h);
358:     PetscMalloc1(jac->n, &jac->beta);
359:     PetscMalloc1(jac->n, &jac->s);
360:     PetscMalloc1(jac->n, &jac->g);
361:     jac->lwork = 12 * jac->n;
362: #if defined(PETSC_USE_COMPLEX)
363:     PetscMalloc1(jac->lwork, &jac->rwork);
364: #endif
365:     PetscMalloc1(jac->lwork, &jac->work);
366:   }

368:   return 0;
369: }

371: static PetscErrorCode SNESReset_Composite(SNES snes)
372: {
373:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
374:   SNES_CompositeLink next = jac->head;

376:   while (next) {
377:     SNESReset(next->snes);
378:     next = next->next;
379:   }
380:   VecDestroy(&jac->Xorig);
381:   if (jac->Xes) VecDestroyVecs(jac->nsnes, &jac->Xes);
382:   if (jac->Fes) VecDestroyVecs(jac->nsnes, &jac->Fes);
383:   PetscFree(jac->fnorms);
384:   PetscFree(jac->h);
385:   PetscFree(jac->s);
386:   PetscFree(jac->g);
387:   PetscFree(jac->beta);
388:   PetscFree(jac->work);
389:   PetscFree(jac->rwork);
390:   return 0;
391: }

393: static PetscErrorCode SNESDestroy_Composite(SNES snes)
394: {
395:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
396:   SNES_CompositeLink next = jac->head, next_tmp;

398:   SNESReset_Composite(snes);
399:   while (next) {
400:     SNESDestroy(&next->snes);
401:     next_tmp = next;
402:     next     = next->next;
403:     PetscFree(next_tmp);
404:   }
405:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL);
406:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL);
407:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL);
408:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL);
409:   PetscFree(snes->data);
410:   return 0;
411: }

413: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems *PetscOptionsObject)
414: {
415:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
416:   PetscInt           nmax = 8, i;
417:   SNES_CompositeLink next;
418:   char              *sneses[8];
419:   PetscReal          dmps[8];
420:   PetscBool          flg;

422:   PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
423:   PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg);
424:   if (flg) SNESCompositeSetType(snes, jac->type);
425:   PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg);
426:   if (flg) {
427:     for (i = 0; i < nmax; i++) {
428:       SNESCompositeAddSNES(snes, sneses[i]);
429:       PetscFree(sneses[i]); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
430:     }
431:   }
432:   PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg);
433:   if (flg) {
434:     for (i = 0; i < nmax; i++) SNESCompositeSetDamping(snes, i, dmps[i]);
435:   }
436:   PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL);
437:   PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL);
438:   PetscOptionsHeadEnd();

440:   next = jac->head;
441:   while (next) {
442:     SNESSetFromOptions(next->snes);
443:     next = next->next;
444:   }
445:   return 0;
446: }

448: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
449: {
450:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
451:   SNES_CompositeLink next = jac->head;
452:   PetscBool          iascii;

454:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
455:   if (iascii) {
456:     PetscViewerASCIIPrintf(viewer, "  type - %s\n", SNESCompositeTypes[jac->type]);
457:     PetscViewerASCIIPrintf(viewer, "  SNESes on composite preconditioner follow\n");
458:     PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n");
459:   }
460:   if (iascii) PetscViewerASCIIPushTab(viewer);
461:   while (next) {
462:     SNESView(next->snes, viewer);
463:     next = next->next;
464:   }
465:   if (iascii) {
466:     PetscViewerASCIIPopTab(viewer);
467:     PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n");
468:   }
469:   return 0;
470: }

472: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
473: {
474:   SNES_Composite *jac = (SNES_Composite *)snes->data;

476:   jac->type = type;
477:   return 0;
478: }

480: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
481: {
482:   SNES_Composite    *jac;
483:   SNES_CompositeLink next, ilink;
484:   PetscInt           cnt = 0;
485:   const char        *prefix;
486:   char               newprefix[20];
487:   DM                 dm;

489:   PetscNew(&ilink);
490:   ilink->next = NULL;
491:   SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes);
492:   SNESGetDM(snes, &dm);
493:   SNESSetDM(ilink->snes, dm);
494:   SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs);
495:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes);
496:   jac  = (SNES_Composite *)snes->data;
497:   next = jac->head;
498:   if (!next) {
499:     jac->head       = ilink;
500:     ilink->previous = NULL;
501:   } else {
502:     cnt++;
503:     while (next->next) {
504:       next = next->next;
505:       cnt++;
506:     }
507:     next->next      = ilink;
508:     ilink->previous = next;
509:   }
510:   SNESGetOptionsPrefix(snes, &prefix);
511:   SNESSetOptionsPrefix(ilink->snes, prefix);
512:   PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%d_", (int)cnt);
513:   SNESAppendOptionsPrefix(ilink->snes, newprefix);
514:   PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1);
515:   SNESSetType(ilink->snes, type);
516:   SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);

518:   ilink->dmp = 1.0;
519:   jac->nsnes++;
520:   return 0;
521: }

523: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
524: {
525:   SNES_Composite    *jac;
526:   SNES_CompositeLink next;
527:   PetscInt           i;

529:   jac  = (SNES_Composite *)snes->data;
530:   next = jac->head;
531:   for (i = 0; i < n; i++) {
533:     next = next->next;
534:   }
535:   *subsnes = next->snes;
536:   return 0;
537: }

539: /*@C
540:    SNESCompositeSetType - Sets the type of composite preconditioner.

542:    Logically Collective

544:    Input Parameters:
545: +  snes - the preconditioner context
546: -  type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`

548:    Options Database Key:
549: .  -snes_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

551:    Level: Developer

553: .seealso: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`
554: @*/
555: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
556: {
559:   PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
560:   return 0;
561: }

563: /*@C
564:    SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`

566:    Collective

568:    Input Parameters:
569: +  snes - the snes context of type `SNESCOMPOSITE`
570: -  type - the type of the new solver

572:    Level: Developer

574: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
575: @*/
576: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
577: {
579:   PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
580:   return 0;
581: }

583: /*@
584:    SNESCompositeGetSNES - Gets one of the `SNES` objects in the composite  `SNESCOMPOSITE`

586:    Not Collective

588:    Input Parameters:
589: +  snes - the snes context
590: -  n - the number of the composed snes requested

592:    Output Parameters:
593: .  subsnes - the `SNES` requested

595:    Level: Developer

597: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`
598: @*/
599: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
600: {
603:   PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
604:   return 0;
605: }

607: /*@
608:    SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`

610:    Logically Collective

612:    Input Parameter:
613:    snes - the snes context

615:    Output Parameter:
616:    n - the number of subsolvers

618:    Level: Developer

620: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
621: @*/
622: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
623: {
624:   SNES_Composite    *jac;
625:   SNES_CompositeLink next;

627:   jac  = (SNES_Composite *)snes->data;
628:   next = jac->head;

630:   *n = 0;
631:   while (next) {
632:     *n   = *n + 1;
633:     next = next->next;
634:   }
635:   return 0;
636: }

638: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
639: {
640:   SNES_Composite    *jac;
641:   SNES_CompositeLink next;
642:   PetscInt           i;

644:   jac  = (SNES_Composite *)snes->data;
645:   next = jac->head;
646:   for (i = 0; i < n; i++) {
648:     next = next->next;
649:   }
650:   next->dmp = dmp;
651:   return 0;
652: }

654: /*@
655:    SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` `SNESCOMPOSITE`

657:    Not Collective

659:    Input Parameters:
660: +  snes - the snes context
661: .  n - the number of the snes requested
662: -  dmp - the damping

664:    Level: intermediate

666: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
667:           `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
668: @*/
669: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
670: {
672:   PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
673:   return 0;
674: }

676: static PetscErrorCode SNESSolve_Composite(SNES snes)
677: {
678:   Vec              F, X, B, Y;
679:   PetscInt         i;
680:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
681:   SNESNormSchedule normtype;
682:   SNES_Composite  *comp = (SNES_Composite *)snes->data;

684:   X = snes->vec_sol;
685:   F = snes->vec_func;
686:   B = snes->vec_rhs;
687:   Y = snes->vec_sol_update;

689:   PetscObjectSAWsTakeAccess((PetscObject)snes);
690:   snes->iter          = 0;
691:   snes->norm          = 0.;
692:   comp->innerFailures = 0;
693:   PetscObjectSAWsGrantAccess((PetscObject)snes);
694:   snes->reason = SNES_CONVERGED_ITERATING;
695:   SNESGetNormSchedule(snes, &normtype);
696:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
697:     if (!snes->vec_func_init_set) {
698:       SNESComputeFunction(snes, X, F);
699:     } else snes->vec_func_init_set = PETSC_FALSE;

701:     if (snes->xl && snes->xu) {
702:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
703:     } else {
704:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
705:     }
706:     SNESCheckFunctionNorm(snes, fnorm);
707:     PetscObjectSAWsTakeAccess((PetscObject)snes);
708:     snes->iter = 0;
709:     snes->norm = fnorm;
710:     PetscObjectSAWsGrantAccess((PetscObject)snes);
711:     SNESLogConvergenceHistory(snes, snes->norm, 0);
712:     SNESMonitor(snes, 0, snes->norm);

714:     /* test convergence */
715:     PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
716:     if (snes->reason) return 0;
717:   } else {
718:     PetscObjectSAWsGrantAccess((PetscObject)snes);
719:     SNESLogConvergenceHistory(snes, snes->norm, 0);
720:     SNESMonitor(snes, 0, snes->norm);
721:   }

723:   for (i = 0; i < snes->max_its; i++) {
724:     /* Call general purpose update function */
725:     PetscTryTypeMethod(snes, update, snes->iter);

727:     /* Copy the state before modification by application of the composite solver;
728:        we will subtract the new state after application */
729:     VecCopy(X, Y);

731:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
732:       SNESCompositeApply_Additive(snes, X, B, F, &fnorm);
733:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
734:       SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm);
735:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
736:       SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm);
737:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
738:     if (snes->reason < 0) break;

740:     /* Compute the solution update for convergence testing */
741:     VecAYPX(Y, -1.0, X);

743:     if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
744:       SNESComputeFunction(snes, X, F);

746:       if (snes->xl && snes->xu) {
747:         VecNormBegin(X, NORM_2, &xnorm);
748:         VecNormBegin(Y, NORM_2, &snorm);
749:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
750:         VecNormEnd(X, NORM_2, &xnorm);
751:         VecNormEnd(Y, NORM_2, &snorm);
752:       } else {
753:         VecNormBegin(F, NORM_2, &fnorm);
754:         VecNormBegin(X, NORM_2, &xnorm);
755:         VecNormBegin(Y, NORM_2, &snorm);

757:         VecNormEnd(F, NORM_2, &fnorm);
758:         VecNormEnd(X, NORM_2, &xnorm);
759:         VecNormEnd(Y, NORM_2, &snorm);
760:       }
761:       SNESCheckFunctionNorm(snes, fnorm);
762:     } else if (normtype == SNES_NORM_ALWAYS) {
763:       VecNormBegin(X, NORM_2, &xnorm);
764:       VecNormBegin(Y, NORM_2, &snorm);
765:       VecNormEnd(X, NORM_2, &xnorm);
766:       VecNormEnd(Y, NORM_2, &snorm);
767:     }
768:     /* Monitor convergence */
769:     PetscObjectSAWsTakeAccess((PetscObject)snes);
770:     snes->iter  = i + 1;
771:     snes->norm  = fnorm;
772:     snes->xnorm = xnorm;
773:     snes->ynorm = snorm;
774:     PetscObjectSAWsGrantAccess((PetscObject)snes);
775:     SNESLogConvergenceHistory(snes, snes->norm, 0);
776:     SNESMonitor(snes, snes->iter, snes->norm);
777:     /* Test for convergence */
778:     if (normtype == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, snorm, fnorm, &snes->reason, snes->cnvP);
779:     if (snes->reason) break;
780:   }
781:   if (normtype == SNES_NORM_ALWAYS) {
782:     if (i == snes->max_its) {
783:       PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
784:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
785:     }
786:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
787:   return 0;
788: }

790: /*MC
791:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

793:    Options Database Keys:
794: +  -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
795: -  -snes_composite_sneses - <snes0,snes1,...> list of SNESes to compose

797:    Level: intermediate

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

803: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
804:           `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`,
805:           `SNESCompositeSetDamping()`
806: M*/

808: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
809: {
810:   SNES_Composite *jac;

812:   PetscNew(&jac);

814:   snes->ops->solve          = SNESSolve_Composite;
815:   snes->ops->setup          = SNESSetUp_Composite;
816:   snes->ops->reset          = SNESReset_Composite;
817:   snes->ops->destroy        = SNESDestroy_Composite;
818:   snes->ops->setfromoptions = SNESSetFromOptions_Composite;
819:   snes->ops->view           = SNESView_Composite;

821:   snes->usesksp = PETSC_FALSE;

823:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

825:   snes->data  = (void *)jac;
826:   jac->type   = SNES_COMPOSITE_ADDITIVEOPTIMAL;
827:   jac->Fes    = NULL;
828:   jac->Xes    = NULL;
829:   jac->fnorms = NULL;
830:   jac->nsnes  = 0;
831:   jac->head   = NULL;
832:   jac->stol   = 0.1;
833:   jac->rtol   = 1.1;

835:   jac->h     = NULL;
836:   jac->s     = NULL;
837:   jac->beta  = NULL;
838:   jac->work  = NULL;
839:   jac->rwork = NULL;

841:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite);
842:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite);
843:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite);
844:   PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite);
845:   return 0;
846: }