Actual source code: snescomposite.c

petsc-3.12.5 2020-03-29
Report Typos and Errors

  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",0};

 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:   PetscErrorCode      ierr;
 51:   SNES_Composite      *jac = (SNES_Composite*)snes->data;
 52:   SNES_CompositeLink  next = jac->head;
 53:   Vec                 FSub;
 54:   SNESConvergedReason reason;

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

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

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

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

178: /* approximately solve the overdetermined system:

180:  2*F(x_i)\cdot F(\x_j)\alpha_i = 0
181:  \alpha_i                      = 1

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

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

202:   if (!next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");

204:   if (snes->normschedule == SNES_NORM_ALWAYS) {
205:     next = jac->head;
206:     SNESSetInitialFunction(next->snes,F);
207:     while (next->next) {
208:       next = next->next;
209:       SNESSetInitialFunction(next->snes,F);
210:     }
211:   }

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

240:   /* all the solutions are collected; combine optimally */
241:   for (i=0;i<jac->n;i++) {
242:     for (j=0;j<i+1;j++) {
243:       VecDotBegin(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
244:     }
245:     VecDotBegin(Fes[i],F,&jac->g[i]);
246:   }

248:   for (i=0;i<jac->n;i++) {
249:     for (j=0;j<i+1;j++) {
250:       VecDotEnd(Fes[i],Fes[j],&jac->h[i + j*jac->n]);
251:       if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j*jac->n]));
252:     }
253:     VecDotEnd(Fes[i],F,&jac->g[i]);
254:   }

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

258:   for (i=0; i<jac->n; i++) {
259:     for (j=i+1;j<jac->n;j++) {
260:       jac->h[i + j*jac->n] = jac->h[j + i*jac->n];
261:     }
262:   }

264:   for (i=0; i<jac->n; i++) {
265:     for (j=0; j<jac->n; j++) {
266:       jac->h[i + j*jac->n] = jac->h[i + j*jac->n] - jac->g[j] - jac->g[i] + ftf;
267:     }
268:     jac->beta[i] = ftf - jac->g[i];
269:   }

271: #if defined(PETSC_MISSING_LAPACK_GELSS)
272:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"SNESCOMPOSITE with ADDITIVEOPTIMAL requires the LAPACK GELSS routine.");
273: #else
274:   jac->info  = 0;
275:   jac->rcond = -1.;
276:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
277: #if defined(PETSC_USE_COMPLEX)
278:   PetscStackCall("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));
279: #else
280:   PetscStackCall("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));
281: #endif
282:   PetscFPTrapPop();
283:   if (jac->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
284:   if (jac->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
285: #endif
286:   tot = 0.;
287:   total = 0.;
288:   for (i=0; i<jac->n; i++) {
289:     if (snes->errorifnotconverged && PetscIsInfOrNanScalar(jac->beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
290:     PetscInfo2(snes,"%D: %g\n",i,(double)PetscRealPart(jac->beta[i]));
291:     tot += jac->beta[i];
292:     total += PetscAbsScalar(jac->beta[i]);
293:   }
294:   VecScale(X,(1. - tot));
295:   VecMAXPY(X,jac->n,jac->beta,Xes);
296:   SNESComputeFunction(snes,X,F);

298:   if (snes->xl && snes->xu) {
299:     SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm);
300:   } else {
301:     VecNorm(F, NORM_2, fnorm);
302:   }

304:   /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
305:   min_fnorm = jac->fnorms[0];
306:   min_i     = 0;
307:   for (i=0; i<jac->n; i++) {
308:     if (jac->fnorms[i] < min_fnorm) {
309:       min_fnorm = jac->fnorms[i];
310:       min_i     = i;
311:     }
312:   }

314:   /* stagnation or divergence restart to the solution of the solver that failed the least */
315:   if (PetscRealPart(total) < jac->stol || min_fnorm*jac->rtol < *fnorm) {
316:     VecCopy(jac->Xes[min_i],X);
317:     VecCopy(jac->Fes[min_i],F);
318:     *fnorm = min_fnorm;
319:   }
320:   return(0);
321: }

323: static PetscErrorCode SNESSetUp_Composite(SNES snes)
324: {
325:   PetscErrorCode     ierr;
326:   DM                 dm;
327:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
328:   SNES_CompositeLink next = jac->head;
329:   PetscInt           n=0,i;
330:   Vec                F;

333:   SNESGetDM(snes,&dm);

335:   if (snes->ops->computevariablebounds) {
336:     /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
337:     if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
338:     if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
339:     (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
340:   }

342:   while (next) {
343:     n++;
344:     SNESSetDM(next->snes,dm);
345:     SNESSetApplicationContext(next->snes, snes->user);
346:     if (snes->xl && snes->xu) {
347:       if (snes->ops->computevariablebounds) {
348:         SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds);
349:       } else {
350:         SNESVISetVariableBounds(next->snes,snes->xl,snes->xu);
351:       }
352:     }

354:     next = next->next;
355:   }
356:   jac->nsnes = n;
357:   SNESGetFunction(snes,&F,NULL,NULL);
358:   if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
359:     VecDuplicateVecs(F,jac->nsnes,&jac->Xes);
360:     PetscMalloc1(n,&jac->Fes);
361:     PetscMalloc1(n,&jac->fnorms);
362:     next = jac->head;
363:     i = 0;
364:     while (next) {
365:       SNESGetFunction(next->snes,&F,NULL,NULL);
366:       jac->Fes[i] = F;
367:       PetscObjectReference((PetscObject)F);
368:       next = next->next;
369:       i++;
370:     }
371:     /* allocate the subspace direct solve area */
372:     jac->nrhs  = 1;
373:     jac->lda   = jac->nsnes;
374:     jac->ldb   = jac->nsnes;
375:     jac->n     = jac->nsnes;

377:     PetscMalloc1(jac->n*jac->n,&jac->h);
378:     PetscMalloc1(jac->n,&jac->beta);
379:     PetscMalloc1(jac->n,&jac->s);
380:     PetscMalloc1(jac->n,&jac->g);
381:     jac->lwork = 12*jac->n;
382: #if defined(PETSC_USE_COMPLEX)
383:     PetscMalloc1(jac->lwork,&jac->rwork);
384: #endif
385:     PetscMalloc1(jac->lwork,&jac->work);
386:   }

388:   return(0);
389: }

391: static PetscErrorCode SNESReset_Composite(SNES snes)
392: {
393:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
394:   PetscErrorCode   ierr;
395:   SNES_CompositeLink next = jac->head;

398:   while (next) {
399:     SNESReset(next->snes);
400:     next = next->next;
401:   }
402:   VecDestroy(&jac->Xorig);
403:   if (jac->Xes) {VecDestroyVecs(jac->nsnes,&jac->Xes);}
404:   if (jac->Fes) {VecDestroyVecs(jac->nsnes,&jac->Fes);}
405:   PetscFree(jac->fnorms);
406:   PetscFree(jac->h);
407:   PetscFree(jac->s);
408:   PetscFree(jac->g);
409:   PetscFree(jac->beta);
410:   PetscFree(jac->work);
411:   PetscFree(jac->rwork);
412:   return(0);
413: }

415: static PetscErrorCode SNESDestroy_Composite(SNES snes)
416: {
417:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
418:   PetscErrorCode     ierr;
419:   SNES_CompositeLink next = jac->head,next_tmp;

422:   SNESReset_Composite(snes);
423:   while (next) {
424:     SNESDestroy(&next->snes);
425:     next_tmp = next;
426:     next     = next->next;
427:     PetscFree(next_tmp);
428:   }
429:   PetscFree(snes->data);
430:   return(0);
431: }

433: static PetscErrorCode SNESSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,SNES snes)
434: {
435:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
436:   PetscErrorCode     ierr;
437:   PetscInt           nmax = 8,i;
438:   SNES_CompositeLink next;
439:   char               *sneses[8];
440:   PetscReal          dmps[8];
441:   PetscBool          flg;

444:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
445:   PetscOptionsEnum("-snes_composite_type","Type of composition","SNESCompositeSetType",SNESCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
446:   if (flg) {
447:     SNESCompositeSetType(snes,jac->type);
448:   }
449:   PetscOptionsStringArray("-snes_composite_sneses","List of composite solvers","SNESCompositeAddSNES",sneses,&nmax,&flg);
450:   if (flg) {
451:     for (i=0; i<nmax; i++) {
452:       SNESCompositeAddSNES(snes,sneses[i]);
453:       PetscFree(sneses[i]);   /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
454:     }
455:   }
456:   PetscOptionsRealArray("-snes_composite_damping","Damping of the additive composite solvers","SNESCompositeSetDamping",dmps,&nmax,&flg);
457:   if (flg) {
458:     for (i=0; i<nmax; i++) {
459:       SNESCompositeSetDamping(snes,i,dmps[i]);
460:     }
461:   }
462:   PetscOptionsReal("-snes_composite_stol","Step tolerance for restart on the additive composite solvers","",jac->stol,&jac->stol,NULL);
463:   PetscOptionsReal("-snes_composite_rtol","Residual tolerance for the additive composite solvers","",jac->rtol,&jac->rtol,NULL);
464:   PetscOptionsTail();

466:   next = jac->head;
467:   while (next) {
468:     SNESSetFromOptions(next->snes);
469:     next = next->next;
470:   }
471:   return(0);
472: }

474: static PetscErrorCode SNESView_Composite(SNES snes,PetscViewer viewer)
475: {
476:   SNES_Composite     *jac = (SNES_Composite*)snes->data;
477:   PetscErrorCode     ierr;
478:   SNES_CompositeLink next = jac->head;
479:   PetscBool          iascii;

482:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
483:   if (iascii) {
484:     PetscViewerASCIIPrintf(viewer,"  type - %s\n",SNESCompositeTypes[jac->type]);
485:     PetscViewerASCIIPrintf(viewer,"  SNESes on composite preconditioner follow\n");
486:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
487:   }
488:   if (iascii) {
489:     PetscViewerASCIIPushTab(viewer);
490:   }
491:   while (next) {
492:     SNESView(next->snes,viewer);
493:     next = next->next;
494:   }
495:   if (iascii) {
496:     PetscViewerASCIIPopTab(viewer);
497:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
498:   }
499:   return(0);
500: }

502: /* ------------------------------------------------------------------------------*/

504: static PetscErrorCode  SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)
505: {
506:   SNES_Composite *jac = (SNES_Composite*)snes->data;

509:   jac->type = type;
510:   return(0);
511: }

513: static PetscErrorCode  SNESCompositeAddSNES_Composite(SNES snes,SNESType type)
514: {
515:   SNES_Composite     *jac;
516:   SNES_CompositeLink next,ilink;
517:   PetscErrorCode     ierr;
518:   PetscInt           cnt = 0;
519:   const char         *prefix;
520:   char               newprefix[20];
521:   DM                 dm;

524:   PetscNewLog(snes,&ilink);
525:   ilink->next = 0;
526:   SNESCreate(PetscObjectComm((PetscObject)snes),&ilink->snes);
527:   PetscLogObjectParent((PetscObject)snes,(PetscObject)ilink->snes);
528:   SNESGetDM(snes,&dm);
529:   SNESSetDM(ilink->snes,dm);
530:   SNESSetTolerances(ilink->snes,ilink->snes->abstol,ilink->snes->rtol,ilink->snes->stol,1,ilink->snes->max_funcs);
531:   PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)ilink->snes);
532:   jac  = (SNES_Composite*)snes->data;
533:   next = jac->head;
534:   if (!next) {
535:     jac->head       = ilink;
536:     ilink->previous = NULL;
537:   } else {
538:     cnt++;
539:     while (next->next) {
540:       next = next->next;
541:       cnt++;
542:     }
543:     next->next      = ilink;
544:     ilink->previous = next;
545:   }
546:   SNESGetOptionsPrefix(snes,&prefix);
547:   SNESSetOptionsPrefix(ilink->snes,prefix);
548:   sprintf(newprefix,"sub_%d_",(int)cnt);
549:   SNESAppendOptionsPrefix(ilink->snes,newprefix);
550:   PetscObjectIncrementTabLevel((PetscObject)ilink->snes,(PetscObject)snes,1);
551:   SNESSetType(ilink->snes,type);
552:   SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY);

554:   ilink->dmp = 1.0;
555:   jac->nsnes++;
556:   return(0);
557: }

559: static PetscErrorCode  SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES *subsnes)
560: {
561:   SNES_Composite     *jac;
562:   SNES_CompositeLink next;
563:   PetscInt           i;

566:   jac  = (SNES_Composite*)snes->data;
567:   next = jac->head;
568:   for (i=0; i<n; i++) {
569:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
570:     next = next->next;
571:   }
572:   *subsnes = next->snes;
573:   return(0);
574: }

576: /* -------------------------------------------------------------------------------- */
577: /*@C
578:    SNESCompositeSetType - Sets the type of composite preconditioner.

580:    Logically Collective on SNES

582:    Input Parameter:
583: +  snes - the preconditioner context
584: -  type - SNES_COMPOSITE_ADDITIVE (default), SNES_COMPOSITE_MULTIPLICATIVE

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

589:    Level: Developer

591: @*/
592: PetscErrorCode  SNESCompositeSetType(SNES snes,SNESCompositeType type)
593: {

599:   PetscTryMethod(snes,"SNESCompositeSetType_C",(SNES,SNESCompositeType),(snes,type));
600:   return(0);
601: }

603: /*@C
604:    SNESCompositeAddSNES - Adds another SNES to the composite SNES.

606:    Collective on SNES

608:    Input Parameters:
609: +  snes - the preconditioner context
610: -  type - the type of the new preconditioner

612:    Level: Developer

614: @*/
615: PetscErrorCode  SNESCompositeAddSNES(SNES snes,SNESType type)
616: {

621:   PetscTryMethod(snes,"SNESCompositeAddSNES_C",(SNES,SNESType),(snes,type));
622:   return(0);
623: }
624: /*@
625:    SNESCompositeGetSNES - Gets one of the SNES objects in the composite SNES.

627:    Not Collective

629:    Input Parameter:
630: +  snes - the preconditioner context
631: -  n - the number of the snes requested

633:    Output Parameters:
634: .  subsnes - the SNES requested

636:    Level: Developer

638: .seealso: SNESCompositeAddSNES()
639: @*/
640: PetscErrorCode  SNESCompositeGetSNES(SNES snes,PetscInt n,SNES *subsnes)
641: {

647:   PetscUseMethod(snes,"SNESCompositeGetSNES_C",(SNES,PetscInt,SNES*),(snes,n,subsnes));
648:   return(0);
649: }

651: /*@
652:    SNESCompositeGetNumber - Get the number of subsolvers in the composite SNES.

654:    Logically Collective on SNES

656:    Input Parameter:
657:    snes - the preconditioner context

659:    Output Parameter:
660:    n - the number of subsolvers

662:    Level: Developer

664: @*/
665: PetscErrorCode  SNESCompositeGetNumber(SNES snes,PetscInt *n)
666: {
667:   SNES_Composite     *jac;
668:   SNES_CompositeLink next;

671:   jac  = (SNES_Composite*)snes->data;
672:   next = jac->head;

674:   *n = 0;
675:   while (next) {
676:     *n = *n + 1;
677:     next = next->next;
678:   }
679:   return(0);
680: }

682: static PetscErrorCode  SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)
683: {
684:   SNES_Composite     *jac;
685:   SNES_CompositeLink next;
686:   PetscInt           i;

689:   jac  = (SNES_Composite*)snes->data;
690:   next = jac->head;
691:   for (i=0; i<n; i++) {
692:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Not enough SNESes in composite preconditioner");
693:     next = next->next;
694:   }
695:   next->dmp = dmp;
696:   return(0);
697: }

699: /*@
700:    SNESCompositeSetDamping - Sets the damping of a subsolver when using additive composite SNES.

702:    Not Collective

704:    Input Parameter:
705: +  snes - the preconditioner context
706: .  n - the number of the snes requested
707: -  dmp - the damping

709:    Level: Developer

711: .seealso: SNESCompositeAddSNES()
712: @*/
713: PetscErrorCode  SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)
714: {

719:   PetscUseMethod(snes,"SNESCompositeSetDamping_C",(SNES,PetscInt,PetscReal),(snes,n,dmp));
720:   return(0);
721: }

723: static PetscErrorCode SNESSolve_Composite(SNES snes)
724: {
725:   Vec              F;
726:   Vec              X;
727:   Vec              B;
728:   PetscInt         i;
729:   PetscReal        fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
730:   PetscErrorCode   ierr;
731:   SNESNormSchedule normtype;
732:   SNES_Composite   *comp = (SNES_Composite*)snes->data;

735:   X = snes->vec_sol;
736:   F = snes->vec_func;
737:   B = snes->vec_rhs;

739:   PetscObjectSAWsTakeAccess((PetscObject)snes);
740:   snes->iter   = 0;
741:   snes->norm   = 0.;
742:   comp->innerFailures = 0;
743:   PetscObjectSAWsGrantAccess((PetscObject)snes);
744:   SNESSetWorkVecs(snes, 1);
745:   snes->reason = SNES_CONVERGED_ITERATING;
746:   SNESGetNormSchedule(snes, &normtype);
747:   if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
748:     if (!snes->vec_func_init_set) {
749:       SNESComputeFunction(snes,X,F);
750:     } else snes->vec_func_init_set = PETSC_FALSE;

752:     if (snes->xl && snes->xu) {
753:       SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
754:     } else {
755:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
756:     }
757:     SNESCheckFunctionNorm(snes,fnorm);
758:     PetscObjectSAWsTakeAccess((PetscObject)snes);
759:     snes->iter = 0;
760:     snes->norm = fnorm;
761:     PetscObjectSAWsGrantAccess((PetscObject)snes);
762:     SNESLogConvergenceHistory(snes,snes->norm,0);
763:     SNESMonitor(snes,0,snes->norm);

765:     /* test convergence */
766:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
767:     if (snes->reason) return(0);
768:   } else {
769:     PetscObjectSAWsGrantAccess((PetscObject)snes);
770:     SNESLogConvergenceHistory(snes,snes->norm,0);
771:     SNESMonitor(snes,0,snes->norm);
772:   }

774:   /* Call general purpose update function */
775:   if (snes->ops->update) {
776:     (*snes->ops->update)(snes, snes->iter);
777:   }

779:   for (i = 0; i < snes->max_its; i++) {
780:     /* Copy the state before modification by application of the composite solver;
781:        we will subtract the new state after application */
782:     VecCopy(X, snes->work[0]);

784:     if (comp->type == SNES_COMPOSITE_ADDITIVE) {
785:       SNESCompositeApply_Additive(snes,X,B,F,&fnorm);
786:     } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
787:       SNESCompositeApply_Multiplicative(snes,X,B,F,&fnorm);
788:     } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
789:       SNESCompositeApply_AdditiveOptimal(snes,X,B,F,&fnorm);
790:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Unsupported SNESComposite type");
791:     if (snes->reason < 0) break;

793:     /* Compute the solution update for convergence testing */
794:     VecAXPY(snes->work[0], -1.0, X);
795:     VecScale(snes->work[0], -1.0);

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

800:       if (snes->xl && snes->xu) {
801:         VecNormBegin(X, NORM_2, &xnorm);
802:         VecNormBegin(snes->work[0], NORM_2, &snorm);
803:         SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
804:         VecNormEnd(X, NORM_2, &xnorm);
805:         VecNormEnd(snes->work[0], NORM_2, &snorm);
806:       } else {
807:         VecNormBegin(F, NORM_2, &fnorm);
808:         VecNormBegin(X, NORM_2, &xnorm);
809:         VecNormBegin(snes->work[0], NORM_2, &snorm);

811:         VecNormEnd(F, NORM_2, &fnorm);
812:         VecNormEnd(X, NORM_2, &xnorm);
813:         VecNormEnd(snes->work[0], NORM_2, &snorm);
814:       }
815:       SNESCheckFunctionNorm(snes,fnorm);
816:     } else if (normtype == SNES_NORM_ALWAYS) {
817:       VecNormBegin(X, NORM_2, &xnorm);
818:       VecNormBegin(snes->work[0], NORM_2, &snorm);
819:       VecNormEnd(X, NORM_2, &xnorm);
820:       VecNormEnd(snes->work[0], NORM_2, &snorm);
821:     }
822:     /* Monitor convergence */
823:     PetscObjectSAWsTakeAccess((PetscObject)snes);
824:     snes->iter = i+1;
825:     snes->norm = fnorm;
826:     snes->xnorm = xnorm;
827:     snes->ynorm = snorm;
828:     PetscObjectSAWsGrantAccess((PetscObject)snes);
829:     SNESLogConvergenceHistory(snes,snes->norm,0);
830:     SNESMonitor(snes,snes->iter,snes->norm);
831:     /* Test for convergence */
832:     if (normtype == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,xnorm,snorm,fnorm,&snes->reason,snes->cnvP);}
833:     if (snes->reason) break;
834:     /* Call general purpose update function */
835:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
836:   }
837:   if (normtype == SNES_NORM_ALWAYS) {
838:     if (i == snes->max_its) {
839:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
840:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
841:     }
842:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
843:   return(0);
844: }

846: /* -------------------------------------------------------------------------------------------*/

848: /*MC
849:      SNESCOMPOSITE - Build a preconditioner by composing together several nonlinear solvers

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

855:    Level: intermediate

857: .seealso:  SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES,
858:            SNESSHELL, SNESCompositeSetType(), SNESCompositeSpecialSetAlpha(), SNESCompositeAddSNES(),
859:            SNESCompositeGetSNES()

861:    References:
862: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 
863:    SIAM Review, 57(4), 2015

865: M*/

867: PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
868: {
870:   SNES_Composite *jac;

873:   PetscNewLog(snes,&jac);

875:   snes->ops->solve           = SNESSolve_Composite;
876:   snes->ops->setup           = SNESSetUp_Composite;
877:   snes->ops->reset           = SNESReset_Composite;
878:   snes->ops->destroy         = SNESDestroy_Composite;
879:   snes->ops->setfromoptions  = SNESSetFromOptions_Composite;
880:   snes->ops->view            = SNESView_Composite;

882:   snes->usesksp        = PETSC_FALSE;

884:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

886:   snes->data = (void*)jac;
887:   jac->type  = SNES_COMPOSITE_ADDITIVEOPTIMAL;
888:   jac->Fes   = NULL;
889:   jac->Xes   = NULL;
890:   jac->fnorms = NULL;
891:   jac->nsnes = 0;
892:   jac->head  = 0;
893:   jac->stol  = 0.1;
894:   jac->rtol  = 1.1;

896:   jac->h     = NULL;
897:   jac->s     = NULL;
898:   jac->beta  = NULL;
899:   jac->work  = NULL;
900:   jac->rwork = NULL;

902:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetType_C",SNESCompositeSetType_Composite);
903:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeAddSNES_C",SNESCompositeAddSNES_Composite);
904:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeGetSNES_C",SNESCompositeGetSNES_Composite);
905:   PetscObjectComposeFunction((PetscObject)snes,"SNESCompositeSetDamping_C",SNESCompositeSetDamping_Composite);
906:   return(0);
907: }