Actual source code: ex31.c
1: static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";
3: /*
5: Concepts: TS
6: Useful command line parameters:
7: -problem <hull1972a1>: choose which problem to solve (see references
8: for complete listing of problems).
9: -ts_type <euler>: specify time-integrator
10: -ts_adapt_type <basic>: specify time-step adapting (none,basic,advanced)
11: -refinement_levels <1>: number of refinement levels for convergence analysis
12: -refinement_factor <2.0>: factor to refine time step size by for convergence analysis
13: -dt <0.01>: specify time step (initial time step for convergence analysis)
15: */
17: /*
18: List of cases and their names in the code:-
19: From Hull, T.E., Enright, W.H., Fellen, B.M., and Sedgwick, A.E.,
20: "Comparing Numerical Methods for Ordinary Differential
21: Equations", SIAM J. Numer. Anal., 9(4), 1972, pp. 603 - 635
22: A1 -> "hull1972a1" (exact solution available)
23: A2 -> "hull1972a2" (exact solution available)
24: A3 -> "hull1972a3" (exact solution available)
25: A4 -> "hull1972a4" (exact solution available)
26: A5 -> "hull1972a5"
27: B1 -> "hull1972b1"
28: B2 -> "hull1972b2"
29: B3 -> "hull1972b3"
30: B4 -> "hull1972b4"
31: B5 -> "hull1972b5"
32: C1 -> "hull1972c1"
33: C2 -> "hull1972c2"
34: C3 -> "hull1972c3"
35: C4 -> "hull1972c4"
37: From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
38: https://arxiv.org/abs/1503.05166, 2016
40: Kulikov2013I -> "kulik2013i"
42: */
44: #include <petscts.h>
46: /* Function declarations */
47: PetscErrorCode (*RHSFunction) (TS,PetscReal,Vec,Vec,void*);
48: PetscErrorCode (*RHSJacobian) (TS,PetscReal,Vec,Mat,Mat,void*);
49: PetscErrorCode (*IFunction) (TS,PetscReal,Vec,Vec,Vec,void*);
50: PetscErrorCode (*IJacobian) (TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
52: /* Returns the size of the system of equations depending on problem specification */
53: PetscInt GetSize(const char *p)
54: {
56: if ((!strcmp(p,"hull1972a1"))
57: ||(!strcmp(p,"hull1972a2"))
58: ||(!strcmp(p,"hull1972a3"))
59: ||(!strcmp(p,"hull1972a4"))
60: ||(!strcmp(p,"hull1972a5"))) PetscFunctionReturn(1);
61: else if (!strcmp(p,"hull1972b1")) PetscFunctionReturn(2);
62: else if ((!strcmp(p,"hull1972b2"))
63: ||(!strcmp(p,"hull1972b3"))
64: ||(!strcmp(p,"hull1972b4"))
65: ||(!strcmp(p,"hull1972b5"))) PetscFunctionReturn(3);
66: else if ((!strcmp(p,"kulik2013i"))) PetscFunctionReturn(4);
67: else if ((!strcmp(p,"hull1972c1"))
68: ||(!strcmp(p,"hull1972c2"))
69: ||(!strcmp(p,"hull1972c3"))) PetscFunctionReturn(10);
70: else if (!strcmp(p,"hull1972c4")) PetscFunctionReturn(51);
71: else PetscFunctionReturn(-1);
72: }
74: /****************************************************************/
76: /* Problem specific functions */
78: /* Hull, 1972, Problem A1 */
80: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
81: {
82: PetscErrorCode ierr;
83: PetscScalar *f;
84: const PetscScalar *y;
87: VecGetArrayRead(Y,&y);
88: VecGetArray(F,&f);
89: f[0] = -y[0];
90: VecRestoreArrayRead(Y,&y);
91: VecRestoreArray(F,&f);
92: return(0);
93: }
95: PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
96: {
97: PetscErrorCode ierr;
98: const PetscScalar *y;
99: PetscInt row = 0,col = 0;
100: PetscScalar value = -1.0;
103: VecGetArrayRead(Y,&y);
104: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
105: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
107: VecRestoreArrayRead(Y,&y);
108: return(0);
109: }
111: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
112: {
113: PetscErrorCode ierr;
114: const PetscScalar *y;
115: PetscScalar *f;
118: VecGetArrayRead(Y,&y);
119: VecGetArray(F,&f);
120: f[0] = -y[0];
121: VecRestoreArrayRead(Y,&y);
122: VecRestoreArray(F,&f);
123: /* Left hand side = ydot - f(y) */
124: VecAYPX(F,-1.0,Ydot);
125: return(0);
126: }
128: PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
129: {
130: PetscErrorCode ierr;
131: const PetscScalar *y;
132: PetscInt row = 0,col = 0;
133: PetscScalar value = a - 1.0;
136: VecGetArrayRead(Y,&y);
137: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
138: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
139: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
140: VecRestoreArrayRead(Y,&y);
141: return(0);
142: }
144: /* Hull, 1972, Problem A2 */
146: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
147: {
148: PetscErrorCode ierr;
149: const PetscScalar *y;
150: PetscScalar *f;
153: VecGetArrayRead(Y,&y);
154: VecGetArray(F,&f);
155: f[0] = -0.5*y[0]*y[0]*y[0];
156: VecRestoreArrayRead(Y,&y);
157: VecRestoreArray(F,&f);
158: return(0);
159: }
161: PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
162: {
163: PetscErrorCode ierr;
164: const PetscScalar *y;
165: PetscInt row = 0,col = 0;
166: PetscScalar value;
169: VecGetArrayRead(Y,&y);
170: value = -0.5*3.0*y[0]*y[0];
171: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
172: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
173: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
174: VecRestoreArrayRead(Y,&y);
175: return(0);
176: }
178: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
179: {
180: PetscErrorCode ierr;
181: PetscScalar *f;
182: const PetscScalar *y;
185: VecGetArrayRead(Y,&y);
186: VecGetArray(F,&f);
187: f[0] = -0.5*y[0]*y[0]*y[0];
188: VecRestoreArrayRead(Y,&y);
189: VecRestoreArray(F,&f);
190: /* Left hand side = ydot - f(y) */
191: VecAYPX(F,-1.0,Ydot);
192: return(0);
193: }
195: PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
196: {
197: PetscErrorCode ierr;
198: const PetscScalar *y;
199: PetscInt row = 0,col = 0;
200: PetscScalar value;
203: VecGetArrayRead(Y,&y);
204: value = a + 0.5*3.0*y[0]*y[0];
205: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
206: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
208: VecRestoreArrayRead(Y,&y);
209: return(0);
210: }
212: /* Hull, 1972, Problem A3 */
214: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
215: {
216: PetscErrorCode ierr;
217: const PetscScalar *y;
218: PetscScalar *f;
221: VecGetArrayRead(Y,&y);
222: VecGetArray(F,&f);
223: f[0] = y[0]*PetscCosReal(t);
224: VecRestoreArrayRead(Y,&y);
225: VecRestoreArray(F,&f);
226: return(0);
227: }
229: PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
230: {
231: PetscErrorCode ierr;
232: const PetscScalar *y;
233: PetscInt row = 0,col = 0;
234: PetscScalar value = PetscCosReal(t);
237: VecGetArrayRead(Y,&y);
238: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
239: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
240: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
241: VecRestoreArrayRead(Y,&y);
242: return(0);
243: }
245: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
246: {
247: PetscErrorCode ierr;
248: PetscScalar *f;
249: const PetscScalar *y;
252: VecGetArrayRead(Y,&y);
253: VecGetArray(F,&f);
254: f[0] = y[0]*PetscCosReal(t);
255: VecRestoreArrayRead(Y,&y);
256: VecRestoreArray(F,&f);
257: /* Left hand side = ydot - f(y) */
258: VecAYPX(F,-1.0,Ydot);
259: return(0);
260: }
262: PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
263: {
264: PetscErrorCode ierr;
265: const PetscScalar *y;
266: PetscInt row = 0,col = 0;
267: PetscScalar value = a - PetscCosReal(t);
270: VecGetArrayRead(Y,&y);
271: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
272: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
273: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
274: VecRestoreArrayRead(Y,&y);
275: return(0);
276: }
278: /* Hull, 1972, Problem A4 */
280: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
281: {
282: PetscErrorCode ierr;
283: PetscScalar *f;
284: const PetscScalar *y;
287: VecGetArrayRead(Y,&y);
288: VecGetArray(F,&f);
289: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
290: VecRestoreArrayRead(Y,&y);
291: VecRestoreArray(F,&f);
292: return(0);
293: }
295: PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
296: {
297: PetscErrorCode ierr;
298: const PetscScalar *y;
299: PetscInt row = 0,col = 0;
300: PetscScalar value;
303: VecGetArrayRead(Y,&y);
304: value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
305: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
306: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
307: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
308: VecRestoreArrayRead(Y,&y);
309: return(0);
310: }
312: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
313: {
314: PetscErrorCode ierr;
315: PetscScalar *f;
316: const PetscScalar *y;
319: VecGetArrayRead(Y,&y);
320: VecGetArray(F,&f);
321: f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
322: VecRestoreArrayRead(Y,&y);
323: VecRestoreArray(F,&f);
324: /* Left hand side = ydot - f(y) */
325: VecAYPX(F,-1.0,Ydot);
326: return(0);
327: }
329: PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
330: {
331: PetscErrorCode ierr;
332: const PetscScalar *y;
333: PetscInt row = 0,col = 0;
334: PetscScalar value;
337: VecGetArrayRead(Y,&y);
338: value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
339: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
340: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
341: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
342: VecRestoreArrayRead(Y,&y);
343: return(0);
344: }
346: /* Hull, 1972, Problem A5 */
348: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
349: {
350: PetscErrorCode ierr;
351: PetscScalar *f;
352: const PetscScalar *y;
355: VecGetArrayRead(Y,&y);
356: VecGetArray(F,&f);
357: f[0] = (y[0]-t)/(y[0]+t);
358: VecRestoreArrayRead(Y,&y);
359: VecRestoreArray(F,&f);
360: return(0);
361: }
363: PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
364: {
365: PetscErrorCode ierr;
366: const PetscScalar *y;
367: PetscInt row = 0,col = 0;
368: PetscScalar value;
371: VecGetArrayRead(Y,&y);
372: value = 2*t/((t+y[0])*(t+y[0]));
373: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
374: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
375: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
376: VecRestoreArrayRead(Y,&y);
377: return(0);
378: }
380: PetscErrorCode IFunction_Hull1972A5(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] = (y[0]-t)/(y[0]+t);
390: VecRestoreArrayRead(Y,&y);
391: VecRestoreArray(F,&f);
392: /* Left hand side = ydot - f(y) */
393: VecAYPX(F,-1.0,Ydot);
394: return(0);
395: }
397: PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
398: {
399: PetscErrorCode ierr;
400: const PetscScalar *y;
401: PetscInt row = 0,col = 0;
402: PetscScalar value;
405: VecGetArrayRead(Y,&y);
406: value = a - 2*t/((t+y[0])*(t+y[0]));
407: MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);
408: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
409: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
410: VecRestoreArrayRead(Y,&y);
411: return(0);
412: }
414: /* Hull, 1972, Problem B1 */
416: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
417: {
418: PetscErrorCode ierr;
419: PetscScalar *f;
420: const PetscScalar *y;
423: VecGetArrayRead(Y,&y);
424: VecGetArray(F,&f);
425: f[0] = 2.0*(y[0] - y[0]*y[1]);
426: f[1] = -(y[1]-y[0]*y[1]);
427: VecRestoreArrayRead(Y,&y);
428: VecRestoreArray(F,&f);
429: return(0);
430: }
432: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
433: {
434: PetscErrorCode ierr;
435: PetscScalar *f;
436: const PetscScalar *y;
439: VecGetArrayRead(Y,&y);
440: VecGetArray(F,&f);
441: f[0] = 2.0*(y[0] - y[0]*y[1]);
442: f[1] = -(y[1]-y[0]*y[1]);
443: VecRestoreArrayRead(Y,&y);
444: VecRestoreArray(F,&f);
445: /* Left hand side = ydot - f(y) */
446: VecAYPX(F,-1.0,Ydot);
447: return(0);
448: }
450: PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
451: {
452: PetscErrorCode ierr;
453: const PetscScalar *y;
454: PetscInt row[2] = {0,1};
455: PetscScalar value[2][2];
458: VecGetArrayRead(Y,&y);
459: value[0][0] = a - 2.0*(1.0-y[1]); value[0][1] = 2.0*y[0];
460: value[1][0] = -y[1]; value[1][1] = a + 1.0 - y[0];
461: MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);
462: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
464: VecRestoreArrayRead(Y,&y);
465: return(0);
466: }
468: /* Hull, 1972, Problem B2 */
470: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
471: {
472: PetscErrorCode ierr;
473: PetscScalar *f;
474: const PetscScalar *y;
477: VecGetArrayRead(Y,&y);
478: VecGetArray(F,&f);
479: f[0] = -y[0] + y[1];
480: f[1] = y[0] - 2.0*y[1] + y[2];
481: f[2] = y[1] - y[2];
482: VecRestoreArrayRead(Y,&y);
483: VecRestoreArray(F,&f);
484: return(0);
485: }
487: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
488: {
489: PetscErrorCode ierr;
490: PetscScalar *f;
491: const PetscScalar *y;
494: VecGetArrayRead(Y,&y);
495: VecGetArray(F,&f);
496: f[0] = -y[0] + y[1];
497: f[1] = y[0] - 2.0*y[1] + y[2];
498: f[2] = y[1] - y[2];
499: VecRestoreArrayRead(Y,&y);
500: VecRestoreArray(F,&f);
501: /* Left hand side = ydot - f(y) */
502: VecAYPX(F,-1.0,Ydot);
503: return(0);
504: }
506: PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
507: {
508: PetscErrorCode ierr;
509: const PetscScalar *y;
510: PetscInt row[3] = {0,1,2};
511: PetscScalar value[3][3];
514: VecGetArrayRead(Y,&y);
515: value[0][0] = a + 1.0; value[0][1] = -1.0; value[0][2] = 0;
516: value[1][0] = -1.0; value[1][1] = a + 2.0; value[1][2] = -1.0;
517: value[2][0] = 0; value[2][1] = -1.0; value[2][2] = a + 1.0;
518: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
519: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
520: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
521: VecRestoreArrayRead(Y,&y);
522: return(0);
523: }
525: /* Hull, 1972, Problem B3 */
527: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
528: {
529: PetscErrorCode ierr;
530: PetscScalar *f;
531: const PetscScalar *y;
534: VecGetArrayRead(Y,&y);
535: VecGetArray(F,&f);
536: f[0] = -y[0];
537: f[1] = y[0] - y[1]*y[1];
538: f[2] = y[1]*y[1];
539: VecRestoreArrayRead(Y,&y);
540: VecRestoreArray(F,&f);
541: return(0);
542: }
544: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
545: {
546: PetscErrorCode ierr;
547: PetscScalar *f;
548: const PetscScalar *y;
551: VecGetArrayRead(Y,&y);
552: VecGetArray(F,&f);
553: f[0] = -y[0];
554: f[1] = y[0] - y[1]*y[1];
555: f[2] = y[1]*y[1];
556: VecRestoreArrayRead(Y,&y);
557: VecRestoreArray(F,&f);
558: /* Left hand side = ydot - f(y) */
559: VecAYPX(F,-1.0,Ydot);
560: return(0);
561: }
563: PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
564: {
565: PetscErrorCode ierr;
566: const PetscScalar *y;
567: PetscInt row[3] = {0,1,2};
568: PetscScalar value[3][3];
571: VecGetArrayRead(Y,&y);
572: value[0][0] = a + 1.0; value[0][1] = 0; value[0][2] = 0;
573: value[1][0] = -1.0; value[1][1] = a + 2.0*y[1]; value[1][2] = 0;
574: value[2][0] = 0; value[2][1] = -2.0*y[1]; value[2][2] = a;
575: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
576: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
577: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
578: VecRestoreArrayRead(Y,&y);
579: return(0);
580: }
582: /* Hull, 1972, Problem B4 */
584: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
585: {
586: PetscErrorCode ierr;
587: PetscScalar *f;
588: const PetscScalar *y;
591: VecGetArrayRead(Y,&y);
592: VecGetArray(F,&f);
593: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
594: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
595: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
596: VecRestoreArrayRead(Y,&y);
597: VecRestoreArray(F,&f);
598: return(0);
599: }
601: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
602: {
603: PetscErrorCode ierr;
604: PetscScalar *f;
605: const PetscScalar *y;
608: VecGetArrayRead(Y,&y);
609: VecGetArray(F,&f);
610: f[0] = -y[1] - y[0]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
611: f[1] = y[0] - y[1]*y[2]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
612: f[2] = y[0]/PetscSqrtScalar(y[0]*y[0]+y[1]*y[1]);
613: VecRestoreArrayRead(Y,&y);
614: VecRestoreArray(F,&f);
615: /* Left hand side = ydot - f(y) */
616: VecAYPX(F,-1.0,Ydot);
617: return(0);
618: }
620: PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
621: {
622: PetscErrorCode ierr;
623: const PetscScalar *y;
624: PetscInt row[3] = {0,1,2};
625: PetscScalar value[3][3],fac,fac2;
628: VecGetArrayRead(Y,&y);
629: fac = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-1.5);
630: fac2 = PetscPowScalar(y[0]*y[0]+y[1]*y[1],-0.5);
631: value[0][0] = a + (y[1]*y[1]*y[2])*fac;
632: value[0][1] = 1.0 - (y[0]*y[1]*y[2])*fac;
633: value[0][2] = y[0]*fac2;
634: value[1][0] = -1.0 - y[0]*y[1]*y[2]*fac;
635: value[1][1] = a + y[0]*y[0]*y[2]*fac;
636: value[1][2] = y[1]*fac2;
637: value[2][0] = -y[1]*y[1]*fac;
638: value[2][1] = y[0]*y[1]*fac;
639: value[2][2] = a;
640: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
641: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
642: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
643: VecRestoreArrayRead(Y,&y);
644: return(0);
645: }
647: /* Hull, 1972, Problem B5 */
649: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
650: {
651: PetscErrorCode ierr;
652: PetscScalar *f;
653: const PetscScalar *y;
656: VecGetArrayRead(Y,&y);
657: VecGetArray(F,&f);
658: f[0] = y[1]*y[2];
659: f[1] = -y[0]*y[2];
660: f[2] = -0.51*y[0]*y[1];
661: VecRestoreArrayRead(Y,&y);
662: VecRestoreArray(F,&f);
663: return(0);
664: }
666: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
667: {
668: PetscErrorCode ierr;
669: PetscScalar *f;
670: const PetscScalar *y;
673: VecGetArrayRead(Y,&y);
674: VecGetArray(F,&f);
675: f[0] = y[1]*y[2];
676: f[1] = -y[0]*y[2];
677: f[2] = -0.51*y[0]*y[1];
678: VecRestoreArrayRead(Y,&y);
679: VecRestoreArray(F,&f);
680: /* Left hand side = ydot - f(y) */
681: VecAYPX(F,-1.0,Ydot);
682: return(0);
683: }
685: PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
686: {
687: PetscErrorCode ierr;
688: const PetscScalar *y;
689: PetscInt row[3] = {0,1,2};
690: PetscScalar value[3][3];
693: VecGetArrayRead(Y,&y);
694: value[0][0] = a; value[0][1] = -y[2]; value[0][2] = -y[1];
695: value[1][0] = y[2]; value[1][1] = a; value[1][2] = y[0];
696: value[2][0] = 0.51*y[1]; value[2][1] = 0.51*y[0]; value[2][2] = a;
697: MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);
698: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
699: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
700: VecRestoreArrayRead(Y,&y);
701: return(0);
702: }
705: /* Kulikov, 2013, Problem I */
707: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
708: {
709: PetscErrorCode ierr;
710: PetscScalar *f;
711: const PetscScalar *y;
714: VecGetArrayRead(Y,&y);
715: VecGetArray(F,&f);
716: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
717: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
718: f[2] = 2.*t*y[3];
719: f[3] = -2.*t*PetscLogScalar(y[0]);
720: VecRestoreArrayRead(Y,&y);
721: VecRestoreArray(F,&f);
722: return(0);
723: }
725: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
726: {
727: PetscErrorCode ierr;
728: const PetscScalar *y;
729: PetscInt row[4] = {0,1,2,3};
730: PetscScalar value[4][4];
731: PetscScalar m1,m2;
733: VecGetArrayRead(Y,&y);
734: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
735: m2=2.*t*PetscPowScalar(y[1],1./5.);
736: value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
737: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
738: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
739: value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
740: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t;
741: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.;
742: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
743: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
744: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
745: VecRestoreArrayRead(Y,&y);
746: return(0);
747: }
749: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
750: {
751: PetscErrorCode ierr;
752: PetscScalar *f;
753: const PetscScalar *y;
756: VecGetArrayRead(Y,&y);
757: VecGetArray(F,&f);
758: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
759: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
760: f[2] = 2.*t*y[3];
761: f[3] = -2.*t*PetscLogScalar(y[0]);
762: VecRestoreArrayRead(Y,&y);
763: VecRestoreArray(F,&f);
764: /* Left hand side = ydot - f(y) */
765: VecAYPX(F,-1.0,Ydot);
766: return(0);
767: }
769: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
770: {
771: PetscErrorCode ierr;
772: const PetscScalar *y;
773: PetscInt row[4] = {0,1,2,3};
774: PetscScalar value[4][4];
775: PetscScalar m1,m2;
778: VecGetArrayRead(Y,&y);
779: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
780: m2=2.*t*PetscPowScalar(y[1],1./5.);
781: value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
782: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
783: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
784: value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2;
785: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t;
786: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a;
787: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
788: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
789: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
790: VecRestoreArrayRead(Y,&y);
791: return(0);
792: }
795: /* Hull, 1972, Problem C1 */
797: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
798: {
799: PetscErrorCode ierr;
800: PetscScalar *f;
801: const PetscScalar *y;
802: PetscInt N,i;
805: VecGetSize (Y,&N);
806: VecGetArrayRead(Y,&y);
807: VecGetArray(F,&f);
808: f[0] = -y[0];
809: for (i = 1; i < N-1; i++) {
810: f[i] = y[i-1] - y[i];
811: }
812: f[N-1] = y[N-2];
813: VecRestoreArrayRead(Y,&y);
814: VecRestoreArray(F,&f);
815: return(0);
816: }
818: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
819: {
820: PetscErrorCode ierr;
821: PetscScalar *f;
822: const PetscScalar *y;
823: PetscInt N,i;
826: VecGetSize (Y,&N);
827: VecGetArrayRead(Y,&y);
828: VecGetArray(F,&f);
829: f[0] = -y[0];
830: for (i = 1; i < N-1; i++) {
831: f[i] = y[i-1] - y[i];
832: }
833: f[N-1] = y[N-2];
834: VecRestoreArrayRead(Y,&y);
835: VecRestoreArray(F,&f);
836: /* Left hand side = ydot - f(y) */
837: VecAYPX(F,-1.0,Ydot);
838: return(0);
839: }
841: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
842: {
843: PetscErrorCode ierr;
844: const PetscScalar *y;
845: PetscInt N,i,col[2];
846: PetscScalar value[2];
849: VecGetSize (Y,&N);
850: VecGetArrayRead(Y,&y);
851: i = 0;
852: value[0] = a+1; col[0] = 0;
853: value[1] = 0; col[1] = 1;
854: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
855: for (i = 0; i < N; i++) {
856: value[0] = -1; col[0] = i-1;
857: value[1] = a+1; col[1] = i;
858: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
859: }
860: i = N-1;
861: value[0] = -1; col[0] = N-2;
862: value[1] = a; col[1] = N-1;
863: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
864: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
865: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
866: VecRestoreArrayRead(Y,&y);
867: return(0);
868: }
870: /* Hull, 1972, Problem C2 */
872: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
873: {
874: PetscErrorCode ierr;
875: const PetscScalar *y;
876: PetscScalar *f;
877: PetscInt N,i;
880: VecGetSize (Y,&N);
881: VecGetArrayRead(Y,&y);
882: VecGetArray(F,&f);
883: f[0] = -y[0];
884: for (i = 1; i < N-1; i++) {
885: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
886: }
887: f[N-1] = (PetscReal)(N-1)*y[N-2];
888: VecRestoreArrayRead(Y,&y);
889: VecRestoreArray(F,&f);
890: return(0);
891: }
893: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
894: {
895: PetscErrorCode ierr;
896: PetscScalar *f;
897: const PetscScalar *y;
898: PetscInt N,i;
901: VecGetSize (Y,&N);
902: VecGetArrayRead(Y,&y);
903: VecGetArray(F,&f);
904: f[0] = -y[0];
905: for (i = 1; i < N-1; i++) {
906: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
907: }
908: f[N-1] = (PetscReal)(N-1)*y[N-2];
909: VecRestoreArrayRead(Y,&y);
910: VecRestoreArray(F,&f);
911: /* Left hand side = ydot - f(y) */
912: VecAYPX(F,-1.0,Ydot);
913: return(0);
914: }
916: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
917: {
918: PetscErrorCode ierr;
919: const PetscScalar *y;
920: PetscInt N,i,col[2];
921: PetscScalar value[2];
924: VecGetSize (Y,&N);
925: VecGetArrayRead(Y,&y);
926: i = 0;
927: value[0] = a+1; col[0] = 0;
928: value[1] = 0; col[1] = 1;
929: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
930: for (i = 0; i < N; i++) {
931: value[0] = -(PetscReal) i; col[0] = i-1;
932: value[1] = a+(PetscReal)(i+1); col[1] = i;
933: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
934: }
935: i = N-1;
936: value[0] = -(PetscReal) (N-1); col[0] = N-2;
937: value[1] = a; col[1] = N-1;
938: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
939: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
940: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
941: VecRestoreArrayRead(Y,&y);
942: return(0);
943: }
945: /* Hull, 1972, Problem C3 and C4 */
947: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
948: {
949: PetscErrorCode ierr;
950: PetscScalar *f;
951: const PetscScalar *y;
952: PetscInt N,i;
955: VecGetSize (Y,&N);
956: VecGetArrayRead(Y,&y);
957: VecGetArray(F,&f);
958: f[0] = -2.0*y[0] + y[1];
959: for (i = 1; i < N-1; i++) {
960: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
961: }
962: f[N-1] = y[N-2] - 2.0*y[N-1];
963: VecRestoreArrayRead(Y,&y);
964: VecRestoreArray(F,&f);
965: return(0);
966: }
968: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
969: {
970: PetscErrorCode ierr;
971: PetscScalar *f;
972: const PetscScalar *y;
973: PetscInt N,i;
976: VecGetSize (Y,&N);
977: VecGetArrayRead(Y,&y);
978: VecGetArray(F,&f);
979: f[0] = -2.0*y[0] + y[1];
980: for (i = 1; i < N-1; i++) {
981: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
982: }
983: f[N-1] = y[N-2] - 2.0*y[N-1];
984: VecRestoreArrayRead(Y,&y);
985: VecRestoreArray(F,&f);
986: /* Left hand side = ydot - f(y) */
987: VecAYPX(F,-1.0,Ydot);
988: return(0);
989: }
991: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
992: {
993: PetscErrorCode ierr;
994: const PetscScalar *y;
995: PetscScalar value[3];
996: PetscInt N,i,col[3];
999: VecGetSize (Y,&N);
1000: VecGetArrayRead(Y,&y);
1001: for (i = 0; i < N; i++) {
1002: if (i == 0) {
1003: value[0] = a+2; col[0] = i;
1004: value[1] = -1; col[1] = i+1;
1005: value[2] = 0; col[2] = i+2;
1006: } else if (i == N-1) {
1007: value[0] = 0; col[0] = i-2;
1008: value[1] = -1; col[1] = i-1;
1009: value[2] = a+2; col[2] = i;
1010: } else {
1011: value[0] = -1; col[0] = i-1;
1012: value[1] = a+2; col[1] = i;
1013: value[2] = -1; col[2] = i+1;
1014: }
1015: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
1016: }
1017: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1018: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
1019: VecRestoreArrayRead(Y,&y);
1020: return(0);
1021: }
1023: /***************************************************************************/
1025: /* Sets the initial solution for the IVP and sets up the function pointers*/
1026: PetscErrorCode Initialize(Vec Y, void* s)
1027: {
1029: char *p = (char*) s;
1030: PetscScalar *y;
1031: PetscReal t0;
1032: PetscInt N = GetSize((const char *)s);
1033: PetscBool flg;
1036: VecZeroEntries(Y);
1037: VecGetArray(Y,&y);
1038: if (!strcmp(p,"hull1972a1")) {
1039: y[0] = 1.0;
1040: RHSFunction = RHSFunction_Hull1972A1;
1041: RHSJacobian = RHSJacobian_Hull1972A1;
1042: IFunction = IFunction_Hull1972A1;
1043: IJacobian = IJacobian_Hull1972A1;
1044: } else if (!strcmp(p,"hull1972a2")) {
1045: y[0] = 1.0;
1046: RHSFunction = RHSFunction_Hull1972A2;
1047: RHSJacobian = RHSJacobian_Hull1972A2;
1048: IFunction = IFunction_Hull1972A2;
1049: IJacobian = IJacobian_Hull1972A2;
1050: } else if (!strcmp(p,"hull1972a3")) {
1051: y[0] = 1.0;
1052: RHSFunction = RHSFunction_Hull1972A3;
1053: RHSJacobian = RHSJacobian_Hull1972A3;
1054: IFunction = IFunction_Hull1972A3;
1055: IJacobian = IJacobian_Hull1972A3;
1056: } else if (!strcmp(p,"hull1972a4")) {
1057: y[0] = 1.0;
1058: RHSFunction = RHSFunction_Hull1972A4;
1059: RHSJacobian = RHSJacobian_Hull1972A4;
1060: IFunction = IFunction_Hull1972A4;
1061: IJacobian = IJacobian_Hull1972A4;
1062: } else if (!strcmp(p,"hull1972a5")) {
1063: y[0] = 4.0;
1064: RHSFunction = RHSFunction_Hull1972A5;
1065: RHSJacobian = RHSJacobian_Hull1972A5;
1066: IFunction = IFunction_Hull1972A5;
1067: IJacobian = IJacobian_Hull1972A5;
1068: } else if (!strcmp(p,"hull1972b1")) {
1069: y[0] = 1.0;
1070: y[1] = 3.0;
1071: RHSFunction = RHSFunction_Hull1972B1;
1072: IFunction = IFunction_Hull1972B1;
1073: IJacobian = IJacobian_Hull1972B1;
1074: } else if (!strcmp(p,"hull1972b2")) {
1075: y[0] = 2.0;
1076: y[1] = 0.0;
1077: y[2] = 1.0;
1078: RHSFunction = RHSFunction_Hull1972B2;
1079: IFunction = IFunction_Hull1972B2;
1080: IJacobian = IJacobian_Hull1972B2;
1081: } else if (!strcmp(p,"hull1972b3")) {
1082: y[0] = 1.0;
1083: y[1] = 0.0;
1084: y[2] = 0.0;
1085: RHSFunction = RHSFunction_Hull1972B3;
1086: IFunction = IFunction_Hull1972B3;
1087: IJacobian = IJacobian_Hull1972B3;
1088: } else if (!strcmp(p,"hull1972b4")) {
1089: y[0] = 3.0;
1090: y[1] = 0.0;
1091: y[2] = 0.0;
1092: RHSFunction = RHSFunction_Hull1972B4;
1093: IFunction = IFunction_Hull1972B4;
1094: IJacobian = IJacobian_Hull1972B4;
1095: } else if (!strcmp(p,"hull1972b5")) {
1096: y[0] = 0.0;
1097: y[1] = 1.0;
1098: y[2] = 1.0;
1099: RHSFunction = RHSFunction_Hull1972B5;
1100: IFunction = IFunction_Hull1972B5;
1101: IJacobian = IJacobian_Hull1972B5;
1102: } else if (!strcmp(p,"kulik2013i")) {
1103: t0=0.;
1104: y[0] = PetscExpReal(PetscSinReal(t0*t0));
1105: y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1106: y[2] = PetscSinReal(t0*t0)+1.0;
1107: y[3] = PetscCosReal(t0*t0);
1108: RHSFunction = RHSFunction_Kulikov2013I;
1109: RHSJacobian = RHSJacobian_Kulikov2013I;
1110: IFunction = IFunction_Kulikov2013I;
1111: IJacobian = IJacobian_Kulikov2013I;
1112: } else if (!strcmp(p,"hull1972c1")) {
1113: y[0] = 1.0;
1114: RHSFunction = RHSFunction_Hull1972C1;
1115: IFunction = IFunction_Hull1972C1;
1116: IJacobian = IJacobian_Hull1972C1;
1117: } else if (!strcmp(p,"hull1972c2")) {
1118: y[0] = 1.0;
1119: RHSFunction = RHSFunction_Hull1972C2;
1120: IFunction = IFunction_Hull1972C2;
1121: IJacobian = IJacobian_Hull1972C2;
1122: } else if ((!strcmp(p,"hull1972c3"))
1123: ||(!strcmp(p,"hull1972c4"))){
1124: y[0] = 1.0;
1125: RHSFunction = RHSFunction_Hull1972C34;
1126: IFunction = IFunction_Hull1972C34;
1127: IJacobian = IJacobian_Hull1972C34;
1128: }
1129: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1130: 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));
1131: VecRestoreArray(Y,&y);
1132: return(0);
1133: }
1135: /* Calculates the exact solution to problems that have one */
1136: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1137: {
1139: char *p = (char*) s;
1140: PetscScalar *y;
1143: if (!strcmp(p,"hull1972a1")) {
1144: VecGetArray(Y,&y);
1145: y[0] = PetscExpReal(-t);
1146: *flag = PETSC_TRUE;
1147: VecRestoreArray(Y,&y);
1148: } else if (!strcmp(p,"hull1972a2")) {
1149: VecGetArray(Y,&y);
1150: y[0] = 1.0/PetscSqrtReal(t+1);
1151: *flag = PETSC_TRUE;
1152: VecRestoreArray(Y,&y);
1153: } else if (!strcmp(p,"hull1972a3")) {
1154: VecGetArray(Y,&y);
1155: y[0] = PetscExpReal(PetscSinReal(t));
1156: *flag = PETSC_TRUE;
1157: VecRestoreArray(Y,&y);
1158: } else if (!strcmp(p,"hull1972a4")) {
1159: VecGetArray(Y,&y);
1160: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1161: *flag = PETSC_TRUE;
1162: VecRestoreArray(Y,&y);
1163: } else if (!strcmp(p,"kulik2013i")) {
1164: VecGetArray(Y,&y);
1165: y[0] = PetscExpReal(PetscSinReal(t*t));
1166: y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1167: y[2] = PetscSinReal(t*t)+1.0;
1168: y[3] = PetscCosReal(t*t);
1169: *flag = PETSC_TRUE;
1170: VecRestoreArray(Y,&y);
1171: } else {
1172: VecSet(Y,0);
1173: *flag = PETSC_FALSE;
1174: }
1175: return(0);
1176: }
1178: /* Solves the specified ODE and computes the error if exact solution is available */
1179: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1180: {
1181: PetscErrorCode ierr; /* Error code */
1182: TS ts; /* time-integrator */
1183: Vec Y; /* Solution vector */
1184: Vec Yex; /* Exact solution */
1185: PetscInt N; /* Size of the system of equations */
1186: TSType time_scheme; /* Type of time-integration scheme */
1187: Mat Jac = NULL; /* Jacobian matrix */
1188: Vec Yerr; /* Auxiliary solution vector */
1189: PetscReal err_norm; /* Estimated error norm */
1190: PetscReal final_time; /* Actual final time from the integrator */
1193: N = GetSize((const char *)&ptype[0]);
1194: if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1195: VecCreate(PETSC_COMM_WORLD,&Y);
1196: VecSetSizes(Y,N,PETSC_DECIDE);
1197: VecSetUp(Y);
1198: VecSet(Y,0);
1200: /* Initialize the problem */
1201: Initialize(Y,&ptype[0]);
1203: /* Create and initialize the time-integrator */
1204: TSCreate(PETSC_COMM_WORLD,&ts);
1205: /* Default time integration options */
1206: TSSetType(ts,TSRK);
1207: TSSetMaxSteps(ts,maxiter);
1208: TSSetMaxTime(ts,tfinal);
1209: TSSetTimeStep(ts,dt);
1210: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1211: /* Read command line options for time integration */
1212: TSSetFromOptions(ts);
1213: /* Set solution vector */
1214: TSSetSolution(ts,Y);
1215: /* Specify left/right-hand side functions */
1216: TSGetType(ts,&time_scheme);
1218: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1219: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1220: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1221: MatCreate(PETSC_COMM_WORLD,&Jac);
1222: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1223: MatSetFromOptions(Jac);
1224: MatSetUp(Jac);
1225: TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1226: } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1227: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1228: /* and its Jacobian function */
1229: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1230: MatCreate(PETSC_COMM_WORLD,&Jac);
1231: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1232: MatSetFromOptions(Jac);
1233: MatSetUp(Jac);
1234: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1235: }
1237: /* Solve */
1238: TSSolve(ts,Y);
1239: TSGetTime(ts,&final_time);
1241: /* Get the estimated error, if available */
1242: VecDuplicate(Y,&Yerr);
1243: VecZeroEntries(Yerr);
1244: TSGetTimeError(ts,0,&Yerr);
1245: VecNorm(Yerr,NORM_2,&err_norm);
1246: VecDestroy(&Yerr);
1247: PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);
1249: /* Exact solution */
1250: VecDuplicate(Y,&Yex);
1251: if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1252: PetscPrintf(PETSC_COMM_WORLD,"Note: There is a difference between the prescribed final time %g and the actual final time, %g.\n",(double)tfinal,(double)final_time);
1253: }
1254: ExactSolution(Yex,&ptype[0],final_time,exact_flag);
1256: /* Calculate Error */
1257: VecAYPX(Yex,-1.0,Y);
1258: VecNorm(Yex,NORM_2,error);
1259: *error = PetscSqrtReal(((*error)*(*error))/N);
1261: /* Clean up and finalize */
1262: MatDestroy(&Jac);
1263: TSDestroy(&ts);
1264: VecDestroy(&Yex);
1265: VecDestroy(&Y);
1267: return(0);
1268: }
1270: int main(int argc, char **argv)
1271: {
1272: PetscErrorCode ierr; /* Error code */
1273: char ptype[256] = "hull1972a1"; /* Problem specification */
1274: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1275: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1276: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1277: PetscReal dt;
1278: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1279: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1280: PetscReal *error; /* Array to store the errors for convergence analysis */
1281: PetscMPIInt size; /* No of processors */
1282: PetscBool flag; /* Flag denoting availability of exact solution */
1283: PetscInt r;
1285: /* Initialize program */
1286: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
1288: /* Check if running with only 1 proc */
1289: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1290: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1292: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1293: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1294: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1295: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1296: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1297: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1298: PetscOptionsEnd();
1300: PetscMalloc1(n_refine,&error);
1301: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1302: error[r] = 0;
1303: if (r > 0) dt /= refine_fac;
1305: 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]));
1306: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1307: if (flag) {
1308: /* If exact solution available for the specified ODE */
1309: if (r > 0) {
1310: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1311: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1312: } else {
1313: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1314: }
1315: }
1316: }
1317: PetscFree(error);
1318: PetscFinalize();
1319: return ierr;
1320: }
1322: /*TEST
1324: test:
1325: suffix: 2
1326: args: -ts_type glee -final_time 5 -ts_adapt_type none
1327: timeoutfactor: 3
1328: requires: !single
1330: test:
1331: suffix: 3
1332: args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor -ts_max_steps 50 -problem hull1972a3 -ts_adapt_glee_use_local 1
1333: timeoutfactor: 3
1334: requires: !single
1336: test:
1337: suffix: 4
1338: args: -ts_type glee -final_time 5 -ts_adapt_type glee -ts_adapt_monitor -ts_max_steps 50 -problem hull1972a3 -ts_max_reject 100 -ts_adapt_glee_use_local 0
1339: timeoutfactor: 3
1340: requires: !single !__float128
1342: TEST*/