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;
55: if (!strcmp(p, "hull1972a1") || !strcmp(p, "hull1972a2") || !strcmp(p, "hull1972a3") || !strcmp(p, "hull1972a4") || !strcmp(p, "hull1972a5")) *sz = 1;
56: else if (!strcmp(p, "hull1972b1")) *sz = 2;
57: else if (!strcmp(p, "hull1972b2") || !strcmp(p, "hull1972b3") || !strcmp(p, "hull1972b4") || !strcmp(p, "hull1972b5")) *sz = 3;
58: else if (!strcmp(p, "kulik2013i")) *sz = 4;
59: else if (!strcmp(p, "hull1972c1") || !strcmp(p, "hull1972c2") || !strcmp(p, "hull1972c3")) *sz = 10;
60: else if (!strcmp(p, "hull1972c4")) *sz = 52;
61: else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Incorrect problem name %s", p);
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /****************************************************************/
67: /* Problem specific functions */
69: /* Hull, 1972, Problem A1 */
71: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
72: {
73: PetscScalar *f;
74: const PetscScalar *y;
76: PetscFunctionBeginUser;
77: PetscCall(VecGetArrayRead(Y, &y));
78: PetscCall(VecGetArray(F, &f));
79: f[0] = -y[0];
80: PetscCall(VecRestoreArrayRead(Y, &y));
81: PetscCall(VecRestoreArray(F, &f));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
86: {
87: const PetscScalar *y;
88: PetscInt row = 0, col = 0;
89: PetscScalar value = -1.0;
91: PetscFunctionBeginUser;
92: PetscCall(VecGetArrayRead(Y, &y));
93: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
94: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
95: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
96: PetscCall(VecRestoreArrayRead(Y, &y));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
101: {
102: const PetscScalar *y;
103: PetscScalar *f;
105: PetscFunctionBeginUser;
106: PetscCall(VecGetArrayRead(Y, &y));
107: PetscCall(VecGetArray(F, &f));
108: f[0] = -y[0];
109: PetscCall(VecRestoreArrayRead(Y, &y));
110: PetscCall(VecRestoreArray(F, &f));
111: /* Left hand side = ydot - f(y) */
112: PetscCall(VecAYPX(F, -1.0, Ydot));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
117: {
118: const PetscScalar *y;
119: PetscInt row = 0, col = 0;
120: PetscScalar value = a - 1.0;
122: PetscFunctionBeginUser;
123: PetscCall(VecGetArrayRead(Y, &y));
124: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
125: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(VecRestoreArrayRead(Y, &y));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /* Hull, 1972, Problem A2 */
133: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
134: {
135: const PetscScalar *y;
136: PetscScalar *f;
138: PetscFunctionBeginUser;
139: PetscCall(VecGetArrayRead(Y, &y));
140: PetscCall(VecGetArray(F, &f));
141: f[0] = -0.5 * y[0] * y[0] * y[0];
142: PetscCall(VecRestoreArrayRead(Y, &y));
143: PetscCall(VecRestoreArray(F, &f));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
148: {
149: const PetscScalar *y;
150: PetscInt row = 0, col = 0;
151: PetscScalar value;
153: PetscFunctionBeginUser;
154: PetscCall(VecGetArrayRead(Y, &y));
155: value = -0.5 * 3.0 * y[0] * y[0];
156: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
157: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
158: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
159: PetscCall(VecRestoreArrayRead(Y, &y));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
164: {
165: PetscScalar *f;
166: const PetscScalar *y;
168: PetscFunctionBeginUser;
169: PetscCall(VecGetArrayRead(Y, &y));
170: PetscCall(VecGetArray(F, &f));
171: f[0] = -0.5 * y[0] * y[0] * y[0];
172: PetscCall(VecRestoreArrayRead(Y, &y));
173: PetscCall(VecRestoreArray(F, &f));
174: /* Left hand side = ydot - f(y) */
175: PetscCall(VecAYPX(F, -1.0, Ydot));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
180: {
181: const PetscScalar *y;
182: PetscInt row = 0, col = 0;
183: PetscScalar value;
185: PetscFunctionBeginUser;
186: PetscCall(VecGetArrayRead(Y, &y));
187: value = a + 0.5 * 3.0 * y[0] * y[0];
188: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
189: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
190: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
191: PetscCall(VecRestoreArrayRead(Y, &y));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: /* Hull, 1972, Problem A3 */
197: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
198: {
199: const PetscScalar *y;
200: PetscScalar *f;
202: PetscFunctionBeginUser;
203: PetscCall(VecGetArrayRead(Y, &y));
204: PetscCall(VecGetArray(F, &f));
205: f[0] = y[0] * PetscCosReal(t);
206: PetscCall(VecRestoreArrayRead(Y, &y));
207: PetscCall(VecRestoreArray(F, &f));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
212: {
213: const PetscScalar *y;
214: PetscInt row = 0, col = 0;
215: PetscScalar value = PetscCosReal(t);
217: PetscFunctionBeginUser;
218: PetscCall(VecGetArrayRead(Y, &y));
219: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
220: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
221: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
222: PetscCall(VecRestoreArrayRead(Y, &y));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
227: {
228: PetscScalar *f;
229: const PetscScalar *y;
231: PetscFunctionBeginUser;
232: PetscCall(VecGetArrayRead(Y, &y));
233: PetscCall(VecGetArray(F, &f));
234: f[0] = y[0] * PetscCosReal(t);
235: PetscCall(VecRestoreArrayRead(Y, &y));
236: PetscCall(VecRestoreArray(F, &f));
237: /* Left hand side = ydot - f(y) */
238: PetscCall(VecAYPX(F, -1.0, Ydot));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
243: {
244: const PetscScalar *y;
245: PetscInt row = 0, col = 0;
246: PetscScalar value = a - PetscCosReal(t);
248: PetscFunctionBeginUser;
249: PetscCall(VecGetArrayRead(Y, &y));
250: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
251: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
252: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
253: PetscCall(VecRestoreArrayRead(Y, &y));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /* Hull, 1972, Problem A4 */
259: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
260: {
261: PetscScalar *f;
262: const PetscScalar *y;
264: PetscFunctionBeginUser;
265: PetscCall(VecGetArrayRead(Y, &y));
266: PetscCall(VecGetArray(F, &f));
267: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
268: PetscCall(VecRestoreArrayRead(Y, &y));
269: PetscCall(VecRestoreArray(F, &f));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
274: {
275: const PetscScalar *y;
276: PetscInt row = 0, col = 0;
277: PetscScalar value;
279: PetscFunctionBeginUser;
280: PetscCall(VecGetArrayRead(Y, &y));
281: value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05;
282: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
283: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
284: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
285: PetscCall(VecRestoreArrayRead(Y, &y));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
290: {
291: PetscScalar *f;
292: const PetscScalar *y;
294: PetscFunctionBeginUser;
295: PetscCall(VecGetArrayRead(Y, &y));
296: PetscCall(VecGetArray(F, &f));
297: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
298: PetscCall(VecRestoreArrayRead(Y, &y));
299: PetscCall(VecRestoreArray(F, &f));
300: /* Left hand side = ydot - f(y) */
301: PetscCall(VecAYPX(F, -1.0, Ydot));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
306: {
307: const PetscScalar *y;
308: PetscInt row = 0, col = 0;
309: PetscScalar value;
311: PetscFunctionBeginUser;
312: PetscCall(VecGetArrayRead(Y, &y));
313: value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05;
314: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
315: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
316: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
317: PetscCall(VecRestoreArrayRead(Y, &y));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /* Hull, 1972, Problem A5 */
323: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
324: {
325: PetscScalar *f;
326: const PetscScalar *y;
328: PetscFunctionBeginUser;
329: PetscCall(VecGetArrayRead(Y, &y));
330: PetscCall(VecGetArray(F, &f));
331: f[0] = (y[0] - t) / (y[0] + t);
332: PetscCall(VecRestoreArrayRead(Y, &y));
333: PetscCall(VecRestoreArray(F, &f));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
338: {
339: const PetscScalar *y;
340: PetscInt row = 0, col = 0;
341: PetscScalar value;
343: PetscFunctionBeginUser;
344: PetscCall(VecGetArrayRead(Y, &y));
345: value = 2 * t / ((t + y[0]) * (t + y[0]));
346: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
347: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
348: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
349: PetscCall(VecRestoreArrayRead(Y, &y));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
354: {
355: PetscScalar *f;
356: const PetscScalar *y;
358: PetscFunctionBeginUser;
359: PetscCall(VecGetArrayRead(Y, &y));
360: PetscCall(VecGetArray(F, &f));
361: f[0] = (y[0] - t) / (y[0] + t);
362: PetscCall(VecRestoreArrayRead(Y, &y));
363: PetscCall(VecRestoreArray(F, &f));
364: /* Left hand side = ydot - f(y) */
365: PetscCall(VecAYPX(F, -1.0, Ydot));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
370: {
371: const PetscScalar *y;
372: PetscInt row = 0, col = 0;
373: PetscScalar value;
375: PetscFunctionBeginUser;
376: PetscCall(VecGetArrayRead(Y, &y));
377: value = a - 2 * t / ((t + y[0]) * (t + y[0]));
378: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
379: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
380: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
381: PetscCall(VecRestoreArrayRead(Y, &y));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /* Hull, 1972, Problem B1 */
387: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
388: {
389: PetscScalar *f;
390: const PetscScalar *y;
392: PetscFunctionBeginUser;
393: PetscCall(VecGetArrayRead(Y, &y));
394: PetscCall(VecGetArray(F, &f));
395: f[0] = 2.0 * (y[0] - y[0] * y[1]);
396: f[1] = -(y[1] - y[0] * y[1]);
397: PetscCall(VecRestoreArrayRead(Y, &y));
398: PetscCall(VecRestoreArray(F, &f));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
403: {
404: PetscScalar *f;
405: const PetscScalar *y;
407: PetscFunctionBeginUser;
408: PetscCall(VecGetArrayRead(Y, &y));
409: PetscCall(VecGetArray(F, &f));
410: f[0] = 2.0 * (y[0] - y[0] * y[1]);
411: f[1] = -(y[1] - y[0] * y[1]);
412: PetscCall(VecRestoreArrayRead(Y, &y));
413: PetscCall(VecRestoreArray(F, &f));
414: /* Left hand side = ydot - f(y) */
415: PetscCall(VecAYPX(F, -1.0, Ydot));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
420: {
421: const PetscScalar *y;
422: PetscInt row[2] = {0, 1};
423: PetscScalar value[2][2];
425: PetscFunctionBeginUser;
426: PetscCall(VecGetArrayRead(Y, &y));
427: value[0][0] = a - 2.0 * (1.0 - y[1]);
428: value[0][1] = 2.0 * y[0];
429: value[1][0] = -y[1];
430: value[1][1] = a + 1.0 - y[0];
431: PetscCall(MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES));
432: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
433: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
434: PetscCall(VecRestoreArrayRead(Y, &y));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /* Hull, 1972, Problem B2 */
440: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
441: {
442: PetscScalar *f;
443: const PetscScalar *y;
445: PetscFunctionBeginUser;
446: PetscCall(VecGetArrayRead(Y, &y));
447: PetscCall(VecGetArray(F, &f));
448: f[0] = -y[0] + y[1];
449: f[1] = y[0] - 2.0 * y[1] + y[2];
450: f[2] = y[1] - y[2];
451: PetscCall(VecRestoreArrayRead(Y, &y));
452: PetscCall(VecRestoreArray(F, &f));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
457: {
458: PetscScalar *f;
459: const PetscScalar *y;
461: PetscFunctionBeginUser;
462: PetscCall(VecGetArrayRead(Y, &y));
463: PetscCall(VecGetArray(F, &f));
464: f[0] = -y[0] + y[1];
465: f[1] = y[0] - 2.0 * y[1] + y[2];
466: f[2] = y[1] - y[2];
467: PetscCall(VecRestoreArrayRead(Y, &y));
468: PetscCall(VecRestoreArray(F, &f));
469: /* Left hand side = ydot - f(y) */
470: PetscCall(VecAYPX(F, -1.0, Ydot));
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
475: {
476: const PetscScalar *y;
477: PetscInt row[3] = {0, 1, 2};
478: PetscScalar value[3][3];
480: PetscFunctionBeginUser;
481: PetscCall(VecGetArrayRead(Y, &y));
482: value[0][0] = a + 1.0;
483: value[0][1] = -1.0;
484: value[0][2] = 0;
485: value[1][0] = -1.0;
486: value[1][1] = a + 2.0;
487: value[1][2] = -1.0;
488: value[2][0] = 0;
489: value[2][1] = -1.0;
490: value[2][2] = a + 1.0;
491: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
492: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
493: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
494: PetscCall(VecRestoreArrayRead(Y, &y));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /* Hull, 1972, Problem B3 */
500: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
501: {
502: PetscScalar *f;
503: const PetscScalar *y;
505: PetscFunctionBeginUser;
506: PetscCall(VecGetArrayRead(Y, &y));
507: PetscCall(VecGetArray(F, &f));
508: f[0] = -y[0];
509: f[1] = y[0] - y[1] * y[1];
510: f[2] = y[1] * y[1];
511: PetscCall(VecRestoreArrayRead(Y, &y));
512: PetscCall(VecRestoreArray(F, &f));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
517: {
518: PetscScalar *f;
519: const PetscScalar *y;
521: PetscFunctionBeginUser;
522: PetscCall(VecGetArrayRead(Y, &y));
523: PetscCall(VecGetArray(F, &f));
524: f[0] = -y[0];
525: f[1] = y[0] - y[1] * y[1];
526: f[2] = y[1] * y[1];
527: PetscCall(VecRestoreArrayRead(Y, &y));
528: PetscCall(VecRestoreArray(F, &f));
529: /* Left hand side = ydot - f(y) */
530: PetscCall(VecAYPX(F, -1.0, Ydot));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
535: {
536: const PetscScalar *y;
537: PetscInt row[3] = {0, 1, 2};
538: PetscScalar value[3][3];
540: PetscFunctionBeginUser;
541: PetscCall(VecGetArrayRead(Y, &y));
542: value[0][0] = a + 1.0;
543: value[0][1] = 0;
544: value[0][2] = 0;
545: value[1][0] = -1.0;
546: value[1][1] = a + 2.0 * y[1];
547: value[1][2] = 0;
548: value[2][0] = 0;
549: value[2][1] = -2.0 * y[1];
550: value[2][2] = a;
551: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
552: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
553: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
554: PetscCall(VecRestoreArrayRead(Y, &y));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /* Hull, 1972, Problem B4 */
560: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
561: {
562: PetscScalar *f;
563: const PetscScalar *y;
565: PetscFunctionBeginUser;
566: PetscCall(VecGetArrayRead(Y, &y));
567: PetscCall(VecGetArray(F, &f));
568: f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
569: f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
570: f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
571: PetscCall(VecRestoreArrayRead(Y, &y));
572: PetscCall(VecRestoreArray(F, &f));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
577: {
578: PetscScalar *f;
579: const PetscScalar *y;
581: PetscFunctionBeginUser;
582: PetscCall(VecGetArrayRead(Y, &y));
583: PetscCall(VecGetArray(F, &f));
584: f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
585: f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
586: f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]);
587: PetscCall(VecRestoreArrayRead(Y, &y));
588: PetscCall(VecRestoreArray(F, &f));
589: /* Left hand side = ydot - f(y) */
590: PetscCall(VecAYPX(F, -1.0, Ydot));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
595: {
596: const PetscScalar *y;
597: PetscInt row[3] = {0, 1, 2};
598: PetscScalar value[3][3], fac, fac2;
600: PetscFunctionBeginUser;
601: PetscCall(VecGetArrayRead(Y, &y));
602: fac = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -1.5);
603: fac2 = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -0.5);
604: value[0][0] = a + (y[1] * y[1] * y[2]) * fac;
605: value[0][1] = 1.0 - (y[0] * y[1] * y[2]) * fac;
606: value[0][2] = y[0] * fac2;
607: value[1][0] = -1.0 - y[0] * y[1] * y[2] * fac;
608: value[1][1] = a + y[0] * y[0] * y[2] * fac;
609: value[1][2] = y[1] * fac2;
610: value[2][0] = -y[1] * y[1] * fac;
611: value[2][1] = y[0] * y[1] * fac;
612: value[2][2] = a;
613: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
614: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
615: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
616: PetscCall(VecRestoreArrayRead(Y, &y));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: /* Hull, 1972, Problem B5 */
622: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
623: {
624: PetscScalar *f;
625: const PetscScalar *y;
627: PetscFunctionBeginUser;
628: PetscCall(VecGetArrayRead(Y, &y));
629: PetscCall(VecGetArray(F, &f));
630: f[0] = y[1] * y[2];
631: f[1] = -y[0] * y[2];
632: f[2] = -0.51 * y[0] * y[1];
633: PetscCall(VecRestoreArrayRead(Y, &y));
634: PetscCall(VecRestoreArray(F, &f));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
639: {
640: PetscScalar *f;
641: const PetscScalar *y;
643: PetscFunctionBeginUser;
644: PetscCall(VecGetArrayRead(Y, &y));
645: PetscCall(VecGetArray(F, &f));
646: f[0] = y[1] * y[2];
647: f[1] = -y[0] * y[2];
648: f[2] = -0.51 * y[0] * y[1];
649: PetscCall(VecRestoreArrayRead(Y, &y));
650: PetscCall(VecRestoreArray(F, &f));
651: /* Left hand side = ydot - f(y) */
652: PetscCall(VecAYPX(F, -1.0, Ydot));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
657: {
658: const PetscScalar *y;
659: PetscInt row[3] = {0, 1, 2};
660: PetscScalar value[3][3];
662: PetscFunctionBeginUser;
663: PetscCall(VecGetArrayRead(Y, &y));
664: value[0][0] = a;
665: value[0][1] = -y[2];
666: value[0][2] = -y[1];
667: value[1][0] = y[2];
668: value[1][1] = a;
669: value[1][2] = y[0];
670: value[2][0] = 0.51 * y[1];
671: value[2][1] = 0.51 * y[0];
672: value[2][2] = a;
673: PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES));
674: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
675: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
676: PetscCall(VecRestoreArrayRead(Y, &y));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: /* Kulikov, 2013, Problem I */
682: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
683: {
684: PetscScalar *f;
685: const PetscScalar *y;
687: PetscFunctionBeginUser;
688: PetscCall(VecGetArrayRead(Y, &y));
689: PetscCall(VecGetArray(F, &f));
690: f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3];
691: f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.));
692: f[2] = 2. * t * y[3];
693: f[3] = -2. * t * PetscLogScalar(y[0]);
694: PetscCall(VecRestoreArrayRead(Y, &y));
695: PetscCall(VecRestoreArray(F, &f));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
700: {
701: const PetscScalar *y;
702: PetscInt row[4] = {0, 1, 2, 3};
703: PetscScalar value[4][4];
704: 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));
1257: PetscFunctionReturn(PETSC_SUCCESS);
1258: }
1260: int main(int argc, char **argv)
1261: {
1262: char ptype[256] = "hull1972a1"; /* Problem specification */
1263: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1264: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1265: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1266: PetscReal dt;
1267: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1268: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1269: PetscReal *error; /* Array to store the errors for convergence analysis */
1270: PetscMPIInt size; /* No of processors */
1271: PetscBool flag; /* Flag denoting availability of exact solution */
1272: PetscInt r, N;
1274: /* Initialize program */
1275: PetscFunctionBeginUser;
1276: PetscCall(GetSize(&ptype[0], &N));
1277: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1279: /* Check if running with only 1 proc */
1280: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1281: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
1283: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
1284: PetscCall(PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL));
1285: PetscCall(PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL));
1286: PetscCall(PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL));
1287: PetscCall(PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL));
1288: PetscCall(PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL));
1289: PetscOptionsEnd();
1291: PetscCall(PetscMalloc1(n_refine, &error));
1292: for (r = 0, dt = dt_initial; r < n_refine; r++) {
1293: error[r] = 0;
1294: if (r > 0) dt /= refine_fac;
1296: 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));
1297: PetscCall(SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag));
1298: if (flag) {
1299: /* If exact solution available for the specified ODE */
1300: if (r > 0) {
1301: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac));
1302: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate));
1303: } else {
1304: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E.\n", (double)error[r]));
1305: }
1306: }
1307: }
1308: PetscCall(PetscFree(error));
1309: PetscCall(PetscFinalize());
1310: return 0;
1311: }
1313: /*TEST
1315: test:
1316: suffix: 2
1317: args: -ts_type glee -final_time 5 -ts_adapt_type none
1318: timeoutfactor: 3
1319: requires: !single
1321: test:
1322: suffix: 3
1323: 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
1324: timeoutfactor: 3
1325: requires: !single
1327: test:
1328: suffix: 4
1329: # the test is so sensitive that I could not even replace VecMAXPY with a sequence of VECAXPY, so I just disable this GEMV optimization
1330: args: -vec_maxpy_use_gemv 0
1331: 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
1332: timeoutfactor: 3
1333: requires: !single !__float128
1335: TEST*/