Actual source code: ms.c
petsc-3.4.5 2014-06-29
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,PetscReal,&t->gamma,nstages,PetscReal,&t->delta,nstages,PetscReal,&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: MatStructure mstruct;
297: PetscInt i;
301: snes->reason = SNES_CONVERGED_ITERATING;
302: PetscObjectAMSTakeAccess((PetscObject)snes);
303: snes->iter = 0;
304: snes->norm = 0.;
305: PetscObjectAMSGrantAccess((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,&mstruct);
316: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,mstruct);
317: }
318: if (ms->norms) {
319: if (!snes->norm_init_set) {
320: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
321: if (PetscIsInfOrNanReal(fnorm)) {
322: snes->reason = SNES_DIVERGED_FNORM_NAN;
323: return(0);
324: }
325: } else {
326: fnorm = snes->norm_init;
327: snes->norm_init_set = PETSC_FALSE;
328: }
329: /* Monitor convergence */
330: PetscObjectAMSTakeAccess((PetscObject)snes);
331: snes->iter = 0;
332: snes->norm = fnorm;
333: PetscObjectAMSGrantAccess((PetscObject)snes);
334: SNESLogConvergenceHistory(snes,snes->norm,0);
335: SNESMonitor(snes,snes->iter,snes->norm);
337: /* set parameter for default relative tolerance convergence test */
338: snes->ttol = fnorm*snes->rtol;
339: /* Test for convergence */
340: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
341: if (snes->reason) return(0);
342: }
344: /* Call general purpose update function */
345: if (snes->ops->update) {
346: (*snes->ops->update)(snes,snes->iter);
347: }
348: for (i = 0; i < snes->max_its; i++) {
349: SNESMSStep_3Sstar(snes,X,F);
351: if (i+1 < snes->max_its || ms->norms) {
352: SNESComputeFunction(snes,X,F);
353: if (snes->domainerror) {
354: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
355: return(0);
356: }
357: }
359: if (ms->norms) {
360: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
361: if (PetscIsInfOrNanReal(fnorm)) {
362: snes->reason = SNES_DIVERGED_FNORM_NAN;
363: return(0);
364: }
366: /* Monitor convergence */
367: PetscObjectAMSTakeAccess((PetscObject)snes);
368: snes->iter = i+1;
369: snes->norm = fnorm;
370: PetscObjectAMSGrantAccess((PetscObject)snes);
371: SNESLogConvergenceHistory(snes,snes->norm,0);
372: SNESMonitor(snes,snes->iter,snes->norm);
374: /* Test for convergence */
375: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
376: if (snes->reason) return(0);
377: }
379: /* Call general purpose update function */
380: if (snes->ops->update) {
381: (*snes->ops->update)(snes, snes->iter);
382: }
383: }
384: if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
385: return(0);
386: }
390: static PetscErrorCode SNESSetUp_MS(SNES snes)
391: {
392: SNES_MS *ms = (SNES_MS*)snes->data;
396: if (!ms->tableau) {SNESMSSetType(snes,SNESMSDefault);}
397: SNESSetWorkVecs(snes,3);
398: SNESSetUpMatrices(snes);
399: return(0);
400: }
404: static PetscErrorCode SNESReset_MS(SNES snes)
405: {
408: return(0);
409: }
413: static PetscErrorCode SNESDestroy_MS(SNES snes)
414: {
418: PetscFree(snes->data);
419: PetscObjectComposeFunction((PetscObject)snes,"",NULL);
420: return(0);
421: }
425: static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
426: {
427: PetscBool iascii;
429: SNES_MS *ms = (SNES_MS*)snes->data;
432: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
433: if (iascii) {
434: SNESMSTableau tab = ms->tableau;
435: PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab ? tab->name : "not yet set");
436: }
437: return(0);
438: }
442: static PetscErrorCode SNESSetFromOptions_MS(SNES snes)
443: {
444: SNES_MS *ms = (SNES_MS*)snes->data;
448: PetscOptionsHead("SNES MS options");
449: {
450: SNESMSTableauLink link;
451: PetscInt count,choice;
452: PetscBool flg;
453: const char **namelist;
454: char mstype[256];
456: PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));
457: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
458: PetscMalloc(count*sizeof(char*),&namelist);
459: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
460: PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);
461: SNESMSSetType(snes,flg ? namelist[choice] : mstype);
462: PetscFree(namelist);
463: PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);
464: PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);
465: }
466: PetscOptionsTail();
467: return(0);
468: }
472: PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
473: {
474: PetscErrorCode ierr;
475: SNES_MS *ms = (SNES_MS*)snes->data;
476: SNESMSTableauLink link;
477: PetscBool match;
480: if (ms->tableau) {
481: PetscStrcmp(ms->tableau->name,mstype,&match);
482: if (match) return(0);
483: }
484: for (link = SNESMSTableauList; link; link=link->next) {
485: PetscStrcmp(link->tab.name,mstype,&match);
486: if (match) {
487: SNESReset_MS(snes);
488: ms->tableau = &link->tab;
489: return(0);
490: }
491: }
492: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
493: return(0);
494: }
498: /*@C
499: SNESMSSetType - Set the type of multistage smoother
501: Logically collective
503: Input Parameter:
504: + snes - nonlinear solver context
505: - mstype - type of multistage method
507: Level: beginner
509: .seealso: SNESMSGetType(), SNESMS
510: @*/
511: PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
512: {
517: PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));
518: return(0);
519: }
521: /* -------------------------------------------------------------------------- */
522: /*MC
523: SNESMS - multi-stage smoothers
525: Options Database:
527: + -snes_ms_type - type of multi-stage smoother
528: - -snes_ms_damping - damping for multi-stage method
530: Notes:
531: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
532: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
534: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
536: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
538: References:
540: Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
542: Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
544: Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
546: Level: beginner
548: .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
550: M*/
553: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
554: {
556: SNES_MS *ms;
559: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
560: SNESMSInitializePackage();
561: #endif
563: snes->ops->setup = SNESSetUp_MS;
564: snes->ops->solve = SNESSolve_MS;
565: snes->ops->destroy = SNESDestroy_MS;
566: snes->ops->setfromoptions = SNESSetFromOptions_MS;
567: snes->ops->view = SNESView_MS;
568: snes->ops->reset = SNESReset_MS;
570: snes->usespc = PETSC_FALSE;
571: snes->usesksp = PETSC_TRUE;
573: PetscNewLog(snes,SNES_MS,&ms);
574: snes->data = (void*)ms;
575: ms->damping = 0.9;
576: ms->norms = PETSC_FALSE;
578: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);
579: return(0);
580: }