Actual source code: ms.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/snesimpl.h>
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;
31: /*@C
32: SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS
34: Not Collective, but should be called by all processes which will need the schemes to be registered
36: Level: advanced
38: .keywords: SNES, SNESMS, register, all
40: .seealso: SNESMSRegisterDestroy()
41: @*/
42: PetscErrorCode SNESMSRegisterAll(void)
43: {
47: if (SNESMSRegisterAllCalled) return(0);
48: SNESMSRegisterAllCalled = PETSC_TRUE;
50: {
51: PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}};
52: PetscReal delta[1] = {0.0};
53: PetscReal betasub[1] = {1.0};
54: SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);
55: }
57: {
58: PetscReal gamma[3][6] = {
59: {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01},
60: {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01},
61: {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01}
62: };
63: PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
64: PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
65: SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);
66: }
67: {
68: PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
69: PetscReal delta[4] = {0,0,0,0};
70: PetscReal betasub[4] = {0.25, 0.5, 0.55, 1.0};
71: SNESMSRegister(SNESMSJAMESON83,4,3,1.0,&gamma[0][0],delta,betasub);
72: }
73: { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
74: PetscReal gamma[3][2] = {{0,0},{0,0},{1,1}};
75: PetscReal delta[2] = {0,0};
76: PetscReal betasub[2] = {0.3333,1.0};
77: SNESMSRegister(SNESMSVLTP21,2,3,1.0,&gamma[0][0],delta,betasub);
78: }
79: { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
80: PetscReal gamma[3][3] = {{0,0,0},{0,0,0},{1,1,1}};
81: PetscReal delta[3] = {0,0,0};
82: PetscReal betasub[3] = {0.1481,0.4000,1.0};
83: SNESMSRegister(SNESMSVLTP31,3,3,1.5,&gamma[0][0],delta,betasub);
84: }
85: { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
86: PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
87: PetscReal delta[4] = {0,0,0,0};
88: PetscReal betasub[4] = {0.0833,0.2069,0.4265,1.0};
89: SNESMSRegister(SNESMSVLTP41,4,3,2.0,&gamma[0][0],delta,betasub);
90: }
91: { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
92: PetscReal gamma[3][5] = {{0,0,0,0,0},{0,0,0,0,0},{1,1,1,1,1}};
93: PetscReal delta[5] = {0,0,0,0,0};
94: PetscReal betasub[5] = {0.0533,0.1263,0.2375,0.4414,1.0};
95: SNESMSRegister(SNESMSVLTP51,5,3,2.5,&gamma[0][0],delta,betasub);
96: }
97: { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
98: PetscReal gamma[3][6] = {{0,0,0,0,0,0},{0,0,0,0,0,0},{1,1,1,1,1,1}};
99: PetscReal delta[6] = {0,0,0,0,0,0};
100: PetscReal betasub[6] = {0.0370,0.0851,0.1521,0.2562,0.4512,1.0};
101: SNESMSRegister(SNESMSVLTP61,6,3,3.0,&gamma[0][0],delta,betasub);
102: }
103: return(0);
104: }
106: /*@C
107: SNESMSRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister().
109: Not Collective
111: Level: advanced
113: .keywords: TSRosW, register, destroy
114: .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegister()
115: @*/
116: PetscErrorCode SNESMSRegisterDestroy(void)
117: {
118: PetscErrorCode ierr;
119: SNESMSTableauLink link;
122: while ((link = SNESMSTableauList)) {
123: SNESMSTableau t = &link->tab;
124: SNESMSTableauList = link->next;
126: PetscFree3(t->gamma,t->delta,t->betasub);
127: PetscFree(t->name);
128: PetscFree(link);
129: }
130: SNESMSRegisterAllCalled = PETSC_FALSE;
131: return(0);
132: }
134: /*@C
135: SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
136: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to SNESCreate_MS()
137: when using static libraries.
139: Level: developer
141: .keywords: SNES, SNESMS, initialize, package
142: .seealso: PetscInitialize()
143: @*/
144: PetscErrorCode SNESMSInitializePackage(void)
145: {
149: if (SNESMSPackageInitialized) return(0);
150: SNESMSPackageInitialized = PETSC_TRUE;
152: SNESMSRegisterAll();
153: PetscRegisterFinalize(SNESMSFinalizePackage);
154: return(0);
155: }
157: /*@C
158: SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
159: called from PetscFinalize().
161: Level: developer
163: .keywords: Petsc, destroy, package
164: .seealso: PetscFinalize()
165: @*/
166: PetscErrorCode SNESMSFinalizePackage(void)
167: {
171: SNESMSPackageInitialized = PETSC_FALSE;
173: SNESMSRegisterDestroy();
174: return(0);
175: }
177: /*@C
178: SNESMSRegister - register a multistage scheme
180: Not Collective, but the same schemes should be registered on all processes on which they will be used
182: Input Parameters:
183: + name - identifier for method
184: . nstages - number of stages
185: . nregisters - number of registers used by low-storage implementation
186: . gamma - coefficients, see Ketcheson's paper
187: . delta - coefficients, see Ketcheson's paper
188: - betasub - subdiagonal of Shu-Osher form
190: Notes:
191: The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
193: Level: advanced
195: .keywords: SNES, register
197: .seealso: SNESMS
198: @*/
199: PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
200: {
201: PetscErrorCode ierr;
202: SNESMSTableauLink link;
203: SNESMSTableau t;
207: if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
208: if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
213: SNESMSInitializePackage();
214: PetscNew(&link);
215: t = &link->tab;
216: PetscStrallocpy(name,&t->name);
217: t->nstages = nstages;
218: t->nregisters = nregisters;
219: t->stability = stability;
221: PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&t->betasub);
222: PetscMemcpy(t->gamma,gamma,nstages*nregisters*sizeof(PetscReal));
223: PetscMemcpy(t->delta,delta,nstages*sizeof(PetscReal));
224: PetscMemcpy(t->betasub,betasub,nstages*sizeof(PetscReal));
226: link->next = SNESMSTableauList;
227: SNESMSTableauList = link;
228: return(0);
229: }
231: /*
232: X - initial state, updated in-place.
233: F - residual, computed at the initial X on input
234: */
235: static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
236: {
237: PetscErrorCode ierr;
238: SNES_MS *ms = (SNES_MS*)snes->data;
239: SNESMSTableau t = ms->tableau;
240: const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
241: Vec S1,S2,S3,Y;
242: PetscInt i,nstages = t->nstages;;
246: Y = snes->work[0];
247: S1 = X;
248: S2 = snes->work[1];
249: S3 = snes->work[2];
250: VecZeroEntries(S2);
251: VecCopy(X,S3);
252: for (i=0; i<nstages; i++) {
253: Vec Ss[4];
254: PetscScalar scoeff[4];
256: Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
258: scoeff[0] = gamma[0*nstages+i]-(PetscReal)1.0;
259: scoeff[1] = gamma[1*nstages+i];
260: scoeff[2] = gamma[2*nstages+i];
261: scoeff[3] = -betasub[i]*ms->damping;
263: VecAXPY(S2,delta[i],S1);
264: if (i>0) {
265: SNESComputeFunction(snes,S1,F);
266: }
267: KSPSolve(snes->ksp,F,Y);
268: VecMAXPY(S1,4,scoeff,Ss);
269: }
270: return(0);
271: }
273: static PetscErrorCode SNESSolve_MS(SNES snes)
274: {
275: SNES_MS *ms = (SNES_MS*)snes->data;
276: Vec X = snes->vec_sol,F = snes->vec_func;
277: PetscReal fnorm;
278: PetscInt i;
282: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
284: PetscCitationsRegister(SNESCitation,&SNEScite);
285: snes->reason = SNES_CONVERGED_ITERATING;
286: PetscObjectSAWsTakeAccess((PetscObject)snes);
287: snes->iter = 0;
288: snes->norm = 0.;
289: PetscObjectSAWsGrantAccess((PetscObject)snes);
290: if (!snes->vec_func_init_set) {
291: SNESComputeFunction(snes,X,F);
292: } else snes->vec_func_init_set = PETSC_FALSE;
294: if (snes->jacobian) { /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
295: SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);
296: }
297: if (ms->norms) {
298: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
299: SNESCheckFunctionNorm(snes,fnorm);
300: /* Monitor convergence */
301: PetscObjectSAWsTakeAccess((PetscObject)snes);
302: snes->iter = 0;
303: snes->norm = fnorm;
304: PetscObjectSAWsGrantAccess((PetscObject)snes);
305: SNESLogConvergenceHistory(snes,snes->norm,0);
306: SNESMonitor(snes,snes->iter,snes->norm);
308: /* Test for convergence */
309: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
310: if (snes->reason) return(0);
311: }
313: /* Call general purpose update function */
314: if (snes->ops->update) {
315: (*snes->ops->update)(snes,snes->iter);
316: }
317: for (i = 0; i < snes->max_its; i++) {
318: SNESMSStep_3Sstar(snes,X,F);
320: if (i+1 < snes->max_its || ms->norms) {
321: SNESComputeFunction(snes,X,F);
322: }
324: if (ms->norms) {
325: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
326: SNESCheckFunctionNorm(snes,fnorm);
328: /* Monitor convergence */
329: PetscObjectSAWsTakeAccess((PetscObject)snes);
330: snes->iter = i+1;
331: snes->norm = fnorm;
332: PetscObjectSAWsGrantAccess((PetscObject)snes);
333: SNESLogConvergenceHistory(snes,snes->norm,0);
334: SNESMonitor(snes,snes->iter,snes->norm);
336: /* Test for convergence */
337: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
338: if (snes->reason) return(0);
339: }
341: /* Call general purpose update function */
342: if (snes->ops->update) {
343: (*snes->ops->update)(snes, snes->iter);
344: }
345: }
346: if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
347: return(0);
348: }
350: static PetscErrorCode SNESSetUp_MS(SNES snes)
351: {
352: SNES_MS *ms = (SNES_MS*)snes->data;
356: if (!ms->tableau) {SNESMSSetType(snes,SNESMSDefault);}
357: SNESSetWorkVecs(snes,3);
358: SNESSetUpMatrices(snes);
359: return(0);
360: }
362: static PetscErrorCode SNESReset_MS(SNES snes)
363: {
366: return(0);
367: }
369: static PetscErrorCode SNESDestroy_MS(SNES snes)
370: {
374: PetscFree(snes->data);
375: PetscObjectComposeFunction((PetscObject)snes,"",NULL);
376: return(0);
377: }
379: static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
380: {
381: PetscBool iascii;
383: SNES_MS *ms = (SNES_MS*)snes->data;
386: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
387: if (iascii) {
388: SNESMSTableau tab = ms->tableau;
389: PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab ? tab->name : "not yet set");
390: }
391: return(0);
392: }
394: static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
395: {
396: SNES_MS *ms = (SNES_MS*)snes->data;
400: PetscOptionsHead(PetscOptionsObject,"SNES MS options");
401: {
402: SNESMSTableauLink link;
403: PetscInt count,choice;
404: PetscBool flg;
405: const char **namelist;
406: char mstype[256];
408: PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));
409: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
410: PetscMalloc1(count,(char***)&namelist);
411: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
412: PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);
413: SNESMSSetType(snes,flg ? namelist[choice] : mstype);
414: PetscFree(namelist);
415: PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);
416: PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);
417: }
418: PetscOptionsTail();
419: return(0);
420: }
422: PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
423: {
424: PetscErrorCode ierr;
425: SNES_MS *ms = (SNES_MS*)snes->data;
426: SNESMSTableauLink link;
427: PetscBool match;
430: if (ms->tableau) {
431: PetscStrcmp(ms->tableau->name,mstype,&match);
432: if (match) return(0);
433: }
434: for (link = SNESMSTableauList; link; link=link->next) {
435: PetscStrcmp(link->tab.name,mstype,&match);
436: if (match) {
437: SNESReset_MS(snes);
438: ms->tableau = &link->tab;
439: return(0);
440: }
441: }
442: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
443: return(0);
444: }
446: /*@C
447: SNESMSSetType - Set the type of multistage smoother
449: Logically collective
451: Input Parameter:
452: + snes - nonlinear solver context
453: - mstype - type of multistage method
455: Level: beginner
457: .seealso: SNESMSGetType(), SNESMS
458: @*/
459: PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
460: {
465: PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));
466: return(0);
467: }
469: /* -------------------------------------------------------------------------- */
470: /*MC
471: SNESMS - multi-stage smoothers
473: Options Database:
475: + -snes_ms_type - type of multi-stage smoother
476: - -snes_ms_damping - damping for multi-stage method
478: Notes:
479: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
480: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
482: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
484: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
486: References:
487: + 1. - Ketcheson (2010) Runge Kutta methods with minimum storage implementations.
488: . 2. - Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
489: - 3. - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
491: Level: beginner
493: .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
495: M*/
496: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
497: {
499: SNES_MS *ms;
502: SNESMSInitializePackage();
504: snes->ops->setup = SNESSetUp_MS;
505: snes->ops->solve = SNESSolve_MS;
506: snes->ops->destroy = SNESDestroy_MS;
507: snes->ops->setfromoptions = SNESSetFromOptions_MS;
508: snes->ops->view = SNESView_MS;
509: snes->ops->reset = SNESReset_MS;
511: snes->usesnpc = PETSC_FALSE;
512: snes->usesksp = PETSC_TRUE;
514: snes->alwayscomputesfinalresidual = PETSC_FALSE;
516: PetscNewLog(snes,&ms);
517: snes->data = (void*)ms;
518: ms->damping = 0.9;
519: ms->norms = PETSC_FALSE;
521: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);
522: return(0);
523: }