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;

 55:   PetscFunctionBegin;
 56:   PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
 57:   if (snes->normschedule == SNES_NORM_ALWAYS) PetscCall(SNESSetInitialFunction(next->snes, F));
 58:   PetscCall(SNESSolve(next->snes, B, X));
 59:   PetscCall(SNESGetConvergedReason(next->snes, &reason));
 60:   if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 61:     jac->innerFailures++;
 62:     if (jac->innerFailures >= snes->maxFailures) {
 63:       snes->reason = SNES_DIVERGED_INNER;
 64:       PetscFunctionReturn(PETSC_SUCCESS);
 65:     }
 66:   }

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

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

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

174: /* approximately solve the overdetermined system:

176:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
177:  \alpha_i                      = 1

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

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

196:   PetscFunctionBegin;
197:   PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");

199:   if (snes->normschedule == SNES_NORM_ALWAYS) {
200:     next = jac->head;
201:     PetscCall(SNESSetInitialFunction(next->snes, F));
202:     while (next->next) {
203:       next = next->next;
204:       PetscCall(SNESSetInitialFunction(next->snes, F));
205:     }
206:   }

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

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

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

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

251:   for (i = 0; i < jac->n; i++) {
252:     for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
253:   }

255:   for (i = 0; i < jac->n; i++) {
256:     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;
257:     jac->beta[i] = ftf - jac->g[i];
258:   }

260:   jac->info  = 0;
261:   jac->rcond = -1.;
262:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
263: #if defined(PETSC_USE_COMPLEX)
264:   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));
265: #else
266:   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));
267: #endif
268:   PetscCall(PetscFPTrapPop());
269:   PetscCheck(jac->info >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "Bad argument to GELSS");
270:   PetscCheck(jac->info <= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD failed to converge");
271:   tot   = 0.;
272:   total = 0.;
273:   for (i = 0; i < jac->n; i++) {
274:     PetscCheck(!snes->errorifnotconverged || !PetscIsInfOrNanScalar(jac->beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
275:     PetscCall(PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i])));
276:     tot += jac->beta[i];
277:     total += PetscAbsScalar(jac->beta[i]);
278:   }
279:   PetscCall(VecScale(X, (1. - tot)));
280:   PetscCall(VecMAXPY(X, jac->n, jac->beta, Xes));
281:   PetscCall(SNESComputeFunction(snes, X, F));

283:   if (snes->xl && snes->xu) {
284:     PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
285:   } else {
286:     PetscCall(VecNorm(F, NORM_2, fnorm));
287:   }

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

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

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

316:   PetscFunctionBegin;
317:   PetscCall(SNESGetDM(snes, &dm));

319:   if (snes->ops->computevariablebounds) {
320:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
321:     if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
322:     if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
323:     PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
324:   }

326:   while (next) {
327:     n++;
328:     PetscCall(SNESSetDM(next->snes, dm));
329:     PetscCall(SNESSetApplicationContext(next->snes, snes->user));
330:     if (snes->xl && snes->xu) {
331:       if (snes->ops->computevariablebounds) {
332:         PetscCall(SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds));
333:       } else {
334:         PetscCall(SNESVISetVariableBounds(next->snes, snes->xl, snes->xu));
335:       }
336:     }

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

361:     PetscCall(PetscMalloc1(jac->n * jac->n, &jac->h));
362:     PetscCall(PetscMalloc1(jac->n, &jac->beta));
363:     PetscCall(PetscMalloc1(jac->n, &jac->s));
364:     PetscCall(PetscMalloc1(jac->n, &jac->g));
365:     jac->lwork = 12 * jac->n;
366: #if defined(PETSC_USE_COMPLEX)
367:     PetscCall(PetscMalloc1(jac->lwork, &jac->rwork));
368: #endif
369:     PetscCall(PetscMalloc1(jac->lwork, &jac->work));
370:   }

372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: static PetscErrorCode SNESReset_Composite(SNES snes)
376: {
377:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
378:   SNES_CompositeLink next = jac->head;

380:   PetscFunctionBegin;
381:   while (next) {
382:     PetscCall(SNESReset(next->snes));
383:     next = next->next;
384:   }
385:   PetscCall(VecDestroy(&jac->Xorig));
386:   if (jac->Xes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Xes));
387:   if (jac->Fes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Fes));
388:   PetscCall(PetscFree(jac->fnorms));
389:   PetscCall(PetscFree(jac->h));
390:   PetscCall(PetscFree(jac->s));
391:   PetscCall(PetscFree(jac->g));
392:   PetscCall(PetscFree(jac->beta));
393:   PetscCall(PetscFree(jac->work));
394:   PetscCall(PetscFree(jac->rwork));
395:   PetscFunctionReturn(PETSC_SUCCESS);
396: }

398: static PetscErrorCode SNESDestroy_Composite(SNES snes)
399: {
400:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
401:   SNES_CompositeLink next = jac->head, next_tmp;

403:   PetscFunctionBegin;
404:   PetscCall(SNESReset_Composite(snes));
405:   while (next) {
406:     PetscCall(SNESDestroy(&next->snes));
407:     next_tmp = next;
408:     next     = next->next;
409:     PetscCall(PetscFree(next_tmp));
410:   }
411:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL));
412:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL));
413:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL));
414:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL));
415:   PetscCall(PetscFree(snes->data));
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems *PetscOptionsObject)
420: {
421:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
422:   PetscInt           nmax = 8, i;
423:   SNES_CompositeLink next;
424:   char              *sneses[8];
425:   PetscReal          dmps[8];
426:   PetscBool          flg;

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

447:   next = jac->head;
448:   while (next) {
449:     PetscCall(SNESSetFromOptions(next->snes));
450:     next = next->next;
451:   }
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

455: static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
456: {
457:   SNES_Composite    *jac  = (SNES_Composite *)snes->data;
458:   SNES_CompositeLink next = jac->head;
459:   PetscBool          iascii;

461:   PetscFunctionBegin;
462:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
463:   if (iascii) {
464:     PetscCall(PetscViewerASCIIPrintf(viewer, "  type - %s\n", SNESCompositeTypes[jac->type]));
465:     PetscCall(PetscViewerASCIIPrintf(viewer, "  SNESes on composite preconditioner follow\n"));
466:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
467:   }
468:   if (iascii) PetscCall(PetscViewerASCIIPushTab(viewer));
469:   while (next) {
470:     PetscCall(SNESView(next->snes, viewer));
471:     next = next->next;
472:   }
473:   if (iascii) {
474:     PetscCall(PetscViewerASCIIPopTab(viewer));
475:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
476:   }
477:   PetscFunctionReturn(PETSC_SUCCESS);
478: }

480: static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
481: {
482:   SNES_Composite *jac = (SNES_Composite *)snes->data;

484:   PetscFunctionBegin;
485:   jac->type = type;
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
490: {
491:   SNES_Composite    *jac;
492:   SNES_CompositeLink next, ilink;
493:   PetscInt           cnt = 0;
494:   const char        *prefix;
495:   char               newprefix[20];
496:   DM                 dm;

498:   PetscFunctionBegin;
499:   PetscCall(PetscNew(&ilink));
500:   ilink->next = NULL;
501:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes));
502:   PetscCall(SNESGetDM(snes, &dm));
503:   PetscCall(SNESSetDM(ilink->snes, dm));
504:   PetscCall(SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs));
505:   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes));
506:   jac  = (SNES_Composite *)snes->data;
507:   next = jac->head;
508:   if (!next) {
509:     jac->head       = ilink;
510:     ilink->previous = NULL;
511:   } else {
512:     cnt++;
513:     while (next->next) {
514:       next = next->next;
515:       cnt++;
516:     }
517:     next->next      = ilink;
518:     ilink->previous = next;
519:   }
520:   PetscCall(SNESGetOptionsPrefix(snes, &prefix));
521:   PetscCall(SNESSetOptionsPrefix(ilink->snes, prefix));
522:   PetscCall(PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%d_", (int)cnt));
523:   PetscCall(SNESAppendOptionsPrefix(ilink->snes, newprefix));
524:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1));
525:   PetscCall(SNESSetType(ilink->snes, type));
526:   PetscCall(SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY));

528:   ilink->dmp = 1.0;
529:   jac->nsnes++;
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
534: {
535:   SNES_Composite    *jac;
536:   SNES_CompositeLink next;
537:   PetscInt           i;

539:   PetscFunctionBegin;
540:   jac  = (SNES_Composite *)snes->data;
541:   next = jac->head;
542:   for (i = 0; i < n; i++) {
543:     PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
544:     next = next->next;
545:   }
546:   *subsnes = next->snes;
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }

550: /*@C
551:    SNESCompositeSetType - Sets the type of composite preconditioner.

553:    Logically Collective

555:    Input Parameters:
556: +  snes - the preconditioner context
557: -  type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`

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

562:    Level: Developer

564: .seealso: `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`
565: @*/
566: PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
567: {
568:   PetscFunctionBegin;
571:   PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
572:   PetscFunctionReturn(PETSC_SUCCESS);
573: }

575: /*@C
576:    SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`

578:    Collective

580:    Input Parameters:
581: +  snes - the snes context of type `SNESCOMPOSITE`
582: -  type - the type of the new solver

584:    Level: Developer

586: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
587: @*/
588: PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
589: {
590:   PetscFunctionBegin;
592:   PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

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

599:    Not Collective

601:    Input Parameters:
602: +  snes - the snes context
603: -  n - the number of the composed snes requested

605:    Output Parameter:
606: .  subsnes - the `SNES` requested

608:    Level: Developer

610: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`
611: @*/
612: PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
613: {
614:   PetscFunctionBegin;
617:   PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
618:   PetscFunctionReturn(PETSC_SUCCESS);
619: }

621: /*@
622:    SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`

624:    Logically Collective

626:    Input Parameter:
627: .   snes - the snes context

629:    Output Parameter:
630: .   n - the number of subsolvers

632:    Level: Developer

634: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
635: @*/
636: PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
637: {
638:   SNES_Composite    *jac;
639:   SNES_CompositeLink next;

641:   PetscFunctionBegin;
642:   jac  = (SNES_Composite *)snes->data;
643:   next = jac->head;

645:   *n = 0;
646:   while (next) {
647:     *n   = *n + 1;
648:     next = next->next;
649:   }
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
654: {
655:   SNES_Composite    *jac;
656:   SNES_CompositeLink next;
657:   PetscInt           i;

659:   PetscFunctionBegin;
660:   jac  = (SNES_Composite *)snes->data;
661:   next = jac->head;
662:   for (i = 0; i < n; i++) {
663:     PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
664:     next = next->next;
665:   }
666:   next->dmp = dmp;
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }

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

673:    Not Collective

675:    Input Parameters:
676: +  snes - the snes context
677: .  n - the number of the sub-`SNES` object requested
678: -  dmp - the damping

680:    Level: intermediate

682: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
683:           `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
684: @*/
685: PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
686: {
687:   PetscFunctionBegin;
689:   PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
690:   PetscFunctionReturn(PETSC_SUCCESS);
691: }

693: static PetscErrorCode SNESSolve_Composite(SNES snes)
694: {
695:   Vec              F, X, B, Y;
696:   PetscInt         i;
697:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
698:   SNESNormSchedule normtype;
699:   SNES_Composite  *comp = (SNES_Composite *)snes->data;

701:   PetscFunctionBegin;
702:   X = snes->vec_sol;
703:   F = snes->vec_func;
704:   B = snes->vec_rhs;
705:   Y = snes->vec_sol_update;

707:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
708:   snes->iter          = 0;
709:   snes->norm          = 0.;
710:   comp->innerFailures = 0;
711:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
712:   snes->reason = SNES_CONVERGED_ITERATING;
713:   PetscCall(SNESGetNormSchedule(snes, &normtype));
714:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
715:     if (!snes->vec_func_init_set) {
716:       PetscCall(SNESComputeFunction(snes, X, F));
717:     } else snes->vec_func_init_set = PETSC_FALSE;

719:     if (snes->xl && snes->xu) {
720:       PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
721:     } else {
722:       PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
723:     }
724:     SNESCheckFunctionNorm(snes, fnorm);
725:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
726:     snes->iter = 0;
727:     snes->norm = fnorm;
728:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
729:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
730:     PetscCall(SNESMonitor(snes, 0, snes->norm));

732:     /* test convergence */
733:     PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
734:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
735:   } else {
736:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
737:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
738:     PetscCall(SNESMonitor(snes, 0, snes->norm));
739:   }

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

745:     /* Copy the state before modification by application of the composite solver;
746:        we will subtract the new state after application */
747:     PetscCall(VecCopy(X, Y));

749:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
750:       PetscCall(SNESCompositeApply_Additive(snes, X, B, F, &fnorm));
751:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
752:       PetscCall(SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm));
753:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
754:       PetscCall(SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm));
755:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
756:     if (snes->reason < 0) break;

758:     /* Compute the solution update for convergence testing */
759:     PetscCall(VecAYPX(Y, -1.0, X));

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

764:       if (snes->xl && snes->xu) {
765:         PetscCall(VecNormBegin(X, NORM_2, &xnorm));
766:         PetscCall(VecNormBegin(Y, NORM_2, &snorm));
767:         PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
768:         PetscCall(VecNormEnd(X, NORM_2, &xnorm));
769:         PetscCall(VecNormEnd(Y, NORM_2, &snorm));
770:       } else {
771:         PetscCall(VecNormBegin(F, NORM_2, &fnorm));
772:         PetscCall(VecNormBegin(X, NORM_2, &xnorm));
773:         PetscCall(VecNormBegin(Y, NORM_2, &snorm));

775:         PetscCall(VecNormEnd(F, NORM_2, &fnorm));
776:         PetscCall(VecNormEnd(X, NORM_2, &xnorm));
777:         PetscCall(VecNormEnd(Y, NORM_2, &snorm));
778:       }
779:       SNESCheckFunctionNorm(snes, fnorm);
780:     } else if (normtype == SNES_NORM_ALWAYS) {
781:       PetscCall(VecNormBegin(X, NORM_2, &xnorm));
782:       PetscCall(VecNormBegin(Y, NORM_2, &snorm));
783:       PetscCall(VecNormEnd(X, NORM_2, &xnorm));
784:       PetscCall(VecNormEnd(Y, NORM_2, &snorm));
785:     }
786:     /* Monitor convergence */
787:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
788:     snes->iter  = i + 1;
789:     snes->norm  = fnorm;
790:     snes->xnorm = xnorm;
791:     snes->ynorm = snorm;
792:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
793:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
794:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
795:     /* Test for convergence */
796:     if (normtype == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, snorm, fnorm, &snes->reason, snes->cnvP);
797:     if (snes->reason) break;
798:   }
799:   if (normtype == SNES_NORM_ALWAYS) {
800:     if (i == snes->max_its) {
801:       PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its));
802:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
803:     }
804:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

808: /*MC
809:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

811:    Options Database Keys:
812: +  -snes_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
813: -  -snes_composite_sneses - <snes0,snes1,...> list of `SNES` to compose

815:    Level: intermediate

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

821: .seealso: `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
822:           `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`,
823:           `SNESCompositeSetDamping()`
824: M*/

826: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
827: {
828:   SNES_Composite *jac;

830:   PetscFunctionBegin;
831:   PetscCall(PetscNew(&jac));

833:   snes->ops->solve          = SNESSolve_Composite;
834:   snes->ops->setup          = SNESSetUp_Composite;
835:   snes->ops->reset          = SNESReset_Composite;
836:   snes->ops->destroy        = SNESDestroy_Composite;
837:   snes->ops->setfromoptions = SNESSetFromOptions_Composite;
838:   snes->ops->view           = SNESView_Composite;

840:   snes->usesksp = PETSC_FALSE;

842:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

844:   snes->data  = (void *)jac;
845:   jac->type   = SNES_COMPOSITE_ADDITIVEOPTIMAL;
846:   jac->Fes    = NULL;
847:   jac->Xes    = NULL;
848:   jac->fnorms = NULL;
849:   jac->nsnes  = 0;
850:   jac->head   = NULL;
851:   jac->stol   = 0.1;
852:   jac->rtol   = 1.1;

854:   jac->h     = NULL;
855:   jac->s     = NULL;
856:   jac->beta  = NULL;
857:   jac->work  = NULL;
858:   jac->rwork = NULL;

860:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite));
861:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite));
862:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite));
863:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite));
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }