Actual source code: ms.c
petsc-3.13.6 2020-09-29
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: .seealso: SNESMSRegisterDestroy()
39: @*/
40: PetscErrorCode SNESMSRegisterAll(void)
41: {
45: if (SNESMSRegisterAllCalled) return(0);
46: SNESMSRegisterAllCalled = PETSC_TRUE;
48: {
49: PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}};
50: PetscReal delta[1] = {0.0};
51: PetscReal betasub[1] = {1.0};
52: SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);
53: }
55: {
56: PetscReal gamma[3][6] = {
57: {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01},
58: {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01},
59: {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01}
60: };
61: PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
62: PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
63: SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);
64: }
65: {
66: PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
67: PetscReal delta[4] = {0,0,0,0};
68: PetscReal betasub[4] = {0.25, 0.5, 0.55, 1.0};
69: SNESMSRegister(SNESMSJAMESON83,4,3,1.0,&gamma[0][0],delta,betasub);
70: }
71: { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
72: PetscReal gamma[3][2] = {{0,0},{0,0},{1,1}};
73: PetscReal delta[2] = {0,0};
74: PetscReal betasub[2] = {0.3333,1.0};
75: SNESMSRegister(SNESMSVLTP21,2,3,1.0,&gamma[0][0],delta,betasub);
76: }
77: { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
78: PetscReal gamma[3][3] = {{0,0,0},{0,0,0},{1,1,1}};
79: PetscReal delta[3] = {0,0,0};
80: PetscReal betasub[3] = {0.1481,0.4000,1.0};
81: SNESMSRegister(SNESMSVLTP31,3,3,1.5,&gamma[0][0],delta,betasub);
82: }
83: { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
84: PetscReal gamma[3][4] = {{0,0,0,0},{0,0,0,0},{1,1,1,1}};
85: PetscReal delta[4] = {0,0,0,0};
86: PetscReal betasub[4] = {0.0833,0.2069,0.4265,1.0};
87: SNESMSRegister(SNESMSVLTP41,4,3,2.0,&gamma[0][0],delta,betasub);
88: }
89: { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
90: PetscReal gamma[3][5] = {{0,0,0,0,0},{0,0,0,0,0},{1,1,1,1,1}};
91: PetscReal delta[5] = {0,0,0,0,0};
92: PetscReal betasub[5] = {0.0533,0.1263,0.2375,0.4414,1.0};
93: SNESMSRegister(SNESMSVLTP51,5,3,2.5,&gamma[0][0],delta,betasub);
94: }
95: { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
96: PetscReal gamma[3][6] = {{0,0,0,0,0,0},{0,0,0,0,0,0},{1,1,1,1,1,1}};
97: PetscReal delta[6] = {0,0,0,0,0,0};
98: PetscReal betasub[6] = {0.0370,0.0851,0.1521,0.2562,0.4512,1.0};
99: SNESMSRegister(SNESMSVLTP61,6,3,3.0,&gamma[0][0],delta,betasub);
100: }
101: return(0);
102: }
104: /*@C
105: SNESMSRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister().
107: Not Collective
109: Level: advanced
111: .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegister()
112: @*/
113: PetscErrorCode SNESMSRegisterDestroy(void)
114: {
115: PetscErrorCode ierr;
116: SNESMSTableauLink link;
119: while ((link = SNESMSTableauList)) {
120: SNESMSTableau t = &link->tab;
121: SNESMSTableauList = link->next;
123: PetscFree3(t->gamma,t->delta,t->betasub);
124: PetscFree(t->name);
125: PetscFree(link);
126: }
127: SNESMSRegisterAllCalled = PETSC_FALSE;
128: return(0);
129: }
131: /*@C
132: SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called
133: from SNESInitializePackage().
135: Level: developer
137: .seealso: PetscInitialize()
138: @*/
139: PetscErrorCode SNESMSInitializePackage(void)
140: {
144: if (SNESMSPackageInitialized) return(0);
145: SNESMSPackageInitialized = PETSC_TRUE;
147: SNESMSRegisterAll();
148: PetscRegisterFinalize(SNESMSFinalizePackage);
149: return(0);
150: }
152: /*@C
153: SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is
154: called from PetscFinalize().
156: Level: developer
158: .seealso: PetscFinalize()
159: @*/
160: PetscErrorCode SNESMSFinalizePackage(void)
161: {
165: SNESMSPackageInitialized = PETSC_FALSE;
167: SNESMSRegisterDestroy();
168: return(0);
169: }
171: /*@C
172: SNESMSRegister - register a multistage scheme
174: Not Collective, but the same schemes should be registered on all processes on which they will be used
176: Input Parameters:
177: + name - identifier for method
178: . nstages - number of stages
179: . nregisters - number of registers used by low-storage implementation
180: . gamma - coefficients, see Ketcheson's paper
181: . delta - coefficients, see Ketcheson's paper
182: - betasub - subdiagonal of Shu-Osher form
184: Notes:
185: The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
187: Level: advanced
189: .seealso: SNESMS
190: @*/
191: PetscErrorCode SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
192: {
193: PetscErrorCode ierr;
194: SNESMSTableauLink link;
195: SNESMSTableau t;
199: if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
200: if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
205: SNESMSInitializePackage();
206: PetscNew(&link);
207: t = &link->tab;
208: PetscStrallocpy(name,&t->name);
209: t->nstages = nstages;
210: t->nregisters = nregisters;
211: t->stability = stability;
213: PetscMalloc3(nstages*nregisters,&t->gamma,nstages,&t->delta,nstages,&t->betasub);
214: PetscArraycpy(t->gamma,gamma,nstages*nregisters);
215: PetscArraycpy(t->delta,delta,nstages);
216: PetscArraycpy(t->betasub,betasub,nstages);
218: link->next = SNESMSTableauList;
219: SNESMSTableauList = link;
220: return(0);
221: }
223: /*
224: X - initial state, updated in-place.
225: F - residual, computed at the initial X on input
226: */
227: static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)
228: {
229: PetscErrorCode ierr;
230: SNES_MS *ms = (SNES_MS*)snes->data;
231: SNESMSTableau t = ms->tableau;
232: const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub;
233: Vec S1,S2,S3,Y;
234: PetscInt i,nstages = t->nstages;
238: Y = snes->work[0];
239: S1 = X;
240: S2 = snes->work[1];
241: S3 = snes->work[2];
242: VecZeroEntries(S2);
243: VecCopy(X,S3);
244: for (i = 0; i < nstages; i++) {
245: Vec Ss[4];
246: PetscScalar scoeff[4];
248: Ss[0] = S1; Ss[1] = S2; Ss[2] = S3; Ss[3] = Y;
250: scoeff[0] = gamma[0*nstages+i] - 1;
251: scoeff[1] = gamma[1*nstages+i];
252: scoeff[2] = gamma[2*nstages+i];
253: scoeff[3] = -betasub[i]*ms->damping;
255: VecAXPY(S2,delta[i],S1);
256: if (i > 0) {
257: SNESComputeFunction(snes,S1,F);
258: }
259: KSPSolve(snes->ksp,F,Y);
260: VecMAXPY(S1,4,scoeff,Ss);
261: }
262: return(0);
263: }
265: static PetscErrorCode SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)
266: {
267: SNES_MS *ms = (SNES_MS*)snes->data;
268: PetscReal fnorm;
272: if (ms->norms) {
273: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
274: SNESCheckFunctionNorm(snes,fnorm);
275: /* Monitor convergence */
276: PetscObjectSAWsTakeAccess((PetscObject)snes);
277: snes->iter = iter;
278: snes->norm = fnorm;
279: PetscObjectSAWsGrantAccess((PetscObject)snes);
280: SNESLogConvergenceHistory(snes,snes->norm,0);
281: SNESMonitor(snes,snes->iter,snes->norm);
282: /* Test for convergence */
283: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
284: } else if (iter > 0) {
285: PetscObjectSAWsTakeAccess((PetscObject)snes);
286: snes->iter = iter;
287: PetscObjectSAWsGrantAccess((PetscObject)snes);
288: }
289: return(0);
290: }
293: static PetscErrorCode SNESSolve_MS(SNES snes)
294: {
295: SNES_MS *ms = (SNES_MS*)snes->data;
296: Vec X = snes->vec_sol,F = snes->vec_func;
297: PetscInt i;
301: 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);
302: PetscCitationsRegister(SNESCitation,&SNEScite);
304: snes->reason = SNES_CONVERGED_ITERATING;
305: PetscObjectSAWsTakeAccess((PetscObject)snes);
306: snes->iter = 0;
307: snes->norm = 0;
308: PetscObjectSAWsGrantAccess((PetscObject)snes);
310: if (!snes->vec_func_init_set) {
311: SNESComputeFunction(snes,X,F);
312: } else snes->vec_func_init_set = PETSC_FALSE;
314: SNESMSStep_Norms(snes,0,F);
315: if (snes->reason) return(0);
317: for (i = 0; i < snes->max_its; i++) {
319: /* Call general purpose update function */
320: if (snes->ops->update) {
321: (*snes->ops->update)(snes,snes->iter);
322: }
324: if (i == 0 && snes->jacobian) {
325: /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
326: SNESComputeJacobian(snes,snes->vec_sol,snes->jacobian,snes->jacobian_pre);
327: SNESCheckJacobianDomainerror(snes);
328: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
329: }
331: SNESMSStep_3Sstar(snes,X,F);
333: if (i < snes->max_its-1 || ms->norms) {
334: SNESComputeFunction(snes,X,F);
335: }
337: SNESMSStep_Norms(snes,i+1,F);
338: if (snes->reason) return(0);
339: }
340: if (!snes->reason) snes->reason = SNES_CONVERGED_ITS;
341: return(0);
342: }
344: static PetscErrorCode SNESSetUp_MS(SNES snes)
345: {
346: SNES_MS *ms = (SNES_MS*)snes->data;
350: if (!ms->tableau) {SNESMSSetType(snes,SNESMSDefault);}
351: SNESSetWorkVecs(snes,3);
352: SNESSetUpMatrices(snes);
353: return(0);
354: }
356: static PetscErrorCode SNESReset_MS(SNES snes)
357: {
360: return(0);
361: }
363: static PetscErrorCode SNESDestroy_MS(SNES snes)
364: {
368: SNESReset_MS(snes);
369: PetscFree(snes->data);
370: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",NULL);
371: return(0);
372: }
374: static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer)
375: {
376: PetscBool iascii;
378: SNES_MS *ms = (SNES_MS*)snes->data;
381: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
382: if (iascii) {
383: SNESMSTableau tab = ms->tableau;
384: PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab ? tab->name : "not yet set");
385: }
386: return(0);
387: }
389: static PetscErrorCode SNESSetFromOptions_MS(PetscOptionItems *PetscOptionsObject,SNES snes)
390: {
391: SNES_MS *ms = (SNES_MS*)snes->data;
395: PetscOptionsHead(PetscOptionsObject,"SNES MS options");
396: {
397: SNESMSTableauLink link;
398: PetscInt count,choice;
399: PetscBool flg;
400: const char **namelist;
401: char mstype[256];
403: PetscStrncpy(mstype,SNESMSDefault,sizeof(mstype));
404: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ;
405: PetscMalloc1(count,(char***)&namelist);
406: for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
407: PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);
408: SNESMSSetType(snes,flg ? namelist[choice] : mstype);
409: PetscFree(namelist);
410: PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,NULL);
411: PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,NULL);
412: }
413: PetscOptionsTail();
414: return(0);
415: }
417: PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype)
418: {
419: PetscErrorCode ierr;
420: SNES_MS *ms = (SNES_MS*)snes->data;
421: SNESMSTableauLink link;
422: PetscBool match;
425: if (ms->tableau) {
426: PetscStrcmp(ms->tableau->name,mstype,&match);
427: if (match) return(0);
428: }
429: for (link = SNESMSTableauList; link; link=link->next) {
430: PetscStrcmp(link->tab.name,mstype,&match);
431: if (match) {
432: SNESReset_MS(snes);
433: ms->tableau = &link->tab;
434: return(0);
435: }
436: }
437: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
438: return(0);
439: }
441: /*@C
442: SNESMSSetType - Set the type of multistage smoother
444: Logically collective
446: Input Parameter:
447: + snes - nonlinear solver context
448: - mstype - type of multistage method
450: Level: beginner
452: .seealso: SNESMSGetType(), SNESMS
453: @*/
454: PetscErrorCode SNESMSSetType(SNES snes,SNESMSType rostype)
455: {
460: PetscTryMethod(snes,"SNESMSSetType_C",(SNES,SNESMSType),(snes,rostype));
461: return(0);
462: }
464: /* -------------------------------------------------------------------------- */
465: /*MC
466: SNESMS - multi-stage smoothers
468: Options Database:
470: + -snes_ms_type - type of multi-stage smoother
471: - -snes_ms_damping - damping for multi-stage method
473: Notes:
474: These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
475: FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
477: Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
479: The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister().
481: References:
482: + 1. - Ketcheson (2010) Runge Kutta methods with minimum storage implementations.
483: . 2. - Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method.
484: - 3. - Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes.
486: Level: beginner
488: .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYSHEV
490: M*/
491: PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
492: {
494: SNES_MS *ms;
497: SNESMSInitializePackage();
499: snes->ops->setup = SNESSetUp_MS;
500: snes->ops->solve = SNESSolve_MS;
501: snes->ops->destroy = SNESDestroy_MS;
502: snes->ops->setfromoptions = SNESSetFromOptions_MS;
503: snes->ops->view = SNESView_MS;
504: snes->ops->reset = SNESReset_MS;
506: snes->usesnpc = PETSC_FALSE;
507: snes->usesksp = PETSC_TRUE;
509: snes->alwayscomputesfinalresidual = PETSC_FALSE;
511: PetscNewLog(snes,&ms);
512: snes->data = (void*)ms;
513: ms->damping = 0.9;
514: ms->norms = PETSC_FALSE;
516: PetscObjectComposeFunction((PetscObject)snes,"SNESMSSetType_C",SNESMSSetType_MS);
517: return(0);
518: }