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: PetscInt GetSize(const char *p)
53: {
54: if ((!strcmp(p,"hull1972a1"))
55: ||(!strcmp(p,"hull1972a2"))
56: ||(!strcmp(p,"hull1972a3"))
57: ||(!strcmp(p,"hull1972a4"))
58: ||(!strcmp(p,"hull1972a5"))) return 1;
59: else if (!strcmp(p,"hull1972b1")) return 2;
60: else if ((!strcmp(p,"hull1972b2"))
61: ||(!strcmp(p,"hull1972b3"))
62: ||(!strcmp(p,"hull1972b4"))
63: ||(!strcmp(p,"hull1972b5"))) return 3;
64: else if ((!strcmp(p,"kulik2013i"))) return 4;
65: else if ((!strcmp(p,"hull1972c1"))
66: ||(!strcmp(p,"hull1972c2"))
67: ||(!strcmp(p,"hull1972c3"))) return 10;
68: else if (!strcmp(p,"hull1972c4")) return 51;
69: else PetscFunctionReturn(-1);
70: }
72: /****************************************************************/
74: /* Problem specific functions */
76: /* Hull, 1972, Problem A1 */
78: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
79: {
80: PetscScalar *f;
81: const PetscScalar *y;
83: VecGetArrayRead(Y,&y);
84: VecGetArray(F,&f);
85: f[0] = -y[0];
86: VecRestoreArrayRead(Y,&y);
87: VecRestoreArray(F,&f);
88: return 0;
89: }
91: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
92: {
93: const PetscScalar *y;
94: PetscInt row = 0,col = 0;
95: PetscScalar value = -1.0;
97: VecGetArrayRead(Y,&y);
98: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
99: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
101: VecRestoreArrayRead(Y,&y);
102: return 0;
103: }
105: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
106: {
107: const PetscScalar *y;
108: PetscScalar *f;
110: VecGetArrayRead(Y,&y);
111: VecGetArray(F,&f);
112: f[0] = -y[0];
113: VecRestoreArrayRead(Y,&y);
114: VecRestoreArray(F,&f);
115: /* Left hand side = ydot - f(y) */
116: VecAYPX(F,-1.0,Ydot);
117: return 0;
118: }
120: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
121: {
122: const PetscScalar *y;
123: PetscInt row = 0,col = 0;
124: PetscScalar value = a - 1.0;
126: VecGetArrayRead(Y,&y);
127: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
128: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
130: VecRestoreArrayRead(Y,&y);
131: return 0;
132: }
134: /* Hull, 1972, Problem A2 */
136: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
137: {
138: const PetscScalar *y;
139: PetscScalar *f;
141: VecGetArrayRead(Y,&y);
142: VecGetArray(F,&f);
143: f[0] = -0.5*y[0]*y[0]*y[0];
144: VecRestoreArrayRead(Y,&y);
145: VecRestoreArray(F,&f);
146: return 0;
147: }
149: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
150: {
151: const PetscScalar *y;
152: PetscInt row = 0,col = 0;
153: 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;
169: VecGetArrayRead(Y,&y);
170: VecGetArray(F,&f);
171: f[0] = -0.5*y[0]*y[0]*y[0];
172: VecRestoreArrayRead(Y,&y);
173: VecRestoreArray(F,&f);
174: /* Left hand side = ydot - f(y) */
175: VecAYPX(F,-1.0,Ydot);
176: return 0;
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: VecGetArrayRead(Y,&y);
186: value = a + 0.5*3.0*y[0]*y[0];
187: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
188: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
189: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
190: VecRestoreArrayRead(Y,&y);
191: return 0;
192: }
194: /* Hull, 1972, Problem A3 */
196: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
197: {
198: const PetscScalar *y;
199: PetscScalar *f;
201: VecGetArrayRead(Y,&y);
202: VecGetArray(F,&f);
203: f[0] = y[0]*PetscCosReal(t);
204: VecRestoreArrayRead(Y,&y);
205: VecRestoreArray(F,&f);
206: return 0;
207: }
209: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
210: {
211: const PetscScalar *y;
212: PetscInt row = 0,col = 0;
213: PetscScalar value = PetscCosReal(t);
215: VecGetArrayRead(Y,&y);
216: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
217: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
218: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
219: VecRestoreArrayRead(Y,&y);
220: return 0;
221: }
223: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
224: {
225: PetscScalar *f;
226: const PetscScalar *y;
228: VecGetArrayRead(Y,&y);
229: VecGetArray(F,&f);
230: f[0] = y[0]*PetscCosReal(t);
231: VecRestoreArrayRead(Y,&y);
232: VecRestoreArray(F,&f);
233: /* Left hand side = ydot - f(y) */
234: VecAYPX(F,-1.0,Ydot);
235: return 0;
236: }
238: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
239: {
240: const PetscScalar *y;
241: PetscInt row = 0,col = 0;
242: PetscScalar value = a - PetscCosReal(t);
244: VecGetArrayRead(Y,&y);
245: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
246: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
247: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
248: VecRestoreArrayRead(Y,&y);
249: return 0;
250: }
252: /* Hull, 1972, Problem A4 */
254: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
255: {
256: PetscScalar *f;
257: const PetscScalar *y;
259: VecGetArrayRead(Y,&y);
260: VecGetArray(F,&f);
261: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
262: VecRestoreArrayRead(Y,&y);
263: VecRestoreArray(F,&f);
264: return 0;
265: }
267: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
268: {
269: const PetscScalar *y;
270: PetscInt row = 0,col = 0;
271: PetscScalar value;
273: VecGetArrayRead(Y,&y);
274: value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
275: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
276: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
278: VecRestoreArrayRead(Y,&y);
279: return 0;
280: }
282: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
283: {
284: PetscScalar *f;
285: const PetscScalar *y;
287: VecGetArrayRead(Y,&y);
288: VecGetArray(F,&f);
289: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
290: VecRestoreArrayRead(Y,&y);
291: VecRestoreArray(F,&f);
292: /* Left hand side = ydot - f(y) */
293: VecAYPX(F,-1.0,Ydot);
294: return 0;
295: }
297: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
298: {
299: const PetscScalar *y;
300: PetscInt row = 0,col = 0;
301: PetscScalar value;
303: VecGetArrayRead(Y,&y);
304: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
305: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
306: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
307: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
308: VecRestoreArrayRead(Y,&y);
309: return 0;
310: }
312: /* Hull, 1972, Problem A5 */
314: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
315: {
316: PetscScalar *f;
317: const PetscScalar *y;
319: VecGetArrayRead(Y,&y);
320: VecGetArray(F,&f);
321: f[0] = (y[0]-t)/(y[0]+t);
322: VecRestoreArrayRead(Y,&y);
323: VecRestoreArray(F,&f);
324: return 0;
325: }
327: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
328: {
329: const PetscScalar *y;
330: PetscInt row = 0,col = 0;
331: PetscScalar value;
333: VecGetArrayRead(Y,&y);
334: value = 2*t/((t+y[0])*(t+y[0]));
335: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
336: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
337: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
338: VecRestoreArrayRead(Y,&y);
339: return 0;
340: }
342: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
343: {
344: PetscScalar *f;
345: const PetscScalar *y;
347: VecGetArrayRead(Y,&y);
348: VecGetArray(F,&f);
349: f[0] = (y[0]-t)/(y[0]+t);
350: VecRestoreArrayRead(Y,&y);
351: VecRestoreArray(F,&f);
352: /* Left hand side = ydot - f(y) */
353: VecAYPX(F,-1.0,Ydot);
354: return 0;
355: }
357: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
358: {
359: const PetscScalar *y;
360: PetscInt row = 0,col = 0;
361: PetscScalar value;
363: VecGetArrayRead(Y,&y);
364: value = a - 2*t/((t+y[0])*(t+y[0]));
365: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
366: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
367: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
368: VecRestoreArrayRead(Y,&y);
369: return 0;
370: }
372: /* Hull, 1972, Problem B1 */
374: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
375: {
376: PetscScalar *f;
377: const PetscScalar *y;
379: VecGetArrayRead(Y,&y);
380: VecGetArray(F,&f);
381: f[0] = 2.0*(y[0] - y[0]*y[1]);
382: f[1] = -(y[1]-y[0]*y[1]);
383: VecRestoreArrayRead(Y,&y);
384: VecRestoreArray(F,&f);
385: return 0;
386: }
388: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
389: {
390: PetscScalar *f;
391: const PetscScalar *y;
393: VecGetArrayRead(Y,&y);
394: VecGetArray(F,&f);
395: f[0] = 2.0*(y[0] - y[0]*y[1]);
396: f[1] = -(y[1]-y[0]*y[1]);
397: VecRestoreArrayRead(Y,&y);
398: VecRestoreArray(F,&f);
399: /* Left hand side = ydot - f(y) */
400: VecAYPX(F,-1.0,Ydot);
401: return 0;
402: }
404: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
405: {
406: const PetscScalar *y;
407: PetscInt row[2] = {0,1};
408: PetscScalar value[2][2];
410: VecGetArrayRead(Y,&y);
411: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
412: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
413: MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
414: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
415: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
416: VecRestoreArrayRead(Y,&y);
417: return 0;
418: }
420: /* Hull, 1972, Problem B2 */
422: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
423: {
424: PetscScalar *f;
425: const PetscScalar *y;
427: VecGetArrayRead(Y,&y);
428: VecGetArray(F,&f);
429: f[0] = -y[0] + y[1];
430: f[1] = y[0] - 2.0*y[1] + y[2];
431: f[2] = y[1] - y[2];
432: VecRestoreArrayRead(Y,&y);
433: VecRestoreArray(F,&f);
434: return 0;
435: }
437: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
438: {
439: PetscScalar *f;
440: const PetscScalar *y;
442: VecGetArrayRead(Y,&y);
443: VecGetArray(F,&f);
444: f[0] = -y[0] + y[1];
445: f[1] = y[0] - 2.0*y[1] + y[2];
446: f[2] = y[1] - y[2];
447: VecRestoreArrayRead(Y,&y);
448: VecRestoreArray(F,&f);
449: /* Left hand side = ydot - f(y) */
450: VecAYPX(F,-1.0,Ydot);
451: return 0;
452: }
454: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
455: {
456: const PetscScalar *y;
457: PetscInt row[3] = {0,1,2};
458: PetscScalar value[3][3];
460: VecGetArrayRead(Y,&y);
461: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
462: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
463: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
464: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
465: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
466: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
467: VecRestoreArrayRead(Y,&y);
468: return 0;
469: }
471: /* Hull, 1972, Problem B3 */
473: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
474: {
475: PetscScalar *f;
476: const PetscScalar *y;
478: VecGetArrayRead(Y,&y);
479: VecGetArray(F,&f);
480: f[0] = -y[0];
481: f[1] = y[0] - y[1]*y[1];
482: f[2] = y[1]*y[1];
483: VecRestoreArrayRead(Y,&y);
484: VecRestoreArray(F,&f);
485: return 0;
486: }
488: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
489: {
490: PetscScalar *f;
491: const PetscScalar *y;
493: VecGetArrayRead(Y,&y);
494: VecGetArray(F,&f);
495: f[0] = -y[0];
496: f[1] = y[0] - y[1]*y[1];
497: f[2] = y[1]*y[1];
498: VecRestoreArrayRead(Y,&y);
499: VecRestoreArray(F,&f);
500: /* Left hand side = ydot - f(y) */
501: VecAYPX(F,-1.0,Ydot);
502: return 0;
503: }
505: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
506: {
507: const PetscScalar *y;
508: PetscInt row[3] = {0,1,2};
509: PetscScalar value[3][3];
511: VecGetArrayRead(Y,&y);
512: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
513: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
514: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
515: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
516: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
517: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
518: VecRestoreArrayRead(Y,&y);
519: return 0;
520: }
522: /* Hull, 1972, Problem B4 */
524: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
525: {
526: PetscScalar *f;
527: const PetscScalar *y;
529: VecGetArrayRead(Y,&y);
530: VecGetArray(F,&f);
531: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
532: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
533: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
534: VecRestoreArrayRead(Y,&y);
535: VecRestoreArray(F,&f);
536: return 0;
537: }
539: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
540: {
541: PetscScalar *f;
542: const PetscScalar *y;
544: VecGetArrayRead(Y,&y);
545: VecGetArray(F,&f);
546: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
547: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
548: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
549: VecRestoreArrayRead(Y,&y);
550: VecRestoreArray(F,&f);
551: /* Left hand side = ydot - f(y) */
552: VecAYPX(F,-1.0,Ydot);
553: return 0;
554: }
556: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
557: {
558: const PetscScalar *y;
559: PetscInt row[3] = {0,1,2};
560: PetscScalar value[3][3],fac,fac2;
562: VecGetArrayRead(Y,&y);
563: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
564: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
565: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
566: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
567: value[0][2] = y[0]*fac2;
568: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
569: value[1][1] = a + y[0]*y[0]*y[2]*fac;
570: value[1][2] = y[1]*fac2;
571: value[2][0] = -y[1]*y[1]*fac;
572: value[2][1] = y[0]*y[1]*fac;
573: value[2][2] = a;
574: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
575: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
576: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
577: VecRestoreArrayRead(Y,&y);
578: return 0;
579: }
581: /* Hull, 1972, Problem B5 */
583: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
584: {
585: PetscScalar *f;
586: const PetscScalar *y;
588: VecGetArrayRead(Y,&y);
589: VecGetArray(F,&f);
590: f[0] = y[1]*y[2];
591: f[1] = -y[0]*y[2];
592: f[2] = -0.51*y[0]*y[1];
593: VecRestoreArrayRead(Y,&y);
594: VecRestoreArray(F,&f);
595: return 0;
596: }
598: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
599: {
600: PetscScalar *f;
601: const PetscScalar *y;
603: VecGetArrayRead(Y,&y);
604: VecGetArray(F,&f);
605: f[0] = y[1]*y[2];
606: f[1] = -y[0]*y[2];
607: f[2] = -0.51*y[0]*y[1];
608: VecRestoreArrayRead(Y,&y);
609: VecRestoreArray(F,&f);
610: /* Left hand side = ydot - f(y) */
611: VecAYPX(F,-1.0,Ydot);
612: return 0;
613: }
615: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
616: {
617: const PetscScalar *y;
618: PetscInt row[3] = {0,1,2};
619: PetscScalar value[3][3];
621: VecGetArrayRead(Y,&y);
622: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
623: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
624: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
625: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
626: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
627: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
628: VecRestoreArrayRead(Y,&y);
629: return 0;
630: }
632: /* Kulikov, 2013, Problem I */
634: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
635: {
636: PetscScalar *f;
637: const PetscScalar *y;
639: VecGetArrayRead(Y,&y);
640: VecGetArray(F,&f);
641: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
642: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
643: f[2] = 2.*t*y[3];
644: f[3] = -2.*t*PetscLogScalar(y[0]);
645: VecRestoreArrayRead(Y,&y);
646: VecRestoreArray(F,&f);
647: return 0;
648: }
650: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
651: {
652: const PetscScalar *y;
653: PetscInt row[4] = {0,1,2,3};
654: PetscScalar value[4][4];
655: PetscScalar m1,m2;
656: VecGetArrayRead(Y,&y);
657: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
658: m2=2.*t*PetscPowScalar(y[1],1./5.);
659: value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
660: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
661: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
662: value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
663: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t;
664: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.;
665: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
666: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
667: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
668: VecRestoreArrayRead(Y,&y);
669: return 0;
670: }
672: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
673: {
674: PetscScalar *f;
675: const PetscScalar *y;
677: VecGetArrayRead(Y,&y);
678: VecGetArray(F,&f);
679: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
680: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
681: f[2] = 2.*t*y[3];
682: f[3] = -2.*t*PetscLogScalar(y[0]);
683: VecRestoreArrayRead(Y,&y);
684: VecRestoreArray(F,&f);
685: /* Left hand side = ydot - f(y) */
686: VecAYPX(F,-1.0,Ydot);
687: return 0;
688: }
690: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
691: {
692: const PetscScalar *y;
693: PetscInt row[4] = {0,1,2,3};
694: PetscScalar value[4][4];
695: PetscScalar m1,m2;
697: VecGetArrayRead(Y,&y);
698: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
699: m2=2.*t*PetscPowScalar(y[1],1./5.);
700: value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
701: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
702: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
703: value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2;
704: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t;
705: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a;
706: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
707: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
708: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
709: VecRestoreArrayRead(Y,&y);
710: return 0;
711: }
713: /* Hull, 1972, Problem C1 */
715: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
716: {
717: PetscScalar *f;
718: const PetscScalar *y;
719: PetscInt N,i;
721: VecGetSize (Y,&N);
722: VecGetArrayRead(Y,&y);
723: VecGetArray(F,&f);
724: f[0] = -y[0];
725: for (i = 1; i < N-1; i++) {
726: f[i] = y[i-1] - y[i];
727: }
728: f[N-1] = y[N-2];
729: VecRestoreArrayRead(Y,&y);
730: VecRestoreArray(F,&f);
731: return 0;
732: }
734: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
735: {
736: PetscScalar *f;
737: const PetscScalar *y;
738: PetscInt N,i;
740: VecGetSize (Y,&N);
741: VecGetArrayRead(Y,&y);
742: VecGetArray(F,&f);
743: f[0] = -y[0];
744: for (i = 1; i < N-1; i++) {
745: f[i] = y[i-1] - y[i];
746: }
747: f[N-1] = y[N-2];
748: VecRestoreArrayRead(Y,&y);
749: VecRestoreArray(F,&f);
750: /* Left hand side = ydot - f(y) */
751: VecAYPX(F,-1.0,Ydot);
752: return 0;
753: }
755: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
756: {
757: const PetscScalar *y;
758: PetscInt N,i,col[2];
759: PetscScalar value[2];
761: VecGetSize (Y,&N);
762: VecGetArrayRead(Y,&y);
763: i = 0;
764: value[0] = a+1; col[0] = 0;
765: value[1] = 0; col[1] = 1;
766: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
767: for (i = 0; i < N; i++) {
768: value[0] = -1; col[0] = i-1;
769: value[1] = a+1; col[1] = i;
770: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
771: }
772: i = N-1;
773: value[0] = -1; col[0] = N-2;
774: value[1] = a; col[1] = N-1;
775: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
776: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
777: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
778: VecRestoreArrayRead(Y,&y);
779: return 0;
780: }
782: /* Hull, 1972, Problem C2 */
784: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
785: {
786: const PetscScalar *y;
787: PetscScalar *f;
788: PetscInt N,i;
790: VecGetSize (Y,&N);
791: VecGetArrayRead(Y,&y);
792: VecGetArray(F,&f);
793: f[0] = -y[0];
794: for (i = 1; i < N-1; i++) {
795: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
796: }
797: f[N-1] = (PetscReal)(N-1)*y[N-2];
798: VecRestoreArrayRead(Y,&y);
799: VecRestoreArray(F,&f);
800: return 0;
801: }
803: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
804: {
805: PetscScalar *f;
806: const PetscScalar *y;
807: PetscInt N,i;
809: VecGetSize (Y,&N);
810: VecGetArrayRead(Y,&y);
811: VecGetArray(F,&f);
812: f[0] = -y[0];
813: for (i = 1; i < N-1; i++) {
814: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
815: }
816: f[N-1] = (PetscReal)(N-1)*y[N-2];
817: VecRestoreArrayRead(Y,&y);
818: VecRestoreArray(F,&f);
819: /* Left hand side = ydot - f(y) */
820: VecAYPX(F,-1.0,Ydot);
821: return 0;
822: }
824: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
825: {
826: const PetscScalar *y;
827: PetscInt N,i,col[2];
828: PetscScalar value[2];
830: VecGetSize (Y,&N);
831: VecGetArrayRead(Y,&y);
832: i = 0;
833: value[0] = a+1; col[0] = 0;
834: value[1] = 0; col[1] = 1;
835: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
836: for (i = 0; i < N; i++) {
837: value[0] = -(PetscReal) i; col[0] = i-1;
838: value[1] = a+(PetscReal)(i+1); col[1] = i;
839: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
840: }
841: i = N-1;
842: value[0] = -(PetscReal) (N-1); col[0] = N-2;
843: value[1] = a; col[1] = N-1;
844: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
845: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
846: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
847: VecRestoreArrayRead(Y,&y);
848: return 0;
849: }
851: /* Hull, 1972, Problem C3 and C4 */
853: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
854: {
855: PetscScalar *f;
856: const PetscScalar *y;
857: PetscInt N,i;
859: VecGetSize (Y,&N);
860: VecGetArrayRead(Y,&y);
861: VecGetArray(F,&f);
862: f[0] = -2.0*y[0] + y[1];
863: for (i = 1; i < N-1; i++) {
864: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
865: }
866: f[N-1] = y[N-2] - 2.0*y[N-1];
867: VecRestoreArrayRead(Y,&y);
868: VecRestoreArray(F,&f);
869: return 0;
870: }
872: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
873: {
874: PetscScalar *f;
875: const PetscScalar *y;
876: PetscInt N,i;
878: VecGetSize (Y,&N);
879: VecGetArrayRead(Y,&y);
880: VecGetArray(F,&f);
881: f[0] = -2.0*y[0] + y[1];
882: for (i = 1; i < N-1; i++) {
883: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
884: }
885: f[N-1] = y[N-2] - 2.0*y[N-1];
886: VecRestoreArrayRead(Y,&y);
887: VecRestoreArray(F,&f);
888: /* Left hand side = ydot - f(y) */
889: VecAYPX(F,-1.0,Ydot);
890: return 0;
891: }
893: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
894: {
895: const PetscScalar *y;
896: PetscScalar value[3];
897: PetscInt N,i,col[3];
899: VecGetSize (Y,&N);
900: VecGetArrayRead(Y,&y);
901: for (i = 0; i < N; i++) {
902: if (i == 0) {
903: value[0] = a+2; col[0] = i;
904: value[1] = -1; col[1] = i+1;
905: value[2] = 0; col[2] = i+2;
906: } else if (i == N-1) {
907: value[0] = 0; col[0] = i-2;
908: value[1] = -1; col[1] = i-1;
909: value[2] = a+2; col[2] = i;
910: } else {
911: value[0] = -1; col[0] = i-1;
912: value[1] = a+2; col[1] = i;
913: value[2] = -1; col[2] = i+1;
914: }
915: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
916: }
917: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
918: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
919: VecRestoreArrayRead(Y,&y);
920: return 0;
921: }
923: /***************************************************************************/
925: /* Sets the initial solution for the IVP and sets up the function pointers*/
926: PetscErrorCode Initialize(Vec Y, void* s)
927: {
928: char *p = (char*) s;
929: PetscScalar *y;
930: PetscReal t0;
931: PetscInt N = GetSize((const char *)s);
932: PetscBool flg;
934: VecZeroEntries(Y);
935: VecGetArray(Y,&y);
936: if (!strcmp(p,"hull1972a1")) {
937: y[0] = 1.0;
938: RHSFunction = RHSFunction_Hull1972A1;
939: RHSJacobian = RHSJacobian_Hull1972A1;
940: IFunction = IFunction_Hull1972A1;
941: IJacobian = IJacobian_Hull1972A1;
942: } else if (!strcmp(p,"hull1972a2")) {
943: y[0] = 1.0;
944: RHSFunction = RHSFunction_Hull1972A2;
945: RHSJacobian = RHSJacobian_Hull1972A2;
946: IFunction = IFunction_Hull1972A2;
947: IJacobian = IJacobian_Hull1972A2;
948: } else if (!strcmp(p,"hull1972a3")) {
949: y[0] = 1.0;
950: RHSFunction = RHSFunction_Hull1972A3;
951: RHSJacobian = RHSJacobian_Hull1972A3;
952: IFunction = IFunction_Hull1972A3;
953: IJacobian = IJacobian_Hull1972A3;
954: } else if (!strcmp(p,"hull1972a4")) {
955: y[0] = 1.0;
956: RHSFunction = RHSFunction_Hull1972A4;
957: RHSJacobian = RHSJacobian_Hull1972A4;
958: IFunction = IFunction_Hull1972A4;
959: IJacobian = IJacobian_Hull1972A4;
960: } else if (!strcmp(p,"hull1972a5")) {
961: y[0] = 4.0;
962: RHSFunction = RHSFunction_Hull1972A5;
963: RHSJacobian = RHSJacobian_Hull1972A5;
964: IFunction = IFunction_Hull1972A5;
965: IJacobian = IJacobian_Hull1972A5;
966: } else if (!strcmp(p,"hull1972b1")) {
967: y[0] = 1.0;
968: y[1] = 3.0;
969: RHSFunction = RHSFunction_Hull1972B1;
970: IFunction = IFunction_Hull1972B1;
971: IJacobian = IJacobian_Hull1972B1;
972: } else if (!strcmp(p,"hull1972b2")) {
973: y[0] = 2.0;
974: y[1] = 0.0;
975: y[2] = 1.0;
976: RHSFunction = RHSFunction_Hull1972B2;
977: IFunction = IFunction_Hull1972B2;
978: IJacobian = IJacobian_Hull1972B2;
979: } else if (!strcmp(p,"hull1972b3")) {
980: y[0] = 1.0;
981: y[1] = 0.0;
982: y[2] = 0.0;
983: RHSFunction = RHSFunction_Hull1972B3;
984: IFunction = IFunction_Hull1972B3;
985: IJacobian = IJacobian_Hull1972B3;
986: } else if (!strcmp(p,"hull1972b4")) {
987: y[0] = 3.0;
988: y[1] = 0.0;
989: y[2] = 0.0;
990: RHSFunction = RHSFunction_Hull1972B4;
991: IFunction = IFunction_Hull1972B4;
992: IJacobian = IJacobian_Hull1972B4;
993: } else if (!strcmp(p,"hull1972b5")) {
994: y[0] = 0.0;
995: y[1] = 1.0;
996: y[2] = 1.0;
997: RHSFunction = RHSFunction_Hull1972B5;
998: IFunction = IFunction_Hull1972B5;
999: IJacobian = IJacobian_Hull1972B5;
1000: } else if (!strcmp(p,"kulik2013i")) {
1001: t0=0.;
1002: y[0] = PetscExpReal(PetscSinReal(t0*t0));
1003: y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1004: y[2] = PetscSinReal(t0*t0)+1.0;
1005: y[3] = PetscCosReal(t0*t0);
1006: RHSFunction = RHSFunction_Kulikov2013I;
1007: RHSJacobian = RHSJacobian_Kulikov2013I;
1008: IFunction = IFunction_Kulikov2013I;
1009: IJacobian = IJacobian_Kulikov2013I;
1010: } else if (!strcmp(p,"hull1972c1")) {
1011: y[0] = 1.0;
1012: RHSFunction = RHSFunction_Hull1972C1;
1013: IFunction = IFunction_Hull1972C1;
1014: IJacobian = IJacobian_Hull1972C1;
1015: } else if (!strcmp(p,"hull1972c2")) {
1016: y[0] = 1.0;
1017: RHSFunction = RHSFunction_Hull1972C2;
1018: IFunction = IFunction_Hull1972C2;
1019: IJacobian = IJacobian_Hull1972C2;
1020: } else if ((!strcmp(p,"hull1972c3"))
1021: ||(!strcmp(p,"hull1972c4"))) {
1022: y[0] = 1.0;
1023: RHSFunction = RHSFunction_Hull1972C34;
1024: IFunction = IFunction_Hull1972C34;
1025: IJacobian = IJacobian_Hull1972C34;
1026: }
1027: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1029: VecRestoreArray(Y,&y);
1030: return 0;
1031: }
1033: /* Calculates the exact solution to problems that have one */
1034: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1035: {
1036: char *p = (char*) s;
1037: PetscScalar *y;
1039: if (!strcmp(p,"hull1972a1")) {
1040: VecGetArray(Y,&y);
1041: y[0] = PetscExpReal(-t);
1042: *flag = PETSC_TRUE;
1043: VecRestoreArray(Y,&y);
1044: } else if (!strcmp(p,"hull1972a2")) {
1045: VecGetArray(Y,&y);
1046: y[0] = 1.0/PetscSqrtReal(t+1);
1047: *flag = PETSC_TRUE;
1048: VecRestoreArray(Y,&y);
1049: } else if (!strcmp(p,"hull1972a3")) {
1050: VecGetArray(Y,&y);
1051: y[0] = PetscExpReal(PetscSinReal(t));
1052: *flag = PETSC_TRUE;
1053: VecRestoreArray(Y,&y);
1054: } else if (!strcmp(p,"hull1972a4")) {
1055: VecGetArray(Y,&y);
1056: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1057: *flag = PETSC_TRUE;
1058: VecRestoreArray(Y,&y);
1059: } else if (!strcmp(p,"kulik2013i")) {
1060: VecGetArray(Y,&y);
1061: y[0] = PetscExpReal(PetscSinReal(t*t));
1062: y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1063: y[2] = PetscSinReal(t*t)+1.0;
1064: y[3] = PetscCosReal(t*t);
1065: *flag = PETSC_TRUE;
1066: VecRestoreArray(Y,&y);
1067: } else {
1068: VecSet(Y,0);
1069: *flag = PETSC_FALSE;
1070: }
1071: return 0;
1072: }
1074: /* Solves the specified ODE and computes the error if exact solution is available */
1075: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1076: {
1077: TS ts; /* time-integrator */
1078: Vec Y; /* Solution vector */
1079: Vec Yex; /* Exact solution */
1080: PetscInt N; /* Size of the system of equations */
1081: TSType time_scheme; /* Type of time-integration scheme */
1082: Mat Jac = NULL; /* Jacobian matrix */
1083: Vec Yerr; /* Auxiliary solution vector */
1084: PetscReal err_norm; /* Estimated error norm */
1085: PetscReal final_time; /* Actual final time from the integrator */
1087: N = GetSize((const char *)&ptype[0]);
1089: VecCreate(PETSC_COMM_WORLD,&Y);
1090: VecSetSizes(Y,N,PETSC_DECIDE);
1091: VecSetUp(Y);
1092: VecSet(Y,0);
1094: /* Initialize the problem */
1095: Initialize(Y,&ptype[0]);
1097: /* Create and initialize the time-integrator */
1098: TSCreate(PETSC_COMM_WORLD,&ts);
1099: /* Default time integration options */
1100: TSSetType(ts,TSRK);
1101: TSSetMaxSteps(ts,maxiter);
1102: TSSetMaxTime(ts,tfinal);
1103: TSSetTimeStep(ts,dt);
1104: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1105: /* Read command line options for time integration */
1106: TSSetFromOptions(ts);
1107: /* Set solution vector */
1108: TSSetSolution(ts,Y);
1109: /* Specify left/right-hand side functions */
1110: TSGetType(ts,&time_scheme);
1112: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1113: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1114: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1115: MatCreate(PETSC_COMM_WORLD,&Jac);
1116: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1117: MatSetFromOptions(Jac);
1118: MatSetUp(Jac);
1119: TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1120: } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1121: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1122: /* and its Jacobian function */
1123: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1124: MatCreate(PETSC_COMM_WORLD,&Jac);
1125: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1126: MatSetFromOptions(Jac);
1127: MatSetUp(Jac);
1128: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1129: }
1131: /* Solve */
1132: TSSolve(ts,Y);
1133: TSGetTime(ts,&final_time);
1135: /* Get the estimated error, if available */
1136: VecDuplicate(Y,&Yerr);
1137: VecZeroEntries(Yerr);
1138: TSGetTimeError(ts,0,&Yerr);
1139: VecNorm(Yerr,NORM_2,&err_norm);
1140: VecDestroy(&Yerr);
1141: PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);
1143: /* Exact solution */
1144: VecDuplicate(Y,&Yex);
1145: if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1146: 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);
1147: }
1148: ExactSolution(Yex,&ptype[0],final_time,exact_flag);
1150: /* Calculate Error */
1151: VecAYPX(Yex,-1.0,Y);
1152: VecNorm(Yex,NORM_2,error);
1153: *error = PetscSqrtReal(((*error)*(*error))/N);
1155: /* Clean up and finalize */
1156: MatDestroy(&Jac);
1157: TSDestroy(&ts);
1158: VecDestroy(&Yex);
1159: VecDestroy(&Y);
1161: return 0;
1162: }
1164: int main(int argc, char **argv)
1165: {
1166: PetscErrorCode ierr; /* Error code */
1167: char ptype[256] = "hull1972a1"; /* Problem specification */
1168: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1169: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1170: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1171: PetscReal dt;
1172: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1173: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1174: PetscReal *error; /* Array to store the errors for convergence analysis */
1175: PetscMPIInt size; /* No of processors */
1176: PetscBool flag; /* Flag denoting availability of exact solution */
1177: PetscInt r;
1179: /* Initialize program */
1180: PetscInitialize(&argc,&argv,(char*)0,help);
1182: /* Check if running with only 1 proc */
1183: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1186: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1187: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1188: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1189: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1190: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1191: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1192: PetscOptionsEnd();
1194: PetscMalloc1(n_refine,&error);
1195: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1196: error[r] = 0;
1197: if (r > 0) dt /= refine_fac;
1199: 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]));
1200: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1201: if (flag) {
1202: /* If exact solution available for the specified ODE */
1203: if (r > 0) {
1204: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1205: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1206: } else {
1207: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1208: }
1209: }
1210: }
1211: PetscFree(error);
1212: PetscFinalize();
1213: return 0;
1214: }
1216: /*TEST
1218: test:
1219: suffix: 2
1220: args: -ts_type glee -final_time 5 -ts_adapt_type none
1221: timeoutfactor: 3
1222: requires: !single
1224: test:
1225: suffix: 3
1226: 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
1227: timeoutfactor: 3
1228: requires: !single
1230: test:
1231: suffix: 4
1232: 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
1233: timeoutfactor: 3
1234: requires: !single !__float128
1236: TEST*/