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: }
704: /* Kulikov, 2013, Problem I */
706: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
707: {
708: PetscErrorCode ierr;
709: PetscScalar *f;
710: const PetscScalar *y;
713: VecGetArrayRead(Y,&y);
714: VecGetArray(F,&f);
715: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
716: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
717: f[2] = 2.*t*y[3];
718: f[3] = -2.*t*PetscLogScalar(y[0]);
719: VecRestoreArrayRead(Y,&y);
720: VecRestoreArray(F,&f);
721: return(0);
722: }
724: PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s)
725: {
726: PetscErrorCode ierr;
727: const PetscScalar *y;
728: PetscInt row[4] = {0,1,2,3};
729: PetscScalar value[4][4];
730: PetscScalar m1,m2;
732: VecGetArrayRead(Y,&y);
733: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
734: m2=2.*t*PetscPowScalar(y[1],1./5.);
735: value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
736: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
737: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
738: value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2;
739: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t;
740: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.;
741: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
742: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
743: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
744: VecRestoreArrayRead(Y,&y);
745: return(0);
746: }
748: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
749: {
750: PetscErrorCode ierr;
751: PetscScalar *f;
752: const PetscScalar *y;
755: VecGetArrayRead(Y,&y);
756: VecGetArray(F,&f);
757: f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3];
758: f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
759: f[2] = 2.*t*y[3];
760: f[3] = -2.*t*PetscLogScalar(y[0]);
761: VecRestoreArrayRead(Y,&y);
762: VecRestoreArray(F,&f);
763: /* Left hand side = ydot - f(y) */
764: VecAYPX(F,-1.0,Ydot);
765: return(0);
766: }
768: PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
769: {
770: PetscErrorCode ierr;
771: const PetscScalar *y;
772: PetscInt row[4] = {0,1,2,3};
773: PetscScalar value[4][4];
774: PetscScalar m1,m2;
777: VecGetArrayRead(Y,&y);
778: m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.));
779: m2=2.*t*PetscPowScalar(y[1],1./5.);
780: value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2;
781: m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.));
782: m2=10.*t*PetscExpScalar(5.0*(y[2]-1.));
783: value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2;
784: value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t;
785: value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a;
786: MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);
787: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
788: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
789: VecRestoreArrayRead(Y,&y);
790: return(0);
791: }
793: /* Hull, 1972, Problem C1 */
795: PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
796: {
797: PetscErrorCode ierr;
798: PetscScalar *f;
799: const PetscScalar *y;
800: PetscInt N,i;
803: VecGetSize (Y,&N);
804: VecGetArrayRead(Y,&y);
805: VecGetArray(F,&f);
806: f[0] = -y[0];
807: for (i = 1; i < N-1; i++) {
808: f[i] = y[i-1] - y[i];
809: }
810: f[N-1] = y[N-2];
811: VecRestoreArrayRead(Y,&y);
812: VecRestoreArray(F,&f);
813: return(0);
814: }
816: PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
817: {
818: PetscErrorCode ierr;
819: PetscScalar *f;
820: const PetscScalar *y;
821: PetscInt N,i;
824: VecGetSize (Y,&N);
825: VecGetArrayRead(Y,&y);
826: VecGetArray(F,&f);
827: f[0] = -y[0];
828: for (i = 1; i < N-1; i++) {
829: f[i] = y[i-1] - y[i];
830: }
831: f[N-1] = y[N-2];
832: VecRestoreArrayRead(Y,&y);
833: VecRestoreArray(F,&f);
834: /* Left hand side = ydot - f(y) */
835: VecAYPX(F,-1.0,Ydot);
836: return(0);
837: }
839: PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
840: {
841: PetscErrorCode ierr;
842: const PetscScalar *y;
843: PetscInt N,i,col[2];
844: PetscScalar value[2];
847: VecGetSize (Y,&N);
848: VecGetArrayRead(Y,&y);
849: i = 0;
850: value[0] = a+1; col[0] = 0;
851: value[1] = 0; col[1] = 1;
852: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
853: for (i = 0; i < N; i++) {
854: value[0] = -1; col[0] = i-1;
855: value[1] = a+1; col[1] = i;
856: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
857: }
858: i = N-1;
859: value[0] = -1; col[0] = N-2;
860: value[1] = a; col[1] = N-1;
861: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
862: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
863: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
864: VecRestoreArrayRead(Y,&y);
865: return(0);
866: }
868: /* Hull, 1972, Problem C2 */
870: PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
871: {
872: PetscErrorCode ierr;
873: const PetscScalar *y;
874: PetscScalar *f;
875: PetscInt N,i;
878: VecGetSize (Y,&N);
879: VecGetArrayRead(Y,&y);
880: VecGetArray(F,&f);
881: f[0] = -y[0];
882: for (i = 1; i < N-1; i++) {
883: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
884: }
885: f[N-1] = (PetscReal)(N-1)*y[N-2];
886: VecRestoreArrayRead(Y,&y);
887: VecRestoreArray(F,&f);
888: return(0);
889: }
891: PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
892: {
893: PetscErrorCode ierr;
894: PetscScalar *f;
895: const PetscScalar *y;
896: PetscInt N,i;
899: VecGetSize (Y,&N);
900: VecGetArrayRead(Y,&y);
901: VecGetArray(F,&f);
902: f[0] = -y[0];
903: for (i = 1; i < N-1; i++) {
904: f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i];
905: }
906: f[N-1] = (PetscReal)(N-1)*y[N-2];
907: VecRestoreArrayRead(Y,&y);
908: VecRestoreArray(F,&f);
909: /* Left hand side = ydot - f(y) */
910: VecAYPX(F,-1.0,Ydot);
911: return(0);
912: }
914: PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
915: {
916: PetscErrorCode ierr;
917: const PetscScalar *y;
918: PetscInt N,i,col[2];
919: PetscScalar value[2];
922: VecGetSize (Y,&N);
923: VecGetArrayRead(Y,&y);
924: i = 0;
925: value[0] = a+1; col[0] = 0;
926: value[1] = 0; col[1] = 1;
927: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
928: for (i = 0; i < N; i++) {
929: value[0] = -(PetscReal) i; col[0] = i-1;
930: value[1] = a+(PetscReal)(i+1); col[1] = i;
931: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
932: }
933: i = N-1;
934: value[0] = -(PetscReal) (N-1); col[0] = N-2;
935: value[1] = a; col[1] = N-1;
936: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
937: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
938: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
939: VecRestoreArrayRead(Y,&y);
940: return(0);
941: }
943: /* Hull, 1972, Problem C3 and C4 */
945: PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s)
946: {
947: PetscErrorCode ierr;
948: PetscScalar *f;
949: const PetscScalar *y;
950: PetscInt N,i;
953: VecGetSize (Y,&N);
954: VecGetArrayRead(Y,&y);
955: VecGetArray(F,&f);
956: f[0] = -2.0*y[0] + y[1];
957: for (i = 1; i < N-1; i++) {
958: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
959: }
960: f[N-1] = y[N-2] - 2.0*y[N-1];
961: VecRestoreArrayRead(Y,&y);
962: VecRestoreArray(F,&f);
963: return(0);
964: }
966: PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
967: {
968: PetscErrorCode ierr;
969: PetscScalar *f;
970: const PetscScalar *y;
971: PetscInt N,i;
974: VecGetSize (Y,&N);
975: VecGetArrayRead(Y,&y);
976: VecGetArray(F,&f);
977: f[0] = -2.0*y[0] + y[1];
978: for (i = 1; i < N-1; i++) {
979: f[i] = y[i-1] - 2.0*y[i] + y[i+1];
980: }
981: f[N-1] = y[N-2] - 2.0*y[N-1];
982: VecRestoreArrayRead(Y,&y);
983: VecRestoreArray(F,&f);
984: /* Left hand side = ydot - f(y) */
985: VecAYPX(F,-1.0,Ydot);
986: return(0);
987: }
989: PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s)
990: {
991: PetscErrorCode ierr;
992: const PetscScalar *y;
993: PetscScalar value[3];
994: PetscInt N,i,col[3];
997: VecGetSize (Y,&N);
998: VecGetArrayRead(Y,&y);
999: for (i = 0; i < N; i++) {
1000: if (i == 0) {
1001: value[0] = a+2; col[0] = i;
1002: value[1] = -1; col[1] = i+1;
1003: value[2] = 0; col[2] = i+2;
1004: } else if (i == N-1) {
1005: value[0] = 0; col[0] = i-2;
1006: value[1] = -1; col[1] = i-1;
1007: value[2] = a+2; col[2] = i;
1008: } else {
1009: value[0] = -1; col[0] = i-1;
1010: value[1] = a+2; col[1] = i;
1011: value[2] = -1; col[2] = i+1;
1012: }
1013: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
1014: }
1015: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1016: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
1017: VecRestoreArrayRead(Y,&y);
1018: return(0);
1019: }
1021: /***************************************************************************/
1023: /* Sets the initial solution for the IVP and sets up the function pointers*/
1024: PetscErrorCode Initialize(Vec Y, void* s)
1025: {
1027: char *p = (char*) s;
1028: PetscScalar *y;
1029: PetscReal t0;
1030: PetscInt N = GetSize((const char *)s);
1031: PetscBool flg;
1034: VecZeroEntries(Y);
1035: VecGetArray(Y,&y);
1036: if (!strcmp(p,"hull1972a1")) {
1037: y[0] = 1.0;
1038: RHSFunction = RHSFunction_Hull1972A1;
1039: RHSJacobian = RHSJacobian_Hull1972A1;
1040: IFunction = IFunction_Hull1972A1;
1041: IJacobian = IJacobian_Hull1972A1;
1042: } else if (!strcmp(p,"hull1972a2")) {
1043: y[0] = 1.0;
1044: RHSFunction = RHSFunction_Hull1972A2;
1045: RHSJacobian = RHSJacobian_Hull1972A2;
1046: IFunction = IFunction_Hull1972A2;
1047: IJacobian = IJacobian_Hull1972A2;
1048: } else if (!strcmp(p,"hull1972a3")) {
1049: y[0] = 1.0;
1050: RHSFunction = RHSFunction_Hull1972A3;
1051: RHSJacobian = RHSJacobian_Hull1972A3;
1052: IFunction = IFunction_Hull1972A3;
1053: IJacobian = IJacobian_Hull1972A3;
1054: } else if (!strcmp(p,"hull1972a4")) {
1055: y[0] = 1.0;
1056: RHSFunction = RHSFunction_Hull1972A4;
1057: RHSJacobian = RHSJacobian_Hull1972A4;
1058: IFunction = IFunction_Hull1972A4;
1059: IJacobian = IJacobian_Hull1972A4;
1060: } else if (!strcmp(p,"hull1972a5")) {
1061: y[0] = 4.0;
1062: RHSFunction = RHSFunction_Hull1972A5;
1063: RHSJacobian = RHSJacobian_Hull1972A5;
1064: IFunction = IFunction_Hull1972A5;
1065: IJacobian = IJacobian_Hull1972A5;
1066: } else if (!strcmp(p,"hull1972b1")) {
1067: y[0] = 1.0;
1068: y[1] = 3.0;
1069: RHSFunction = RHSFunction_Hull1972B1;
1070: IFunction = IFunction_Hull1972B1;
1071: IJacobian = IJacobian_Hull1972B1;
1072: } else if (!strcmp(p,"hull1972b2")) {
1073: y[0] = 2.0;
1074: y[1] = 0.0;
1075: y[2] = 1.0;
1076: RHSFunction = RHSFunction_Hull1972B2;
1077: IFunction = IFunction_Hull1972B2;
1078: IJacobian = IJacobian_Hull1972B2;
1079: } else if (!strcmp(p,"hull1972b3")) {
1080: y[0] = 1.0;
1081: y[1] = 0.0;
1082: y[2] = 0.0;
1083: RHSFunction = RHSFunction_Hull1972B3;
1084: IFunction = IFunction_Hull1972B3;
1085: IJacobian = IJacobian_Hull1972B3;
1086: } else if (!strcmp(p,"hull1972b4")) {
1087: y[0] = 3.0;
1088: y[1] = 0.0;
1089: y[2] = 0.0;
1090: RHSFunction = RHSFunction_Hull1972B4;
1091: IFunction = IFunction_Hull1972B4;
1092: IJacobian = IJacobian_Hull1972B4;
1093: } else if (!strcmp(p,"hull1972b5")) {
1094: y[0] = 0.0;
1095: y[1] = 1.0;
1096: y[2] = 1.0;
1097: RHSFunction = RHSFunction_Hull1972B5;
1098: IFunction = IFunction_Hull1972B5;
1099: IJacobian = IJacobian_Hull1972B5;
1100: } else if (!strcmp(p,"kulik2013i")) {
1101: t0=0.;
1102: y[0] = PetscExpReal(PetscSinReal(t0*t0));
1103: y[1] = PetscExpReal(5.*PetscSinReal(t0*t0));
1104: y[2] = PetscSinReal(t0*t0)+1.0;
1105: y[3] = PetscCosReal(t0*t0);
1106: RHSFunction = RHSFunction_Kulikov2013I;
1107: RHSJacobian = RHSJacobian_Kulikov2013I;
1108: IFunction = IFunction_Kulikov2013I;
1109: IJacobian = IJacobian_Kulikov2013I;
1110: } else if (!strcmp(p,"hull1972c1")) {
1111: y[0] = 1.0;
1112: RHSFunction = RHSFunction_Hull1972C1;
1113: IFunction = IFunction_Hull1972C1;
1114: IJacobian = IJacobian_Hull1972C1;
1115: } else if (!strcmp(p,"hull1972c2")) {
1116: y[0] = 1.0;
1117: RHSFunction = RHSFunction_Hull1972C2;
1118: IFunction = IFunction_Hull1972C2;
1119: IJacobian = IJacobian_Hull1972C2;
1120: } else if ((!strcmp(p,"hull1972c3"))
1121: ||(!strcmp(p,"hull1972c4"))) {
1122: y[0] = 1.0;
1123: RHSFunction = RHSFunction_Hull1972C34;
1124: IFunction = IFunction_Hull1972C34;
1125: IJacobian = IJacobian_Hull1972C34;
1126: }
1127: PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);
1128: 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));
1129: VecRestoreArray(Y,&y);
1130: return(0);
1131: }
1133: /* Calculates the exact solution to problems that have one */
1134: PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1135: {
1137: char *p = (char*) s;
1138: PetscScalar *y;
1141: if (!strcmp(p,"hull1972a1")) {
1142: VecGetArray(Y,&y);
1143: y[0] = PetscExpReal(-t);
1144: *flag = PETSC_TRUE;
1145: VecRestoreArray(Y,&y);
1146: } else if (!strcmp(p,"hull1972a2")) {
1147: VecGetArray(Y,&y);
1148: y[0] = 1.0/PetscSqrtReal(t+1);
1149: *flag = PETSC_TRUE;
1150: VecRestoreArray(Y,&y);
1151: } else if (!strcmp(p,"hull1972a3")) {
1152: VecGetArray(Y,&y);
1153: y[0] = PetscExpReal(PetscSinReal(t));
1154: *flag = PETSC_TRUE;
1155: VecRestoreArray(Y,&y);
1156: } else if (!strcmp(p,"hull1972a4")) {
1157: VecGetArray(Y,&y);
1158: y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1159: *flag = PETSC_TRUE;
1160: VecRestoreArray(Y,&y);
1161: } else if (!strcmp(p,"kulik2013i")) {
1162: VecGetArray(Y,&y);
1163: y[0] = PetscExpReal(PetscSinReal(t*t));
1164: y[1] = PetscExpReal(5.*PetscSinReal(t*t));
1165: y[2] = PetscSinReal(t*t)+1.0;
1166: y[3] = PetscCosReal(t*t);
1167: *flag = PETSC_TRUE;
1168: VecRestoreArray(Y,&y);
1169: } else {
1170: VecSet(Y,0);
1171: *flag = PETSC_FALSE;
1172: }
1173: return(0);
1174: }
1176: /* Solves the specified ODE and computes the error if exact solution is available */
1177: PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag)
1178: {
1179: PetscErrorCode ierr; /* Error code */
1180: TS ts; /* time-integrator */
1181: Vec Y; /* Solution vector */
1182: Vec Yex; /* Exact solution */
1183: PetscInt N; /* Size of the system of equations */
1184: TSType time_scheme; /* Type of time-integration scheme */
1185: Mat Jac = NULL; /* Jacobian matrix */
1186: Vec Yerr; /* Auxiliary solution vector */
1187: PetscReal err_norm; /* Estimated error norm */
1188: PetscReal final_time; /* Actual final time from the integrator */
1191: N = GetSize((const char *)&ptype[0]);
1192: if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1193: VecCreate(PETSC_COMM_WORLD,&Y);
1194: VecSetSizes(Y,N,PETSC_DECIDE);
1195: VecSetUp(Y);
1196: VecSet(Y,0);
1198: /* Initialize the problem */
1199: Initialize(Y,&ptype[0]);
1201: /* Create and initialize the time-integrator */
1202: TSCreate(PETSC_COMM_WORLD,&ts);
1203: /* Default time integration options */
1204: TSSetType(ts,TSRK);
1205: TSSetMaxSteps(ts,maxiter);
1206: TSSetMaxTime(ts,tfinal);
1207: TSSetTimeStep(ts,dt);
1208: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
1209: /* Read command line options for time integration */
1210: TSSetFromOptions(ts);
1211: /* Set solution vector */
1212: TSSetSolution(ts,Y);
1213: /* Specify left/right-hand side functions */
1214: TSGetType(ts,&time_scheme);
1216: if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) {
1217: /* Explicit time-integration -> specify right-hand side function ydot = f(y) */
1218: TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);
1219: MatCreate(PETSC_COMM_WORLD,&Jac);
1220: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1221: MatSetFromOptions(Jac);
1222: MatSetUp(Jac);
1223: TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);
1224: } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) {
1225: /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */
1226: /* and its Jacobian function */
1227: TSSetIFunction(ts,NULL,IFunction,&ptype[0]);
1228: MatCreate(PETSC_COMM_WORLD,&Jac);
1229: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
1230: MatSetFromOptions(Jac);
1231: MatSetUp(Jac);
1232: TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);
1233: }
1235: /* Solve */
1236: TSSolve(ts,Y);
1237: TSGetTime(ts,&final_time);
1239: /* Get the estimated error, if available */
1240: VecDuplicate(Y,&Yerr);
1241: VecZeroEntries(Yerr);
1242: TSGetTimeError(ts,0,&Yerr);
1243: VecNorm(Yerr,NORM_2,&err_norm);
1244: VecDestroy(&Yerr);
1245: PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);
1247: /* Exact solution */
1248: VecDuplicate(Y,&Yex);
1249: if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1250: 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);
1251: }
1252: ExactSolution(Yex,&ptype[0],final_time,exact_flag);
1254: /* Calculate Error */
1255: VecAYPX(Yex,-1.0,Y);
1256: VecNorm(Yex,NORM_2,error);
1257: *error = PetscSqrtReal(((*error)*(*error))/N);
1259: /* Clean up and finalize */
1260: MatDestroy(&Jac);
1261: TSDestroy(&ts);
1262: VecDestroy(&Yex);
1263: VecDestroy(&Y);
1265: return(0);
1266: }
1268: int main(int argc, char **argv)
1269: {
1270: PetscErrorCode ierr; /* Error code */
1271: char ptype[256] = "hull1972a1"; /* Problem specification */
1272: PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */
1273: PetscReal refine_fac = 2.0; /* Refinement factor for dt */
1274: PetscReal dt_initial = 0.01; /* Initial default value of dt */
1275: PetscReal dt;
1276: PetscReal tfinal = 20.0; /* Final time for the time-integration */
1277: PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */
1278: PetscReal *error; /* Array to store the errors for convergence analysis */
1279: PetscMPIInt size; /* No of processors */
1280: PetscBool flag; /* Flag denoting availability of exact solution */
1281: PetscInt r;
1283: /* Initialize program */
1284: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
1286: /* Check if running with only 1 proc */
1287: MPI_Comm_size(PETSC_COMM_WORLD,&size);
1288: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1290: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);
1291: PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);
1292: PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);
1293: PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);
1294: PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);
1295: PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);
1296: PetscOptionsEnd();
1298: PetscMalloc1(n_refine,&error);
1299: for (r = 0,dt = dt_initial; r < n_refine; r++) {
1300: error[r] = 0;
1301: if (r > 0) dt /= refine_fac;
1303: PetscPrintf(PETSC_COMM_WORLD,"Solving ODE \"%s\" with dt %f, final time %f and system size %D.\n",ptype,(double)dt,(double)tfinal,GetSize(&ptype[0]));
1304: SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);
1305: if (flag) {
1306: /* If exact solution available for the specified ODE */
1307: if (r > 0) {
1308: PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac));
1309: PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);
1310: } else {
1311: PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);
1312: }
1313: }
1314: }
1315: PetscFree(error);
1316: PetscFinalize();
1317: return ierr;
1318: }
1320: /*TEST
1322: test:
1323: suffix: 2
1324: args: -ts_type glee -final_time 5 -ts_adapt_type none
1325: timeoutfactor: 3
1326: requires: !single
1328: test:
1329: suffix: 3
1330: 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
1331: timeoutfactor: 3
1332: requires: !single
1334: test:
1335: suffix: 4
1336: 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
1337: timeoutfactor: 3
1338: requires: !single !__float128
1340: TEST*/