Actual source code: ms.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  1: #include <petsc-private/snesimpl.h>   /*I "petscsnes.h" I*/

  3: static SNESMSType SNESMSDefault = SNESMSM62;
  4: static PetscBool  SNESMSRegisterAllCalled;
  5: static PetscBool  SNESMSPackageInitialized;

  7: typedef struct _SNESMSTableau *SNESMSTableau;
  8: struct _SNESMSTableau {
  9:   char      *name;
 10:   PetscInt  nstages;            /* Number of stages */
 11:   PetscInt  nregisters;         /* Number of registers */
 12:   PetscReal stability;          /* Scaled stability region */
 13:   PetscReal *gamma;             /* Coefficients of 3S* method */
 14:   PetscReal *delta;             /* Coefficients of 3S* method */
 15:   PetscReal *betasub;           /* Subdiagonal of beta in Shu-Osher form */
 16: };

 18: typedef struct _SNESMSTableauLink *SNESMSTableauLink;
 19: struct _SNESMSTableauLink {
 20:   struct _SNESMSTableau tab;
 21:   SNESMSTableauLink     next;
 22: };
 23: static SNESMSTableauLink SNESMSTableauList;

 25: typedef struct {
 26:   SNESMSTableau tableau;        /* Tableau in low-storage form */
 27:   PetscReal     damping;        /* Damping parameter, like length of (pseudo) time step */
 28:   PetscBool     norms;          /* Compute norms, usually only for monitoring purposes */
 29: } SNES_MS;

 33: /*@C
 34:   SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS

 36:   Not Collective, but should be called by all processes which will need the schemes to be registered

 38:   Level: advanced

 40: .keywords: SNES, SNESMS, register, all

 42: .seealso:  SNESMSRegisterDestroy()
 43: @*/
 44: PetscErrorCode SNESMSRegisterAll(void)
 45: {

 49:   if (SNESMSRegisterAllCalled) return(0);
 50:   SNESMSRegisterAllCalled = PETSC_TRUE;

 52:   {
 53:     PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}};
 54:     PetscReal delta[1]    = {0.0};
 55:     PetscReal betasub[1]  = {1.0};
 56:     SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);
 57:   }

 59:   {
 60:     PetscReal gamma[3][6] = {
 61:       {0.0000000000000000E+00,  -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02,  -1.4204296130641869E-01},
 62:       {1.0000000000000000E+00,  1.1111025767083920E+00,  5.6150921583923230E-01,  7.4151723494934041E-01,  3.1714538168600587E-01,  4.6479276238548706E-01},
 63:       {0.0000000000000000E+00,  0.0000000000000000E+00,  0.0000000000000000E+00,  6.7968174970583317E-01,  -4.1755042846051737E-03, -1.9115668129923846E-01}
 64:     };
 65:     PetscReal delta[6]   = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
 66:     PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
 67:     SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);
 68:   }
 69:   {
 70:     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
 71:     PetscReal delta[4]    = {0,0,0,0};
 72:     PetscReal betasub[4]  = {0.25, 0.5, 0.55, 1.0};
 73:     SNESMSRegister(SNESMSJAMESON83,4,3,1.0,&gamma[0][0],delta,betasub);
 74:   }
 75:   {                             /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
 76:     PetscReal gamma[3][2] = {{0,0},{0,0},{1,1}};
 77:     PetscReal delta[2]    = {0,0};
 78:     PetscReal betasub[2]  = {0.3333,1.0};
 79:     SNESMSRegister(SNESMSVLTP21,2,3,1.0,&gamma[0][0],delta,betasub);
 80:   }
 81:   {                             /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
 82:     PetscReal gamma[3][3] = {{0,0,0},{0,0,0},{1,1,1}};
 83:     PetscReal delta[3]    = {0,0,0};
 84:     PetscReal betasub[3]  = {0.1481,0.4000,1.0};
 85:     SNESMSRegister(SNESMSVLTP31,3,3,1.5,&gamma[0][0],delta,betasub);
 86:   }
 87:   {                             /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
 88:     PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
 89:     PetscReal delta[4]    = {0,0,0,0};
 90:     PetscReal betasub[4]  = {0.0833,0.2069,0.4265,1.0};
 91:     SNESMSRegister(SNESMSVLTP41,4,3,2.0,&gamma[0][0],delta,betasub);
 92:   }
 93:   {                             /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
 94:     PetscReal gamma[3][5] = {{0,0,0,0,0},{0,0,0,0,0},{1,1,1,1,1}};
 95:     PetscReal delta[5]    = {0,0,0,0,0};
 96:     PetscReal betasub[5]  = {0.0533,0.1263,0.2375,0.4414,1.0};
 97:     SNESMSRegister(SNESMSVLTP51,5,3,2.5,&gamma[0][0],delta,betasub);
 98:   }
 99:   {                             /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
100:     PetscReal gamma[3][6] = {{0,0,0,0,0,0},{0,0,0,0,0,0},{1,1,1,1,1,1}};
101:     PetscReal delta[6]    = {0,0,0,0,0,0};
102:     PetscReal betasub[6]  = {0.0370,0.0851,0.1521,0.2562,0.4512,1.0};
103:     SNESMSRegister(SNESMSVLTP61,6,3,3.0,&gamma[0][0],delta,betasub);
104:   }
105:   return(0);
106: }

110: /*@C
111:    SNESMSRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister().

113:    Not Collective

115:    Level: advanced

117: .keywords: TSRosW, register, destroy
118: .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegister()
119: @*/
120: PetscErrorCode SNESMSRegisterDestroy(void)
121: {
122:   PetscErrorCode    ierr;
123:   SNESMSTableauLink link;

126:   while ((link = SNESMSTableauList)) {
127:     SNESMSTableau t = &link->tab;
128:     SNESMSTableauList = link->next;

130:     PetscFree3(t->gamma,t->delta,t->betasub);
131:     PetscFree(t->name);
132:     PetscFree(link);
133:   }
134:   SNESMSRegisterAllCalled = PETSC_FALSE;
135:   return(0);
136: }

140: /*@C
141:   SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
142:   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to SNESCreate_MS()
143:   when using static libraries.

145:   Level: developer

147: .keywords: SNES, SNESMS, initialize, package
148: .seealso: PetscInitialize()
149: @*/
150: PetscErrorCode SNESMSInitializePackage(void)
151: {

155:   if (SNESMSPackageInitialized) return(0);
156:   SNESMSPackageInitialized = PETSC_TRUE;

158:   SNESMSRegisterAll();
159:   PetscRegisterFinalize(SNESMSFinalizePackage);
160:   return(0);
161: }

165: /*@C
166:   SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
167:   called from PetscFinalize().

169:   Level: developer

171: .keywords: Petsc, destroy, package
172: .seealso: PetscFinalize()
173: @*/
174: PetscErrorCode SNESMSFinalizePackage(void)
175: {

179:   SNESMSPackageInitialized = PETSC_FALSE;

181:   SNESMSRegisterDestroy();
182:   return(0);
183: }

187: /*@C
188:    SNESMSRegister - register a multistage scheme

190:    Not Collective, but the same schemes should be registered on all processes on which they will be used

192:    Input Parameters:
193: +  name - identifier for method
194: .  nstages - number of stages
195: .  nregisters - number of registers used by low-storage implementation
196: .  gamma - coefficients, see Ketcheson's paper
197: .  delta - coefficients, see Ketcheson's paper
198: -  betasub - subdiagonal of Shu-Osher form

200:    Notes:
201:    The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.

203:    Level: advanced

205: .keywords: SNES, register

207: .seealso: SNESMS
208: @*/
209: PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
210: {
211:   PetscErrorCode    ierr;
212:   SNESMSTableauLink link;
213:   SNESMSTableau     t;

217:   if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
218:   if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");

223:   PetscMalloc(sizeof(*link),&link);
224:   PetscMemzero(link,sizeof(*link));
225:   t             = &link->tab;
226:   PetscStrallocpy(name,&t->name);
227:   t->nstages    = nstages;
228:   t->nregisters = nregisters;
229:   t->stability  = stability;

231:   PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&t->betasub);
232:   PetscMemcpy(t->gamma,gamma,nstages*nregisters*sizeof(PetscReal));
233:   PetscMemcpy(t->delta,delta,nstages*sizeof(PetscReal));
234:   PetscMemcpy(t->betasub,betasub,nstages*sizeof(PetscReal));

236:   link->next        = SNESMSTableauList;
237:   SNESMSTableauList = link;
238:   return(0);
239: }

243: /*
244:   X - initial state, updated in-place.
245:   F - residual, computed at the initial X on input
246: */
247: static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
248: {
249:   PetscErrorCode  ierr;
250:   SNES_MS         *ms    = (SNES_MS*)snes->data;
251:   SNESMSTableau   t      = ms->tableau;
252:   const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
253:   Vec             S1,S2,S3,Y;
254:   PetscInt        i,nstages = t->nstages;;


258:   Y    = snes->work[0];
259:   S1   = X;
260:   S2   = snes->work[1];
261:   S3   = snes->work[2];
262:   VecZeroEntries(S2);
263:   VecCopy(X,S3);
264:   for (i=0; i<nstages; i++) {
265:     Vec         Ss[4];
266:     PetscScalar scoeff[4];

268:     Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;

270:     scoeff[0] = gamma[0*nstages+i]-(PetscReal)1.0;
271:     scoeff[1] = gamma[1*nstages+i];
272:     scoeff[2] = gamma[2*nstages+i];
273:     scoeff[3] = -betasub[i]*ms->damping;

275:     VecAXPY(S2,delta[i],S1);
276:     if (i>0) {
277:       SNESComputeFunction(snes,S1,F);
278:       if (snes->domainerror) {
279:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
280:         return(0);
281:       }
282:     }
283:     KSPSolve(snes->ksp,F,Y);
284:     VecMAXPY(S1,4,scoeff,Ss);
285:   }
286:   return(0);
287: }

291: static PetscErrorCode SNESSolve_MS(SNES snes)
292: {
293:   SNES_MS        *ms = (SNES_MS*)snes->data;
294:   Vec            X   = snes->vec_sol,F = snes->vec_func;
295:   PetscReal      fnorm;
296:   PetscInt       i;

300:   PetscCitationsRegister(SNESCitation,&SNEScite);
301:   snes->reason = SNES_CONVERGED_ITERATING;
302:   PetscObjectSAWsTakeAccess((PetscObject)snes);
303:   snes->iter   = 0;
304:   snes->norm   = 0.;
305:   PetscObjectSAWsGrantAccess((PetscObject)snes);
306:   if (!snes->vec_func_init_set) {
307:     SNESComputeFunction(snes,X,F);
308:     if (snes->domainerror) {
309:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
310:       return(0);
311:     }
312:   } else snes->vec_func_init_set = PETSC_FALSE;

314:   if (snes->jacobian) {         /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
315:     SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);
316:   }
317:   if (ms->norms) {
318:     VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F||  */
319:     if (PetscIsInfOrNanReal(fnorm)) {
320:       snes->reason = SNES_DIVERGED_FNORM_NAN;
321:       return(0);
322:     }
323:     /* Monitor convergence */
324:     PetscObjectSAWsTakeAccess((PetscObject)snes);
325:     snes->iter = 0;
326:     snes->norm = fnorm;
327:     PetscObjectSAWsGrantAccess((PetscObject)snes);
328:     SNESLogConvergenceHistory(snes,snes->norm,0);
329:     SNESMonitor(snes,snes->iter,snes->norm);

331:     /* Test for convergence */
332:     (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
333:     if (snes->reason) return(0);
334:   }

336:   /* Call general purpose update function */
337:   if (snes->ops->update) {
338:     (*snes->ops->update)(snes,snes->iter);
339:   }
340:   for (i = 0; i < snes->max_its; i++) {
341:     SNESMSStep_3Sstar(snes,X,F);

343:     if (i+1 < snes->max_its || ms->norms) {
344:       SNESComputeFunction(snes,X,F);
345:       if (snes->domainerror) {
346:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
347:         return(0);
348:       }
349:     }

351:     if (ms->norms) {
352:       VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F||  */
353:       if (PetscIsInfOrNanReal(fnorm)) {
354:         snes->reason = SNES_DIVERGED_FNORM_NAN;
355:         return(0);
356:       }

358:       /* Monitor convergence */
359:       PetscObjectSAWsTakeAccess((PetscObject)snes);
360:       snes->iter = i+1;
361:       snes->norm = fnorm;
362:       PetscObjectSAWsGrantAccess((PetscObject)snes);
363:       SNESLogConvergenceHistory(snes,snes->norm,0);
364:       SNESMonitor(snes,snes->iter,snes->norm);

366:       /* Test for convergence */
367:       (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
368:       if (snes->reason) return(0);
369:     }

371:     /* Call general purpose update function */
372:     if (snes->ops->update) {
373:       (*snes->ops->update)(snes, snes->iter);
374:     }
375:   }
376:   if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
377:   return(0);
378: }

382: static PetscErrorCode SNESSetUp_MS(SNES snes)
383: {
384:   SNES_MS        *ms = (SNES_MS*)snes->data;

388:   if (!ms->tableau) {SNESMSSetType(snes,SNESMSDefault);}
389:   SNESSetWorkVecs(snes,3);
390:   SNESSetUpMatrices(snes);
391:   return(0);
392: }

396: static PetscErrorCode SNESReset_MS(SNES snes)
397: {

400:   return(0);
401: }

405: static PetscErrorCode SNESDestroy_MS(SNES snes)
406: {

410:   PetscFree(snes->data);
411:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
412:   return(0);
413: }

417: static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
418: {
419:   PetscBool      iascii;
421:   SNES_MS        *ms = (SNES_MS*)snes->data;

424:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
425:   if (iascii) {
426:     SNESMSTableau tab = ms->tableau;
427:     PetscViewerASCIIPrintf(viewer,"  multi-stage method type: %s\n",tab ? tab->name : "not yet set");
428:   }
429:   return(0);
430: }

434: static PetscErrorCode SNESSetFromOptions_MS(SNES snes)
435: {
436:   SNES_MS        *ms = (SNES_MS*)snes->data;

440:   PetscOptionsHead("SNES MS options");
441:   {
442:     SNESMSTableauLink link;
443:     PetscInt          count,choice;
444:     PetscBool         flg;
445:     const char        **namelist;
446:     char              mstype[256];

448:     PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));
449:     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
450:     PetscMalloc1(count,&namelist);
451:     for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
452:     PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);
453:     SNESMSSetType(snes,flg ? namelist[choice] : mstype);
454:     PetscFree(namelist);
455:     PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);
456:     PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);
457:   }
458:   PetscOptionsTail();
459:   return(0);
460: }

464: PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
465: {
466:   PetscErrorCode    ierr;
467:   SNES_MS           *ms = (SNES_MS*)snes->data;
468:   SNESMSTableauLink link;
469:   PetscBool         match;

472:   if (ms->tableau) {
473:     PetscStrcmp(ms->tableau->name,mstype,&match);
474:     if (match) return(0);
475:   }
476:   for (link = SNESMSTableauList; link; link=link->next) {
477:     PetscStrcmp(link->tab.name,mstype,&match);
478:     if (match) {
479:       SNESReset_MS(snes);
480:       ms->tableau = &link->tab;
481:       return(0);
482:     }
483:   }
484:   SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
485:   return(0);
486: }

490: /*@C
491:   SNESMSSetType - Set the type of multistage smoother

493:   Logically collective

495:   Input Parameter:
496: +  snes - nonlinear solver context
497: -  mstype - type of multistage method

499:   Level: beginner

501: .seealso: SNESMSGetType(), SNESMS
502: @*/
503: PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
504: {

509:   PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));
510:   return(0);
511: }

513: /* -------------------------------------------------------------------------- */
514: /*MC
515:       SNESMS - multi-stage smoothers

517:       Options Database:

519: +     -snes_ms_type - type of multi-stage smoother
520: -     -snes_ms_damping - damping for multi-stage method

522:       Notes:
523:       These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
524:       FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).

526:       Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.

528:       The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().

530:       References:

532:       Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.

534:       Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.

536:       Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.

538:       Level: beginner

540: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV

542: M*/
545: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
546: {
548:   SNES_MS        *ms;

551:   SNESMSInitializePackage();

553:   snes->ops->setup          = SNESSetUp_MS;
554:   snes->ops->solve          = SNESSolve_MS;
555:   snes->ops->destroy        = SNESDestroy_MS;
556:   snes->ops->setfromoptions = SNESSetFromOptions_MS;
557:   snes->ops->view           = SNESView_MS;
558:   snes->ops->reset          = SNESReset_MS;

560:   snes->usespc  = PETSC_FALSE;
561:   snes->usesksp = PETSC_TRUE;

563:   PetscNewLog(snes,&ms);
564:   snes->data  = (void*)ms;
565:   ms->damping = 0.9;
566:   ms->norms   = PETSC_FALSE;

568:   PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);
569:   return(0);
570: }