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: Concepts: TS
6: Useful command line parameters:
7: -problem <hull1972a1>: choose which problem to solve (see references
8: for complete listing of problems).
9: -ts_type <euler>: specify time-integrator
10: -ts_adapt_type <basic>: specify time-step adapting (none,basic,advanced)
11: -refinement_levels <1>: number of refinement levels for convergence analysis
12: -refinement_factor <2.0>: factor to refine time step size by for convergence analysis
13: -dt <0.01>: specify time step (initial time step for convergence analysis)
15: */
17: /*
18: List of cases and their names in the code:-
19: From Hull, T.E., Enright, W.H., Fellen, B.M., and Sedgwick, A.E.,
20: "Comparing Numerical Methods for Ordinary Differential
21: Equations", SIAM J. Numer. Anal., 9(4), 1972, pp. 603 - 635
22: A1 -> "hull1972a1" (exact solution available)
23: A2 -> "hull1972a2" (exact solution available)
24: A3 -> "hull1972a3" (exact solution available)
25: A4 -> "hull1972a4" (exact solution available)
26: A5 -> "hull1972a5"
27: B1 -> "hull1972b1"
28: B2 -> "hull1972b2"
29: B3 -> "hull1972b3"
30: B4 -> "hull1972b4"
31: B5 -> "hull1972b5"
32: C1 -> "hull1972c1"
33: C2 -> "hull1972c2"
34: C3 -> "hull1972c3"
35: C4 -> "hull1972c4"
37: From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
38: https://arxiv.org/abs/1503.05166, 2016
40: Kulikov2013I -> "kulik2013i"
42: */
44: #include <petscts.h>
46: /* Function declarations */
47: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
48: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
49: PetscErrorCode (*IFunction) (TS,PetscReal,Vec,Vec,Vec,void*);
50: PetscErrorCode (*IJacobian) (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
52: /* Returns the size of the system of equations depending on problem specification */
53: PetscInt GetSize(const char *p)
54: {
55: if ((!strcmp(p,"hull1972a1"))
56: ||(!strcmp(p,"hull1972a2"))
57: ||(!strcmp(p,"hull1972a3"))
58: ||(!strcmp(p,"hull1972a4"))
59: ||(!strcmp(p,"hull1972a5"))) return 1;
60: else if (!strcmp(p,"hull1972b1")) return 2;
61: else if ((!strcmp(p,"hull1972b2"))
62: ||(!strcmp(p,"hull1972b3"))
63: ||(!strcmp(p,"hull1972b4"))
64: ||(!strcmp(p,"hull1972b5"))) return 3;
65: else if ((!strcmp(p,"kulik2013i"))) return 4;
66: else if ((!strcmp(p,"hull1972c1"))
67: ||(!strcmp(p,"hull1972c2"))
68: ||(!strcmp(p,"hull1972c3"))) return 10;
69: else if (!strcmp(p,"hull1972c4")) return 51;
70: else PetscFunctionReturn(-1);
71: }
73: /****************************************************************/
75: /* Problem specific functions */
77: /* Hull, 1972, Problem A1 */
79: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
80: {
81: PetscScalar *f;
82: const PetscScalar *y;
84: VecGetArrayRead(Y,&y);
85: VecGetArray(F,&f);
86: f[0] = -y[0];
87: VecRestoreArrayRead(Y,&y);
88: VecRestoreArray(F,&f);
89: return 0;
90: }
92: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
93: {
94: const PetscScalar *y;
95: PetscInt row = 0,col = 0;
96: PetscScalar value = -1.0;
98: VecGetArrayRead(Y,&y);
99: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
102: VecRestoreArrayRead(Y,&y);
103: return 0;
104: }
106: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
107: {
108: const PetscScalar *y;
109: PetscScalar *f;
111: VecGetArrayRead(Y,&y);
112: VecGetArray(F,&f);
113: f[0] = -y[0];
114: VecRestoreArrayRead(Y,&y);
115: VecRestoreArray(F,&f);
116: /* Left hand side = ydot - f(y) */
117: VecAYPX(F,-1.0,Ydot);
118: return 0;
119: }
121: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
122: {
123: const PetscScalar *y;
124: PetscInt row = 0,col = 0;
125: PetscScalar value = a - 1.0;
127: VecGetArrayRead(Y,&y);
128: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
129: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
131: VecRestoreArrayRead(Y,&y);
132: return 0;
133: }
135: /* Hull, 1972, Problem A2 */
137: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
138: {
139: const PetscScalar *y;
140: PetscScalar *f;
142: VecGetArrayRead(Y,&y);
143: VecGetArray(F,&f);
144: f[0] = -0.5*y[0]*y[0]*y[0];
145: VecRestoreArrayRead(Y,&y);
146: VecRestoreArray(F,&f);
147: return 0;
148: }
150: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
151: {
152: const PetscScalar *y;
153: PetscInt row = 0,col = 0;
154: PetscScalar value;
156: VecGetArrayRead(Y,&y);
157: value = -0.5*3.0*y[0]*y[0];
158: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
159: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
160: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
161: VecRestoreArrayRead(Y,&y);
162: return 0;
163: }
165: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
166: {
167: PetscScalar *f;
168: 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;
186: VecGetArrayRead(Y,&y);
187: value = a + 0.5*3.0*y[0]*y[0];
188: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
189: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
190: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
191: VecRestoreArrayRead(Y,&y);
192: return 0;
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: VecGetArrayRead(Y,&y);
203: VecGetArray(F,&f);
204: f[0] = y[0]*PetscCosReal(t);
205: VecRestoreArrayRead(Y,&y);
206: VecRestoreArray(F,&f);
207: return 0;
208: }
210: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
211: {
212: const PetscScalar *y;
213: PetscInt row = 0,col = 0;
214: PetscScalar value = PetscCosReal(t);
216: VecGetArrayRead(Y,&y);
217: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
218: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
219: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
220: VecRestoreArrayRead(Y,&y);
221: return 0;
222: }
224: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
225: {
226: PetscScalar *f;
227: const PetscScalar *y;
229: VecGetArrayRead(Y,&y);
230: VecGetArray(F,&f);
231: f[0] = y[0]*PetscCosReal(t);
232: VecRestoreArrayRead(Y,&y);
233: VecRestoreArray(F,&f);
234: /* Left hand side = ydot - f(y) */
235: VecAYPX(F,-1.0,Ydot);
236: return 0;
237: }
239: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
240: {
241: const PetscScalar *y;
242: PetscInt row = 0,col = 0;
243: PetscScalar value = a - PetscCosReal(t);
245: VecGetArrayRead(Y,&y);
246: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
247: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
248: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
249: VecRestoreArrayRead(Y,&y);
250: return 0;
251: }
253: /* Hull, 1972, Problem A4 */
255: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
256: {
257: PetscScalar *f;
258: const PetscScalar *y;
260: VecGetArrayRead(Y,&y);
261: VecGetArray(F,&f);
262: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
263: VecRestoreArrayRead(Y,&y);
264: VecRestoreArray(F,&f);
265: return 0;
266: }
268: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
269: {
270: const PetscScalar *y;
271: PetscInt row = 0,col = 0;
272: PetscScalar value;
274: VecGetArrayRead(Y,&y);
275: value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
276: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
277: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
278: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
279: VecRestoreArrayRead(Y,&y);
280: return 0;
281: }
283: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
284: {
285: PetscScalar *f;
286: const PetscScalar *y;
288: VecGetArrayRead(Y,&y);
289: VecGetArray(F,&f);
290: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
291: VecRestoreArrayRead(Y,&y);
292: VecRestoreArray(F,&f);
293: /* Left hand side = ydot - f(y) */
294: VecAYPX(F,-1.0,Ydot);
295: return 0;
296: }
298: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
299: {
300: const PetscScalar *y;
301: PetscInt row = 0,col = 0;
302: PetscScalar value;
304: VecGetArrayRead(Y,&y);
305: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
306: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
307: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
309: VecRestoreArrayRead(Y,&y);
310: return 0;
311: }
313: /* Hull, 1972, Problem A5 */
315: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
316: {
317: PetscScalar *f;
318: const PetscScalar *y;
320: VecGetArrayRead(Y,&y);
321: VecGetArray(F,&f);
322: f[0] = (y[0]-t)/(y[0]+t);
323: VecRestoreArrayRead(Y,&y);
324: VecRestoreArray(F,&f);
325: return 0;
326: }
328: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
329: {
330: const PetscScalar *y;
331: PetscInt row = 0,col = 0;
332: PetscScalar value;
334: VecGetArrayRead(Y,&y);
335: value = 2*t/((t+y[0])*(t+y[0]));
336: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
337: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
338: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
339: VecRestoreArrayRead(Y,&y);
340: return 0;
341: }
343: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
344: {
345: PetscScalar *f;
346: const PetscScalar *y;
348: VecGetArrayRead(Y,&y);
349: VecGetArray(F,&f);
350: f[0] = (y[0]-t)/(y[0]+t);
351: VecRestoreArrayRead(Y,&y);
352: VecRestoreArray(F,&f);
353: /* Left hand side = ydot - f(y) */
354: VecAYPX(F,-1.0,Ydot);
355: return 0;
356: }
358: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
359: {
360: const PetscScalar *y;
361: PetscInt row = 0,col = 0;
362: PetscScalar value;
364: VecGetArrayRead(Y,&y);
365: value = a - 2*t/((t+y[0])*(t+y[0]));
366: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
367: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
368: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
369: VecRestoreArrayRead(Y,&y);
370: return 0;
371: }
373: /* Hull, 1972, Problem B1 */
375: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
376: {
377: PetscScalar *f;
378: const PetscScalar *y;
380: VecGetArrayRead(Y,&y);
381: VecGetArray(F,&f);
382: f[0] = 2.0*(y[0] - y[0]*y[1]);
383: f[1] = -(y[1]-y[0]*y[1]);
384: VecRestoreArrayRead(Y,&y);
385: VecRestoreArray(F,&f);
386: return 0;
387: }
389: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
390: {
391: PetscScalar *f;
392: 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: /* Left hand side = ydot - f(y) */
401: VecAYPX(F,-1.0,Ydot);
402: return 0;
403: }
405: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
406: {
407: const PetscScalar *y;
408: PetscInt row[2] = {0,1};
409: PetscScalar value[2][2];
411: VecGetArrayRead(Y,&y);
412: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
413: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
414: MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
415: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
416: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
417: VecRestoreArrayRead(Y,&y);
418: return 0;
419: }
421: /* Hull, 1972, Problem B2 */
423: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
424: {
425: PetscScalar *f;
426: const PetscScalar *y;
428: VecGetArrayRead(Y,&y);
429: VecGetArray(F,&f);
430: f[0] = -y[0] + y[1];
431: f[1] = y[0] - 2.0*y[1] + y[2];
432: f[2] = y[1] - y[2];
433: VecRestoreArrayRead(Y,&y);
434: VecRestoreArray(F,&f);
435: return 0;
436: }
438: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
439: {
440: PetscScalar *f;
441: const PetscScalar *y;
443: VecGetArrayRead(Y,&y);
444: VecGetArray(F,&f);
445: f[0] = -y[0] + y[1];
446: f[1] = y[0] - 2.0*y[1] + y[2];
447: f[2] = y[1] - y[2];
448: VecRestoreArrayRead(Y,&y);
449: VecRestoreArray(F,&f);
450: /* Left hand side = ydot - f(y) */
451: VecAYPX(F,-1.0,Ydot);
452: return 0;
453: }
455: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
456: {
457: const PetscScalar *y;
458: PetscInt row[3] = {0,1,2};
459: PetscScalar value[3][3];
461: VecGetArrayRead(Y,&y);
462: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
463: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
464: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
465: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
466: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
467: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
468: VecRestoreArrayRead(Y,&y);
469: return 0;
470: }
472: /* Hull, 1972, Problem B3 */
474: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
475: {
476: PetscScalar *f;
477: const PetscScalar *y;
479: VecGetArrayRead(Y,&y);
480: VecGetArray(F,&f);
481: f[0] = -y[0];
482: f[1] = y[0] - y[1]*y[1];
483: f[2] = y[1]*y[1];
484: VecRestoreArrayRead(Y,&y);
485: VecRestoreArray(F,&f);
486: return 0;
487: }
489: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
490: {
491: PetscScalar *f;
492: const PetscScalar *y;
494: VecGetArrayRead(Y,&y);
495: VecGetArray(F,&f);
496: f[0] = -y[0];
497: f[1] = y[0] - y[1]*y[1];
498: f[2] = y[1]*y[1];
499: VecRestoreArrayRead(Y,&y);
500: VecRestoreArray(F,&f);
501: /* Left hand side = ydot - f(y) */
502: VecAYPX(F,-1.0,Ydot);
503: return 0;
504: }
506: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
507: {
508: const PetscScalar *y;
509: PetscInt row[3] = {0,1,2};
510: PetscScalar value[3][3];
512: VecGetArrayRead(Y,&y);
513: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
514: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
515: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
516: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
517: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
518: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
519: VecRestoreArrayRead(Y,&y);
520: return 0;
521: }
523: /* Hull, 1972, Problem B4 */
525: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
526: {
527: PetscScalar *f;
528: const PetscScalar *y;
530: VecGetArrayRead(Y,&y);
531: VecGetArray(F,&f);
532: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
533: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
534: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
535: VecRestoreArrayRead(Y,&y);
536: VecRestoreArray(F,&f);
537: return 0;
538: }
540: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
541: {
542: PetscScalar *f;
543: const PetscScalar *y;
545: VecGetArrayRead(Y,&y);
546: VecGetArray(F,&f);
547: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
548: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
549: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
550: VecRestoreArrayRead(Y,&y);
551: VecRestoreArray(F,&f);
552: /* Left hand side = ydot - f(y) */
553: VecAYPX(F,-1.0,Ydot);
554: return 0;
555: }
557: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
558: {
559: const PetscScalar *y;
560: PetscInt row[3] = {0,1,2};
561: PetscScalar value[3][3],fac,fac2;
563: VecGetArrayRead(Y,&y);
564: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
565: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
566: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
567: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
568: value[0][2] = y[0]*fac2;
569: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
570: value[1][1] = a + y[0]*y[0]*y[2]*fac;
571: value[1][2] = y[1]*fac2;
572: value[2][0] = -y[1]*y[1]*fac;
573: value[2][1] = y[0]*y[1]*fac;
574: value[2][2] = a;
575: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
576: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
577: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
578: VecRestoreArrayRead(Y,&y);
579: return 0;
580: }
582: /* Hull, 1972, Problem B5 */
584: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
585: {
586: PetscScalar *f;
587: const PetscScalar *y;
589: VecGetArrayRead(Y,&y);
590: VecGetArray(F,&f);
591: f[0] = y[1]*y[2];
592: f[1] = -y[0]*y[2];
593: f[2] = -0.51*y[0]*y[1];
594: VecRestoreArrayRead(Y,&y);
595: VecRestoreArray(F,&f);
596: return 0;
597: }
599: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
600: {
601: PetscScalar *f;
602: const PetscScalar *y;
604: VecGetArrayRead(Y,&y);
605: VecGetArray(F,&f);
606: f[0] = y[1]*y[2];
607: f[1] = -y[0]*y[2];
608: f[2] = -0.51*y[0]*y[1];
609: VecRestoreArrayRead(Y,&y);
610: VecRestoreArray(F,&f);
611: /* Left hand side = ydot - f(y) */
612: VecAYPX(F,-1.0,Ydot);
613: return 0;
614: }
616: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
617: {
618: const PetscScalar *y;
619: PetscInt row[3] = {0,1,2};
620: PetscScalar value[3][3];
622: VecGetArrayRead(Y,&y);
623: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
624: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
625: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
626: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
627: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
628: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
629: VecRestoreArrayRead(Y,&y);
630: return 0;
631: }
633: /* Kulikov, 2013, Problem I */
635: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
636: {
637: PetscScalar *f;
638: const PetscScalar *y;
640: VecGetArrayRead(Y,&y);
641: VecGetArray(F,&f);
642: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
643: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
644: f[2] = 2.*t*y[3];
645: f[3] = -2.*t*PetscLogScalar(y[0]);
646: VecRestoreArrayRead(Y,&y);
647: VecRestoreArray(F,&f);
648: return 0;
649: }
651: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
652: {
653: const PetscScalar *y;
654: PetscInt row[4] = {0,1,2,3};
655: PetscScalar value[4][4];
656: PetscScalar m1,m2;
657: VecGetArrayRead(Y,&y);
658: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
659: m2=2.*t*PetscPowScalar(y[1],1./5.);
660: value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
661: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
662: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
663: value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
664: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t;
665: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.;
666: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
667: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
668: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
669: VecRestoreArrayRead(Y,&y);
670: return 0;
671: }
673: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
674: {
675: PetscScalar *f;
676: const PetscScalar *y;
678: VecGetArrayRead(Y,&y);
679: VecGetArray(F,&f);
680: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
681: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
682: f[2] = 2.*t*y[3];
683: f[3] = -2.*t*PetscLogScalar(y[0]);
684: VecRestoreArrayRead(Y,&y);
685: VecRestoreArray(F,&f);
686: /* Left hand side = ydot - f(y) */
687: VecAYPX(F,-1.0,Ydot);
688: return 0;
689: }
691: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
692: {
693: const PetscScalar *y;
694: PetscInt row[4] = {0,1,2,3};
695: PetscScalar value[4][4];
696: PetscScalar m1,m2;
698: VecGetArrayRead(Y,&y);
699: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
700: m2=2.*t*PetscPowScalar(y[1],1./5.);
701: value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
702: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
703: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
704: value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2;
705: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t;
706: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a;
707: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
708: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
709: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
710: VecRestoreArrayRead(Y,&y);
711: return 0;
712: }
714: /* Hull, 1972, Problem C1 */
716: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
717: {
718: PetscScalar *f;
719: const PetscScalar *y;
720: PetscInt N,i;
722: VecGetSize (Y,&N);
723: VecGetArrayRead(Y,&y);
724: VecGetArray(F,&f);
725: f[0] = -y[0];
726: for (i = 1; i < N-1; i++) {
727: f[i] = y[i-1] - y[i];
728: }
729: f[N-1] = y[N-2];
730: VecRestoreArrayRead(Y,&y);
731: VecRestoreArray(F,&f);
732: return 0;
733: }
735: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
736: {
737: PetscScalar *f;
738: const PetscScalar *y;
739: PetscInt N,i;
741: VecGetSize (Y,&N);
742: VecGetArrayRead(Y,&y);
743: VecGetArray(F,&f);
744: f[0] = -y[0];
745: for (i = 1; i < N-1; i++) {
746: f[i] = y[i-1] - y[i];
747: }
748: f[N-1] = y[N-2];
749: VecRestoreArrayRead(Y,&y);
750: VecRestoreArray(F,&f);
751: /* Left hand side = ydot - f(y) */
752: VecAYPX(F,-1.0,Ydot);
753: return 0;
754: }
756: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
757: {
758: const PetscScalar *y;
759: PetscInt N,i,col[2];
760: PetscScalar value[2];
762: VecGetSize (Y,&N);
763: VecGetArrayRead(Y,&y);
764: i = 0;
765: value[0] = a+1; col[0] = 0;
766: value[1] = 0; col[1] = 1;
767: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
768: for (i = 0; i < N; i++) {
769: value[0] = -1; col[0] = i-1;
770: value[1] = a+1; col[1] = i;
771: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
772: }
773: i = N-1;
774: value[0] = -1; col[0] = N-2;
775: value[1] = a; col[1] = N-1;
776: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
777: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
778: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
779: VecRestoreArrayRead(Y,&y);
780: return 0;
781: }
783: /* Hull, 1972, Problem C2 */
785: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
786: {
787: const PetscScalar *y;
788: PetscScalar *f;
789: PetscInt N,i;
791: VecGetSize (Y,&N);
792: VecGetArrayRead(Y,&y);
793: VecGetArray(F,&f);
794: f[0] = -y[0];
795: for (i = 1; i < N-1; i++) {
796: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
797: }
798: f[N-1] = (PetscReal)(N-1)*y[N-2];
799: VecRestoreArrayRead(Y,&y);
800: VecRestoreArray(F,&f);
801: return 0;
802: }
804: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
805: {
806: PetscScalar *f;
807: const PetscScalar *y;
808: PetscInt N,i;
810: VecGetSize (Y,&N);
811: VecGetArrayRead(Y,&y);
812: VecGetArray(F,&f);
813: f[0] = -y[0];
814: for (i = 1; i < N-1; i++) {
815: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
816: }
817: f[N-1] = (PetscReal)(N-1)*y[N-2];
818: VecRestoreArrayRead(Y,&y);
819: VecRestoreArray(F,&f);
820: /* Left hand side = ydot - f(y) */
821: VecAYPX(F,-1.0,Ydot);
822: return 0;
823: }
825: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
826: {
827: const PetscScalar *y;
828: PetscInt N,i,col[2];
829: PetscScalar value[2];
831: VecGetSize (Y,&N);
832: VecGetArrayRead(Y,&y);
833: i = 0;
834: value[0] = a+1; col[0] = 0;
835: value[1] = 0; col[1] = 1;
836: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
837: for (i = 0; i < N; i++) {
838: value[0] = -(PetscReal) i; col[0] = i-1;
839: value[1] = a+(PetscReal)(i+1); col[1] = i;
840: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
841: }
842: i = N-1;
843: value[0] = -(PetscReal) (N-1); col[0] = N-2;
844: value[1] = a; col[1] = N-1;
845: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
846: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
847: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
848: VecRestoreArrayRead(Y,&y);
849: return 0;
850: }
852: /* Hull, 1972, Problem C3 and C4 */
854: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
855: {
856: PetscScalar *f;
857: const PetscScalar *y;
858: PetscInt N,i;
860: VecGetSize (Y,&N);
861: VecGetArrayRead(Y,&y);
862: VecGetArray(F,&f);
863: f[0] = -2.0*y[0] + y[1];
864: for (i = 1; i < N-1; i++) {
865: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
866: }
867: f[N-1] = y[N-2] - 2.0*y[N-1];
868: VecRestoreArrayRead(Y,&y);
869: VecRestoreArray(F,&f);
870: return 0;
871: }
873: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
874: {
875: PetscScalar *f;
876: const PetscScalar *y;
877: PetscInt N,i;
879: VecGetSize (Y,&N);
880: VecGetArrayRead(Y,&y);
881: VecGetArray(F,&f);
882: f[0] = -2.0*y[0] + y[1];
883: for (i = 1; i < N-1; i++) {
884: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
885: }
886: f[N-1] = y[N-2] - 2.0*y[N-1];
887: VecRestoreArrayRead(Y,&y);
888: VecRestoreArray(F,&f);
889: /* Left hand side = ydot - f(y) */
890: VecAYPX(F,-1.0,Ydot);
891: return 0;
892: }
894: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
895: {
896: const PetscScalar *y;
897: PetscScalar value[3];
898: PetscInt N,i,col[3];
900: VecGetSize (Y,&N);
901: VecGetArrayRead(Y,&y);
902: for (i = 0; i < N; i++) {
903: if (i == 0) {
904: value[0] = a+2; col[0] = i;
905: value[1] = -1; col[1] = i+1;
906: value[2] = 0; col[2] = i+2;
907: } else if (i == N-1) {
908: value[0] = 0; col[0] = i-2;
909: value[1] = -1; col[1] = i-1;
910: value[2] = a+2; col[2] = i;
911: } else {
912: value[0] = -1; col[0] = i-1;
913: value[1] = a+2; col[1] = i;
914: value[2] = -1; col[2] = i+1;
915: }
916: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
917: }
918: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
919: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
920: VecRestoreArrayRead(Y,&y);
921: return 0;
922: }
924: /***************************************************************************/
926: /* Sets the initial solution for the IVP and sets up the function pointers*/
927: PetscErrorCode Initialize(Vec Y, void* s)
928: {
929: char *p = (char*) s;
930: PetscScalar *y;
931: PetscReal t0;
932: PetscInt N = GetSize((const char *)s);
933: PetscBool flg;
935: VecZeroEntries(Y);
936: VecGetArray(Y,&y);
937: if (!strcmp(p,"hull1972a1")) {
938: y[0] = 1.0;
939: RHSFunction = RHSFunction_Hull1972A1;
940: RHSJacobian = RHSJacobian_Hull1972A1;
941: IFunction = IFunction_Hull1972A1;
942: IJacobian = IJacobian_Hull1972A1;
943: } else if (!strcmp(p,"hull1972a2")) {
944: y[0] = 1.0;
945: RHSFunction = RHSFunction_Hull1972A2;
946: RHSJacobian = RHSJacobian_Hull1972A2;
947: IFunction = IFunction_Hull1972A2;
948: IJacobian = IJacobian_Hull1972A2;
949: } else if (!strcmp(p,"hull1972a3")) {
950: y[0] = 1.0;
951: RHSFunction = RHSFunction_Hull1972A3;
952: RHSJacobian = RHSJacobian_Hull1972A3;
953: IFunction = IFunction_Hull1972A3;
954: IJacobian = IJacobian_Hull1972A3;
955: } else if (!strcmp(p,"hull1972a4")) {
956: y[0] = 1.0;
957: RHSFunction = RHSFunction_Hull1972A4;
958: RHSJacobian = RHSJacobian_Hull1972A4;
959: IFunction = IFunction_Hull1972A4;
960: IJacobian = IJacobian_Hull1972A4;
961: } else if (!strcmp(p,"hull1972a5")) {
962: y[0] = 4.0;
963: RHSFunction = RHSFunction_Hull1972A5;
964: RHSJacobian = RHSJacobian_Hull1972A5;
965: IFunction = IFunction_Hull1972A5;
966: IJacobian = IJacobian_Hull1972A5;
967: } else if (!strcmp(p,"hull1972b1")) {
968: y[0] = 1.0;
969: y[1] = 3.0;
970: RHSFunction = RHSFunction_Hull1972B1;
971: IFunction = IFunction_Hull1972B1;
972: IJacobian = IJacobian_Hull1972B1;
973: } else if (!strcmp(p,"hull1972b2")) {
974: y[0] = 2.0;
975: y[1] = 0.0;
976: y[2] = 1.0;
977: RHSFunction = RHSFunction_Hull1972B2;
978: IFunction = IFunction_Hull1972B2;
979: IJacobian = IJacobian_Hull1972B2;
980: } else if (!strcmp(p,"hull1972b3")) {
981: y[0] = 1.0;
982: y[1] = 0.0;
983: y[2] = 0.0;
984: RHSFunction = RHSFunction_Hull1972B3;
985: IFunction = IFunction_Hull1972B3;
986: IJacobian = IJacobian_Hull1972B3;
987: } else if (!strcmp(p,"hull1972b4")) {
988: y[0] = 3.0;
989: y[1] = 0.0;
990: y[2] = 0.0;
991: RHSFunction = RHSFunction_Hull1972B4;
992: IFunction = IFunction_Hull1972B4;
993: IJacobian = IJacobian_Hull1972B4;
994: } else if (!strcmp(p,"hull1972b5")) {
995: y[0] = 0.0;
996: y[1] = 1.0;
997: y[2] = 1.0;
998: RHSFunction = RHSFunction_Hull1972B5;
999: IFunction = IFunction_Hull1972B5;
1000: IJacobian = IJacobian_Hull1972B5;
1001: } else if (!strcmp(p,"kulik2013i")) {
1002: t0=0.;
1003: y[0] = PetscExpReal(PetscSinReal(t0*t0));
1004: y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1005: y[2] = PetscSinReal(t0*t0)+1.0;
1006: y[3] = PetscCosReal(t0*t0);
1007: RHSFunction = RHSFunction_Kulikov2013I;
1008: RHSJacobian = RHSJacobian_Kulikov2013I;
1009: IFunction = IFunction_Kulikov2013I;
1010: IJacobian = IJacobian_Kulikov2013I;
1011: } else if (!strcmp(p,"hull1972c1")) {
1012: y[0] = 1.0;
1013: RHSFunction = RHSFunction_Hull1972C1;
1014: IFunction = IFunction_Hull1972C1;
1015: IJacobian = IJacobian_Hull1972C1;
1016: } else if (!strcmp(p,"hull1972c2")) {
1017: y[0] = 1.0;
1018: RHSFunction = RHSFunction_Hull1972C2;
1019: IFunction = IFunction_Hull1972C2;
1020: IJacobian = IJacobian_Hull1972C2;
1021: } else if ((!strcmp(p,"hull1972c3"))
1022: ||(!strcmp(p,"hull1972c4"))) {
1023: y[0] = 1.0;
1024: RHSFunction = RHSFunction_Hull1972C34;
1025: IFunction = IFunction_Hull1972C34;
1026: IJacobian = IJacobian_Hull1972C34;
1027: }
1028: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1030: VecRestoreArray(Y,&y);
1031: return 0;
1032: }
1034: /* Calculates the exact solution to problems that have one */
1035: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1036: {
1037: char *p = (char*) s;
1038: PetscScalar *y;
1040: if (!strcmp(p,"hull1972a1")) {
1041: VecGetArray(Y,&y);
1042: y[0] = PetscExpReal(-t);
1043: *flag = PETSC_TRUE;
1044: VecRestoreArray(Y,&y);
1045: } else if (!strcmp(p,"hull1972a2")) {
1046: VecGetArray(Y,&y);
1047: y[0] = 1.0/PetscSqrtReal(t+1);
1048: *flag = PETSC_TRUE;
1049: VecRestoreArray(Y,&y);
1050: } else if (!strcmp(p,"hull1972a3")) {
1051: VecGetArray(Y,&y);
1052: y[0] = PetscExpReal(PetscSinReal(t));
1053: *flag = PETSC_TRUE;
1054: VecRestoreArray(Y,&y);
1055: } else if (!strcmp(p,"hull1972a4")) {
1056: VecGetArray(Y,&y);
1057: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1058: *flag = PETSC_TRUE;
1059: VecRestoreArray(Y,&y);
1060: } else if (!strcmp(p,"kulik2013i")) {
1061: VecGetArray(Y,&y);
1062: y[0] = PetscExpReal(PetscSinReal(t*t));
1063: y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1064: y[2] = PetscSinReal(t*t)+1.0;
1065: y[3] = PetscCosReal(t*t);
1066: *flag = PETSC_TRUE;
1067: VecRestoreArray(Y,&y);
1068: } else {
1069: VecSet(Y,0);
1070: *flag = PETSC_FALSE;
1071: }
1072: return 0;
1073: }
1075: /* Solves the specified ODE and computes the error if exact solution is available */
1076: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1077: {
1078: TS ts; /* time-integrator */
1079: Vec Y; /* Solution vector */
1080: Vec Yex; /* Exact solution */
1081: PetscInt N; /* Size of the system of equations */
1082: TSType time_scheme; /* Type of time-integration scheme */
1083: Mat Jac = NULL; /* Jacobian matrix */
1084: Vec Yerr; /* Auxiliary solution vector */
1085: PetscReal err_norm; /* Estimated error norm */
1086: PetscReal final_time; /* Actual final time from the integrator */
1088: N = GetSize((const char *)&ptype[0]);
1090: VecCreate(PETSC_COMM_WORLD,&Y);
1091: VecSetSizes(Y,N,PETSC_DECIDE);
1092: VecSetUp(Y);
1093: VecSet(Y,0);
1095: /* Initialize the problem */
1096: Initialize(Y,&ptype[0]);
1098: /* Create and initialize the time-integrator */
1099: TSCreate(PETSC_COMM_WORLD,&ts);
1100: /* Default time integration options */
1101: TSSetType(ts,TSRK);
1102: TSSetMaxSteps(ts,maxiter);
1103: TSSetMaxTime(ts,tfinal);
1104: TSSetTimeStep(ts,dt);
1105: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1106: /* Read command line options for time integration */
1107: TSSetFromOptions(ts);
1108: /* Set solution vector */
1109: TSSetSolution(ts,Y);
1110: /* Specify left/right-hand side functions */
1111: TSGetType(ts,&time_scheme);
1113: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1114: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1115: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1116: MatCreate(PETSC_COMM_WORLD,&Jac);
1117: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1118: MatSetFromOptions(Jac);
1119: MatSetUp(Jac);
1120: TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1121: } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1122: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1123: /* and its Jacobian function */
1124: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1125: MatCreate(PETSC_COMM_WORLD,&Jac);
1126: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1127: MatSetFromOptions(Jac);
1128: MatSetUp(Jac);
1129: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1130: }
1132: /* Solve */
1133: TSSolve(ts,Y);
1134: TSGetTime(ts,&final_time);
1136: /* Get the estimated error, if available */
1137: VecDuplicate(Y,&Yerr);
1138: VecZeroEntries(Yerr);
1139: TSGetTimeError(ts,0,&Yerr);
1140: VecNorm(Yerr,NORM_2,&err_norm);
1141: VecDestroy(&Yerr);
1142: PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);
1144: /* Exact solution */
1145: VecDuplicate(Y,&Yex);
1146: if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1147: 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);
1148: }
1149: ExactSolution(Yex,&ptype[0],final_time,exact_flag);
1151: /* Calculate Error */
1152: VecAYPX(Yex,-1.0,Y);
1153: VecNorm(Yex,NORM_2,error);
1154: *error = PetscSqrtReal(((*error)*(*error))/N);
1156: /* Clean up and finalize */
1157: MatDestroy(&Jac);
1158: TSDestroy(&ts);
1159: VecDestroy(&Yex);
1160: VecDestroy(&Y);
1162: return 0;
1163: }
1165: int main(int argc, char **argv)
1166: {
1167: PetscErrorCode ierr; /* Error code */
1168: char ptype[256] = "hull1972a1"; /* Problem specification */
1169: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1170: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1171: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1172: PetscReal dt;
1173: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1174: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1175: PetscReal *error; /* Array to store the errors for convergence analysis */
1176: PetscMPIInt size; /* No of processors */
1177: PetscBool flag; /* Flag denoting availability of exact solution */
1178: PetscInt r;
1180: /* Initialize program */
1181: PetscInitialize(&argc,&argv,(char*)0,help);
1183: /* Check if running with only 1 proc */
1184: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1187: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1188: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1189: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1190: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1191: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1192: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1193: PetscOptionsEnd();
1195: PetscMalloc1(n_refine,&error);
1196: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1197: error[r] = 0;
1198: if (r > 0) dt /= refine_fac;
1200: PetscPrintf(PETSC_COMM_WORLD,"Solving ODE \"%s\" with dt %f, final time %f and system size %D.\n",ptype,(double)dt,(double)tfinal,GetSize(&ptype[0]));
1201: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1202: if (flag) {
1203: /* If exact solution available for the specified ODE */
1204: if (r > 0) {
1205: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1206: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1207: } else {
1208: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1209: }
1210: }
1211: }
1212: PetscFree(error);
1213: PetscFinalize();
1214: return 0;
1215: }
1217: /*TEST
1219: test:
1220: suffix: 2
1221: args: -ts_type glee -final_time 5 -ts_adapt_type none
1222: timeoutfactor: 3
1223: requires: !single
1225: test:
1226: suffix: 3
1227: 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
1228: timeoutfactor: 3
1229: requires: !single
1231: test:
1232: suffix: 4
1233: 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
1234: timeoutfactor: 3
1235: requires: !single !__float128
1237: TEST*/