Actual source code: ex31.c
petsc-3.7.5 2017-01-01
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: {
81: PetscErrorCode ierr;
82: PetscScalar *f;
83: const PetscScalar *y;
86: VecGetArrayRead(Y,&y);
87: VecGetArray(F,&f);
88: f[0] = -y[0];
89: VecRestoreArrayRead(Y,&y);
90: VecRestoreArray(F,&f);
91: return(0);
92: }
96: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
97: {
98: PetscErrorCode ierr;
99: const PetscScalar *y;
100: PetscScalar *f;
103: VecGetArrayRead(Y,&y);
104: VecGetArray(F,&f);
105: f[0] = -y[0];
106: VecRestoreArrayRead(Y,&y);
107: VecRestoreArray(F,&f);
108: /* Left hand side = ydot - f(y) */
109: VecAYPX(F,-1.0,Ydot);
110: return(0);
111: }
115: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
116: {
117: PetscErrorCode ierr;
118: const PetscScalar *y;
119: PetscInt row = 0,col = 0;
120: PetscScalar value = a - 1.0;
123: VecGetArrayRead(Y,&y);
124: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
125: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
127: VecRestoreArrayRead(Y,&y);
128: return(0);
129: }
131: /* Hull, 1972, Problem A2 */
135: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
136: {
137: PetscErrorCode ierr;
138: const PetscScalar *y;
139: PetscScalar *f;
142: VecGetArrayRead(Y,&y);
143: VecGetArray(F,&f);
144: f[0] = -0.5*y[0]*y[0]*y[0];
145: VecRestoreArrayRead(Y,&y);
146: VecRestoreArray(F,&f);
147: return(0);
148: }
152: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
153: {
154: PetscErrorCode ierr;
155: PetscScalar *f;
156: const PetscScalar *y;
159: VecGetArrayRead(Y,&y);
160: VecGetArray(F,&f);
161: f[0] = -0.5*y[0]*y[0]*y[0];
162: VecRestoreArrayRead(Y,&y);
163: VecRestoreArray(F,&f);
164: /* Left hand side = ydot - f(y) */
165: VecAYPX(F,-1.0,Ydot);
166: return(0);
167: }
171: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
172: {
173: PetscErrorCode ierr;
174: const PetscScalar *y;
175: PetscInt row = 0,col = 0;
176: PetscScalar value;
179: VecGetArrayRead(Y,&y);
180: value = a + 0.5*3.0*y[0]*y[0];
181: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
182: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
183: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
184: VecRestoreArrayRead(Y,&y);
185: return(0);
186: }
188: /* Hull, 1972, Problem A3 */
192: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
193: {
194: PetscErrorCode ierr;
195: const PetscScalar *y;
196: PetscScalar *f;
199: VecGetArrayRead(Y,&y);
200: VecGetArray(F,&f);
201: f[0] = y[0]*PetscCosReal(t);
202: VecRestoreArrayRead(Y,&y);
203: VecRestoreArray(F,&f);
204: return(0);
205: }
209: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
210: {
211: PetscErrorCode ierr;
212: PetscScalar *f;
213: const PetscScalar *y;
216: VecGetArrayRead(Y,&y);
217: VecGetArray(F,&f);
218: f[0] = y[0]*PetscCosReal(t);
219: VecRestoreArrayRead(Y,&y);
220: VecRestoreArray(F,&f);
221: /* Left hand side = ydot - f(y) */
222: VecAYPX(F,-1.0,Ydot);
223: return(0);
224: }
228: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
229: {
230: PetscErrorCode ierr;
231: const PetscScalar *y;
232: PetscInt row = 0,col = 0;
233: PetscScalar value = a - PetscCosReal(t);
236: VecGetArrayRead(Y,&y);
237: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
238: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
239: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
240: VecRestoreArrayRead(Y,&y);
241: return(0);
242: }
244: /* Hull, 1972, Problem A4 */
248: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, 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] = (0.25*y[0])*(1.0-0.05*y[0]);
258: VecRestoreArrayRead(Y,&y);
259: VecRestoreArray(F,&f);
260: return(0);
261: }
265: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
266: {
267: PetscErrorCode ierr;
268: PetscScalar *f;
269: const PetscScalar *y;
272: VecGetArrayRead(Y,&y);
273: VecGetArray(F,&f);
274: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
275: VecRestoreArrayRead(Y,&y);
276: VecRestoreArray(F,&f);
277: /* Left hand side = ydot - f(y) */
278: VecAYPX(F,-1.0,Ydot);
279: return(0);
280: }
284: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
285: {
286: PetscErrorCode ierr;
287: const PetscScalar *y;
288: PetscInt row = 0,col = 0;
289: PetscScalar value;
292: VecGetArrayRead(Y,&y);
293: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
294: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
295: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
296: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
297: VecRestoreArrayRead(Y,&y);
298: return(0);
299: }
301: /* Hull, 1972, Problem A5 */
305: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
306: {
307: PetscErrorCode ierr;
308: PetscScalar *f;
309: const PetscScalar *y;
312: VecGetArrayRead(Y,&y);
313: VecGetArray(F,&f);
314: f[0] = (y[0]-t)/(y[0]+t);
315: VecRestoreArrayRead(Y,&y);
316: VecRestoreArray(F,&f);
317: return(0);
318: }
322: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
323: {
324: PetscErrorCode ierr;
325: PetscScalar *f;
326: const PetscScalar *y;
329: VecGetArrayRead(Y,&y);
330: VecGetArray(F,&f);
331: f[0] = (y[0]-t)/(y[0]+t);
332: VecRestoreArrayRead(Y,&y);
333: VecRestoreArray(F,&f);
334: /* Left hand side = ydot - f(y) */
335: VecAYPX(F,-1.0,Ydot);
336: return(0);
337: }
341: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
342: {
343: PetscErrorCode ierr;
344: const PetscScalar *y;
345: PetscInt row = 0,col = 0;
346: PetscScalar value;
349: VecGetArrayRead(Y,&y);
350: value = a - 2*t/((t+y[0])*(t+y[0]));
351: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
352: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
353: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
354: VecRestoreArrayRead(Y,&y);
355: return(0);
356: }
358: /* Hull, 1972, Problem B1 */
362: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
363: {
364: PetscErrorCode ierr;
365: PetscScalar *f;
366: const PetscScalar *y;
369: VecGetArrayRead(Y,&y);
370: VecGetArray(F,&f);
371: f[0] = 2.0*(y[0] - y[0]*y[1]);
372: f[1] = -(y[1]-y[0]*y[1]);
373: VecRestoreArrayRead(Y,&y);
374: VecRestoreArray(F,&f);
375: return(0);
376: }
380: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
381: {
382: PetscErrorCode ierr;
383: PetscScalar *f;
384: const PetscScalar *y;
387: VecGetArrayRead(Y,&y);
388: VecGetArray(F,&f);
389: f[0] = 2.0*(y[0] - y[0]*y[1]);
390: f[1] = -(y[1]-y[0]*y[1]);
391: VecRestoreArrayRead(Y,&y);
392: VecRestoreArray(F,&f);
393: /* Left hand side = ydot - f(y) */
394: VecAYPX(F,-1.0,Ydot);
395: return(0);
396: }
400: PetscErrorCode IJacobian_Hull1972B1(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[2] = {0,1};
405: PetscScalar value[2][2];
408: VecGetArrayRead(Y,&y);
409: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
410: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
411: MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
412: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
413: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
414: VecRestoreArrayRead(Y,&y);
415: return(0);
416: }
418: /* Hull, 1972, Problem B2 */
422: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
423: {
424: PetscErrorCode ierr;
425: PetscScalar *f;
426: const PetscScalar *y;
429: VecGetArrayRead(Y,&y);
430: VecGetArray(F,&f);
431: f[0] = -y[0] + y[1];
432: f[1] = y[0] - 2.0*y[1] + y[2];
433: f[2] = y[1] - y[2];
434: VecRestoreArrayRead(Y,&y);
435: VecRestoreArray(F,&f);
436: return(0);
437: }
441: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
442: {
443: PetscErrorCode ierr;
444: PetscScalar *f;
445: const PetscScalar *y;
448: VecGetArrayRead(Y,&y);
449: VecGetArray(F,&f);
450: f[0] = -y[0] + y[1];
451: f[1] = y[0] - 2.0*y[1] + y[2];
452: f[2] = y[1] - y[2];
453: VecRestoreArrayRead(Y,&y);
454: VecRestoreArray(F,&f);
455: /* Left hand side = ydot - f(y) */
456: VecAYPX(F,-1.0,Ydot);
457: return(0);
458: }
462: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
463: {
464: PetscErrorCode ierr;
465: const PetscScalar *y;
466: PetscInt row[3] = {0,1,2};
467: PetscScalar value[3][3];
470: VecGetArrayRead(Y,&y);
471: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
472: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
473: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
474: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
475: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
476: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
477: VecRestoreArrayRead(Y,&y);
478: return(0);
479: }
481: /* Hull, 1972, Problem B3 */
485: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
486: {
487: PetscErrorCode ierr;
488: PetscScalar *f;
489: const PetscScalar *y;
492: VecGetArrayRead(Y,&y);
493: VecGetArray(F,&f);
494: f[0] = -y[0];
495: f[1] = y[0] - y[1]*y[1];
496: f[2] = y[1]*y[1];
497: VecRestoreArrayRead(Y,&y);
498: VecRestoreArray(F,&f);
499: return(0);
500: }
504: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
505: {
506: PetscErrorCode ierr;
507: PetscScalar *f;
508: const PetscScalar *y;
511: VecGetArrayRead(Y,&y);
512: VecGetArray(F,&f);
513: f[0] = -y[0];
514: f[1] = y[0] - y[1]*y[1];
515: f[2] = y[1]*y[1];
516: VecRestoreArrayRead(Y,&y);
517: VecRestoreArray(F,&f);
518: /* Left hand side = ydot - f(y) */
519: VecAYPX(F,-1.0,Ydot);
520: return(0);
521: }
525: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
526: {
527: PetscErrorCode ierr;
528: const PetscScalar *y;
529: PetscInt row[3] = {0,1,2};
530: PetscScalar value[3][3];
533: VecGetArrayRead(Y,&y);
534: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
535: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
536: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
537: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
538: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
539: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
540: VecRestoreArrayRead(Y,&y);
541: return(0);
542: }
544: /* Hull, 1972, Problem B4 */
548: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
549: {
550: PetscErrorCode ierr;
551: PetscScalar *f;
552: const PetscScalar *y;
555: VecGetArrayRead(Y,&y);
556: VecGetArray(F,&f);
557: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
558: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
559: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
560: VecRestoreArrayRead(Y,&y);
561: VecRestoreArray(F,&f);
562: return(0);
563: }
567: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
568: {
569: PetscErrorCode ierr;
570: PetscScalar *f;
571: const PetscScalar *y;
574: VecGetArrayRead(Y,&y);
575: VecGetArray(F,&f);
576: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
577: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
578: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
579: VecRestoreArrayRead(Y,&y);
580: VecRestoreArray(F,&f);
581: /* Left hand side = ydot - f(y) */
582: VecAYPX(F,-1.0,Ydot);
583: return(0);
584: }
588: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
589: {
590: PetscErrorCode ierr;
591: const PetscScalar *y;
592: PetscInt row[3] = {0,1,2};
593: PetscScalar value[3][3],fac,fac2;
596: VecGetArrayRead(Y,&y);
597: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
598: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
599: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
600: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
601: value[0][2] = y[0]*fac2;
602: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
603: value[1][1] = a + y[0]*y[0]*y[2]*fac;
604: value[1][2] = y[1]*fac2;
605: value[2][0] = -y[1]*y[1]*fac;
606: value[2][1] = y[0]*y[1]*fac;
607: value[2][2] = a;
608: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
609: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
610: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
611: VecRestoreArrayRead(Y,&y);
612: return(0);
613: }
615: /* Hull, 1972, Problem B5 */
619: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
620: {
621: PetscErrorCode ierr;
622: PetscScalar *f;
623: const PetscScalar *y;
626: VecGetArrayRead(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: VecRestoreArrayRead(Y,&y);
632: VecRestoreArray(F,&f);
633: return(0);
634: }
638: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
639: {
640: PetscErrorCode ierr;
641: PetscScalar *f;
642: const PetscScalar *y;
645: VecGetArrayRead(Y,&y);
646: VecGetArray(F,&f);
647: f[0] = y[1]*y[2];
648: f[1] = -y[0]*y[2];
649: f[2] = -0.51*y[0]*y[1];
650: VecRestoreArrayRead(Y,&y);
651: VecRestoreArray(F,&f);
652: /* Left hand side = ydot - f(y) */
653: VecAYPX(F,-1.0,Ydot);
654: return(0);
655: }
659: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
660: {
661: PetscErrorCode ierr;
662: const PetscScalar *y;
663: PetscInt row[3] = {0,1,2};
664: PetscScalar value[3][3];
667: VecGetArrayRead(Y,&y);
668: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
669: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
670: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
671: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
672: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
673: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
674: VecRestoreArrayRead(Y,&y);
675: return(0);
676: }
678: /* Hull, 1972, Problem C1 */
682: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
683: {
684: PetscErrorCode ierr;
685: PetscScalar *f;
686: const PetscScalar *y;
687: PetscInt N,i;
690: VecGetSize (Y,&N);
691: VecGetArrayRead(Y,&y);
692: VecGetArray(F,&f);
693: f[0] = -y[0];
694: for (i = 1; i < N-1; i++) {
695: f[i] = y[i-1] - y[i];
696: }
697: f[N-1] = y[N-2];
698: VecRestoreArrayRead(Y,&y);
699: VecRestoreArray(F,&f);
700: return(0);
701: }
705: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
706: {
707: PetscErrorCode ierr;
708: PetscScalar *f;
709: const PetscScalar *y;
710: PetscInt N,i;
713: VecGetSize (Y,&N);
714: VecGetArrayRead(Y,&y);
715: VecGetArray(F,&f);
716: f[0] = -y[0];
717: for (i = 1; i < N-1; i++) {
718: f[i] = y[i-1] - y[i];
719: }
720: f[N-1] = y[N-2];
721: VecRestoreArrayRead(Y,&y);
722: VecRestoreArray(F,&f);
723: /* Left hand side = ydot - f(y) */
724: VecAYPX(F,-1.0,Ydot);
725: return(0);
726: }
730: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
731: {
732: PetscErrorCode ierr;
733: const PetscScalar *y;
734: PetscInt N,i,col[2];
735: PetscScalar value[2];
738: VecGetSize (Y,&N);
739: VecGetArrayRead(Y,&y);
740: i = 0;
741: value[0] = a+1; col[0] = 0;
742: value[1] = 0; col[1] = 1;
743: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
744: for (i = 0; i < N; i++) {
745: value[0] = -1; col[0] = i-1;
746: value[1] = a+1; col[1] = i;
747: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
748: }
749: i = N-1;
750: value[0] = -1; col[0] = N-2;
751: value[1] = a; col[1] = N-1;
752: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
753: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
754: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
755: VecRestoreArrayRead(Y,&y);
756: return(0);
757: }
759: /* Hull, 1972, Problem C2 */
763: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
764: {
765: PetscErrorCode ierr;
766: const PetscScalar *y;
767: PetscScalar *f;
768: PetscInt N,i;
771: VecGetSize (Y,&N);
772: VecGetArrayRead(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: VecRestoreArrayRead(Y,&y);
780: VecRestoreArray(F,&f);
781: return(0);
782: }
786: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
787: {
788: PetscErrorCode ierr;
789: PetscScalar *f;
790: const PetscScalar *y;
791: PetscInt N,i;
794: VecGetSize (Y,&N);
795: VecGetArrayRead(Y,&y);
796: VecGetArray(F,&f);
797: f[0] = -y[0];
798: for (i = 1; i < N-1; i++) {
799: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
800: }
801: f[N-1] = (PetscReal)(N-1)*y[N-2];
802: VecRestoreArrayRead(Y,&y);
803: VecRestoreArray(F,&f);
804: /* Left hand side = ydot - f(y) */
805: VecAYPX(F,-1.0,Ydot);
806: return(0);
807: }
811: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
812: {
813: PetscErrorCode ierr;
814: const PetscScalar *y;
815: PetscInt N,i,col[2];
816: PetscScalar value[2];
819: VecGetSize (Y,&N);
820: VecGetArrayRead(Y,&y);
821: i = 0;
822: value[0] = a+1; col[0] = 0;
823: value[1] = 0; col[1] = 1;
824: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
825: for (i = 0; i < N; i++) {
826: value[0] = -(PetscReal) i; col[0] = i-1;
827: value[1] = a+(PetscReal)(i+1); col[1] = i;
828: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
829: }
830: i = N-1;
831: value[0] = -(PetscReal) (N-1); col[0] = N-2;
832: value[1] = a; col[1] = N-1;
833: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
834: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
835: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
836: VecRestoreArrayRead(Y,&y);
837: return(0);
838: }
840: /* Hull, 1972, Problem C3 and C4 */
844: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
845: {
846: PetscErrorCode ierr;
847: PetscScalar *f;
848: const PetscScalar *y;
849: PetscInt N,i;
852: VecGetSize (Y,&N);
853: VecGetArrayRead(Y,&y);
854: VecGetArray(F,&f);
855: f[0] = -2.0*y[0] + y[1];
856: for (i = 1; i < N-1; i++) {
857: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
858: }
859: f[N-1] = y[N-2] - 2.0*y[N-1];
860: VecRestoreArrayRead(Y,&y);
861: VecRestoreArray(F,&f);
862: return(0);
863: }
867: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
868: {
869: PetscErrorCode ierr;
870: PetscScalar *f;
871: const PetscScalar *y;
872: PetscInt N,i;
875: VecGetSize (Y,&N);
876: VecGetArrayRead(Y,&y);
877: VecGetArray(F,&f);
878: f[0] = -2.0*y[0] + y[1];
879: for (i = 1; i < N-1; i++) {
880: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
881: }
882: f[N-1] = y[N-2] - 2.0*y[N-1];
883: VecRestoreArrayRead(Y,&y);
884: VecRestoreArray(F,&f);
885: /* Left hand side = ydot - f(y) */
886: VecAYPX(F,-1.0,Ydot);
887: return(0);
888: }
892: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
893: {
894: PetscErrorCode ierr;
895: const PetscScalar *y;
896: PetscScalar value[3];
897: PetscInt N,i,col[3];
900: VecGetSize (Y,&N);
901: VecGetArrayRead(Y,&y);
902: for (i = 0; i < N; i++) {
903: if (i == 0) {
904: value[0] = a+2; col[0] = i;
905: value[1] = -1; col[1] = i+1;
906: value[2] = 0; col[2] = i+2;
907: } else if (i == N-1) {
908: value[0] = 0; col[0] = i-2;
909: value[1] = -1; col[1] = i-1;
910: value[2] = a+2; col[2] = i;
911: } else {
912: value[0] = -1; col[0] = i-1;
913: value[1] = a+2; col[1] = i;
914: value[2] = -1; col[2] = i+1;
915: }
916: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
917: }
918: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
919: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
920: VecRestoreArrayRead(Y,&y);
921: return(0);
922: }
924: /***************************************************************************/
928: /* Sets the initial solution for the IVP and sets up the function pointers*/
929: PetscErrorCode Initialize(Vec Y, void* s)
930: {
932: char *p = (char*) s;
933: PetscScalar *y;
934: PetscInt N = GetSize((const char *)s);
935: PetscBool flg;
938: VecZeroEntries(Y);
939: VecGetArray(Y,&y);
940: if (!strcmp(p,"hull1972a1")) {
941: y[0] = 1.0;
942: RHSFunction = RHSFunction_Hull1972A1;
943: IFunction = IFunction_Hull1972A1;
944: IJacobian = IJacobian_Hull1972A1;
945: } else if (!strcmp(p,"hull1972a2")) {
946: y[0] = 1.0;
947: RHSFunction = RHSFunction_Hull1972A2;
948: IFunction = IFunction_Hull1972A2;
949: IJacobian = IJacobian_Hull1972A2;
950: } else if (!strcmp(p,"hull1972a3")) {
951: y[0] = 1.0;
952: RHSFunction = RHSFunction_Hull1972A3;
953: IFunction = IFunction_Hull1972A3;
954: IJacobian = IJacobian_Hull1972A3;
955: } else if (!strcmp(p,"hull1972a4")) {
956: y[0] = 1.0;
957: RHSFunction = RHSFunction_Hull1972A4;
958: IFunction = IFunction_Hull1972A4;
959: IJacobian = IJacobian_Hull1972A4;
960: } else if (!strcmp(p,"hull1972a5")) {
961: y[0] = 4.0;
962: RHSFunction = RHSFunction_Hull1972A5;
963: IFunction = IFunction_Hull1972A5;
964: IJacobian = IJacobian_Hull1972A5;
965: } else if (!strcmp(p,"hull1972b1")) {
966: y[0] = 1.0;
967: y[1] = 3.0;
968: RHSFunction = RHSFunction_Hull1972B1;
969: IFunction = IFunction_Hull1972B1;
970: IJacobian = IJacobian_Hull1972B1;
971: } else if (!strcmp(p,"hull1972b2")) {
972: y[0] = 2.0;
973: y[1] = 0.0;
974: y[2] = 1.0;
975: RHSFunction = RHSFunction_Hull1972B2;
976: IFunction = IFunction_Hull1972B2;
977: IJacobian = IJacobian_Hull1972B2;
978: } else if (!strcmp(p,"hull1972b3")) {
979: y[0] = 1.0;
980: y[1] = 0.0;
981: y[2] = 0.0;
982: RHSFunction = RHSFunction_Hull1972B3;
983: IFunction = IFunction_Hull1972B3;
984: IJacobian = IJacobian_Hull1972B3;
985: } else if (!strcmp(p,"hull1972b4")) {
986: y[0] = 3.0;
987: y[1] = 0.0;
988: y[2] = 0.0;
989: RHSFunction = RHSFunction_Hull1972B4;
990: IFunction = IFunction_Hull1972B4;
991: IJacobian = IJacobian_Hull1972B4;
992: } else if (!strcmp(p,"hull1972b5")) {
993: y[0] = 0.0;
994: y[1] = 1.0;
995: y[2] = 1.0;
996: RHSFunction = RHSFunction_Hull1972B5;
997: IFunction = IFunction_Hull1972B5;
998: IJacobian = IJacobian_Hull1972B5;
999: } else if (!strcmp(p,"hull1972c1")) {
1000: y[0] = 1.0;
1001: RHSFunction = RHSFunction_Hull1972C1;
1002: IFunction = IFunction_Hull1972C1;
1003: IJacobian = IJacobian_Hull1972C1;
1004: } else if (!strcmp(p,"hull1972c2")) {
1005: y[0] = 1.0;
1006: RHSFunction = RHSFunction_Hull1972C2;
1007: IFunction = IFunction_Hull1972C2;
1008: IJacobian = IJacobian_Hull1972C2;
1009: } else if ((!strcmp(p,"hull1972c3"))
1010: ||(!strcmp(p,"hull1972c4"))){
1011: y[0] = 1.0;
1012: RHSFunction = RHSFunction_Hull1972C34;
1013: IFunction = IFunction_Hull1972C34;
1014: IJacobian = IJacobian_Hull1972C34;
1015: }
1016: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1017: 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));
1018: VecRestoreArray(Y,&y);
1019: return(0);
1020: }
1024: /* Calculates the exact solution to problems that have one */
1025: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1026: {
1028: char *p = (char*) s;
1029: PetscScalar *y;
1032: if (!strcmp(p,"hull1972a1")) {
1033: VecGetArray(Y,&y);
1034: y[0] = PetscExpReal(-t);
1035: *flag = PETSC_TRUE;
1036: VecRestoreArray(Y,&y);
1037: } else if (!strcmp(p,"hull1972a2")) {
1038: VecGetArray(Y,&y);
1039: y[0] = 1.0/PetscSqrtReal(t+1);
1040: *flag = PETSC_TRUE;
1041: VecRestoreArray(Y,&y);
1042: } else if (!strcmp(p,"hull1972a3")) {
1043: VecGetArray(Y,&y);
1044: y[0] = PetscExpReal(PetscSinReal(t));
1045: *flag = PETSC_TRUE;
1046: VecRestoreArray(Y,&y);
1047: } else if (!strcmp(p,"hull1972a4")) {
1048: VecGetArray(Y,&y);
1049: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1050: *flag = PETSC_TRUE;
1051: VecRestoreArray(Y,&y);
1052: } else {
1053: VecSet(Y,0);
1054: *flag = PETSC_FALSE;
1055: }
1056: return(0);
1057: }
1061: /* Solves the specified ODE and computes the error if exact solution is available */
1062: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1063: {
1064: PetscErrorCode ierr; /* Error code */
1065: TS ts; /* time-integrator */
1066: Vec Y; /* Solution vector */
1067: Vec Yex; /* Exact solution */
1068: PetscInt N; /* Size of the system of equations */
1069: TSType time_scheme; /* Type of time-integration scheme */
1070: Mat Jac = NULL; /* Jacobian matrix */
1073: N = GetSize((const char *)&ptype[0]);
1074: if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1075: VecCreate(PETSC_COMM_WORLD,&Y);
1076: VecSetSizes(Y,N,PETSC_DECIDE);
1077: VecSetUp(Y);
1078: VecSet(Y,0);
1080: /* Initialize the problem */
1081: Initialize(Y,&ptype[0]);
1083: /* Create and initialize the time-integrator */
1084: TSCreate(PETSC_COMM_WORLD,&ts);
1085: /* Default time integration options */
1086: TSSetType(ts,TSRK);
1087: TSSetDuration(ts,maxiter,tfinal);
1088: TSSetInitialTimeStep(ts,0.0,dt);
1089: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1090: /* Read command line options for time integration */
1091: TSSetFromOptions(ts);
1092: /* Set solution vector */
1093: TSSetSolution(ts,Y);
1094: /* Specify left/right-hand side functions */
1095: TSGetType(ts,&time_scheme);
1096: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP))) {
1097: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1098: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1099: } else if ((!strcmp(time_scheme,TSTHETA)) ||
1100: (!strcmp(time_scheme,TSBEULER)) ||
1101: (!strcmp(time_scheme,TSCN)) ||
1102: (!strcmp(time_scheme,TSALPHA)) ||
1103: (!strcmp(time_scheme,TSARKIMEX))) {
1104: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1105: /* and its Jacobian function */
1106: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1107: MatCreate(PETSC_COMM_WORLD,&Jac);
1108: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1109: MatSetFromOptions(Jac);
1110: MatSetUp(Jac);
1111: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1112: }
1114: /* Solve */
1115: TSSolve(ts,Y);
1117: /* Exact solution */
1118: VecDuplicate(Y,&Yex);
1119: ExactSolution(Yex,&ptype[0],tfinal,exact_flag);
1121: /* Calculate Error */
1122: VecAYPX(Yex,-1.0,Y);
1123: VecNorm(Yex,NORM_2,error);
1124: *error = PetscSqrtReal(((*error)*(*error))/N);
1126: /* Clean up and finalize */
1127: MatDestroy(&Jac);
1128: TSDestroy(&ts);
1129: VecDestroy(&Yex);
1130: VecDestroy(&Y);
1132: return(0);
1133: }
1137: int main(int argc, char **argv)
1138: {
1139: PetscErrorCode ierr; /* Error code */
1140: char ptype[256] = "hull1972a1"; /* Problem specification */
1141: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1142: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1143: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1144: PetscReal dt;
1145: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1146: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1147: PetscReal *error; /* Array to store the errors for convergence analysis */
1148: PetscMPIInt size; /* No of processors */
1149: PetscBool flag; /* Flag denoting availability of exact solution */
1150: PetscInt r;
1152: /* Initialize program */
1153: PetscInitialize(&argc,&argv,(char*)0,help);
1155: /* Check if running with only 1 proc */
1156: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1157: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1159: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
1160: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1161: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1162: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1163: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1164: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1165: PetscOptionsEnd();
1167: PetscMalloc1(n_refine,&error);
1168: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1169: error[r] = 0;
1170: if (r > 0) dt /= refine_fac;
1172: 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]));
1173: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1174: if (flag) {
1175: /* If exact solution available for the specified ODE */
1176: if (r > 0) {
1177: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1178: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f\n.",(double)error[r],(double)conv_rate);
1179: } else {
1180: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1181: }
1182: }
1183: }
1184: PetscFree(error);
1186: /* Exit */
1187: PetscFinalize();
1188: return(0);
1189: }