Actual source code: ex31.c
petsc-3.8.4 2018-03-24
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>
45: #if defined(PETSC_HAVE_STRING_H)
46: #include <string.h> /* strcmp */
47: #endif
49: /* Function declarations */
50: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
51: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
52: PetscErrorCode (*IFunction) (TS,PetscReal,Vec,Vec,Vec,void*);
53: PetscErrorCode (*IJacobian) (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
55: /* Returns the size of the system of equations depending on problem specification */
56: PetscInt GetSize(const char *p)
57: {
59: if ((!strcmp(p,"hull1972a1"))
60: ||(!strcmp(p,"hull1972a2"))
61: ||(!strcmp(p,"hull1972a3"))
62: ||(!strcmp(p,"hull1972a4"))
63: ||(!strcmp(p,"hull1972a5")) ) PetscFunctionReturn(1);
64: else if (!strcmp(p,"hull1972b1") ) PetscFunctionReturn(2);
65: else if ((!strcmp(p,"hull1972b2"))
66: ||(!strcmp(p,"hull1972b3"))
67: ||(!strcmp(p,"hull1972b4"))
68: ||(!strcmp(p,"hull1972b5")) ) PetscFunctionReturn(3);
69: else if ((!strcmp(p,"kulik2013i")) ) PetscFunctionReturn(4);
70: else if ((!strcmp(p,"hull1972c1"))
71: ||(!strcmp(p,"hull1972c2"))
72: ||(!strcmp(p,"hull1972c3")) ) PetscFunctionReturn(10);
73: else if (!strcmp(p,"hull1972c4") ) PetscFunctionReturn(51);
74: else PetscFunctionReturn(-1);
75: }
77: /****************************************************************/
79: /* Problem specific functions */
81: /* Hull, 1972, Problem A1 */
83: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
84: {
85: PetscErrorCode ierr;
86: PetscScalar *f;
87: const PetscScalar *y;
90: VecGetArrayRead(Y,&y);
91: VecGetArray(F,&f);
92: f[0] = -y[0];
93: VecRestoreArrayRead(Y,&y);
94: VecRestoreArray(F,&f);
95: return(0);
96: }
98: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
99: {
100: PetscErrorCode ierr;
101: const PetscScalar *y;
102: PetscInt row = 0,col = 0;
103: PetscScalar value = -1.0;
106: VecGetArrayRead(Y,&y);
107: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
108: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
110: VecRestoreArrayRead(Y,&y);
111: return(0);
112: }
114: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
115: {
116: PetscErrorCode ierr;
117: const PetscScalar *y;
118: PetscScalar *f;
121: VecGetArrayRead(Y,&y);
122: VecGetArray(F,&f);
123: f[0] = -y[0];
124: VecRestoreArrayRead(Y,&y);
125: VecRestoreArray(F,&f);
126: /* Left hand side = ydot - f(y) */
127: VecAYPX(F,-1.0,Ydot);
128: return(0);
129: }
131: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
132: {
133: PetscErrorCode ierr;
134: const PetscScalar *y;
135: PetscInt row = 0,col = 0;
136: PetscScalar value = a - 1.0;
139: VecGetArrayRead(Y,&y);
140: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
141: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
142: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
143: VecRestoreArrayRead(Y,&y);
144: return(0);
145: }
147: /* Hull, 1972, Problem A2 */
149: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
150: {
151: PetscErrorCode ierr;
152: const PetscScalar *y;
153: PetscScalar *f;
156: VecGetArrayRead(Y,&y);
157: VecGetArray(F,&f);
158: f[0] = -0.5*y[0]*y[0]*y[0];
159: VecRestoreArrayRead(Y,&y);
160: VecRestoreArray(F,&f);
161: return(0);
162: }
164: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
165: {
166: PetscErrorCode ierr;
167: const PetscScalar *y;
168: PetscInt row = 0,col = 0;
169: PetscScalar value;
172: VecGetArrayRead(Y,&y);
173: value = -0.5*3.0*y[0]*y[0];
174: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
175: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
176: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
177: VecRestoreArrayRead(Y,&y);
178: return(0);
179: }
181: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
182: {
183: PetscErrorCode ierr;
184: PetscScalar *f;
185: const PetscScalar *y;
188: VecGetArrayRead(Y,&y);
189: VecGetArray(F,&f);
190: f[0] = -0.5*y[0]*y[0]*y[0];
191: VecRestoreArrayRead(Y,&y);
192: VecRestoreArray(F,&f);
193: /* Left hand side = ydot - f(y) */
194: VecAYPX(F,-1.0,Ydot);
195: return(0);
196: }
198: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
199: {
200: PetscErrorCode ierr;
201: const PetscScalar *y;
202: PetscInt row = 0,col = 0;
203: PetscScalar value;
206: VecGetArrayRead(Y,&y);
207: value = a + 0.5*3.0*y[0]*y[0];
208: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
209: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
210: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
211: VecRestoreArrayRead(Y,&y);
212: return(0);
213: }
215: /* Hull, 1972, Problem A3 */
217: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
218: {
219: PetscErrorCode ierr;
220: const PetscScalar *y;
221: PetscScalar *f;
224: VecGetArrayRead(Y,&y);
225: VecGetArray(F,&f);
226: f[0] = y[0]*PetscCosReal(t);
227: VecRestoreArrayRead(Y,&y);
228: VecRestoreArray(F,&f);
229: return(0);
230: }
232: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
233: {
234: PetscErrorCode ierr;
235: const PetscScalar *y;
236: PetscInt row = 0,col = 0;
237: PetscScalar value = PetscCosReal(t);
240: VecGetArrayRead(Y,&y);
241: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
242: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
243: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
244: VecRestoreArrayRead(Y,&y);
245: return(0);
246: }
248: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
249: {
250: PetscErrorCode ierr;
251: PetscScalar *f;
252: const PetscScalar *y;
255: VecGetArrayRead(Y,&y);
256: VecGetArray(F,&f);
257: f[0] = y[0]*PetscCosReal(t);
258: VecRestoreArrayRead(Y,&y);
259: VecRestoreArray(F,&f);
260: /* Left hand side = ydot - f(y) */
261: VecAYPX(F,-1.0,Ydot);
262: return(0);
263: }
265: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
266: {
267: PetscErrorCode ierr;
268: const PetscScalar *y;
269: PetscInt row = 0,col = 0;
270: PetscScalar value = a - PetscCosReal(t);
273: VecGetArrayRead(Y,&y);
274: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
275: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
276: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
277: VecRestoreArrayRead(Y,&y);
278: return(0);
279: }
281: /* Hull, 1972, Problem A4 */
283: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
284: {
285: PetscErrorCode ierr;
286: PetscScalar *f;
287: const PetscScalar *y;
290: VecGetArrayRead(Y,&y);
291: VecGetArray(F,&f);
292: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
293: VecRestoreArrayRead(Y,&y);
294: VecRestoreArray(F,&f);
295: return(0);
296: }
298: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
299: {
300: PetscErrorCode ierr;
301: const PetscScalar *y;
302: PetscInt row = 0,col = 0;
303: PetscScalar value;
306: VecGetArrayRead(Y,&y);
307: value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
308: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
309: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
310: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
311: VecRestoreArrayRead(Y,&y);
312: return(0);
313: }
315: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
316: {
317: PetscErrorCode ierr;
318: PetscScalar *f;
319: const PetscScalar *y;
322: VecGetArrayRead(Y,&y);
323: VecGetArray(F,&f);
324: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
325: VecRestoreArrayRead(Y,&y);
326: VecRestoreArray(F,&f);
327: /* Left hand side = ydot - f(y) */
328: VecAYPX(F,-1.0,Ydot);
329: return(0);
330: }
332: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
333: {
334: PetscErrorCode ierr;
335: const PetscScalar *y;
336: PetscInt row = 0,col = 0;
337: PetscScalar value;
340: VecGetArrayRead(Y,&y);
341: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
342: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
343: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
344: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
345: VecRestoreArrayRead(Y,&y);
346: return(0);
347: }
349: /* Hull, 1972, Problem A5 */
351: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
352: {
353: PetscErrorCode ierr;
354: PetscScalar *f;
355: const PetscScalar *y;
358: VecGetArrayRead(Y,&y);
359: VecGetArray(F,&f);
360: f[0] = (y[0]-t)/(y[0]+t);
361: VecRestoreArrayRead(Y,&y);
362: VecRestoreArray(F,&f);
363: return(0);
364: }
366: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
367: {
368: PetscErrorCode ierr;
369: const PetscScalar *y;
370: PetscInt row = 0,col = 0;
371: PetscScalar value;
374: VecGetArrayRead(Y,&y);
375: value = 2*t/((t+y[0])*(t+y[0]));
376: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
377: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
378: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
379: VecRestoreArrayRead(Y,&y);
380: return(0);
381: }
383: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
384: {
385: PetscErrorCode ierr;
386: PetscScalar *f;
387: const PetscScalar *y;
390: VecGetArrayRead(Y,&y);
391: VecGetArray(F,&f);
392: f[0] = (y[0]-t)/(y[0]+t);
393: VecRestoreArrayRead(Y,&y);
394: VecRestoreArray(F,&f);
395: /* Left hand side = ydot - f(y) */
396: VecAYPX(F,-1.0,Ydot);
397: return(0);
398: }
400: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
401: {
402: PetscErrorCode ierr;
403: const PetscScalar *y;
404: PetscInt row = 0,col = 0;
405: PetscScalar value;
408: VecGetArrayRead(Y,&y);
409: value = a - 2*t/((t+y[0])*(t+y[0]));
410: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
411: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
412: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
413: VecRestoreArrayRead(Y,&y);
414: return(0);
415: }
417: /* Hull, 1972, Problem B1 */
419: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
420: {
421: PetscErrorCode ierr;
422: PetscScalar *f;
423: const PetscScalar *y;
426: VecGetArrayRead(Y,&y);
427: VecGetArray(F,&f);
428: f[0] = 2.0*(y[0] - y[0]*y[1]);
429: f[1] = -(y[1]-y[0]*y[1]);
430: VecRestoreArrayRead(Y,&y);
431: VecRestoreArray(F,&f);
432: return(0);
433: }
435: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
436: {
437: PetscErrorCode ierr;
438: PetscScalar *f;
439: const PetscScalar *y;
442: VecGetArrayRead(Y,&y);
443: VecGetArray(F,&f);
444: f[0] = 2.0*(y[0] - y[0]*y[1]);
445: f[1] = -(y[1]-y[0]*y[1]);
446: VecRestoreArrayRead(Y,&y);
447: VecRestoreArray(F,&f);
448: /* Left hand side = ydot - f(y) */
449: VecAYPX(F,-1.0,Ydot);
450: return(0);
451: }
453: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
454: {
455: PetscErrorCode ierr;
456: const PetscScalar *y;
457: PetscInt row[2] = {0,1};
458: PetscScalar value[2][2];
461: VecGetArrayRead(Y,&y);
462: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
463: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
464: MatSetValues(A,2,&row[0],2,&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 B2 */
473: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
474: {
475: PetscErrorCode ierr;
476: PetscScalar *f;
477: const PetscScalar *y;
480: VecGetArrayRead(Y,&y);
481: VecGetArray(F,&f);
482: f[0] = -y[0] + y[1];
483: f[1] = y[0] - 2.0*y[1] + y[2];
484: f[2] = y[1] - y[2];
485: VecRestoreArrayRead(Y,&y);
486: VecRestoreArray(F,&f);
487: return(0);
488: }
490: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
491: {
492: PetscErrorCode ierr;
493: PetscScalar *f;
494: const PetscScalar *y;
497: VecGetArrayRead(Y,&y);
498: VecGetArray(F,&f);
499: f[0] = -y[0] + y[1];
500: f[1] = y[0] - 2.0*y[1] + y[2];
501: f[2] = y[1] - y[2];
502: VecRestoreArrayRead(Y,&y);
503: VecRestoreArray(F,&f);
504: /* Left hand side = ydot - f(y) */
505: VecAYPX(F,-1.0,Ydot);
506: return(0);
507: }
509: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
510: {
511: PetscErrorCode ierr;
512: const PetscScalar *y;
513: PetscInt row[3] = {0,1,2};
514: PetscScalar value[3][3];
517: VecGetArrayRead(Y,&y);
518: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
519: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
520: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
521: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
522: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
523: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
524: VecRestoreArrayRead(Y,&y);
525: return(0);
526: }
528: /* Hull, 1972, Problem B3 */
530: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
531: {
532: PetscErrorCode ierr;
533: PetscScalar *f;
534: const PetscScalar *y;
537: VecGetArrayRead(Y,&y);
538: VecGetArray(F,&f);
539: f[0] = -y[0];
540: f[1] = y[0] - y[1]*y[1];
541: f[2] = y[1]*y[1];
542: VecRestoreArrayRead(Y,&y);
543: VecRestoreArray(F,&f);
544: return(0);
545: }
547: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
548: {
549: PetscErrorCode ierr;
550: PetscScalar *f;
551: const PetscScalar *y;
554: VecGetArrayRead(Y,&y);
555: VecGetArray(F,&f);
556: f[0] = -y[0];
557: f[1] = y[0] - y[1]*y[1];
558: f[2] = y[1]*y[1];
559: VecRestoreArrayRead(Y,&y);
560: VecRestoreArray(F,&f);
561: /* Left hand side = ydot - f(y) */
562: VecAYPX(F,-1.0,Ydot);
563: return(0);
564: }
566: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
567: {
568: PetscErrorCode ierr;
569: const PetscScalar *y;
570: PetscInt row[3] = {0,1,2};
571: PetscScalar value[3][3];
574: VecGetArrayRead(Y,&y);
575: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
576: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
577: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
578: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
579: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
580: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
581: VecRestoreArrayRead(Y,&y);
582: return(0);
583: }
585: /* Hull, 1972, Problem B4 */
587: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
588: {
589: PetscErrorCode ierr;
590: PetscScalar *f;
591: const PetscScalar *y;
594: VecGetArrayRead(Y,&y);
595: VecGetArray(F,&f);
596: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
597: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
598: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
599: VecRestoreArrayRead(Y,&y);
600: VecRestoreArray(F,&f);
601: return(0);
602: }
604: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
605: {
606: PetscErrorCode ierr;
607: PetscScalar *f;
608: const PetscScalar *y;
611: VecGetArrayRead(Y,&y);
612: VecGetArray(F,&f);
613: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
614: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
615: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
616: VecRestoreArrayRead(Y,&y);
617: VecRestoreArray(F,&f);
618: /* Left hand side = ydot - f(y) */
619: VecAYPX(F,-1.0,Ydot);
620: return(0);
621: }
623: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
624: {
625: PetscErrorCode ierr;
626: const PetscScalar *y;
627: PetscInt row[3] = {0,1,2};
628: PetscScalar value[3][3],fac,fac2;
631: VecGetArrayRead(Y,&y);
632: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
633: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
634: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
635: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
636: value[0][2] = y[0]*fac2;
637: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
638: value[1][1] = a + y[0]*y[0]*y[2]*fac;
639: value[1][2] = y[1]*fac2;
640: value[2][0] = -y[1]*y[1]*fac;
641: value[2][1] = y[0]*y[1]*fac;
642: value[2][2] = a;
643: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
644: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
645: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
646: VecRestoreArrayRead(Y,&y);
647: return(0);
648: }
650: /* Hull, 1972, Problem B5 */
652: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
653: {
654: PetscErrorCode ierr;
655: PetscScalar *f;
656: const PetscScalar *y;
659: VecGetArrayRead(Y,&y);
660: VecGetArray(F,&f);
661: f[0] = y[1]*y[2];
662: f[1] = -y[0]*y[2];
663: f[2] = -0.51*y[0]*y[1];
664: VecRestoreArrayRead(Y,&y);
665: VecRestoreArray(F,&f);
666: return(0);
667: }
669: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
670: {
671: PetscErrorCode ierr;
672: PetscScalar *f;
673: const PetscScalar *y;
676: VecGetArrayRead(Y,&y);
677: VecGetArray(F,&f);
678: f[0] = y[1]*y[2];
679: f[1] = -y[0]*y[2];
680: f[2] = -0.51*y[0]*y[1];
681: VecRestoreArrayRead(Y,&y);
682: VecRestoreArray(F,&f);
683: /* Left hand side = ydot - f(y) */
684: VecAYPX(F,-1.0,Ydot);
685: return(0);
686: }
688: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
689: {
690: PetscErrorCode ierr;
691: const PetscScalar *y;
692: PetscInt row[3] = {0,1,2};
693: PetscScalar value[3][3];
696: VecGetArrayRead(Y,&y);
697: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
698: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
699: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
700: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
701: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
702: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
703: VecRestoreArrayRead(Y,&y);
704: return(0);
705: }
708: /* Kulikov, 2013, Problem I */
710: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
711: {
712: PetscErrorCode ierr;
713: PetscScalar *f;
714: const PetscScalar *y;
717: VecGetArrayRead(Y,&y);
718: VecGetArray(F,&f);
719: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
720: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
721: f[2] = 2.*t*y[3];
722: f[3] = -2.*t*PetscLogScalar(y[0]);
723: VecRestoreArrayRead(Y,&y);
724: VecRestoreArray(F,&f);
725: return(0);
726: }
728: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
729: {
730: PetscErrorCode ierr;
731: const PetscScalar *y;
732: PetscInt row[4] = {0,1,2,3};
733: PetscScalar value[4][4];
734: PetscScalar m1,m2;
736: VecGetArrayRead(Y,&y);
737: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
738: m2=2.*t*PetscPowScalar(y[1],1./5.);
739: value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
740: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
741: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
742: value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
743: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t;
744: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.;
745: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
746: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
747: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
748: VecRestoreArrayRead(Y,&y);
749: return(0);
750: }
752: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
753: {
754: PetscErrorCode ierr;
755: PetscScalar *f;
756: const PetscScalar *y;
759: VecGetArrayRead(Y,&y);
760: VecGetArray(F,&f);
761: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
762: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
763: f[2] = 2.*t*y[3];
764: f[3] = -2.*t*PetscLogScalar(y[0]);
765: VecRestoreArrayRead(Y,&y);
766: VecRestoreArray(F,&f);
767: /* Left hand side = ydot - f(y) */
768: VecAYPX(F,-1.0,Ydot);
769: return(0);
770: }
772: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
773: {
774: PetscErrorCode ierr;
775: const PetscScalar *y;
776: PetscInt row[4] = {0,1,2,3};
777: PetscScalar value[4][4];
778: PetscScalar m1,m2;
781: VecGetArrayRead(Y,&y);
782: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
783: m2=2.*t*PetscPowScalar(y[1],1./5.);
784: value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
785: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
786: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
787: value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2;
788: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t;
789: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a;
790: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
791: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
792: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
793: VecRestoreArrayRead(Y,&y);
794: return(0);
795: }
798: /* Hull, 1972, Problem C1 */
800: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
801: {
802: PetscErrorCode ierr;
803: PetscScalar *f;
804: const PetscScalar *y;
805: PetscInt N,i;
808: VecGetSize (Y,&N);
809: VecGetArrayRead(Y,&y);
810: VecGetArray(F,&f);
811: f[0] = -y[0];
812: for (i = 1; i < N-1; i++) {
813: f[i] = y[i-1] - y[i];
814: }
815: f[N-1] = y[N-2];
816: VecRestoreArrayRead(Y,&y);
817: VecRestoreArray(F,&f);
818: return(0);
819: }
821: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
822: {
823: PetscErrorCode ierr;
824: PetscScalar *f;
825: const PetscScalar *y;
826: PetscInt N,i;
829: VecGetSize (Y,&N);
830: VecGetArrayRead(Y,&y);
831: VecGetArray(F,&f);
832: f[0] = -y[0];
833: for (i = 1; i < N-1; i++) {
834: f[i] = y[i-1] - y[i];
835: }
836: f[N-1] = y[N-2];
837: VecRestoreArrayRead(Y,&y);
838: VecRestoreArray(F,&f);
839: /* Left hand side = ydot - f(y) */
840: VecAYPX(F,-1.0,Ydot);
841: return(0);
842: }
844: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
845: {
846: PetscErrorCode ierr;
847: const PetscScalar *y;
848: PetscInt N,i,col[2];
849: PetscScalar value[2];
852: VecGetSize (Y,&N);
853: VecGetArrayRead(Y,&y);
854: i = 0;
855: value[0] = a+1; col[0] = 0;
856: value[1] = 0; col[1] = 1;
857: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
858: for (i = 0; i < N; i++) {
859: value[0] = -1; col[0] = i-1;
860: value[1] = a+1; col[1] = i;
861: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
862: }
863: i = N-1;
864: value[0] = -1; col[0] = N-2;
865: value[1] = a; col[1] = N-1;
866: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
867: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
868: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
869: VecRestoreArrayRead(Y,&y);
870: return(0);
871: }
873: /* Hull, 1972, Problem C2 */
875: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
876: {
877: PetscErrorCode ierr;
878: const PetscScalar *y;
879: PetscScalar *f;
880: PetscInt N,i;
883: VecGetSize (Y,&N);
884: VecGetArrayRead(Y,&y);
885: VecGetArray(F,&f);
886: f[0] = -y[0];
887: for (i = 1; i < N-1; i++) {
888: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
889: }
890: f[N-1] = (PetscReal)(N-1)*y[N-2];
891: VecRestoreArrayRead(Y,&y);
892: VecRestoreArray(F,&f);
893: return(0);
894: }
896: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
897: {
898: PetscErrorCode ierr;
899: PetscScalar *f;
900: const PetscScalar *y;
901: PetscInt N,i;
904: VecGetSize (Y,&N);
905: VecGetArrayRead(Y,&y);
906: VecGetArray(F,&f);
907: f[0] = -y[0];
908: for (i = 1; i < N-1; i++) {
909: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
910: }
911: f[N-1] = (PetscReal)(N-1)*y[N-2];
912: VecRestoreArrayRead(Y,&y);
913: VecRestoreArray(F,&f);
914: /* Left hand side = ydot - f(y) */
915: VecAYPX(F,-1.0,Ydot);
916: return(0);
917: }
919: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
920: {
921: PetscErrorCode ierr;
922: const PetscScalar *y;
923: PetscInt N,i,col[2];
924: PetscScalar value[2];
927: VecGetSize (Y,&N);
928: VecGetArrayRead(Y,&y);
929: i = 0;
930: value[0] = a+1; col[0] = 0;
931: value[1] = 0; col[1] = 1;
932: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
933: for (i = 0; i < N; i++) {
934: value[0] = -(PetscReal) i; col[0] = i-1;
935: value[1] = a+(PetscReal)(i+1); col[1] = i;
936: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
937: }
938: i = N-1;
939: value[0] = -(PetscReal) (N-1); col[0] = N-2;
940: value[1] = a; col[1] = N-1;
941: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
942: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
943: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
944: VecRestoreArrayRead(Y,&y);
945: return(0);
946: }
948: /* Hull, 1972, Problem C3 and C4 */
950: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
951: {
952: PetscErrorCode ierr;
953: PetscScalar *f;
954: const PetscScalar *y;
955: PetscInt N,i;
958: VecGetSize (Y,&N);
959: VecGetArrayRead(Y,&y);
960: VecGetArray(F,&f);
961: f[0] = -2.0*y[0] + y[1];
962: for (i = 1; i < N-1; i++) {
963: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
964: }
965: f[N-1] = y[N-2] - 2.0*y[N-1];
966: VecRestoreArrayRead(Y,&y);
967: VecRestoreArray(F,&f);
968: return(0);
969: }
971: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
972: {
973: PetscErrorCode ierr;
974: PetscScalar *f;
975: const PetscScalar *y;
976: PetscInt N,i;
979: VecGetSize (Y,&N);
980: VecGetArrayRead(Y,&y);
981: VecGetArray(F,&f);
982: f[0] = -2.0*y[0] + y[1];
983: for (i = 1; i < N-1; i++) {
984: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
985: }
986: f[N-1] = y[N-2] - 2.0*y[N-1];
987: VecRestoreArrayRead(Y,&y);
988: VecRestoreArray(F,&f);
989: /* Left hand side = ydot - f(y) */
990: VecAYPX(F,-1.0,Ydot);
991: return(0);
992: }
994: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
995: {
996: PetscErrorCode ierr;
997: const PetscScalar *y;
998: PetscScalar value[3];
999: PetscInt N,i,col[3];
1002: VecGetSize (Y,&N);
1003: VecGetArrayRead(Y,&y);
1004: for (i = 0; i < N; i++) {
1005: if (i == 0) {
1006: value[0] = a+2; col[0] = i;
1007: value[1] = -1; col[1] = i+1;
1008: value[2] = 0; col[2] = i+2;
1009: } else if (i == N-1) {
1010: value[0] = 0; col[0] = i-2;
1011: value[1] = -1; col[1] = i-1;
1012: value[2] = a+2; col[2] = i;
1013: } else {
1014: value[0] = -1; col[0] = i-1;
1015: value[1] = a+2; col[1] = i;
1016: value[2] = -1; col[2] = i+1;
1017: }
1018: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
1019: }
1020: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1021: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
1022: VecRestoreArrayRead(Y,&y);
1023: return(0);
1024: }
1026: /***************************************************************************/
1028: /* Sets the initial solution for the IVP and sets up the function pointers*/
1029: PetscErrorCode Initialize(Vec Y, void* s)
1030: {
1032: char *p = (char*) s;
1033: PetscScalar *y;
1034: PetscReal t0;
1035: PetscInt N = GetSize((const char *)s);
1036: PetscBool flg;
1039: VecZeroEntries(Y);
1040: VecGetArray(Y,&y);
1041: if (!strcmp(p,"hull1972a1")) {
1042: y[0] = 1.0;
1043: RHSFunction = RHSFunction_Hull1972A1;
1044: RHSJacobian = RHSJacobian_Hull1972A1;
1045: IFunction = IFunction_Hull1972A1;
1046: IJacobian = IJacobian_Hull1972A1;
1047: } else if (!strcmp(p,"hull1972a2")) {
1048: y[0] = 1.0;
1049: RHSFunction = RHSFunction_Hull1972A2;
1050: RHSJacobian = RHSJacobian_Hull1972A2;
1051: IFunction = IFunction_Hull1972A2;
1052: IJacobian = IJacobian_Hull1972A2;
1053: } else if (!strcmp(p,"hull1972a3")) {
1054: y[0] = 1.0;
1055: RHSFunction = RHSFunction_Hull1972A3;
1056: RHSJacobian = RHSJacobian_Hull1972A3;
1057: IFunction = IFunction_Hull1972A3;
1058: IJacobian = IJacobian_Hull1972A3;
1059: } else if (!strcmp(p,"hull1972a4")) {
1060: y[0] = 1.0;
1061: RHSFunction = RHSFunction_Hull1972A4;
1062: RHSJacobian = RHSJacobian_Hull1972A4;
1063: IFunction = IFunction_Hull1972A4;
1064: IJacobian = IJacobian_Hull1972A4;
1065: } else if (!strcmp(p,"hull1972a5")) {
1066: y[0] = 4.0;
1067: RHSFunction = RHSFunction_Hull1972A5;
1068: RHSJacobian = RHSJacobian_Hull1972A5;
1069: IFunction = IFunction_Hull1972A5;
1070: IJacobian = IJacobian_Hull1972A5;
1071: } else if (!strcmp(p,"hull1972b1")) {
1072: y[0] = 1.0;
1073: y[1] = 3.0;
1074: RHSFunction = RHSFunction_Hull1972B1;
1075: IFunction = IFunction_Hull1972B1;
1076: IJacobian = IJacobian_Hull1972B1;
1077: } else if (!strcmp(p,"hull1972b2")) {
1078: y[0] = 2.0;
1079: y[1] = 0.0;
1080: y[2] = 1.0;
1081: RHSFunction = RHSFunction_Hull1972B2;
1082: IFunction = IFunction_Hull1972B2;
1083: IJacobian = IJacobian_Hull1972B2;
1084: } else if (!strcmp(p,"hull1972b3")) {
1085: y[0] = 1.0;
1086: y[1] = 0.0;
1087: y[2] = 0.0;
1088: RHSFunction = RHSFunction_Hull1972B3;
1089: IFunction = IFunction_Hull1972B3;
1090: IJacobian = IJacobian_Hull1972B3;
1091: } else if (!strcmp(p,"hull1972b4")) {
1092: y[0] = 3.0;
1093: y[1] = 0.0;
1094: y[2] = 0.0;
1095: RHSFunction = RHSFunction_Hull1972B4;
1096: IFunction = IFunction_Hull1972B4;
1097: IJacobian = IJacobian_Hull1972B4;
1098: } else if (!strcmp(p,"hull1972b5")) {
1099: y[0] = 0.0;
1100: y[1] = 1.0;
1101: y[2] = 1.0;
1102: RHSFunction = RHSFunction_Hull1972B5;
1103: IFunction = IFunction_Hull1972B5;
1104: IJacobian = IJacobian_Hull1972B5;
1105: } else if (!strcmp(p,"kulik2013i")) {
1106: t0=0.;
1107: y[0] = PetscExpReal(PetscSinReal(t0*t0));
1108: y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1109: y[2] = PetscSinReal(t0*t0)+1.0;
1110: y[3] = PetscCosReal(t0*t0);
1111: RHSFunction = RHSFunction_Kulikov2013I;
1112: RHSJacobian = RHSJacobian_Kulikov2013I;
1113: IFunction = IFunction_Kulikov2013I;
1114: IJacobian = IJacobian_Kulikov2013I;
1115: } else if (!strcmp(p,"hull1972c1")) {
1116: y[0] = 1.0;
1117: RHSFunction = RHSFunction_Hull1972C1;
1118: IFunction = IFunction_Hull1972C1;
1119: IJacobian = IJacobian_Hull1972C1;
1120: } else if (!strcmp(p,"hull1972c2")) {
1121: y[0] = 1.0;
1122: RHSFunction = RHSFunction_Hull1972C2;
1123: IFunction = IFunction_Hull1972C2;
1124: IJacobian = IJacobian_Hull1972C2;
1125: } else if ((!strcmp(p,"hull1972c3"))
1126: ||(!strcmp(p,"hull1972c4"))){
1127: y[0] = 1.0;
1128: RHSFunction = RHSFunction_Hull1972C34;
1129: IFunction = IFunction_Hull1972C34;
1130: IJacobian = IJacobian_Hull1972C34;
1131: }
1132: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1133: if ((N != GetSize((const char*)s)) && flg) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Number of initial values %D does not match problem size %D.\n",N,GetSize((const char*)s));
1134: VecRestoreArray(Y,&y);
1135: return(0);
1136: }
1138: /* Calculates the exact solution to problems that have one */
1139: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1140: {
1142: char *p = (char*) s;
1143: PetscScalar *y;
1146: if (!strcmp(p,"hull1972a1")) {
1147: VecGetArray(Y,&y);
1148: y[0] = PetscExpReal(-t);
1149: *flag = PETSC_TRUE;
1150: VecRestoreArray(Y,&y);
1151: } else if (!strcmp(p,"hull1972a2")) {
1152: VecGetArray(Y,&y);
1153: y[0] = 1.0/PetscSqrtReal(t+1);
1154: *flag = PETSC_TRUE;
1155: VecRestoreArray(Y,&y);
1156: } else if (!strcmp(p,"hull1972a3")) {
1157: VecGetArray(Y,&y);
1158: y[0] = PetscExpReal(PetscSinReal(t));
1159: *flag = PETSC_TRUE;
1160: VecRestoreArray(Y,&y);
1161: } else if (!strcmp(p,"hull1972a4")) {
1162: VecGetArray(Y,&y);
1163: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1164: *flag = PETSC_TRUE;
1165: VecRestoreArray(Y,&y);
1166: } else if (!strcmp(p,"kulik2013i")) {
1167: VecGetArray(Y,&y);
1168: y[0] = PetscExpReal(PetscSinReal(t*t));
1169: y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1170: y[2] = PetscSinReal(t*t)+1.0;
1171: y[3] = PetscCosReal(t*t);
1172: *flag = PETSC_TRUE;
1173: VecRestoreArray(Y,&y);
1174: } else {
1175: VecSet(Y,0);
1176: *flag = PETSC_FALSE;
1177: }
1178: return(0);
1179: }
1181: /* Solves the specified ODE and computes the error if exact solution is available */
1182: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1183: {
1184: PetscErrorCode ierr; /* Error code */
1185: TS ts; /* time-integrator */
1186: Vec Y; /* Solution vector */
1187: Vec Yex; /* Exact solution */
1188: PetscInt N; /* Size of the system of equations */
1189: TSType time_scheme; /* Type of time-integration scheme */
1190: Mat Jac = NULL; /* Jacobian matrix */
1191: Vec Yerr; /* Auxiliary solution vector */
1192: PetscReal err_norm; /* Estimated error norm */
1195: N = GetSize((const char *)&ptype[0]);
1196: if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1197: VecCreate(PETSC_COMM_WORLD,&Y);
1198: VecSetSizes(Y,N,PETSC_DECIDE);
1199: VecSetUp(Y);
1200: VecSet(Y,0);
1202: /* Initialize the problem */
1203: Initialize(Y,&ptype[0]);
1205: /* Create and initialize the time-integrator */
1206: TSCreate(PETSC_COMM_WORLD,&ts);
1207: /* Default time integration options */
1208: TSSetType(ts,TSRK);
1209: TSSetMaxSteps(ts,maxiter);
1210: TSSetMaxTime(ts,tfinal);
1211: TSSetTimeStep(ts,dt);
1212: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1213: /* Read command line options for time integration */
1214: TSSetFromOptions(ts);
1215: /* Set solution vector */
1216: TSSetSolution(ts,Y);
1217: /* Specify left/right-hand side functions */
1218: TSGetType(ts,&time_scheme);
1220: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1221: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1222: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1223: MatCreate(PETSC_COMM_WORLD,&Jac);
1224: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1225: MatSetFromOptions(Jac);
1226: MatSetUp(Jac);
1227: TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1228: } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1229: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1230: /* and its Jacobian function */
1231: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1232: MatCreate(PETSC_COMM_WORLD,&Jac);
1233: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1234: MatSetFromOptions(Jac);
1235: MatSetUp(Jac);
1236: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1237: }
1239: /* Solve */
1240: TSSolve(ts,Y);
1242: /* Get the estimated error, if available */
1243: VecDuplicate(Y,&Yerr);
1244: VecZeroEntries(Yerr);
1245: TSGetTimeError(ts,0,&Yerr);
1246: VecNorm(Yerr,NORM_2,&err_norm);
1247: VecDestroy(&Yerr);
1248: PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);
1250: /* Exact solution */
1251: VecDuplicate(Y,&Yex);
1252: ExactSolution(Yex,&ptype[0],tfinal,exact_flag);
1254: /* Calculate Error */
1255: VecAYPX(Yex,-1.0,Y);
1256: VecNorm(Yex,NORM_2,error);
1257: *error = PetscSqrtReal(((*error)*(*error))/N);
1259: /* Clean up and finalize */
1260: MatDestroy(&Jac);
1261: TSDestroy(&ts);
1262: VecDestroy(&Yex);
1263: VecDestroy(&Y);
1265: return(0);
1266: }
1268: int main(int argc, char **argv)
1269: {
1270: PetscErrorCode ierr; /* Error code */
1271: char ptype[256] = "hull1972a1"; /* Problem specification */
1272: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1273: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1274: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1275: PetscReal dt;
1276: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1277: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1278: PetscReal *error; /* Array to store the errors for convergence analysis */
1279: PetscMPIInt size; /* No of processors */
1280: PetscBool flag; /* Flag denoting availability of exact solution */
1281: PetscInt r;
1283: /* Initialize program */
1284: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
1286: /* Check if running with only 1 proc */
1287: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1288: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1290: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
1291: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1292: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1293: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1294: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1295: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1296: PetscOptionsEnd();
1298: PetscMalloc1(n_refine,&error);
1299: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1300: error[r] = 0;
1301: if (r > 0) dt /= refine_fac;
1303: 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]));
1304: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1305: if (flag) {
1306: /* If exact solution available for the specified ODE */
1307: if (r > 0) {
1308: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1309: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1310: } else {
1311: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1312: }
1313: }
1314: }
1315: PetscFree(error);
1316: PetscFinalize();
1317: return ierr;
1318: }