1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2
3 static SNESMSType SNESMSDefault = SNESMSM62;
4 static PetscBool SNESMSRegisterAllCalled;
5 static PetscBool SNESMSPackageInitialized;
6
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 };
17
18 typedef struct _SNESMSTableauLink *SNESMSTableauLink;
19 struct _SNESMSTableauLink {
20 struct _SNESMSTableau tab;
21 SNESMSTableauLink next;
22 };
23 static SNESMSTableauLink SNESMSTableauList;
24
25 typedef struct {
26 SNESMSTableau tableau; /* Tableau in low-storage form */
27 PetscReal damping; /* Damping parameter, like length of (pseudo) time step */
28 } SNES_MS;
29
30 /*@C
31 SNESMSRegisterAll - Registers all of the multi-stage methods in `SNESMS`
32
33 Logically Collective
34
35 Level: developer
36
37 .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegisterDestroy()`
38 @*/
SNESMSRegisterAll(void)39 PetscErrorCode SNESMSRegisterAll(void)
40 {
41 PetscFunctionBegin;
42 if (SNESMSRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
43 SNESMSRegisterAllCalled = PETSC_TRUE;
44
45 {
46 PetscReal alpha[1] = {1.0};
47 PetscCall(SNESMSRegister(SNESMSEULER, 1, 1, 1.0, NULL, NULL, alpha));
48 }
49
50 {
51 PetscReal gamma[3][6] = {
52 {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01},
53 {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01 },
54 {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01}
55 };
56 PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00};
57 PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01};
58 PetscCall(SNESMSRegister(SNESMSM62, 6, 3, 1.0, &gamma[0][0], delta, betasub));
59 }
60
61 { /* Jameson (1983) */
62 PetscReal alpha[4] = {0.25, 0.5, 0.55, 1.0};
63 PetscCall(SNESMSRegister(SNESMSJAMESON83, 4, 1, 1.0, NULL, NULL, alpha));
64 }
65
66 { /* Van Leer, Tai, and Powell (1989) 1 stage, order 1 */
67 PetscReal alpha[1] = {1.0};
68 PetscCall(SNESMSRegister(SNESMSVLTP11, 1, 1, 0.5, NULL, NULL, alpha));
69 }
70 { /* Van Leer, Tai, and Powell (1989) 2 stage, order 1 */
71 PetscReal alpha[2] = {0.3333, 1.0};
72 PetscCall(SNESMSRegister(SNESMSVLTP21, 2, 1, 1.0, NULL, NULL, alpha));
73 }
74 { /* Van Leer, Tai, and Powell (1989) 3 stage, order 1 */
75 PetscReal alpha[3] = {0.1481, 0.4000, 1.0};
76 PetscCall(SNESMSRegister(SNESMSVLTP31, 3, 1, 1.5, NULL, NULL, alpha));
77 }
78 { /* Van Leer, Tai, and Powell (1989) 4 stage, order 1 */
79 PetscReal alpha[4] = {0.0833, 0.2069, 0.4265, 1.0};
80 PetscCall(SNESMSRegister(SNESMSVLTP41, 4, 1, 2.0, NULL, NULL, alpha));
81 }
82 { /* Van Leer, Tai, and Powell (1989) 5 stage, order 1 */
83 PetscReal alpha[5] = {0.0533, 0.1263, 0.2375, 0.4414, 1.0};
84 PetscCall(SNESMSRegister(SNESMSVLTP51, 5, 1, 2.5, NULL, NULL, alpha));
85 }
86 { /* Van Leer, Tai, and Powell (1989) 6 stage, order 1 */
87 PetscReal alpha[6] = {0.0370, 0.0851, 0.1521, 0.2562, 0.4512, 1.0};
88 PetscCall(SNESMSRegister(SNESMSVLTP61, 6, 1, 3.0, NULL, NULL, alpha));
89 }
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
93 /*@C
94 SNESMSRegisterDestroy - Frees the list of schemes that were registered by `SNESMSRegister()`.
95
96 Logically Collective
97
98 Level: developer
99
100 .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`
101 @*/
SNESMSRegisterDestroy(void)102 PetscErrorCode SNESMSRegisterDestroy(void)
103 {
104 SNESMSTableauLink link;
105
106 PetscFunctionBegin;
107 while ((link = SNESMSTableauList)) {
108 SNESMSTableau t = &link->tab;
109 SNESMSTableauList = link->next;
110
111 PetscCall(PetscFree(t->name));
112 PetscCall(PetscFree(t->gamma));
113 PetscCall(PetscFree(t->delta));
114 PetscCall(PetscFree(t->betasub));
115 PetscCall(PetscFree(link));
116 }
117 SNESMSRegisterAllCalled = PETSC_FALSE;
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
121 /*@C
122 SNESMSInitializePackage - This function initializes everything in the `SNESMS` package. It is called
123 from `SNESInitializePackage()`.
124
125 Level: developer
126
127 .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`, `PetscInitialize()`
128 @*/
SNESMSInitializePackage(void)129 PetscErrorCode SNESMSInitializePackage(void)
130 {
131 PetscFunctionBegin;
132 if (SNESMSPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
133 SNESMSPackageInitialized = PETSC_TRUE;
134
135 PetscCall(SNESMSRegisterAll());
136 PetscCall(PetscRegisterFinalize(SNESMSFinalizePackage));
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
140 /*@C
141 SNESMSFinalizePackage - This function destroys everything in the `SNESMS` package. It is
142 called from `PetscFinalize()`.
143
144 Level: developer
145
146 .seealso: [](ch_snes), `SNES`, `SNESMS`, `SNESMSRegister()`, `SNESMSRegisterAll()`, `SNESMSInitializePackage()`, `PetscFinalize()`
147 @*/
SNESMSFinalizePackage(void)148 PetscErrorCode SNESMSFinalizePackage(void)
149 {
150 PetscFunctionBegin;
151 SNESMSPackageInitialized = PETSC_FALSE;
152
153 PetscCall(SNESMSRegisterDestroy());
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
157 /*@C
158 SNESMSRegister - register a multistage scheme for `SNESMS`
159
160 Logically Collective, No Fortran Support
161
162 Input Parameters:
163 + name - identifier for method
164 . nstages - number of stages
165 . nregisters - number of registers used by low-storage implementation
166 . stability - scaled stability region
167 . gamma - coefficients, see Ketcheson's paper {cite}`ketcheson2010runge`
168 . delta - coefficients, see Ketcheson's paper {cite}`ketcheson2010runge`
169 - betasub - subdiagonal of Shu-Osher form
170
171 Level: advanced
172
173 Notes:
174 The notation is described in {cite}`ketcheson2010runge` Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.
175
176 Many multistage schemes are of the form
177
178 $$
179 \begin{align*}
180 X_0 = X^{(n)} \\
181 X_k = X_0 + \alpha_k * F(X_{k-1}), k = 1,\ldots,s \\
182 X^{(n+1)} = X_s
183 \end{align*}
184 $$
185
186 These methods can be registered with
187 .vb
188 SNESMSRegister("name",s,1,stability,NULL,NULL,alpha);
189 .ve
190
191 .seealso: [](ch_snes), `SNES`, `SNESMS`
192 @*/
SNESMSRegister(SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])193 PetscErrorCode SNESMSRegister(SNESMSType name, PetscInt nstages, PetscInt nregisters, PetscReal stability, const PetscReal gamma[], const PetscReal delta[], const PetscReal betasub[])
194 {
195 SNESMSTableauLink link;
196 SNESMSTableau t;
197
198 PetscFunctionBegin;
199 PetscAssertPointer(name, 1);
200 PetscCheck(nstages >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have at least one stage");
201 if (gamma || delta) {
202 PetscCheck(nregisters == 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 3-register form");
203 PetscAssertPointer(gamma, 5);
204 PetscAssertPointer(delta, 6);
205 } else {
206 PetscCheck(nregisters == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support for methods written in 1-register form");
207 }
208 PetscAssertPointer(betasub, 7);
209
210 PetscCall(SNESMSInitializePackage());
211 PetscCall(PetscNew(&link));
212 t = &link->tab;
213 PetscCall(PetscStrallocpy(name, &t->name));
214 t->nstages = nstages;
215 t->nregisters = nregisters;
216 t->stability = stability;
217
218 if (gamma && delta) {
219 PetscCall(PetscMalloc1(nstages * nregisters, &t->gamma));
220 PetscCall(PetscMalloc1(nstages, &t->delta));
221 PetscCall(PetscArraycpy(t->gamma, gamma, nstages * nregisters));
222 PetscCall(PetscArraycpy(t->delta, delta, nstages));
223 }
224 PetscCall(PetscMalloc1(nstages, &t->betasub));
225 PetscCall(PetscArraycpy(t->betasub, betasub, nstages));
226
227 link->next = SNESMSTableauList;
228 SNESMSTableauList = link;
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 /*
233 X - initial state, updated in-place.
234 F - residual, computed at the initial X on input
235 */
SNESMSStep_3Sstar(SNES snes,Vec X,Vec F)236 static PetscErrorCode SNESMSStep_3Sstar(SNES snes, Vec X, Vec F)
237 {
238 SNES_MS *ms = (SNES_MS *)snes->data;
239 SNESMSTableau t = ms->tableau;
240 const PetscInt nstages = t->nstages;
241 const PetscReal *gamma = t->gamma, *delta = t->delta, *betasub = t->betasub;
242 Vec S1 = X, S2 = snes->work[1], S3 = snes->work[2], Y = snes->work[0], S1copy = snes->work[3];
243
244 PetscFunctionBegin;
245 PetscCall(VecZeroEntries(S2));
246 PetscCall(VecCopy(X, S3));
247 for (PetscInt i = 0; i < nstages; i++) {
248 Vec Ss[] = {S1copy, S2, S3, Y};
249 const PetscScalar scoeff[] = {gamma[0 * nstages + i] - 1, gamma[1 * nstages + i], gamma[2 * nstages + i], -betasub[i] * ms->damping};
250
251 PetscCall(VecAXPY(S2, delta[i], S1));
252 if (i > 0) PetscCall(SNESComputeFunction(snes, S1, F));
253 PetscCall(KSPSolve(snes->ksp, F, Y));
254 PetscCall(VecCopy(S1, S1copy));
255 PetscCall(VecMAXPY(S1, 4, scoeff, Ss));
256 }
257 PetscFunctionReturn(PETSC_SUCCESS);
258 }
259
260 /*
261 X - initial state, updated in-place.
262 F - residual, computed at the initial X on input
263 */
SNESMSStep_Basic(SNES snes,Vec X,Vec F)264 static PetscErrorCode SNESMSStep_Basic(SNES snes, Vec X, Vec F)
265 {
266 SNES_MS *ms = (SNES_MS *)snes->data;
267 SNESMSTableau tab = ms->tableau;
268 const PetscReal *alpha = tab->betasub, h = ms->damping;
269 PetscInt i, nstages = tab->nstages;
270 Vec X0 = snes->work[0];
271
272 PetscFunctionBegin;
273 PetscCall(VecCopy(X, X0));
274 for (i = 0; i < nstages; i++) {
275 if (i > 0) PetscCall(SNESComputeFunction(snes, X, F));
276 PetscCall(KSPSolve(snes->ksp, F, X));
277 PetscCall(VecAYPX(X, -alpha[i] * h, X0));
278 }
279 PetscFunctionReturn(PETSC_SUCCESS);
280 }
281
SNESMSStep_Step(SNES snes,Vec X,Vec F)282 static PetscErrorCode SNESMSStep_Step(SNES snes, Vec X, Vec F)
283 {
284 SNES_MS *ms = (SNES_MS *)snes->data;
285 SNESMSTableau tab = ms->tableau;
286
287 PetscFunctionBegin;
288 if (tab->gamma && tab->delta) {
289 PetscCall(SNESMSStep_3Sstar(snes, X, F));
290 } else {
291 PetscCall(SNESMSStep_Basic(snes, X, F));
292 }
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
SNESMSStep_Norms(SNES snes,PetscInt iter,Vec F)296 static PetscErrorCode SNESMSStep_Norms(SNES snes, PetscInt iter, Vec F)
297 {
298 PetscReal fnorm;
299
300 PetscFunctionBegin;
301 if (SNESNeedNorm_Private(snes, iter)) {
302 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
303 SNESCheckFunctionDomainError(snes, fnorm);
304 /* Monitor convergence */
305 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
306 snes->iter = iter;
307 snes->norm = fnorm;
308 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
309 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
310 /* Test for convergence */
311 PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
312 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
313 } else if (iter > 0) {
314 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
315 snes->iter = iter;
316 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
317 }
318 PetscFunctionReturn(PETSC_SUCCESS);
319 }
320
SNESSolve_MS(SNES snes)321 static PetscErrorCode SNESSolve_MS(SNES snes)
322 {
323 Vec X = snes->vec_sol, F = snes->vec_func;
324 PetscInt i;
325
326 PetscFunctionBegin;
327 PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
328 PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
329
330 snes->reason = SNES_CONVERGED_ITERATING;
331 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
332 snes->iter = 0;
333 snes->norm = 0;
334 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
335
336 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
337 else snes->vec_func_init_set = PETSC_FALSE;
338
339 PetscCall(SNESMSStep_Norms(snes, 0, F));
340 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
341
342 for (i = 0; i < snes->max_its; i++) {
343 /* Call general purpose update function */
344 PetscTryTypeMethod(snes, update, snes->iter);
345
346 if (i == 0 && snes->jacobian) {
347 /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */
348 PetscCall(SNESComputeJacobian(snes, snes->vec_sol, snes->jacobian, snes->jacobian_pre));
349 SNESCheckJacobianDomainError(snes);
350 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
351 }
352
353 PetscCall(SNESMSStep_Step(snes, X, F));
354
355 if (i < snes->max_its - 1 || SNESNeedNorm_Private(snes, i + 1)) PetscCall(SNESComputeFunction(snes, X, F));
356
357 PetscCall(SNESMSStep_Norms(snes, i + 1, F));
358 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
359 }
360 PetscFunctionReturn(PETSC_SUCCESS);
361 }
362
SNESSetUp_MS(SNES snes)363 static PetscErrorCode SNESSetUp_MS(SNES snes)
364 {
365 SNES_MS *ms = (SNES_MS *)snes->data;
366 SNESMSTableau tab = ms->tableau;
367 PetscInt nwork = tab->nregisters + 1; // +1 because VecMAXPY() in SNESMSStep_3Sstar()
368 // needs an extra work vec
369
370 PetscFunctionBegin;
371 PetscCall(SNESSetWorkVecs(snes, nwork));
372 PetscCall(SNESSetUpMatrices(snes));
373 PetscFunctionReturn(PETSC_SUCCESS);
374 }
375
SNESDestroy_MS(SNES snes)376 static PetscErrorCode SNESDestroy_MS(SNES snes)
377 {
378 PetscFunctionBegin;
379 PetscCall(PetscFree(snes->data));
380 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", NULL));
381 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", NULL));
382 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", NULL));
383 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", NULL));
384 PetscFunctionReturn(PETSC_SUCCESS);
385 }
386
SNESView_MS(SNES snes,PetscViewer viewer)387 static PetscErrorCode SNESView_MS(SNES snes, PetscViewer viewer)
388 {
389 PetscBool isascii;
390 SNES_MS *ms = (SNES_MS *)snes->data;
391 SNESMSTableau tab = ms->tableau;
392
393 PetscFunctionBegin;
394 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
395 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " multi-stage method type: %s\n", tab->name));
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
SNESSetFromOptions_MS(SNES snes,PetscOptionItems PetscOptionsObject)399 static PetscErrorCode SNESSetFromOptions_MS(SNES snes, PetscOptionItems PetscOptionsObject)
400 {
401 PetscFunctionBegin;
402 PetscOptionsHeadBegin(PetscOptionsObject, "SNES MS options");
403 {
404 SNESMSTableauLink link;
405 PetscInt count, choice;
406 PetscBool flg;
407 const char **namelist;
408 SNESMSType mstype;
409 PetscReal damping;
410 PetscBool norms = PETSC_FALSE;
411
412 PetscCall(SNESMSGetType(snes, &mstype));
413 for (link = SNESMSTableauList, count = 0; link; link = link->next, count++);
414 PetscCall(PetscMalloc1(count, (char ***)&namelist));
415 for (link = SNESMSTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
416 PetscCall(PetscOptionsEList("-snes_ms_type", "Multistage smoother type", "SNESMSSetType", (const char *const *)namelist, count, mstype, &choice, &flg));
417 if (flg) PetscCall(SNESMSSetType(snes, namelist[choice]));
418 PetscCall(PetscFree(namelist));
419 PetscCall(SNESMSGetDamping(snes, &damping));
420 PetscCall(PetscOptionsReal("-snes_ms_damping", "Damping for multistage method", "SNESMSSetDamping", damping, &damping, &flg));
421 if (flg) PetscCall(SNESMSSetDamping(snes, damping));
422
423 /* deprecated option */
424 PetscCall(PetscOptionsDeprecated("-snes_ms_norms", NULL, "3.20", "Use -snes_norm_schedule always"));
425 PetscCall(PetscOptionsBool("-snes_ms_norms", NULL, NULL, norms, &norms, NULL));
426 if (norms) PetscCall(SNESSetNormSchedule(snes, SNES_NORM_ALWAYS));
427 }
428 PetscOptionsHeadEnd();
429 PetscFunctionReturn(PETSC_SUCCESS);
430 }
431
SNESMSGetType_MS(SNES snes,SNESMSType * mstype)432 static PetscErrorCode SNESMSGetType_MS(SNES snes, SNESMSType *mstype)
433 {
434 SNES_MS *ms = (SNES_MS *)snes->data;
435 SNESMSTableau tab = ms->tableau;
436
437 PetscFunctionBegin;
438 *mstype = tab->name;
439 PetscFunctionReturn(PETSC_SUCCESS);
440 }
441
SNESMSSetType_MS(SNES snes,SNESMSType mstype)442 static PetscErrorCode SNESMSSetType_MS(SNES snes, SNESMSType mstype)
443 {
444 SNES_MS *ms = (SNES_MS *)snes->data;
445 SNESMSTableauLink link;
446 PetscBool match;
447
448 PetscFunctionBegin;
449 if (ms->tableau) {
450 PetscCall(PetscStrcmp(ms->tableau->name, mstype, &match));
451 if (match) PetscFunctionReturn(PETSC_SUCCESS);
452 }
453 for (link = SNESMSTableauList; link; link = link->next) {
454 PetscCall(PetscStrcmp(link->tab.name, mstype, &match));
455 if (match) {
456 ms->tableau = &link->tab;
457 if (snes->setupcalled) PetscCall(SNESSetUp_MS(snes));
458 PetscFunctionReturn(PETSC_SUCCESS);
459 }
460 }
461 SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mstype);
462 }
463
464 /*@
465 SNESMSGetType - Get the type of multistage smoother `SNESMS`
466
467 Not Collective
468
469 Input Parameter:
470 . snes - nonlinear solver context
471
472 Output Parameter:
473 . mstype - type of multistage method
474
475 Level: advanced
476
477 .seealso: [](ch_snes), `SNESMS`, `SNESMSSetType()`, `SNESMSType`
478 @*/
SNESMSGetType(SNES snes,SNESMSType * mstype)479 PetscErrorCode SNESMSGetType(SNES snes, SNESMSType *mstype)
480 {
481 PetscFunctionBegin;
482 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
483 PetscAssertPointer(mstype, 2);
484 PetscUseMethod(snes, "SNESMSGetType_C", (SNES, SNESMSType *), (snes, mstype));
485 PetscFunctionReturn(PETSC_SUCCESS);
486 }
487
488 /*@
489 SNESMSSetType - Set the type of multistage smoother `SNESMS`
490
491 Logically Collective
492
493 Input Parameters:
494 + snes - nonlinear solver context
495 - mstype - type of multistage method
496
497 Level: advanced
498
499 .seealso: [](ch_snes), `SNESMS`, `SNESMSGetType()`, `SNESMSType`
500 @*/
SNESMSSetType(SNES snes,SNESMSType mstype)501 PetscErrorCode SNESMSSetType(SNES snes, SNESMSType mstype)
502 {
503 PetscFunctionBegin;
504 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
505 PetscAssertPointer(mstype, 2);
506 PetscTryMethod(snes, "SNESMSSetType_C", (SNES, SNESMSType), (snes, mstype));
507 PetscFunctionReturn(PETSC_SUCCESS);
508 }
509
SNESMSGetDamping_MS(SNES snes,PetscReal * damping)510 static PetscErrorCode SNESMSGetDamping_MS(SNES snes, PetscReal *damping)
511 {
512 SNES_MS *ms = (SNES_MS *)snes->data;
513
514 PetscFunctionBegin;
515 *damping = ms->damping;
516 PetscFunctionReturn(PETSC_SUCCESS);
517 }
518
SNESMSSetDamping_MS(SNES snes,PetscReal damping)519 static PetscErrorCode SNESMSSetDamping_MS(SNES snes, PetscReal damping)
520 {
521 SNES_MS *ms = (SNES_MS *)snes->data;
522
523 PetscFunctionBegin;
524 ms->damping = damping;
525 PetscFunctionReturn(PETSC_SUCCESS);
526 }
527
528 /*@
529 SNESMSGetDamping - Get the damping parameter of `SNESMS` multistage scheme
530
531 Not Collective
532
533 Input Parameter:
534 . snes - nonlinear solver context
535
536 Output Parameter:
537 . damping - damping parameter
538
539 Level: advanced
540
541 .seealso: [](ch_snes), `SNESMSSetDamping()`, `SNESMS`
542 @*/
SNESMSGetDamping(SNES snes,PetscReal * damping)543 PetscErrorCode SNESMSGetDamping(SNES snes, PetscReal *damping)
544 {
545 PetscFunctionBegin;
546 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
547 PetscAssertPointer(damping, 2);
548 PetscUseMethod(snes, "SNESMSGetDamping_C", (SNES, PetscReal *), (snes, damping));
549 PetscFunctionReturn(PETSC_SUCCESS);
550 }
551
552 /*@
553 SNESMSSetDamping - Set the damping parameter for a `SNESMS` multistage scheme
554
555 Logically Collective
556
557 Input Parameters:
558 + snes - nonlinear solver context
559 - damping - damping parameter
560
561 Level: advanced
562
563 .seealso: [](ch_snes), `SNESMSGetDamping()`, `SNESMS`
564 @*/
SNESMSSetDamping(SNES snes,PetscReal damping)565 PetscErrorCode SNESMSSetDamping(SNES snes, PetscReal damping)
566 {
567 PetscFunctionBegin;
568 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
569 PetscValidLogicalCollectiveReal(snes, damping, 2);
570 PetscTryMethod(snes, "SNESMSSetDamping_C", (SNES, PetscReal), (snes, damping));
571 PetscFunctionReturn(PETSC_SUCCESS);
572 }
573
574 /*MC
575 SNESMS - multi-stage smoothers
576
577 Options Database Keys:
578 + -snes_ms_type - type of multi-stage smoother
579 - -snes_ms_damping - damping for multi-stage method
580
581 Level: advanced
582
583 Notes:
584 These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for
585 FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebyshev).
586
587 Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds.
588
589 The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with `SNESMSRegister()`.
590
591 See {cite}`ketcheson2010runge`, {cite}`jameson1983`, {cite}`pierce1997preconditioned`, and {cite}`va1981design`
592
593 .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESMS`, `SNESFAS`, `KSPCHEBYSHEV`, `SNESMSSetDamping()`, `SNESMSGetDamping()`,
594 `SNESMSSetType()`, `SNESMSGetType()`
595 M*/
596
SNESCreate_MS(SNES snes)597 PETSC_EXTERN PetscErrorCode SNESCreate_MS(SNES snes)
598 {
599 SNES_MS *ms;
600
601 PetscFunctionBegin;
602 PetscCall(SNESMSInitializePackage());
603
604 snes->ops->setup = SNESSetUp_MS;
605 snes->ops->solve = SNESSolve_MS;
606 snes->ops->destroy = SNESDestroy_MS;
607 snes->ops->setfromoptions = SNESSetFromOptions_MS;
608 snes->ops->view = SNESView_MS;
609
610 snes->usesnpc = PETSC_FALSE;
611 snes->usesksp = PETSC_TRUE;
612
613 snes->alwayscomputesfinalresidual = PETSC_FALSE;
614
615 PetscCall(SNESParametersInitialize(snes));
616
617 PetscCall(PetscNew(&ms));
618 snes->data = (void *)ms;
619 ms->damping = 0.9;
620
621 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetType_C", SNESMSGetType_MS));
622 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetType_C", SNESMSSetType_MS));
623 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSGetDamping_C", SNESMSGetDamping_MS));
624 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMSSetDamping_C", SNESMSSetDamping_MS));
625
626 PetscCall(SNESMSSetType(snes, SNESMSDefault));
627 PetscFunctionReturn(PETSC_SUCCESS);
628 }
629