Actual source code: ex31.c
petsc-3.5.4 2015-05-23
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"
36: */
38: #include <petscts.h>
39: #if defined(PETSC_HAVE_STRING_H)
40: #include <string.h> /* strcmp */
41: #endif
43: /* Function declarations */
44: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
45: PetscErrorCode (*IFunction) (TS,PetscReal,Vec,Vec,Vec,void*);
46: PetscErrorCode (*IJacobian) (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
50: /* Returns the size of the system of equations depending on problem specification */
51: PetscInt GetSize(const char *p)
52: {
54: if ((!strcmp(p,"hull1972a1"))
55: ||(!strcmp(p,"hull1972a2"))
56: ||(!strcmp(p,"hull1972a3"))
57: ||(!strcmp(p,"hull1972a4"))
58: ||(!strcmp(p,"hull1972a5")) ) PetscFunctionReturn(1);
59: else if (!strcmp(p,"hull1972b1") ) PetscFunctionReturn(2);
60: else if ((!strcmp(p,"hull1972b2"))
61: ||(!strcmp(p,"hull1972b3"))
62: ||(!strcmp(p,"hull1972b4"))
63: ||(!strcmp(p,"hull1972b5")) ) PetscFunctionReturn(3);
64: else if ((!strcmp(p,"hull1972c1"))
65: ||(!strcmp(p,"hull1972c2"))
66: ||(!strcmp(p,"hull1972c3")) ) PetscFunctionReturn(10);
67: else if (!strcmp(p,"hull1972c4") ) PetscFunctionReturn(51);
68: else PetscFunctionReturn(-1);
69: }
71: /****************************************************************/
73: /* Problem specific functions */
75: /* Hull, 1972, Problem A1 */
79: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
80: {
82: PetscScalar *y,*f;
85: VecGetArray(Y,&y);
86: VecGetArray(F,&f);
87: f[0] = -y[0];
88: VecRestoreArray(Y,&y);
89: VecRestoreArray(F,&f);
90: return(0);
91: }
95: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
96: {
98: PetscScalar *y,*f;
101: VecGetArray(Y,&y);
102: VecGetArray(F,&f);
103: f[0] = -y[0];
104: VecRestoreArray(Y,&y);
105: VecRestoreArray(F,&f);
106: /* Left hand side = ydot - f(y) */
107: VecAYPX(F,-1.0,Ydot);
108: return(0);
109: }
113: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
114: {
116: PetscScalar *y;
117: PetscInt row = 0,col = 0;
118: PetscScalar value = a - 1.0;
121: VecGetArray(Y,&y);
122: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
123: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
125: VecRestoreArray(Y,&y);
126: return(0);
127: }
129: /* Hull, 1972, Problem A2 */
133: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
134: {
136: PetscScalar *y,*f;
139: VecGetArray(Y,&y);
140: VecGetArray(F,&f);
141: f[0] = -0.5*y[0]*y[0]*y[0];
142: VecRestoreArray(Y,&y);
143: VecRestoreArray(F,&f);
144: return(0);
145: }
149: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
150: {
152: PetscScalar *y,*f;
155: VecGetArray(Y,&y);
156: VecGetArray(F,&f);
157: f[0] = -0.5*y[0]*y[0]*y[0];
158: VecRestoreArray(Y,&y);
159: VecRestoreArray(F,&f);
160: /* Left hand side = ydot - f(y) */
161: VecAYPX(F,-1.0,Ydot);
162: return(0);
163: }
167: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
168: {
170: PetscScalar *y;
171: PetscInt row = 0,col = 0;
172: PetscScalar value;
175: VecGetArray(Y,&y);
176: value = a + 0.5*3.0*y[0]*y[0];
177: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
178: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
179: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
180: VecRestoreArray(Y,&y);
181: return(0);
182: }
184: /* Hull, 1972, Problem A3 */
188: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
189: {
191: PetscScalar *y,*f;
194: VecGetArray(Y,&y);
195: VecGetArray(F,&f);
196: f[0] = y[0]*PetscCosReal(t);
197: VecRestoreArray(Y,&y);
198: VecRestoreArray(F,&f);
199: return(0);
200: }
204: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
205: {
207: PetscScalar *y,*f;
210: VecGetArray(Y,&y);
211: VecGetArray(F,&f);
212: f[0] = y[0]*PetscCosReal(t);
213: VecRestoreArray(Y,&y);
214: VecRestoreArray(F,&f);
215: /* Left hand side = ydot - f(y) */
216: VecAYPX(F,-1.0,Ydot);
217: return(0);
218: }
222: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
223: {
225: PetscScalar *y;
226: PetscInt row = 0,col = 0;
227: PetscScalar value = a - PetscCosReal(t);
230: VecGetArray(Y,&y);
231: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
232: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
233: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
234: VecRestoreArray(Y,&y);
235: return(0);
236: }
238: /* Hull, 1972, Problem A4 */
242: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
243: {
245: PetscScalar *y,*f;
248: VecGetArray(Y,&y);
249: VecGetArray(F,&f);
250: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
251: VecRestoreArray(Y,&y);
252: VecRestoreArray(F,&f);
253: return(0);
254: }
258: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
259: {
261: PetscScalar *y,*f;
264: VecGetArray(Y,&y);
265: VecGetArray(F,&f);
266: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
267: VecRestoreArray(Y,&y);
268: VecRestoreArray(F,&f);
269: /* Left hand side = ydot - f(y) */
270: VecAYPX(F,-1.0,Ydot);
271: return(0);
272: }
276: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
277: {
279: PetscScalar *y;
280: PetscInt row = 0,col = 0;
281: PetscScalar value;
284: VecGetArray(Y,&y);
285: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
286: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
287: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
288: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
290: VecRestoreArray(Y,&y);
291: return(0);
292: }
294: /* Hull, 1972, Problem A5 */
298: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
299: {
301: PetscScalar *y,*f;
304: VecGetArray(Y,&y);
305: VecGetArray(F,&f);
306: f[0] = (y[0]-t)/(y[0]+t);
307: VecRestoreArray(Y,&y);
308: VecRestoreArray(F,&f);
309: return(0);
310: }
314: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
315: {
317: PetscScalar *y,*f;
320: VecGetArray(Y,&y);
321: VecGetArray(F,&f);
322: f[0] = (y[0]-t)/(y[0]+t);
323: VecRestoreArray(Y,&y);
324: VecRestoreArray(F,&f);
325: /* Left hand side = ydot - f(y) */
326: VecAYPX(F,-1.0,Ydot);
327: return(0);
328: }
332: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
333: {
335: PetscScalar *y;
336: PetscInt row = 0,col = 0;
337: PetscScalar value;
340: VecGetArray(Y,&y);
341: value = a - 2*t/((t+y[0])*(t+y[0]));
342: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
343: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
344: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
345: VecRestoreArray(Y,&y);
346: return(0);
347: }
349: /* Hull, 1972, Problem B1 */
353: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
354: {
356: PetscScalar *y,*f;
359: VecGetArray(Y,&y);
360: VecGetArray(F,&f);
361: f[0] = 2.0*(y[0] - y[0]*y[1]);
362: f[1] = -(y[1]-y[0]*y[1]);
363: VecRestoreArray(Y,&y);
364: VecRestoreArray(F,&f);
365: return(0);
366: }
370: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
371: {
373: PetscScalar *y,*f;
376: VecGetArray(Y,&y);
377: VecGetArray(F,&f);
378: f[0] = 2.0*(y[0] - y[0]*y[1]);
379: f[1] = -(y[1]-y[0]*y[1]);
380: VecRestoreArray(Y,&y);
381: VecRestoreArray(F,&f);
382: /* Left hand side = ydot - f(y) */
383: VecAYPX(F,-1.0,Ydot);
384: return(0);
385: }
389: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
390: {
392: PetscScalar *y;
393: PetscInt row[2] = {0,1};
394: PetscScalar value[2][2];
397: VecGetArray(Y,&y);
398: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
399: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
400: MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
401: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
402: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
403: VecRestoreArray(Y,&y);
404: return(0);
405: }
407: /* Hull, 1972, Problem B2 */
411: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
412: {
414: PetscScalar *y,*f;
417: VecGetArray(Y,&y);
418: VecGetArray(F,&f);
419: f[0] = -y[0] + y[1];
420: f[1] = y[0] - 2.0*y[1] + y[2];
421: f[2] = y[1] - y[2];
422: VecRestoreArray(Y,&y);
423: VecRestoreArray(F,&f);
424: return(0);
425: }
429: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
430: {
431: PetscErrorCode ierr;
432: PetscScalar *y,*f;
435: VecGetArray(Y,&y);
436: VecGetArray(F,&f);
437: f[0] = -y[0] + y[1];
438: f[1] = y[0] - 2.0*y[1] + y[2];
439: f[2] = y[1] - y[2];
440: VecRestoreArray(Y,&y);
441: VecRestoreArray(F,&f);
442: /* Left hand side = ydot - f(y) */
443: VecAYPX(F,-1.0,Ydot);
444: return(0);
445: }
449: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
450: {
452: PetscScalar *y;
453: PetscInt row[3] = {0,1,2};
454: PetscScalar value[3][3];
457: VecGetArray(Y,&y);
458: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
459: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
460: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
461: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
462: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
464: VecRestoreArray(Y,&y);
465: return(0);
466: }
468: /* Hull, 1972, Problem B3 */
472: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
473: {
475: PetscScalar *y,*f;
478: VecGetArray(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: VecRestoreArray(Y,&y);
484: VecRestoreArray(F,&f);
485: return(0);
486: }
490: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
491: {
493: PetscScalar *y,*f;
496: VecGetArray(Y,&y);
497: VecGetArray(F,&f);
498: f[0] = -y[0];
499: f[1] = y[0] - y[1]*y[1];
500: f[2] = y[1]*y[1];
501: VecRestoreArray(Y,&y);
502: VecRestoreArray(F,&f);
503: /* Left hand side = ydot - f(y) */
504: VecAYPX(F,-1.0,Ydot);
505: return(0);
506: }
510: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
511: {
513: PetscScalar *y;
514: PetscInt row[3] = {0,1,2};
515: PetscScalar value[3][3];
518: VecGetArray(Y,&y);
519: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
520: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
521: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
522: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
523: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
524: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
525: VecRestoreArray(Y,&y);
526: return(0);
527: }
529: /* Hull, 1972, Problem B4 */
533: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
534: {
536: PetscScalar *y,*f;
539: VecGetArray(Y,&y);
540: VecGetArray(F,&f);
541: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
542: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
543: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
544: VecRestoreArray(Y,&y);
545: VecRestoreArray(F,&f);
546: return(0);
547: }
551: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
552: {
553: PetscErrorCode ierr;
554: PetscScalar *y,*f;
557: VecGetArray(Y,&y);
558: VecGetArray(F,&f);
559: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
560: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
561: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
562: VecRestoreArray(Y,&y);
563: VecRestoreArray(F,&f);
564: /* Left hand side = ydot - f(y) */
565: VecAYPX(F,-1.0,Ydot);
566: return(0);
567: }
571: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
572: {
574: PetscScalar *y;
575: PetscInt row[3] = {0,1,2};
576: PetscScalar value[3][3],fac,fac2;
579: VecGetArray(Y,&y);
580: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
581: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
582: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
583: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
584: value[0][2] = y[0]*fac2;
585: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
586: value[1][1] = a + y[0]*y[0]*y[2]*fac;
587: value[1][2] = y[1]*fac2;
588: value[2][0] = -y[1]*y[1]*fac;
589: value[2][1] = y[0]*y[1]*fac;
590: value[2][2] = a;
591: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
592: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
593: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
594: VecRestoreArray(Y,&y);
595: return(0);
596: }
598: /* Hull, 1972, Problem B5 */
602: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
603: {
605: PetscScalar *y,*f;
608: VecGetArray(Y,&y);
609: VecGetArray(F,&f);
610: f[0] = y[1]*y[2];
611: f[1] = -y[0]*y[2];
612: f[2] = -0.51*y[0]*y[1];
613: VecRestoreArray(Y,&y);
614: VecRestoreArray(F,&f);
615: return(0);
616: }
620: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
621: {
623: PetscScalar *y,*f;
626: VecGetArray(Y,&y);
627: VecGetArray(F,&f);
628: f[0] = y[1]*y[2];
629: f[1] = -y[0]*y[2];
630: f[2] = -0.51*y[0]*y[1];
631: VecRestoreArray(Y,&y);
632: VecRestoreArray(F,&f);
633: /* Left hand side = ydot - f(y) */
634: VecAYPX(F,-1.0,Ydot);
635: return(0);
636: }
640: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
641: {
643: PetscScalar *y;
644: PetscInt row[3] = {0,1,2};
645: PetscScalar value[3][3];
648: VecGetArray(Y,&y);
649: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
650: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
651: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
652: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
653: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
654: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
655: VecRestoreArray(Y,&y);
656: return(0);
657: }
659: /* Hull, 1972, Problem C1 */
663: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
664: {
666: PetscScalar *y,*f;
667: PetscInt N,i;
670: VecGetSize (Y,&N);
671: VecGetArray(Y,&y);
672: VecGetArray(F,&f);
673: f[0] = -y[0];
674: for (i = 1; i < N-1; i++) {
675: f[i] = y[i-1] - y[i];
676: }
677: f[N-1] = y[N-2];
678: VecRestoreArray(Y,&y);
679: VecRestoreArray(F,&f);
680: return(0);
681: }
685: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
686: {
688: PetscScalar *y,*f;
689: PetscInt N,i;
692: VecGetSize (Y,&N);
693: VecGetArray(Y,&y);
694: VecGetArray(F,&f);
695: f[0] = -y[0];
696: for (i = 1; i < N-1; i++) {
697: f[i] = y[i-1] - y[i];
698: }
699: f[N-1] = y[N-2];
700: VecRestoreArray(Y,&y);
701: VecRestoreArray(F,&f);
702: /* Left hand side = ydot - f(y) */
703: VecAYPX(F,-1.0,Ydot);
704: return(0);
705: }
709: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
710: {
711: PetscErrorCode ierr;
712: PetscScalar *y;
713: PetscInt N,i,col[2];
714: PetscScalar value[2];
717: VecGetSize (Y,&N);
718: VecGetArray(Y,&y);
719: i = 0;
720: value[0] = a+1; col[0] = 0;
721: value[1] = 0; col[1] = 1;
722: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
723: for (i = 0; i < N; i++) {
724: value[0] = -1; col[0] = i-1;
725: value[1] = a+1; col[1] = i;
726: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
727: }
728: i = N-1;
729: value[0] = -1; col[0] = N-2;
730: value[1] = a; col[1] = N-1;
731: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
732: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
733: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
734: VecRestoreArray(Y,&y);
735: return(0);
736: }
738: /* Hull, 1972, Problem C2 */
742: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
743: {
744: PetscErrorCode ierr;
745: PetscScalar *y,*f;
746: PetscInt N,i;
749: VecGetSize (Y,&N);
750: VecGetArray(Y,&y);
751: VecGetArray(F,&f);
752: f[0] = -y[0];
753: for (i = 1; i < N-1; i++) {
754: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
755: }
756: f[N-1] = (PetscReal)(N-1)*y[N-2];
757: VecRestoreArray(Y,&y);
758: VecRestoreArray(F,&f);
759: return(0);
760: }
764: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
765: {
766: PetscErrorCode ierr;
767: PetscScalar *y,*f;
768: PetscInt N,i;
771: VecGetSize (Y,&N);
772: VecGetArray(Y,&y);
773: VecGetArray(F,&f);
774: f[0] = -y[0];
775: for (i = 1; i < N-1; i++) {
776: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
777: }
778: f[N-1] = (PetscReal)(N-1)*y[N-2];
779: VecRestoreArray(Y,&y);
780: VecRestoreArray(F,&f);
781: /* Left hand side = ydot - f(y) */
782: VecAYPX(F,-1.0,Ydot);
783: return(0);
784: }
788: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
789: {
790: PetscErrorCode ierr;
791: PetscScalar *y;
792: PetscInt N,i,col[2];
793: PetscScalar value[2];
796: VecGetSize (Y,&N);
797: VecGetArray(Y,&y);
798: i = 0;
799: value[0] = a+1; col[0] = 0;
800: value[1] = 0; col[1] = 1;
801: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
802: for (i = 0; i < N; i++) {
803: value[0] = -(PetscReal) i; col[0] = i-1;
804: value[1] = a+(PetscReal)(i+1); col[1] = i;
805: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
806: }
807: i = N-1;
808: value[0] = -(PetscReal) (N-1); col[0] = N-2;
809: value[1] = a; col[1] = N-1;
810: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
811: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
812: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
813: VecRestoreArray(Y,&y);
814: return(0);
815: }
817: /* Hull, 1972, Problem C3 and C4 */
821: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
822: {
824: PetscScalar *y,*f;
825: PetscInt N,i;
828: VecGetSize (Y,&N);
829: VecGetArray(Y,&y);
830: VecGetArray(F,&f);
831: f[0] = -2.0*y[0] + y[1];
832: for (i = 1; i < N-1; i++) {
833: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
834: }
835: f[N-1] = y[N-2] - 2.0*y[N-1];
836: VecRestoreArray(Y,&y);
837: VecRestoreArray(F,&f);
838: return(0);
839: }
843: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
844: {
846: PetscScalar *y,*f;
847: PetscInt N,i;
850: VecGetSize (Y,&N);
851: VecGetArray(Y,&y);
852: VecGetArray(F,&f);
853: f[0] = -2.0*y[0] + y[1];
854: for (i = 1; i < N-1; i++) {
855: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
856: }
857: f[N-1] = y[N-2] - 2.0*y[N-1];
858: VecRestoreArray(Y,&y);
859: VecRestoreArray(F,&f);
860: /* Left hand side = ydot - f(y) */
861: VecAYPX(F,-1.0,Ydot);
862: return(0);
863: }
867: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
868: {
870: PetscScalar *y,value[3];
871: PetscInt N,i,col[3];
874: VecGetSize (Y,&N);
875: VecGetArray(Y,&y);
876: for (i = 0; i < N; i++) {
877: if (i == 0) {
878: value[0] = a+2; col[0] = i;
879: value[1] = -1; col[1] = i+1;
880: value[2] = 0; col[2] = i+2;
881: } else if (i == N-1) {
882: value[0] = 0; col[0] = i-2;
883: value[1] = -1; col[1] = i-1;
884: value[2] = a+2; col[2] = i;
885: } else {
886: value[0] = -1; col[0] = i-1;
887: value[1] = a+2; col[1] = i;
888: value[2] = -1; col[2] = i+1;
889: }
890: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
891: }
892: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
893: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
894: VecRestoreArray(Y,&y);
895: return(0);
896: }
898: /***************************************************************************/
902: /* Sets the initial solution for the IVP and sets up the function pointers*/
903: PetscErrorCode Initialize(Vec Y, void* s)
904: {
906: char *p = (char*) s;
907: PetscScalar *y;
908: PetscInt N = GetSize((const char *)s);
909: PetscBool flg;
912: VecZeroEntries(Y);
913: VecGetArray(Y,&y);
914: if (!strcmp(p,"hull1972a1")) {
915: y[0] = 1.0;
916: RHSFunction = RHSFunction_Hull1972A1;
917: IFunction = IFunction_Hull1972A1;
918: IJacobian = IJacobian_Hull1972A1;
919: } else if (!strcmp(p,"hull1972a2")) {
920: y[0] = 1.0;
921: RHSFunction = RHSFunction_Hull1972A2;
922: IFunction = IFunction_Hull1972A2;
923: IJacobian = IJacobian_Hull1972A2;
924: } else if (!strcmp(p,"hull1972a3")) {
925: y[0] = 1.0;
926: RHSFunction = RHSFunction_Hull1972A3;
927: IFunction = IFunction_Hull1972A3;
928: IJacobian = IJacobian_Hull1972A3;
929: } else if (!strcmp(p,"hull1972a4")) {
930: y[0] = 1.0;
931: RHSFunction = RHSFunction_Hull1972A4;
932: IFunction = IFunction_Hull1972A4;
933: IJacobian = IJacobian_Hull1972A4;
934: } else if (!strcmp(p,"hull1972a5")) {
935: y[0] = 4.0;
936: RHSFunction = RHSFunction_Hull1972A5;
937: IFunction = IFunction_Hull1972A5;
938: IJacobian = IJacobian_Hull1972A5;
939: } else if (!strcmp(p,"hull1972b1")) {
940: y[0] = 1.0;
941: y[1] = 3.0;
942: RHSFunction = RHSFunction_Hull1972B1;
943: IFunction = IFunction_Hull1972B1;
944: IJacobian = IJacobian_Hull1972B1;
945: } else if (!strcmp(p,"hull1972b2")) {
946: y[0] = 2.0;
947: y[1] = 0.0;
948: y[2] = 1.0;
949: RHSFunction = RHSFunction_Hull1972B2;
950: IFunction = IFunction_Hull1972B2;
951: IJacobian = IJacobian_Hull1972B2;
952: } else if (!strcmp(p,"hull1972b3")) {
953: y[0] = 1.0;
954: y[1] = 0.0;
955: y[2] = 0.0;
956: RHSFunction = RHSFunction_Hull1972B3;
957: IFunction = IFunction_Hull1972B3;
958: IJacobian = IJacobian_Hull1972B3;
959: } else if (!strcmp(p,"hull1972b4")) {
960: y[0] = 3.0;
961: y[1] = 0.0;
962: y[2] = 0.0;
963: RHSFunction = RHSFunction_Hull1972B4;
964: IFunction = IFunction_Hull1972B4;
965: IJacobian = IJacobian_Hull1972B4;
966: } else if (!strcmp(p,"hull1972b5")) {
967: y[0] = 0.0;
968: y[1] = 1.0;
969: y[2] = 1.0;
970: RHSFunction = RHSFunction_Hull1972B5;
971: IFunction = IFunction_Hull1972B5;
972: IJacobian = IJacobian_Hull1972B5;
973: } else if (!strcmp(p,"hull1972c1")) {
974: y[0] = 1.0;
975: RHSFunction = RHSFunction_Hull1972C1;
976: IFunction = IFunction_Hull1972C1;
977: IJacobian = IJacobian_Hull1972C1;
978: } else if (!strcmp(p,"hull1972c2")) {
979: y[0] = 1.0;
980: RHSFunction = RHSFunction_Hull1972C2;
981: IFunction = IFunction_Hull1972C2;
982: IJacobian = IJacobian_Hull1972C2;
983: } else if ((!strcmp(p,"hull1972c3"))
984: ||(!strcmp(p,"hull1972c4"))){
985: y[0] = 1.0;
986: RHSFunction = RHSFunction_Hull1972C34;
987: IFunction = IFunction_Hull1972C34;
988: IJacobian = IJacobian_Hull1972C34;
989: }
990: PetscOptionsGetScalarArray(NULL,"-yinit",y,&N,&flg);
991: if ((N != GetSize((const char*)s)) && flg) {
992: 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));
993: }
994: VecRestoreArray(Y,&y);
995: return(0);
996: }
1000: /* Calculates the exact solution to problems that have one */
1001: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1002: {
1004: char *p = (char*) s;
1005: PetscScalar *y;
1008: if (!strcmp(p,"hull1972a1")) {
1009: VecGetArray(Y,&y);
1010: y[0] = PetscExpReal(-t);
1011: *flag = PETSC_TRUE;
1012: VecRestoreArray(Y,&y);
1013: } else if (!strcmp(p,"hull1972a2")) {
1014: VecGetArray(Y,&y);
1015: y[0] = 1.0/PetscSqrtReal(t+1);
1016: *flag = PETSC_TRUE;
1017: VecRestoreArray(Y,&y);
1018: } else if (!strcmp(p,"hull1972a3")) {
1019: VecGetArray(Y,&y);
1020: y[0] = PetscExpReal(PetscSinReal(t));
1021: *flag = PETSC_TRUE;
1022: VecRestoreArray(Y,&y);
1023: } else if (!strcmp(p,"hull1972a4")) {
1024: VecGetArray(Y,&y);
1025: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1026: *flag = PETSC_TRUE;
1027: VecRestoreArray(Y,&y);
1028: } else {
1029: VecSet(Y,0);
1030: *flag = PETSC_FALSE;
1031: }
1032: return(0);
1033: }
1037: /* Solves the specified ODE and computes the error if exact solution is available */
1038: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1039: {
1040: PetscErrorCode ierr; /* Error code */
1041: TS ts; /* time-integrator */
1042: Vec Y; /* Solution vector */
1043: Vec Yex; /* Exact solution */
1044: PetscInt N; /* Size of the system of equations */
1045: TSType time_scheme; /* Type of time-integration scheme */
1046: Mat Jac = NULL; /* Jacobian matrix */
1049: N = GetSize((const char *)&ptype[0]);
1050: if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1051: VecCreate(PETSC_COMM_WORLD,&Y);
1052: VecSetSizes(Y,N,PETSC_DECIDE);
1053: VecSetUp(Y);
1054: VecSet(Y,0);
1056: /* Initialize the problem */
1057: Initialize(Y,&ptype[0]);
1059: /* Create and initialize the time-integrator */
1060: TSCreate(PETSC_COMM_WORLD,&ts);
1061: /* Default time integration options */
1062: TSSetType(ts,TSEULER);
1063: TSSetDuration(ts,maxiter,tfinal);
1064: TSSetInitialTimeStep(ts,0.0,dt);
1065: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1066: /* Read command line options for time integration */
1067: TSSetFromOptions(ts);
1068: /* Set solution vector */
1069: TSSetSolution(ts,Y);
1070: /* Specify left/right-hand side functions */
1071: TSGetType(ts,&time_scheme);
1072: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP))) {
1073: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1074: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1075: } else if ((!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSARKIMEX))) {
1076: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1077: /* and its Jacobian function */
1078: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1079: MatCreate(PETSC_COMM_WORLD,&Jac);
1080: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1081: MatSetFromOptions(Jac);
1082: MatSetUp(Jac);
1083: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1084: }
1086: /* Solve */
1087: TSSolve(ts,Y);
1089: /* Exact solution */
1090: VecDuplicate(Y,&Yex);
1091: ExactSolution(Yex,&ptype[0],tfinal,exact_flag);
1093: /* Calculate Error */
1094: VecAYPX(Yex,-1.0,Y);
1095: VecNorm(Yex,NORM_2,error);
1096: *error = PetscSqrtReal(((*error)*(*error))/N);
1098: /* Clean up and finalize */
1099: MatDestroy(&Jac);
1100: TSDestroy(&ts);
1101: VecDestroy(&Yex);
1102: VecDestroy(&Y);
1104: return(0);
1105: }
1109: int main(int argc, char **argv)
1110: {
1111: PetscErrorCode ierr; /* Error code */
1112: char ptype[256] = "hull1972a1"; /* Problem specification */
1113: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1114: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1115: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1116: PetscReal dt;
1117: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1118: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1119: PetscReal *error; /* Array to store the errors for convergence analysis */
1120: PetscMPIInt size; /* No of processors */
1121: PetscBool flag; /* Flag denoting availability of exact solution */
1122: PetscInt r;
1124: /* Initialize program */
1125: PetscInitialize(&argc,&argv,(char*)0,help);
1127: /* Check if running with only 1 proc */
1128: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1129: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1131: PetscOptionsString("-problem","Problem specification","<hull1972a1>",
1132: ptype,ptype,sizeof(ptype),NULL);
1133: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis",
1134: "<1>",n_refine,&n_refine,NULL);
1135: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",
1136: refine_fac,&refine_fac,NULL);
1137: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)",
1138: "<0.01>",dt_initial,&dt_initial,NULL);
1139: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",
1140: tfinal,&tfinal,NULL);
1142: PetscMalloc1(n_refine,&error);
1143: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1144: error[r] = 0;
1145: if (r > 0) dt /= refine_fac;
1147: 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]));
1148: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1149: if (flag) {
1150: /* If exact solution available for the specified ODE */
1151: if (r > 0) {
1152: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1153: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f\n.",(double)error[r],(double)conv_rate);
1154: } else {
1155: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1156: }
1157: }
1158: }
1159: PetscFree(error);
1161: /* Exit */
1162: PetscFinalize();
1163: return(0);
1164: }