Actual source code: ex31.c
1: static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";
3: /*
5: Useful command line parameters:
6: -problem <hull1972a1>: choose which problem to solve (see references
7: for complete listing of problems).
8: -ts_type <euler>: specify time-integrator
9: -ts_adapt_type <basic>: specify time-step adapting (none,basic,advanced)
10: -refinement_levels <1>: number of refinement levels for convergence analysis
11: -refinement_factor <2.0>: factor to refine time step size by for convergence analysis
12: -dt <0.01>: specify time step (initial time step for convergence analysis)
14: */
16: /*
17: List of cases and their names in the code:-
18: From Hull, T.E., Enright, W.H., Fellen, B.M., and Sedgwick, A.E.,
19: "Comparing Numerical Methods for Ordinary Differential
20: Equations", SIAM J. Numer. Anal., 9(4), 1972, pp. 603 - 635
21: A1 -> "hull1972a1" (exact solution available)
22: A2 -> "hull1972a2" (exact solution available)
23: A3 -> "hull1972a3" (exact solution available)
24: A4 -> "hull1972a4" (exact solution available)
25: A5 -> "hull1972a5"
26: B1 -> "hull1972b1"
27: B2 -> "hull1972b2"
28: B3 -> "hull1972b3"
29: B4 -> "hull1972b4"
30: B5 -> "hull1972b5"
31: C1 -> "hull1972c1"
32: C2 -> "hull1972c2"
33: C3 -> "hull1972c3"
34: C4 -> "hull1972c4"
36: From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
37: https://arxiv.org/abs/1503.05166, 2016
39: Kulikov2013I -> "kulik2013i"
41: */
43: #include <petscts.h>
45: /* Function declarations */
46: PetscErrorCode (*RHSFunction)(TS, PetscReal, Vec, Vec, void *);
47: PetscErrorCode (*RHSJacobian)(TS, PetscReal, Vec, Mat, Mat, void *);
48: PetscErrorCode (*IFunction)(TS, PetscReal, Vec, Vec, Vec, void *);
49: PetscErrorCode (*IJacobian)(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
51: /* Returns the size of the system of equations depending on problem specification */
52: PetscErrorCode GetSize(const char *p, PetscInt *sz)
53: {
54: PetscFunctionBeginUser;
56: if (!strcmp(p, "hull1972a1") || !strcmp(p, "hull1972a2") || !strcmp(p, "hull1972a3") || !strcmp(p, "hull1972a4") || !strcmp(p, "hull1972a5")) *sz = 1;
57: else if (!strcmp(p, "hull1972b1")) *sz = 2;
58: else if (!strcmp(p, "hull1972b2") || !strcmp(p, "hull1972b3") || !strcmp(p, "hull1972b4") || !strcmp(p, "hull1972b5")) *sz = 3;
59: else if (!strcmp(p, "kulik2013i")) *sz = 4;
60: else if (!strcmp(p, "hull1972c1") || !strcmp(p, "hull1972c2") || !strcmp(p, "hull1972c3")) *sz = 10;
61: else if (!strcmp(p, "hull1972c4")) *sz = 52;
62: else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Incorrect problem name %s", p);
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /****************************************************************/
68: /* Problem specific functions */
70: /* Hull, 1972, Problem A1 */
72: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
73: {
74: PetscScalar *f;
75: const PetscScalar *y;
77: PetscFunctionBeginUser;
78: PetscCall(VecGetArrayRead(Y, &y));
79: PetscCall(VecGetArray(F, &f));
80: f[0] = -y[0];
81: PetscCall(VecRestoreArrayRead(Y, &y));
82: PetscCall(VecRestoreArray(F, &f));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
87: {
88: const PetscScalar *y;
89: PetscInt row = 0, col = 0;
90: PetscScalar value = -1.0;
92: PetscFunctionBeginUser;
93: PetscCall(VecGetArrayRead(Y, &y));
94: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
95: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
96: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
97: PetscCall(VecRestoreArrayRead(Y, &y));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
102: {
103: const PetscScalar *y;
104: PetscScalar *f;
106: PetscFunctionBeginUser;
107: PetscCall(VecGetArrayRead(Y, &y));
108: PetscCall(VecGetArray(F, &f));
109: f[0] = -y[0];
110: PetscCall(VecRestoreArrayRead(Y, &y));
111: PetscCall(VecRestoreArray(F, &f));
112: /* Left hand side = ydot - f(y) */
113: PetscCall(VecAYPX(F, -1.0, Ydot));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
118: {
119: const PetscScalar *y;
120: PetscInt row = 0, col = 0;
121: PetscScalar value = a - 1.0;
123: PetscFunctionBeginUser;
124: PetscCall(VecGetArrayRead(Y, &y));
125: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
126: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
128: PetscCall(VecRestoreArrayRead(Y, &y));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /* Hull, 1972, Problem A2 */
134: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
135: {
136: const PetscScalar *y;
137: PetscScalar *f;
139: PetscFunctionBeginUser;
140: PetscCall(VecGetArrayRead(Y, &y));
141: PetscCall(VecGetArray(F, &f));
142: f[0] = -0.5 * y[0] * y[0] * y[0];
143: PetscCall(VecRestoreArrayRead(Y, &y));
144: PetscCall(VecRestoreArray(F, &f));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
149: {
150: const PetscScalar *y;
151: PetscInt row = 0, col = 0;
152: PetscScalar value;
154: PetscFunctionBeginUser;
155: PetscCall(VecGetArrayRead(Y, &y));
156: value = -0.5 * 3.0 * y[0] * y[0];
157: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
158: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
159: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
160: PetscCall(VecRestoreArrayRead(Y, &y));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
165: {
166: PetscScalar *f;
167: const PetscScalar *y;
169: PetscFunctionBeginUser;
170: PetscCall(VecGetArrayRead(Y, &y));
171: PetscCall(VecGetArray(F, &f));
172: f[0] = -0.5 * y[0] * y[0] * y[0];
173: PetscCall(VecRestoreArrayRead(Y, &y));
174: PetscCall(VecRestoreArray(F, &f));
175: /* Left hand side = ydot - f(y) */
176: PetscCall(VecAYPX(F, -1.0, Ydot));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
181: {
182: const PetscScalar *y;
183: PetscInt row = 0, col = 0;
184: PetscScalar value;
186: PetscFunctionBeginUser;
187: PetscCall(VecGetArrayRead(Y, &y));
188: value = a + 0.5 * 3.0 * y[0] * y[0];
189: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
190: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
191: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
192: PetscCall(VecRestoreArrayRead(Y, &y));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: /* Hull, 1972, Problem A3 */
198: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
199: {
200: const PetscScalar *y;
201: PetscScalar *f;
203: PetscFunctionBeginUser;
204: PetscCall(VecGetArrayRead(Y, &y));
205: PetscCall(VecGetArray(F, &f));
206: f[0] = y[0] * PetscCosReal(t);
207: PetscCall(VecRestoreArrayRead(Y, &y));
208: PetscCall(VecRestoreArray(F, &f));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
213: {
214: const PetscScalar *y;
215: PetscInt row = 0, col = 0;
216: PetscScalar value = PetscCosReal(t);
218: PetscFunctionBeginUser;
219: PetscCall(VecGetArrayRead(Y, &y));
220: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
221: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
222: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
223: PetscCall(VecRestoreArrayRead(Y, &y));
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
228: {
229: PetscScalar *f;
230: const PetscScalar *y;
232: PetscFunctionBeginUser;
233: PetscCall(VecGetArrayRead(Y, &y));
234: PetscCall(VecGetArray(F, &f));
235: f[0] = y[0] * PetscCosReal(t);
236: PetscCall(VecRestoreArrayRead(Y, &y));
237: PetscCall(VecRestoreArray(F, &f));
238: /* Left hand side = ydot - f(y) */
239: PetscCall(VecAYPX(F, -1.0, Ydot));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
244: {
245: const PetscScalar *y;
246: PetscInt row = 0, col = 0;
247: PetscScalar value = a - PetscCosReal(t);
249: PetscFunctionBeginUser;
250: PetscCall(VecGetArrayRead(Y, &y));
251: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
252: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
253: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
254: PetscCall(VecRestoreArrayRead(Y, &y));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /* Hull, 1972, Problem A4 */
260: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
261: {
262: PetscScalar *f;
263: const PetscScalar *y;
265: PetscFunctionBeginUser;
266: PetscCall(VecGetArrayRead(Y, &y));
267: PetscCall(VecGetArray(F, &f));
268: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
269: PetscCall(VecRestoreArrayRead(Y, &y));
270: PetscCall(VecRestoreArray(F, &f));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
275: {
276: const PetscScalar *y;
277: PetscInt row = 0, col = 0;
278: PetscScalar value;
280: PetscFunctionBeginUser;
281: PetscCall(VecGetArrayRead(Y, &y));
282: value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05;
283: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
284: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
285: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
286: PetscCall(VecRestoreArrayRead(Y, &y));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
291: {
292: PetscScalar *f;
293: const PetscScalar *y;
295: PetscFunctionBeginUser;
296: PetscCall(VecGetArrayRead(Y, &y));
297: PetscCall(VecGetArray(F, &f));
298: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
299: PetscCall(VecRestoreArrayRead(Y, &y));
300: PetscCall(VecRestoreArray(F, &f));
301: /* Left hand side = ydot - f(y) */
302: PetscCall(VecAYPX(F, -1.0, Ydot));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
307: {
308: const PetscScalar *y;
309: PetscInt row = 0, col = 0;
310: PetscScalar value;
312: PetscFunctionBeginUser;
313: PetscCall(VecGetArrayRead(Y, &y));
314: value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05;
315: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
316: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
317: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
318: PetscCall(VecRestoreArrayRead(Y, &y));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /* Hull, 1972, Problem A5 */
324: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
325: {
326: PetscScalar *f;
327: const PetscScalar *y;
329: PetscFunctionBeginUser;
330: PetscCall(VecGetArrayRead(Y, &y));
331: PetscCall(VecGetArray(F, &f));
332: f[0] = (y[0] - t) / (y[0] + t);
333: PetscCall(VecRestoreArrayRead(Y, &y));
334: PetscCall(VecRestoreArray(F, &f));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
339: {
340: const PetscScalar *y;
341: PetscInt row = 0, col = 0;
342: PetscScalar value;
344: PetscFunctionBeginUser;
345: PetscCall(VecGetArrayRead(Y, &y));
346: value = 2 * t / ((t + y[0]) * (t + y[0]));
347: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
348: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
349: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
350: PetscCall(VecRestoreArrayRead(Y, &y));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
355: {
356: PetscScalar *f;
357: const PetscScalar *y;
359: PetscFunctionBeginUser;
360: PetscCall(VecGetArrayRead(Y, &y));
361: PetscCall(VecGetArray(F, &f));
362: f[0] = (y[0] - t) / (y[0] + t);
363: PetscCall(VecRestoreArrayRead(Y, &y));
364: PetscCall(VecRestoreArray(F, &f));
365: /* Left hand side = ydot - f(y) */
366: PetscCall(VecAYPX(F, -1.0, Ydot));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
371: {
372: const PetscScalar *y;
373: PetscInt row = 0, col = 0;
374: PetscScalar value;
376: PetscFunctionBeginUser;
377: PetscCall(VecGetArrayRead(Y, &y));
378: value = a - 2 * t / ((t + y[0]) * (t + y[0]));
379: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
380: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
381: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
382: PetscCall(VecRestoreArrayRead(Y, &y));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /* Hull, 1972, Problem B1 */
388: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
389: {
390: PetscScalar *f;
391: const PetscScalar *y;
393: PetscFunctionBeginUser;
394: PetscCall(VecGetArrayRead(Y, &y));
395: PetscCall(VecGetArray(F, &f));
396: f[0] = 2.0 * (y[0] - y[0] * y[1]);
397: f[1] = -(y[1] - y[0] * y[1]);
398: PetscCall(VecRestoreArrayRead(Y, &y));
399: PetscCall(VecRestoreArray(F, &f));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
404: {
405: PetscScalar *f;
406: const PetscScalar *y;
408: PetscFunctionBeginUser;
409: PetscCall(VecGetArrayRead(Y, &y));
410: PetscCall(VecGetArray(F, &f));
411: f[0] = 2.0 * (y[0] - y[0] * y[1]);
412: f[1] = -(y[1] - y[0] * y[1]);
413: PetscCall(VecRestoreArrayRead(Y, &y));
414: PetscCall(VecRestoreArray(F, &f));
415: /* Left hand side = ydot - f(y) */
416: PetscCall(VecAYPX(F, -1.0, Ydot));
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
421: {
422: const PetscScalar *y;
423: PetscInt row[2] = {0, 1};
424: PetscScalar value[2][2];
426: PetscFunctionBeginUser;
427: PetscCall(VecGetArrayRead(Y, &y));
428: value[0][0] = a - 2.0 * (1.0 - y[1]);
429: value[0][1] = 2.0 * y[0];
430: value[1][0] = -y[1];
431: value[1][1] = a + 1.0 - y[0];
432: PetscCall(MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES));
433: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
434: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
435: PetscCall(VecRestoreArrayRead(Y, &y));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /* Hull, 1972, Problem B2 */
441: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
442: {
443: PetscScalar *f;
444: const PetscScalar *y;
446: PetscFunctionBeginUser;
447: PetscCall(VecGetArrayRead(Y, &y));
448: PetscCall(VecGetArray(F, &f));
449: f[0] = -y[0] + y[1];
450: f[1] = y[0] - 2.0 * y[1] + y[2];
451: f[2] = y[1] - y[2];
452: PetscCall(VecRestoreArrayRead(Y, &y));
453: PetscCall(VecRestoreArray(F, &f));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
458: {
459: PetscScalar *f;
460: const PetscScalar *y;
462: PetscFunctionBeginUser;
463: PetscCall(VecGetArrayRead(Y, &y));
464: PetscCall(VecGetArray(F, &f));
465: f[0] = -y[0] + y[1];
466: f[1] = y[0] - 2.0 * y[1] + y[2];
467: f[2] = y[1] - y[2];
468: PetscCall(VecRestoreArrayRead(Y, &y));
469: PetscCall(VecRestoreArray(F, &f));
470: /* Left hand side = ydot - f(y) */
471: PetscCall(VecAYPX(F, -1.0, Ydot));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
476: {
477: const PetscScalar *y;
478: PetscInt row[3] = {0, 1, 2};
479: PetscScalar value[3][3];
481: PetscFunctionBeginUser;
482: PetscCall(VecGetArrayRead(Y, &y));
483: value[0][0] = a + 1.0;
484: value[0][1] = -1.0;
485: value[0][2] = 0;
486: value[1][0] = -1.0;
487: value[1][1] = a + 2.0;
488: value[1][2] = -1.0;
489: value[2][0] = 0;
490: value[2][1] = -1.0;
491: value[2][2] = a + 1.0;
492: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
493: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
494: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
495: PetscCall(VecRestoreArrayRead(Y, &y));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /* Hull, 1972, Problem B3 */
501: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
502: {
503: PetscScalar *f;
504: const PetscScalar *y;
506: PetscFunctionBeginUser;
507: PetscCall(VecGetArrayRead(Y, &y));
508: PetscCall(VecGetArray(F, &f));
509: f[0] = -y[0];
510: f[1] = y[0] - y[1] * y[1];
511: f[2] = y[1] * y[1];
512: PetscCall(VecRestoreArrayRead(Y, &y));
513: PetscCall(VecRestoreArray(F, &f));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
518: {
519: PetscScalar *f;
520: const PetscScalar *y;
522: PetscFunctionBeginUser;
523: PetscCall(VecGetArrayRead(Y, &y));
524: PetscCall(VecGetArray(F, &f));
525: f[0] = -y[0];
526: f[1] = y[0] - y[1] * y[1];
527: f[2] = y[1] * y[1];
528: PetscCall(VecRestoreArrayRead(Y, &y));
529: PetscCall(VecRestoreArray(F, &f));
530: /* Left hand side = ydot - f(y) */
531: PetscCall(VecAYPX(F, -1.0, Ydot));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
536: {
537: const PetscScalar *y;
538: PetscInt row[3] = {0, 1, 2};
539: PetscScalar value[3][3];
541: PetscFunctionBeginUser;
542: PetscCall(VecGetArrayRead(Y, &y));
543: value[0][0] = a + 1.0;
544: value[0][1] = 0;
545: value[0][2] = 0;
546: value[1][0] = -1.0;
547: value[1][1] = a + 2.0 * y[1];
548: value[1][2] = 0;
549: value[2][0] = 0;
550: value[2][1] = -2.0 * y[1];
551: value[2][2] = a;
552: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
553: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
554: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
555: PetscCall(VecRestoreArrayRead(Y, &y));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /* Hull, 1972, Problem B4 */
561: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
562: {
563: PetscScalar *f;
564: const PetscScalar *y;
566: PetscFunctionBeginUser;
567: PetscCall(VecGetArrayRead(Y, &y));
568: PetscCall(VecGetArray(F, &f));
569: f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
570: f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
571: f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
572: PetscCall(VecRestoreArrayRead(Y, &y));
573: PetscCall(VecRestoreArray(F, &f));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
578: {
579: PetscScalar *f;
580: const PetscScalar *y;
582: PetscFunctionBeginUser;
583: PetscCall(VecGetArrayRead(Y, &y));
584: PetscCall(VecGetArray(F, &f));
585: f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
586: f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
587: f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
588: PetscCall(VecRestoreArrayRead(Y, &y));
589: PetscCall(VecRestoreArray(F, &f));
590: /* Left hand side = ydot - f(y) */
591: PetscCall(VecAYPX(F, -1.0, Ydot));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
596: {
597: const PetscScalar *y;
598: PetscInt row[3] = {0, 1, 2};
599: PetscScalar value[3][3], fac, fac2;
601: PetscFunctionBeginUser;
602: PetscCall(VecGetArrayRead(Y, &y));
603: fac = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -1.5);
604: fac2 = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -0.5);
605: value[0][0] = a + (y[1] * y[1] * y[2]) * fac;
606: value[0][1] = 1.0 - (y[0] * y[1] * y[2]) * fac;
607: value[0][2] = y[0] * fac2;
608: value[1][0] = -1.0 - y[0] * y[1] * y[2] * fac;
609: value[1][1] = a + y[0] * y[0] * y[2] * fac;
610: value[1][2] = y[1] * fac2;
611: value[2][0] = -y[1] * y[1] * fac;
612: value[2][1] = y[0] * y[1] * fac;
613: value[2][2] = a;
614: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
615: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
616: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
617: PetscCall(VecRestoreArrayRead(Y, &y));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /* Hull, 1972, Problem B5 */
623: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
624: {
625: PetscScalar *f;
626: const PetscScalar *y;
628: PetscFunctionBeginUser;
629: PetscCall(VecGetArrayRead(Y, &y));
630: PetscCall(VecGetArray(F, &f));
631: f[0] = y[1] * y[2];
632: f[1] = -y[0] * y[2];
633: f[2] = -0.51 * y[0] * y[1];
634: PetscCall(VecRestoreArrayRead(Y, &y));
635: PetscCall(VecRestoreArray(F, &f));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
640: {
641: PetscScalar *f;
642: const PetscScalar *y;
644: PetscFunctionBeginUser;
645: PetscCall(VecGetArrayRead(Y, &y));
646: PetscCall(VecGetArray(F, &f));
647: f[0] = y[1] * y[2];
648: f[1] = -y[0] * y[2];
649: f[2] = -0.51 * y[0] * y[1];
650: PetscCall(VecRestoreArrayRead(Y, &y));
651: PetscCall(VecRestoreArray(F, &f));
652: /* Left hand side = ydot - f(y) */
653: PetscCall(VecAYPX(F, -1.0, Ydot));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
658: {
659: const PetscScalar *y;
660: PetscInt row[3] = {0, 1, 2};
661: PetscScalar value[3][3];
663: PetscFunctionBeginUser;
664: PetscCall(VecGetArrayRead(Y, &y));
665: value[0][0] = a;
666: value[0][1] = -y[2];
667: value[0][2] = -y[1];
668: value[1][0] = y[2];
669: value[1][1] = a;
670: value[1][2] = y[0];
671: value[2][0] = 0.51 * y[1];
672: value[2][1] = 0.51 * y[0];
673: value[2][2] = a;
674: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
675: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
676: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
677: PetscCall(VecRestoreArrayRead(Y, &y));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /* Kulikov, 2013, Problem I */
683: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
684: {
685: PetscScalar *f;
686: const PetscScalar *y;
688: PetscFunctionBeginUser;
689: PetscCall(VecGetArrayRead(Y, &y));
690: PetscCall(VecGetArray(F, &f));
691: f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
692: f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
693: f[2] = 2. * t * y[3];
694: f[3] = -2. * t * PetscLogScalar(y[0]);
695: PetscCall(VecRestoreArrayRead(Y, &y));
696: PetscCall(VecRestoreArray(F, &f));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
701: {
702: const PetscScalar *y;
703: PetscInt row[4] = {0, 1, 2, 3};
704: PetscScalar value[4][4];
705: PetscScalar m1, m2;
706: PetscFunctionBeginUser;
707: PetscCall(VecGetArrayRead(Y, &y));
708: m1 = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
709: m2 = 2. * t * PetscPowScalar(y[1], 1. / 5.);
710: value[0][0] = 0.;
711: value[0][1] = m1;
712: value[0][2] = 0.;
713: value[0][3] = m2;
714: m1 = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
715: m2 = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
716: value[1][0] = 0.;
717: value[1][1] = 0.;
718: value[1][2] = m1;
719: value[1][3] = m2;
720: value[2][0] = 0.;
721: value[2][1] = 0.;
722: value[2][2] = 0.;
723: value[2][3] = 2 * t;
724: value[3][0] = -2. * t / y[0];
725: value[3][1] = 0.;
726: value[3][2] = 0.;
727: value[3][3] = 0.;
728: PetscCall(MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES));
729: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
730: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
731: PetscCall(VecRestoreArrayRead(Y, &y));
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
736: {
737: PetscScalar *f;
738: const PetscScalar *y;
740: PetscFunctionBeginUser;
741: PetscCall(VecGetArrayRead(Y, &y));
742: PetscCall(VecGetArray(F, &f));
743: f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
744: f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
745: f[2] = 2. * t * y[3];
746: f[3] = -2. * t * PetscLogScalar(y[0]);
747: PetscCall(VecRestoreArrayRead(Y, &y));
748: PetscCall(VecRestoreArray(F, &f));
749: /* Left hand side = ydot - f(y) */
750: PetscCall(VecAYPX(F, -1.0, Ydot));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
755: {
756: const PetscScalar *y;
757: PetscInt row[4] = {0, 1, 2, 3};
758: PetscScalar value[4][4];
759: PetscScalar m1, m2;
761: PetscFunctionBeginUser;
762: PetscCall(VecGetArrayRead(Y, &y));
763: m1 = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.));
764: m2 = 2. * t * PetscPowScalar(y[1], 1. / 5.);
765: value[0][0] = a;
766: value[0][1] = m1;
767: value[0][2] = 0.;
768: value[0][3] = m2;
769: m1 = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
770: m2 = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.));
771: value[1][0] = 0.;
772: value[1][1] = a;
773: value[1][2] = m1;
774: value[1][3] = m2;
775: value[2][0] = 0.;
776: value[2][1] = 0.;
777: value[2][2] = a;
778: value[2][3] = 2 * t;
779: value[3][0] = -2. * t / y[0];
780: value[3][1] = 0.;
781: value[3][2] = 0.;
782: value[3][3] = a;
783: PetscCall(MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES));
784: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
785: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
786: PetscCall(VecRestoreArrayRead(Y, &y));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /* Hull, 1972, Problem C1 */
792: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
793: {
794: PetscScalar *f;
795: const PetscScalar *y;
796: PetscInt N, i;
798: PetscFunctionBeginUser;
799: PetscCall(VecGetSize(Y, &N));
800: PetscCall(VecGetArrayRead(Y, &y));
801: PetscCall(VecGetArray(F, &f));
802: f[0] = -y[0];
803: for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
804: f[N - 1] = y[N - 2];
805: PetscCall(VecRestoreArrayRead(Y, &y));
806: PetscCall(VecRestoreArray(F, &f));
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
811: {
812: PetscScalar *f;
813: const PetscScalar *y;
814: PetscInt N, i;
816: PetscFunctionBeginUser;
817: PetscCall(VecGetSize(Y, &N));
818: PetscCall(VecGetArrayRead(Y, &y));
819: PetscCall(VecGetArray(F, &f));
820: f[0] = -y[0];
821: for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - y[i];
822: f[N - 1] = y[N - 2];
823: PetscCall(VecRestoreArrayRead(Y, &y));
824: PetscCall(VecRestoreArray(F, &f));
825: /* Left hand side = ydot - f(y) */
826: PetscCall(VecAYPX(F, -1.0, Ydot));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
831: {
832: const PetscScalar *y;
833: PetscInt N, i, col[2];
834: PetscScalar value[2];
836: PetscFunctionBeginUser;
837: PetscCall(VecGetSize(Y, &N));
838: PetscCall(VecGetArrayRead(Y, &y));
839: i = 0;
840: value[0] = a + 1;
841: col[0] = 0;
842: value[1] = 0;
843: col[1] = 1;
844: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
845: for (i = 0; i < N; i++) {
846: value[0] = -1;
847: col[0] = i - 1;
848: value[1] = a + 1;
849: col[1] = i;
850: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
851: }
852: i = N - 1;
853: value[0] = -1;
854: col[0] = N - 2;
855: value[1] = a;
856: col[1] = N - 1;
857: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
858: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
859: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
860: PetscCall(VecRestoreArrayRead(Y, &y));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /* Hull, 1972, Problem C2 */
866: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
867: {
868: const PetscScalar *y;
869: PetscScalar *f;
870: PetscInt N, i;
872: PetscFunctionBeginUser;
873: PetscCall(VecGetSize(Y, &N));
874: PetscCall(VecGetArrayRead(Y, &y));
875: PetscCall(VecGetArray(F, &f));
876: f[0] = -y[0];
877: for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
878: f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
879: PetscCall(VecRestoreArrayRead(Y, &y));
880: PetscCall(VecRestoreArray(F, &f));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
885: {
886: PetscScalar *f;
887: const PetscScalar *y;
888: PetscInt N, i;
890: PetscFunctionBeginUser;
891: PetscCall(VecGetSize(Y, &N));
892: PetscCall(VecGetArrayRead(Y, &y));
893: PetscCall(VecGetArray(F, &f));
894: f[0] = -y[0];
895: for (i = 1; i < N - 1; i++) f[i] = (PetscReal)i * y[i - 1] - (PetscReal)(i + 1) * y[i];
896: f[N - 1] = (PetscReal)(N - 1) * y[N - 2];
897: PetscCall(VecRestoreArrayRead(Y, &y));
898: PetscCall(VecRestoreArray(F, &f));
899: /* Left hand side = ydot - f(y) */
900: PetscCall(VecAYPX(F, -1.0, Ydot));
901: PetscFunctionReturn(PETSC_SUCCESS);
902: }
904: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
905: {
906: const PetscScalar *y;
907: PetscInt N, i, col[2];
908: PetscScalar value[2];
910: PetscFunctionBeginUser;
911: PetscCall(VecGetSize(Y, &N));
912: PetscCall(VecGetArrayRead(Y, &y));
913: i = 0;
914: value[0] = a + 1;
915: col[0] = 0;
916: value[1] = 0;
917: col[1] = 1;
918: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
919: for (i = 0; i < N; i++) {
920: value[0] = -(PetscReal)i;
921: col[0] = i - 1;
922: value[1] = a + (PetscReal)(i + 1);
923: col[1] = i;
924: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
925: }
926: i = N - 1;
927: value[0] = -(PetscReal)(N - 1);
928: col[0] = N - 2;
929: value[1] = a;
930: col[1] = N - 1;
931: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
932: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
933: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
934: PetscCall(VecRestoreArrayRead(Y, &y));
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /* Hull, 1972, Problem C3 and C4 */
940: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
941: {
942: PetscScalar *f;
943: const PetscScalar *y;
944: PetscInt N, i;
946: PetscFunctionBeginUser;
947: PetscCall(VecGetSize(Y, &N));
948: PetscCall(VecGetArrayRead(Y, &y));
949: PetscCall(VecGetArray(F, &f));
950: f[0] = -2.0 * y[0] + y[1];
951: for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
952: f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
953: PetscCall(VecRestoreArrayRead(Y, &y));
954: PetscCall(VecRestoreArray(F, &f));
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
959: {
960: PetscScalar *f;
961: const PetscScalar *y;
962: PetscInt N, i;
964: PetscFunctionBeginUser;
965: PetscCall(VecGetSize(Y, &N));
966: PetscCall(VecGetArrayRead(Y, &y));
967: PetscCall(VecGetArray(F, &f));
968: f[0] = -2.0 * y[0] + y[1];
969: for (i = 1; i < N - 1; i++) f[i] = y[i - 1] - 2.0 * y[i] + y[i + 1];
970: f[N - 1] = y[N - 2] - 2.0 * y[N - 1];
971: PetscCall(VecRestoreArrayRead(Y, &y));
972: PetscCall(VecRestoreArray(F, &f));
973: /* Left hand side = ydot - f(y) */
974: PetscCall(VecAYPX(F, -1.0, Ydot));
975: PetscFunctionReturn(PETSC_SUCCESS);
976: }
978: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
979: {
980: const PetscScalar *y;
981: PetscScalar value[3];
982: PetscInt N, i, col[3];
984: PetscFunctionBeginUser;
985: PetscCall(VecGetSize(Y, &N));
986: PetscCall(VecGetArrayRead(Y, &y));
987: for (i = 0; i < N; i++) {
988: if (i == 0) {
989: value[0] = a + 2;
990: col[0] = i;
991: value[1] = -1;
992: col[1] = i + 1;
993: value[2] = 0;
994: col[2] = i + 2;
995: } else if (i == N - 1) {
996: value[0] = 0;
997: col[0] = i - 2;
998: value[1] = -1;
999: col[1] = i - 1;
1000: value[2] = a + 2;
1001: col[2] = i;
1002: } else {
1003: value[0] = -1;
1004: col[0] = i - 1;
1005: value[1] = a + 2;
1006: col[1] = i;
1007: value[2] = -1;
1008: col[2] = i + 1;
1009: }
1010: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
1011: }
1012: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1013: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1014: PetscCall(VecRestoreArrayRead(Y, &y));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /***************************************************************************/
1020: /* Sets the initial solution for the IVP and sets up the function pointers*/
1021: PetscErrorCode Initialize(Vec Y, void *s)
1022: {
1023: char *p = (char *)s;
1024: PetscScalar *y;
1025: PetscReal t0;
1026: PetscInt N;
1027: PetscBool flg;
1029: PetscFunctionBeginUser;
1030: PetscCall(GetSize((const char *)s, &N));
1031: PetscCall(VecZeroEntries(Y));
1032: PetscCall(VecGetArray(Y, &y));
1033: if (!strcmp(p, "hull1972a1")) {
1034: y[0] = 1.0;
1035: RHSFunction = RHSFunction_Hull1972A1;
1036: RHSJacobian = RHSJacobian_Hull1972A1;
1037: IFunction = IFunction_Hull1972A1;
1038: IJacobian = IJacobian_Hull1972A1;
1039: } else if (!strcmp(p, "hull1972a2")) {
1040: y[0] = 1.0;
1041: RHSFunction = RHSFunction_Hull1972A2;
1042: RHSJacobian = RHSJacobian_Hull1972A2;
1043: IFunction = IFunction_Hull1972A2;
1044: IJacobian = IJacobian_Hull1972A2;
1045: } else if (!strcmp(p, "hull1972a3")) {
1046: y[0] = 1.0;
1047: RHSFunction = RHSFunction_Hull1972A3;
1048: RHSJacobian = RHSJacobian_Hull1972A3;
1049: IFunction = IFunction_Hull1972A3;
1050: IJacobian = IJacobian_Hull1972A3;
1051: } else if (!strcmp(p, "hull1972a4")) {
1052: y[0] = 1.0;
1053: RHSFunction = RHSFunction_Hull1972A4;
1054: RHSJacobian = RHSJacobian_Hull1972A4;
1055: IFunction = IFunction_Hull1972A4;
1056: IJacobian = IJacobian_Hull1972A4;
1057: } else if (!strcmp(p, "hull1972a5")) {
1058: y[0] = 4.0;
1059: RHSFunction = RHSFunction_Hull1972A5;
1060: RHSJacobian = RHSJacobian_Hull1972A5;
1061: IFunction = IFunction_Hull1972A5;
1062: IJacobian = IJacobian_Hull1972A5;
1063: } else if (!strcmp(p, "hull1972b1")) {
1064: y[0] = 1.0;
1065: y[1] = 3.0;
1066: RHSFunction = RHSFunction_Hull1972B1;
1067: IFunction = IFunction_Hull1972B1;
1068: IJacobian = IJacobian_Hull1972B1;
1069: } else if (!strcmp(p, "hull1972b2")) {
1070: y[0] = 2.0;
1071: y[1] = 0.0;
1072: y[2] = 1.0;
1073: RHSFunction = RHSFunction_Hull1972B2;
1074: IFunction = IFunction_Hull1972B2;
1075: IJacobian = IJacobian_Hull1972B2;
1076: } else if (!strcmp(p, "hull1972b3")) {
1077: y[0] = 1.0;
1078: y[1] = 0.0;
1079: y[2] = 0.0;
1080: RHSFunction = RHSFunction_Hull1972B3;
1081: IFunction = IFunction_Hull1972B3;
1082: IJacobian = IJacobian_Hull1972B3;
1083: } else if (!strcmp(p, "hull1972b4")) {
1084: y[0] = 3.0;
1085: y[1] = 0.0;
1086: y[2] = 0.0;
1087: RHSFunction = RHSFunction_Hull1972B4;
1088: IFunction = IFunction_Hull1972B4;
1089: IJacobian = IJacobian_Hull1972B4;
1090: } else if (!strcmp(p, "hull1972b5")) {
1091: y[0] = 0.0;
1092: y[1] = 1.0;
1093: y[2] = 1.0;
1094: RHSFunction = RHSFunction_Hull1972B5;
1095: IFunction = IFunction_Hull1972B5;
1096: IJacobian = IJacobian_Hull1972B5;
1097: } else if (!strcmp(p, "kulik2013i")) {
1098: t0 = 0.;
1099: y[0] = PetscExpReal(PetscSinReal(t0 * t0));
1100: y[1] = PetscExpReal(5. * PetscSinReal(t0 * t0));
1101: y[2] = PetscSinReal(t0 * t0) + 1.0;
1102: y[3] = PetscCosReal(t0 * t0);
1103: RHSFunction = RHSFunction_Kulikov2013I;
1104: RHSJacobian = RHSJacobian_Kulikov2013I;
1105: IFunction = IFunction_Kulikov2013I;
1106: IJacobian = IJacobian_Kulikov2013I;
1107: } else if (!strcmp(p, "hull1972c1")) {
1108: y[0] = 1.0;
1109: RHSFunction = RHSFunction_Hull1972C1;
1110: IFunction = IFunction_Hull1972C1;
1111: IJacobian = IJacobian_Hull1972C1;
1112: } else if (!strcmp(p, "hull1972c2")) {
1113: y[0] = 1.0;
1114: RHSFunction = RHSFunction_Hull1972C2;
1115: IFunction = IFunction_Hull1972C2;
1116: IJacobian = IJacobian_Hull1972C2;
1117: } else if (!strcmp(p, "hull1972c3") || !strcmp(p, "hull1972c4")) {
1118: y[0] = 1.0;
1119: RHSFunction = RHSFunction_Hull1972C34;
1120: IFunction = IFunction_Hull1972C34;
1121: IJacobian = IJacobian_Hull1972C34;
1122: }
1123: PetscCall(PetscOptionsGetScalarArray(NULL, NULL, "-yinit", y, &N, &flg));
1124: PetscCall(VecRestoreArray(Y, &y));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: /* Calculates the exact solution to problems that have one */
1129: PetscErrorCode ExactSolution(Vec Y, void *s, PetscReal t, PetscBool *flag)
1130: {
1131: char *p = (char *)s;
1132: PetscScalar *y;
1134: PetscFunctionBeginUser;
1135: if (!strcmp(p, "hull1972a1")) {
1136: PetscCall(VecGetArray(Y, &y));
1137: y[0] = PetscExpReal(-t);
1138: *flag = PETSC_TRUE;
1139: PetscCall(VecRestoreArray(Y, &y));
1140: } else if (!strcmp(p, "hull1972a2")) {
1141: PetscCall(VecGetArray(Y, &y));
1142: y[0] = 1.0 / PetscSqrtReal(t + 1);
1143: *flag = PETSC_TRUE;
1144: PetscCall(VecRestoreArray(Y, &y));
1145: } else if (!strcmp(p, "hull1972a3")) {
1146: PetscCall(VecGetArray(Y, &y));
1147: y[0] = PetscExpReal(PetscSinReal(t));
1148: *flag = PETSC_TRUE;
1149: PetscCall(VecRestoreArray(Y, &y));
1150: } else if (!strcmp(p, "hull1972a4")) {
1151: PetscCall(VecGetArray(Y, &y));
1152: y[0] = 20.0 / (1 + 19.0 * PetscExpReal(-t / 4.0));
1153: *flag = PETSC_TRUE;
1154: PetscCall(VecRestoreArray(Y, &y));
1155: } else if (!strcmp(p, "kulik2013i")) {
1156: PetscCall(VecGetArray(Y, &y));
1157: y[0] = PetscExpReal(PetscSinReal(t * t));
1158: y[1] = PetscExpReal(5. * PetscSinReal(t * t));
1159: y[2] = PetscSinReal(t * t) + 1.0;
1160: y[3] = PetscCosReal(t * t);
1161: *flag = PETSC_TRUE;
1162: PetscCall(VecRestoreArray(Y, &y));
1163: } else {
1164: PetscCall(VecSet(Y, 0));
1165: *flag = PETSC_FALSE;
1166: }
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: /* Solves the specified ODE and computes the error if exact solution is available */
1171: PetscErrorCode SolveODE(char *ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1172: {
1173: TS ts; /* time-integrator */
1174: Vec Y; /* Solution vector */
1175: Vec Yex; /* Exact solution */
1176: PetscInt N; /* Size of the system of equations */
1177: TSType time_scheme; /* Type of time-integration scheme */
1178: Mat Jac = NULL; /* Jacobian matrix */
1179: Vec Yerr; /* Auxiliary solution vector */
1180: PetscReal err_norm; /* Estimated error norm */
1181: PetscReal final_time; /* Actual final time from the integrator */
1183: PetscFunctionBeginUser;
1184: PetscCall(GetSize((const char *)&ptype[0], &N));
1185: PetscCheck(N >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Illegal problem specification.");
1186: PetscCall(VecCreate(PETSC_COMM_WORLD, &Y));
1187: PetscCall(VecSetSizes(Y, N, PETSC_DECIDE));
1188: PetscCall(VecSetUp(Y));
1189: PetscCall(VecSet(Y, 0));
1191: /* Initialize the problem */
1192: PetscCall(Initialize(Y, &ptype[0]));
1194: /* Create and initialize the time-integrator */
1195: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1196: /* Default time integration options */
1197: PetscCall(TSSetType(ts, TSRK));
1198: PetscCall(TSSetMaxSteps(ts, maxiter));
1199: PetscCall(TSSetMaxTime(ts, tfinal));
1200: PetscCall(TSSetTimeStep(ts, dt));
1201: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1202: /* Read command line options for time integration */
1203: PetscCall(TSSetFromOptions(ts));
1204: /* Set solution vector */
1205: PetscCall(TSSetSolution(ts, Y));
1206: /* Specify left/right-hand side functions */
1207: PetscCall(TSGetType(ts, &time_scheme));
1209: if ((!strcmp(time_scheme, TSEULER)) || (!strcmp(time_scheme, TSRK)) || (!strcmp(time_scheme, TSSSP) || (!strcmp(time_scheme, TSGLEE)))) {
1210: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1211: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &ptype[0]));
1212: PetscCall(MatCreate(PETSC_COMM_WORLD, &Jac));
1213: PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
1214: PetscCall(MatSetFromOptions(Jac));
1215: PetscCall(MatSetUp(Jac));
1216: PetscCall(TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &ptype[0]));
1217: } else if ((!strcmp(time_scheme, TSTHETA)) || (!strcmp(time_scheme, TSBEULER)) || (!strcmp(time_scheme, TSCN)) || (!strcmp(time_scheme, TSALPHA)) || (!strcmp(time_scheme, TSARKIMEX))) {
1218: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1219: /* and its Jacobian function */
1220: PetscCall(TSSetIFunction(ts, NULL, IFunction, &ptype[0]));
1221: PetscCall(MatCreate(PETSC_COMM_WORLD, &Jac));
1222: PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
1223: PetscCall(MatSetFromOptions(Jac));
1224: PetscCall(MatSetUp(Jac));
1225: PetscCall(TSSetIJacobian(ts, Jac, Jac, IJacobian, &ptype[0]));
1226: }
1228: /* Solve */
1229: PetscCall(TSSolve(ts, Y));
1230: PetscCall(TSGetTime(ts, &final_time));
1232: /* Get the estimated error, if available */
1233: PetscCall(VecDuplicate(Y, &Yerr));
1234: PetscCall(VecZeroEntries(Yerr));
1235: PetscCall(TSGetTimeError(ts, 0, &Yerr));
1236: PetscCall(VecNorm(Yerr, NORM_2, &err_norm));
1237: PetscCall(VecDestroy(&Yerr));
1238: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Estimated Error = %e.\n", (double)err_norm));
1240: /* Exact solution */
1241: PetscCall(VecDuplicate(Y, &Yex));
1242: if (PetscAbsReal(final_time - tfinal) > 2. * PETSC_MACHINE_EPSILON) {
1243: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference between the prescribed final time %g and the actual final time, %g.\n", (double)tfinal, (double)final_time));
1244: }
1245: PetscCall(ExactSolution(Yex, &ptype[0], final_time, exact_flag));
1247: /* Calculate Error */
1248: PetscCall(VecAYPX(Yex, -1.0, Y));
1249: PetscCall(VecNorm(Yex, NORM_2, error));
1250: *error = PetscSqrtReal(((*error) * (*error)) / N);
1252: /* Clean up and finalize */
1253: PetscCall(MatDestroy(&Jac));
1254: PetscCall(TSDestroy(&ts));
1255: PetscCall(VecDestroy(&Yex));
1256: PetscCall(VecDestroy(&Y));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: int main(int argc, char **argv)
1262: {
1263: char ptype[256] = "hull1972a1"; /* Problem specification */
1264: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1265: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1266: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1267: PetscReal dt;
1268: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1269: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1270: PetscReal *error; /* Array to store the errors for convergence analysis */
1271: PetscMPIInt size; /* No of processors */
1272: PetscBool flag; /* Flag denoting availability of exact solution */
1273: PetscInt r, N;
1275: /* Initialize program */
1276: PetscFunctionBeginUser;
1277: PetscCall(GetSize(&ptype[0], &N));
1278: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1280: /* Check if running with only 1 proc */
1281: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1282: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
1284: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
1285: PetscCall(PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL));
1286: PetscCall(PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL));
1287: PetscCall(PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL));
1288: PetscCall(PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL));
1289: PetscCall(PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL));
1290: PetscOptionsEnd();
1292: PetscCall(PetscMalloc1(n_refine, &error));
1293: for (r = 0, dt = dt_initial; r < n_refine; r++) {
1294: error[r] = 0;
1295: if (r > 0) dt /= refine_fac;
1297: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving ODE \"%s\" with dt %f, final time %f and system size %" PetscInt_FMT ".\n", ptype, (double)dt, (double)tfinal, N));
1298: PetscCall(SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag));
1299: if (flag) {
1300: /* If exact solution available for the specified ODE */
1301: if (r > 0) {
1302: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac));
1303: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate));
1304: } else {
1305: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E.\n", (double)error[r]));
1306: }
1307: }
1308: }
1309: PetscCall(PetscFree(error));
1310: PetscCall(PetscFinalize());
1311: return 0;
1312: }
1314: /*TEST
1316: test:
1317: suffix: 2
1318: args: -ts_type glee -final_time 5 -ts_adapt_type none
1319: timeoutfactor: 3
1320: requires: !single
1322: test:
1323: suffix: 3
1324: args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor -ts_max_steps 50 -problem hull1972a3 -ts_adapt_glee_use_local 1
1325: timeoutfactor: 3
1326: requires: !single
1328: test:
1329: suffix: 4
1330: args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor -ts_max_steps 50 -problem hull1972a3 -ts_max_reject 100 -ts_adapt_glee_use_local 0
1331: timeoutfactor: 3
1332: requires: !single !__float128
1334: TEST*/