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: {
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: return 0;
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;
78: VecGetArrayRead(Y, &y);
79: VecGetArray(F, &f);
80: f[0] = -y[0];
81: VecRestoreArrayRead(Y, &y);
82: VecRestoreArray(F, &f);
83: return 0;
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;
93: VecGetArrayRead(Y, &y);
94: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
95: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
97: VecRestoreArrayRead(Y, &y);
98: return 0;
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;
107: VecGetArrayRead(Y, &y);
108: VecGetArray(F, &f);
109: f[0] = -y[0];
110: VecRestoreArrayRead(Y, &y);
111: VecRestoreArray(F, &f);
112: /* Left hand side = ydot - f(y) */
113: VecAYPX(F, -1.0, Ydot);
114: return 0;
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;
124: VecGetArrayRead(Y, &y);
125: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
126: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
128: VecRestoreArrayRead(Y, &y);
129: return 0;
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;
140: VecGetArrayRead(Y, &y);
141: VecGetArray(F, &f);
142: f[0] = -0.5 * y[0] * y[0] * y[0];
143: VecRestoreArrayRead(Y, &y);
144: VecRestoreArray(F, &f);
145: return 0;
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;
155: VecGetArrayRead(Y, &y);
156: value = -0.5 * 3.0 * y[0] * y[0];
157: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
158: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
159: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
160: VecRestoreArrayRead(Y, &y);
161: return 0;
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;
170: VecGetArrayRead(Y, &y);
171: VecGetArray(F, &f);
172: f[0] = -0.5 * y[0] * y[0] * y[0];
173: VecRestoreArrayRead(Y, &y);
174: VecRestoreArray(F, &f);
175: /* Left hand side = ydot - f(y) */
176: VecAYPX(F, -1.0, Ydot);
177: return 0;
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;
187: VecGetArrayRead(Y, &y);
188: value = a + 0.5 * 3.0 * y[0] * y[0];
189: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
190: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
191: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
192: VecRestoreArrayRead(Y, &y);
193: return 0;
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;
204: VecGetArrayRead(Y, &y);
205: VecGetArray(F, &f);
206: f[0] = y[0] * PetscCosReal(t);
207: VecRestoreArrayRead(Y, &y);
208: VecRestoreArray(F, &f);
209: return 0;
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);
219: VecGetArrayRead(Y, &y);
220: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
221: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
222: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
223: VecRestoreArrayRead(Y, &y);
224: return 0;
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;
233: VecGetArrayRead(Y, &y);
234: VecGetArray(F, &f);
235: f[0] = y[0] * PetscCosReal(t);
236: VecRestoreArrayRead(Y, &y);
237: VecRestoreArray(F, &f);
238: /* Left hand side = ydot - f(y) */
239: VecAYPX(F, -1.0, Ydot);
240: return 0;
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);
250: VecGetArrayRead(Y, &y);
251: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
252: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
253: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
254: VecRestoreArrayRead(Y, &y);
255: return 0;
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;
266: VecGetArrayRead(Y, &y);
267: VecGetArray(F, &f);
268: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
269: VecRestoreArrayRead(Y, &y);
270: VecRestoreArray(F, &f);
271: return 0;
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;
281: VecGetArrayRead(Y, &y);
282: value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05;
283: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
284: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
285: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
286: VecRestoreArrayRead(Y, &y);
287: return 0;
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;
296: VecGetArrayRead(Y, &y);
297: VecGetArray(F, &f);
298: f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]);
299: VecRestoreArrayRead(Y, &y);
300: VecRestoreArray(F, &f);
301: /* Left hand side = ydot - f(y) */
302: VecAYPX(F, -1.0, Ydot);
303: return 0;
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;
313: VecGetArrayRead(Y, &y);
314: value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05;
315: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
316: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
317: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
318: VecRestoreArrayRead(Y, &y);
319: return 0;
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;
330: VecGetArrayRead(Y, &y);
331: VecGetArray(F, &f);
332: f[0] = (y[0] - t) / (y[0] + t);
333: VecRestoreArrayRead(Y, &y);
334: VecRestoreArray(F, &f);
335: return 0;
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;
345: VecGetArrayRead(Y, &y);
346: value = 2 * t / ((t + y[0]) * (t + y[0]));
347: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
348: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
349: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
350: VecRestoreArrayRead(Y, &y);
351: return 0;
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;
360: VecGetArrayRead(Y, &y);
361: VecGetArray(F, &f);
362: f[0] = (y[0] - t) / (y[0] + t);
363: VecRestoreArrayRead(Y, &y);
364: VecRestoreArray(F, &f);
365: /* Left hand side = ydot - f(y) */
366: VecAYPX(F, -1.0, Ydot);
367: return 0;
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;
377: VecGetArrayRead(Y, &y);
378: value = a - 2 * t / ((t + y[0]) * (t + y[0]));
379: MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES);
380: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
381: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
382: VecRestoreArrayRead(Y, &y);
383: return 0;
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;
394: VecGetArrayRead(Y, &y);
395: VecGetArray(F, &f);
396: f[0] = 2.0 * (y[0] - y[0] * y[1]);
397: f[1] = -(y[1] - y[0] * y[1]);
398: VecRestoreArrayRead(Y, &y);
399: VecRestoreArray(F, &f);
400: return 0;
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;
409: VecGetArrayRead(Y, &y);
410: VecGetArray(F, &f);
411: f[0] = 2.0 * (y[0] - y[0] * y[1]);
412: f[1] = -(y[1] - y[0] * y[1]);
413: VecRestoreArrayRead(Y, &y);
414: VecRestoreArray(F, &f);
415: /* Left hand side = ydot - f(y) */
416: VecAYPX(F, -1.0, Ydot);
417: return 0;
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];
427: 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: MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES);
433: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
434: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
435: VecRestoreArrayRead(Y, &y);
436: return 0;
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;
447: VecGetArrayRead(Y, &y);
448: 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: VecRestoreArrayRead(Y, &y);
453: VecRestoreArray(F, &f);
454: return 0;
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;
463: VecGetArrayRead(Y, &y);
464: 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: VecRestoreArrayRead(Y, &y);
469: VecRestoreArray(F, &f);
470: /* Left hand side = ydot - f(y) */
471: VecAYPX(F, -1.0, Ydot);
472: return 0;
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];
482: 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: MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
493: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
494: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
495: VecRestoreArrayRead(Y, &y);
496: return 0;
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;
507: VecGetArrayRead(Y, &y);
508: 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: VecRestoreArrayRead(Y, &y);
513: VecRestoreArray(F, &f);
514: return 0;
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;
523: VecGetArrayRead(Y, &y);
524: 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: VecRestoreArrayRead(Y, &y);
529: VecRestoreArray(F, &f);
530: /* Left hand side = ydot - f(y) */
531: VecAYPX(F, -1.0, Ydot);
532: return 0;
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];
542: 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: MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
553: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
554: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
555: VecRestoreArrayRead(Y, &y);
556: return 0;
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;
567: VecGetArrayRead(Y, &y);
568: 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: VecRestoreArrayRead(Y, &y);
573: VecRestoreArray(F, &f);
574: return 0;
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;
583: VecGetArrayRead(Y, &y);
584: 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: VecRestoreArrayRead(Y, &y);
589: VecRestoreArray(F, &f);
590: /* Left hand side = ydot - f(y) */
591: VecAYPX(F, -1.0, Ydot);
592: return 0;
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;
602: 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: MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
615: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
616: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
617: VecRestoreArrayRead(Y, &y);
618: return 0;
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;
629: VecGetArrayRead(Y, &y);
630: 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: VecRestoreArrayRead(Y, &y);
635: VecRestoreArray(F, &f);
636: return 0;
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;
645: VecGetArrayRead(Y, &y);
646: 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: VecRestoreArrayRead(Y, &y);
651: VecRestoreArray(F, &f);
652: /* Left hand side = ydot - f(y) */
653: VecAYPX(F, -1.0, Ydot);
654: return 0;
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];
664: 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: MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES);
675: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
676: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
677: VecRestoreArrayRead(Y, &y);
678: return 0;
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;
689: VecGetArrayRead(Y, &y);
690: 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: VecRestoreArrayRead(Y, &y);
696: VecRestoreArray(F, &f);
697: return 0;
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;
707: 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: MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES);
729: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
730: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
731: VecRestoreArrayRead(Y, &y);
732: return 0;
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;
741: VecGetArrayRead(Y, &y);
742: 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: VecRestoreArrayRead(Y, &y);
748: VecRestoreArray(F, &f);
749: /* Left hand side = ydot - f(y) */
750: VecAYPX(F, -1.0, Ydot);
751: return 0;
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;
762: 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: MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES);
784: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
785: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
786: VecRestoreArrayRead(Y, &y);
787: return 0;
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;
799: VecGetSize(Y, &N);
800: VecGetArrayRead(Y, &y);
801: 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: VecRestoreArrayRead(Y, &y);
806: VecRestoreArray(F, &f);
807: return 0;
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;
817: VecGetSize(Y, &N);
818: VecGetArrayRead(Y, &y);
819: 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: VecRestoreArrayRead(Y, &y);
824: VecRestoreArray(F, &f);
825: /* Left hand side = ydot - f(y) */
826: VecAYPX(F, -1.0, Ydot);
827: return 0;
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];
837: VecGetSize(Y, &N);
838: VecGetArrayRead(Y, &y);
839: i = 0;
840: value[0] = a + 1;
841: col[0] = 0;
842: value[1] = 0;
843: col[1] = 1;
844: 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: 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: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
858: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
859: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
860: VecRestoreArrayRead(Y, &y);
861: return 0;
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;
873: VecGetSize(Y, &N);
874: VecGetArrayRead(Y, &y);
875: 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: VecRestoreArrayRead(Y, &y);
880: VecRestoreArray(F, &f);
881: return 0;
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;
891: VecGetSize(Y, &N);
892: VecGetArrayRead(Y, &y);
893: 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: VecRestoreArrayRead(Y, &y);
898: VecRestoreArray(F, &f);
899: /* Left hand side = ydot - f(y) */
900: VecAYPX(F, -1.0, Ydot);
901: return 0;
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];
911: VecGetSize(Y, &N);
912: VecGetArrayRead(Y, &y);
913: i = 0;
914: value[0] = a + 1;
915: col[0] = 0;
916: value[1] = 0;
917: col[1] = 1;
918: 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: 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: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
932: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
933: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
934: VecRestoreArrayRead(Y, &y);
935: return 0;
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;
947: VecGetSize(Y, &N);
948: VecGetArrayRead(Y, &y);
949: 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: VecRestoreArrayRead(Y, &y);
954: VecRestoreArray(F, &f);
955: return 0;
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;
965: VecGetSize(Y, &N);
966: VecGetArrayRead(Y, &y);
967: 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: VecRestoreArrayRead(Y, &y);
972: VecRestoreArray(F, &f);
973: /* Left hand side = ydot - f(y) */
974: VecAYPX(F, -1.0, Ydot);
975: return 0;
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];
985: VecGetSize(Y, &N);
986: 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: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
1011: }
1012: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1013: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1014: VecRestoreArrayRead(Y, &y);
1015: return 0;
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;
1030: GetSize((const char *)s, &N);
1031: VecZeroEntries(Y);
1032: 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: PetscOptionsGetScalarArray(NULL, NULL, "-yinit", y, &N, &flg);
1124: VecRestoreArray(Y, &y);
1125: return 0;
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;
1135: if (!strcmp(p, "hull1972a1")) {
1136: VecGetArray(Y, &y);
1137: y[0] = PetscExpReal(-t);
1138: *flag = PETSC_TRUE;
1139: VecRestoreArray(Y, &y);
1140: } else if (!strcmp(p, "hull1972a2")) {
1141: VecGetArray(Y, &y);
1142: y[0] = 1.0 / PetscSqrtReal(t + 1);
1143: *flag = PETSC_TRUE;
1144: VecRestoreArray(Y, &y);
1145: } else if (!strcmp(p, "hull1972a3")) {
1146: VecGetArray(Y, &y);
1147: y[0] = PetscExpReal(PetscSinReal(t));
1148: *flag = PETSC_TRUE;
1149: VecRestoreArray(Y, &y);
1150: } else if (!strcmp(p, "hull1972a4")) {
1151: VecGetArray(Y, &y);
1152: y[0] = 20.0 / (1 + 19.0 * PetscExpReal(-t / 4.0));
1153: *flag = PETSC_TRUE;
1154: VecRestoreArray(Y, &y);
1155: } else if (!strcmp(p, "kulik2013i")) {
1156: 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: VecRestoreArray(Y, &y);
1163: } else {
1164: VecSet(Y, 0);
1165: *flag = PETSC_FALSE;
1166: }
1167: return 0;
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 */
1184: GetSize((const char *)&ptype[0], &N);
1186: VecCreate(PETSC_COMM_WORLD, &Y);
1187: VecSetSizes(Y, N, PETSC_DECIDE);
1188: VecSetUp(Y);
1189: VecSet(Y, 0);
1191: /* Initialize the problem */
1192: Initialize(Y, &ptype[0]);
1194: /* Create and initialize the time-integrator */
1195: TSCreate(PETSC_COMM_WORLD, &ts);
1196: /* Default time integration options */
1197: TSSetType(ts, TSRK);
1198: TSSetMaxSteps(ts, maxiter);
1199: TSSetMaxTime(ts, tfinal);
1200: TSSetTimeStep(ts, dt);
1201: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1202: /* Read command line options for time integration */
1203: TSSetFromOptions(ts);
1204: /* Set solution vector */
1205: TSSetSolution(ts, Y);
1206: /* Specify left/right-hand side functions */
1207: 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: TSSetRHSFunction(ts, NULL, RHSFunction, &ptype[0]);
1212: MatCreate(PETSC_COMM_WORLD, &Jac);
1213: MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N);
1214: MatSetFromOptions(Jac);
1215: MatSetUp(Jac);
1216: 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: TSSetIFunction(ts, NULL, IFunction, &ptype[0]);
1221: MatCreate(PETSC_COMM_WORLD, &Jac);
1222: MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N);
1223: MatSetFromOptions(Jac);
1224: MatSetUp(Jac);
1225: TSSetIJacobian(ts, Jac, Jac, IJacobian, &ptype[0]);
1226: }
1228: /* Solve */
1229: TSSolve(ts, Y);
1230: TSGetTime(ts, &final_time);
1232: /* Get the estimated error, if available */
1233: VecDuplicate(Y, &Yerr);
1234: VecZeroEntries(Yerr);
1235: TSGetTimeError(ts, 0, &Yerr);
1236: VecNorm(Yerr, NORM_2, &err_norm);
1237: VecDestroy(&Yerr);
1238: PetscPrintf(PETSC_COMM_WORLD, "Estimated Error = %e.\n", (double)err_norm);
1240: /* Exact solution */
1241: VecDuplicate(Y, &Yex);
1242: if (PetscAbsReal(final_time - tfinal) > 2. * PETSC_MACHINE_EPSILON) {
1243: 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: ExactSolution(Yex, &ptype[0], final_time, exact_flag);
1247: /* Calculate Error */
1248: VecAYPX(Yex, -1.0, Y);
1249: VecNorm(Yex, NORM_2, error);
1250: *error = PetscSqrtReal(((*error) * (*error)) / N);
1252: /* Clean up and finalize */
1253: MatDestroy(&Jac);
1254: TSDestroy(&ts);
1255: VecDestroy(&Yex);
1256: VecDestroy(&Y);
1258: return 0;
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 */
1277: GetSize(&ptype[0], &N);
1278: PetscInitialize(&argc, &argv, (char *)0, help);
1280: /* Check if running with only 1 proc */
1281: MPI_Comm_size(PETSC_COMM_WORLD, &size);
1284: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
1285: PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL);
1286: PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL);
1287: PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL);
1288: PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL);
1289: PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL);
1290: PetscOptionsEnd();
1292: 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: 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: 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: PetscPrintf(PETSC_COMM_WORLD, "Error = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate);
1304: } else {
1305: PetscPrintf(PETSC_COMM_WORLD, "Error = %E.\n", (double)error[r]);
1306: }
1307: }
1308: }
1309: PetscFree(error);
1310: 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*/