1 static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";
2
3 /*
4
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)
14
15 */
16
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
37 From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
38 https://arxiv.org/abs/1503.05166, 2016
39
40 Kulikov2013I -> "kulik2013i"
41
42 */
43
44 #include <petscts.h>
45
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*);
51
52 /* Returns the size of the system of equations depending on problem specification */
GetSize(const char * p)53 PetscInt GetSize(const char *p)
54 {
55 PetscFunctionBegin;
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 }
73
74 /****************************************************************/
75
76 /* Problem specific functions */
77
78 /* Hull, 1972, Problem A1 */
79
RHSFunction_Hull1972A1(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
85
86 PetscFunctionBegin;
87 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
88 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
89 f[0] = -y[0];
90 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
91 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
92 PetscFunctionReturn(0);
93 }
94
RHSJacobian_Hull1972A1(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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;
101
102 PetscFunctionBegin;
103 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
104 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
105 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
106 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
107 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
108 PetscFunctionReturn(0);
109 }
110
IFunction_Hull1972A1(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
116
117 PetscFunctionBegin;
118 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
119 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
120 f[0] = -y[0];
121 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
122 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
123 /* Left hand side = ydot - f(y) */
124 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
125 PetscFunctionReturn(0);
126 }
127
IJacobian_Hull1972A1(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
134
135 PetscFunctionBegin;
136 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
137 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
138 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
139 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
140 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
141 PetscFunctionReturn(0);
142 }
143
144 /* Hull, 1972, Problem A2 */
145
RHSFunction_Hull1972A2(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
151
152 PetscFunctionBegin;
153 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
154 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
155 f[0] = -0.5*y[0]*y[0]*y[0];
156 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
157 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
158 PetscFunctionReturn(0);
159 }
160
RHSJacobian_Hull1972A2(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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;
167
168 PetscFunctionBegin;
169 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
170 value = -0.5*3.0*y[0]*y[0];
171 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
172 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
173 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
174 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
175 PetscFunctionReturn(0);
176 }
177
IFunction_Hull1972A2(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
183
184 PetscFunctionBegin;
185 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
186 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
187 f[0] = -0.5*y[0]*y[0]*y[0];
188 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
189 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
190 /* Left hand side = ydot - f(y) */
191 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
192 PetscFunctionReturn(0);
193 }
194
IJacobian_Hull1972A2(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
201
202 PetscFunctionBegin;
203 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
204 value = a + 0.5*3.0*y[0]*y[0];
205 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
206 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
207 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
208 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
209 PetscFunctionReturn(0);
210 }
211
212 /* Hull, 1972, Problem A3 */
213
RHSFunction_Hull1972A3(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
219
220 PetscFunctionBegin;
221 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
222 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
223 f[0] = y[0]*PetscCosReal(t);
224 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
225 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
226 PetscFunctionReturn(0);
227 }
228
RHSJacobian_Hull1972A3(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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);
235
236 PetscFunctionBegin;
237 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
238 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
239 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
240 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
241 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
242 PetscFunctionReturn(0);
243 }
244
IFunction_Hull1972A3(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
250
251 PetscFunctionBegin;
252 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
253 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
254 f[0] = y[0]*PetscCosReal(t);
255 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
256 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
257 /* Left hand side = ydot - f(y) */
258 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
259 PetscFunctionReturn(0);
260 }
261
IJacobian_Hull1972A3(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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);
268
269 PetscFunctionBegin;
270 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
271 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
272 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
273 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
274 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
275 PetscFunctionReturn(0);
276 }
277
278 /* Hull, 1972, Problem A4 */
279
RHSFunction_Hull1972A4(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
285
286 PetscFunctionBegin;
287 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
288 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
289 f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
290 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
291 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
292 PetscFunctionReturn(0);
293 }
294
RHSJacobian_Hull1972A4(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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;
301
302 PetscFunctionBegin;
303 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
304 value = 0.25*(1.0-0.05*y[0]) - (0.25*y[0])*0.05;
305 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
306 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
307 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
308 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
309 PetscFunctionReturn(0);
310 }
311
IFunction_Hull1972A4(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
317
318 PetscFunctionBegin;
319 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
320 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
321 f[0] = (0.25*y[0])*(1.0-0.05*y[0]);
322 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
323 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
324 /* Left hand side = ydot - f(y) */
325 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
326 PetscFunctionReturn(0);
327 }
328
IJacobian_Hull1972A4(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
335
336 PetscFunctionBegin;
337 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
338 value = a - 0.25*(1.0-0.05*y[0]) + (0.25*y[0])*0.05;
339 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
340 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
341 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
342 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
343 PetscFunctionReturn(0);
344 }
345
346 /* Hull, 1972, Problem A5 */
347
RHSFunction_Hull1972A5(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
353
354 PetscFunctionBegin;
355 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
356 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
357 f[0] = (y[0]-t)/(y[0]+t);
358 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
359 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
360 PetscFunctionReturn(0);
361 }
362
RHSJacobian_Hull1972A5(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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;
369
370 PetscFunctionBegin;
371 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
372 value = 2*t/((t+y[0])*(t+y[0]));
373 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
374 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
375 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
376 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
377 PetscFunctionReturn(0);
378 }
379
IFunction_Hull1972A5(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
385
386 PetscFunctionBegin;
387 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
388 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
389 f[0] = (y[0]-t)/(y[0]+t);
390 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
391 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
392 /* Left hand side = ydot - f(y) */
393 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
394 PetscFunctionReturn(0);
395 }
396
IJacobian_Hull1972A5(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
403
404 PetscFunctionBegin;
405 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
406 value = a - 2*t/((t+y[0])*(t+y[0]));
407 ierr = MatSetValues(A,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
408 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
409 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
410 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
411 PetscFunctionReturn(0);
412 }
413
414 /* Hull, 1972, Problem B1 */
415
RHSFunction_Hull1972B1(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
421
422 PetscFunctionBegin;
423 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
424 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
425 f[0] = 2.0*(y[0] - y[0]*y[1]);
426 f[1] = -(y[1]-y[0]*y[1]);
427 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
428 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
429 PetscFunctionReturn(0);
430 }
431
IFunction_Hull1972B1(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
437
438 PetscFunctionBegin;
439 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
440 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
441 f[0] = 2.0*(y[0] - y[0]*y[1]);
442 f[1] = -(y[1]-y[0]*y[1]);
443 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
444 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
445 /* Left hand side = ydot - f(y) */
446 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
447 PetscFunctionReturn(0);
448 }
449
IJacobian_Hull1972B1(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
456
457 PetscFunctionBegin;
458 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,2,&row[0],2,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
462 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
463 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
464 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
465 PetscFunctionReturn(0);
466 }
467
468 /* Hull, 1972, Problem B2 */
469
RHSFunction_Hull1972B2(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
475
476 PetscFunctionBegin;
477 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
478 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
483 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
484 PetscFunctionReturn(0);
485 }
486
IFunction_Hull1972B2(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
492
493 PetscFunctionBegin;
494 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
495 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
500 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
501 /* Left hand side = ydot - f(y) */
502 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
503 PetscFunctionReturn(0);
504 }
505
IJacobian_Hull1972B2(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
512
513 PetscFunctionBegin;
514 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
519 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
520 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
521 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
522 PetscFunctionReturn(0);
523 }
524
525 /* Hull, 1972, Problem B3 */
526
RHSFunction_Hull1972B3(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
532
533 PetscFunctionBegin;
534 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
535 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
536 f[0] = -y[0];
537 f[1] = y[0] - y[1]*y[1];
538 f[2] = y[1]*y[1];
539 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
540 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
541 PetscFunctionReturn(0);
542 }
543
IFunction_Hull1972B3(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
549
550 PetscFunctionBegin;
551 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
552 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
553 f[0] = -y[0];
554 f[1] = y[0] - y[1]*y[1];
555 f[2] = y[1]*y[1];
556 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
557 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
558 /* Left hand side = ydot - f(y) */
559 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
560 PetscFunctionReturn(0);
561 }
562
IJacobian_Hull1972B3(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
569
570 PetscFunctionBegin;
571 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
576 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
577 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
578 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
579 PetscFunctionReturn(0);
580 }
581
582 /* Hull, 1972, Problem B4 */
583
RHSFunction_Hull1972B4(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
589
590 PetscFunctionBegin;
591 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
592 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
597 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
598 PetscFunctionReturn(0);
599 }
600
IFunction_Hull1972B4(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
606
607 PetscFunctionBegin;
608 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
609 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
614 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
615 /* Left hand side = ydot - f(y) */
616 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
617 PetscFunctionReturn(0);
618 }
619
IJacobian_Hull1972B4(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
626
627 PetscFunctionBegin;
628 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
641 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
642 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
643 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
644 PetscFunctionReturn(0);
645 }
646
647 /* Hull, 1972, Problem B5 */
648
RHSFunction_Hull1972B5(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
654
655 PetscFunctionBegin;
656 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
657 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
658 f[0] = y[1]*y[2];
659 f[1] = -y[0]*y[2];
660 f[2] = -0.51*y[0]*y[1];
661 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
662 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
663 PetscFunctionReturn(0);
664 }
665
IFunction_Hull1972B5(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
671
672 PetscFunctionBegin;
673 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
674 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
675 f[0] = y[1]*y[2];
676 f[1] = -y[0]*y[2];
677 f[2] = -0.51*y[0]*y[1];
678 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
679 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
680 /* Left hand side = ydot - f(y) */
681 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
682 PetscFunctionReturn(0);
683 }
684
IJacobian_Hull1972B5(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
691
692 PetscFunctionBegin;
693 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,3,&row[0],3,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
698 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
699 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
700 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
701 PetscFunctionReturn(0);
702 }
703
704
705 /* Kulikov, 2013, Problem I */
706
RHSFunction_Kulikov2013I(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
712
713 PetscFunctionBegin;
714 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
715 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
721 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
722 PetscFunctionReturn(0);
723 }
724
RHSJacobian_Kulikov2013I(TS ts,PetscReal t,Vec Y,Mat A,Mat B,void * s)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;
732 PetscFunctionBegin;
733 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
743 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
744 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
745 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
746 PetscFunctionReturn(0);
747 }
748
IFunction_Kulikov2013I(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
754
755 PetscFunctionBegin;
756 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
757 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
763 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
764 /* Left hand side = ydot - f(y) */
765 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
766 PetscFunctionReturn(0);
767 }
768
IJacobian_Kulikov2013I(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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;
776
777 PetscFunctionBegin;
778 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr);
788 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
789 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
790 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
791 PetscFunctionReturn(0);
792 }
793
794
795 /* Hull, 1972, Problem C1 */
796
RHSFunction_Hull1972C1(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
803
804 PetscFunctionBegin;
805 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
806 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
807 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
814 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
815 PetscFunctionReturn(0);
816 }
817
IFunction_Hull1972C1(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
824
825 PetscFunctionBegin;
826 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
827 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
828 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
835 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
836 /* Left hand side = ydot - f(y) */
837 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
838 PetscFunctionReturn(0);
839 }
840
IJacobian_Hull1972C1(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
847
848 PetscFunctionBegin;
849 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
850 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
851 i = 0;
852 value[0] = a+1; col[0] = 0;
853 value[1] = 0; col[1] = 1;
854 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
855 for (i = 0; i < N; i++) {
856 value[0] = -1; col[0] = i-1;
857 value[1] = a+1; col[1] = i;
858 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
859 }
860 i = N-1;
861 value[0] = -1; col[0] = N-2;
862 value[1] = a; col[1] = N-1;
863 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
864 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
865 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
866 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
867 PetscFunctionReturn(0);
868 }
869
870 /* Hull, 1972, Problem C2 */
871
RHSFunction_Hull1972C2(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
878
879 PetscFunctionBegin;
880 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
881 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
882 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
889 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
890 PetscFunctionReturn(0);
891 }
892
IFunction_Hull1972C2(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
899
900 PetscFunctionBegin;
901 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
902 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
903 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
910 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
911 /* Left hand side = ydot - f(y) */
912 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
913 PetscFunctionReturn(0);
914 }
915
IJacobian_Hull1972C2(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
922
923 PetscFunctionBegin;
924 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
925 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
926 i = 0;
927 value[0] = a+1; col[0] = 0;
928 value[1] = 0; col[1] = 1;
929 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
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 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
934 }
935 i = N-1;
936 value[0] = -(PetscReal) (N-1); col[0] = N-2;
937 value[1] = a; col[1] = N-1;
938 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
939 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
940 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
941 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
942 PetscFunctionReturn(0);
943 }
944
945 /* Hull, 1972, Problem C3 and C4 */
946
RHSFunction_Hull1972C34(TS ts,PetscReal t,Vec Y,Vec F,void * s)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;
953
954 PetscFunctionBegin;
955 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
956 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
957 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
964 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
965 PetscFunctionReturn(0);
966 }
967
IFunction_Hull1972C34(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,void * s)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;
974
975 PetscFunctionBegin;
976 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
977 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
978 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
985 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
986 /* Left hand side = ydot - f(y) */
987 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr);
988 PetscFunctionReturn(0);
989 }
990
IJacobian_Hull1972C34(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,void * s)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];
997
998 PetscFunctionBegin;
999 ierr = VecGetSize (Y,&N);CHKERRQ(ierr);
1000 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr);
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 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
1016 }
1017 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1018 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1019 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr);
1020 PetscFunctionReturn(0);
1021 }
1022
1023 /***************************************************************************/
1024
1025 /* Sets the initial solution for the IVP and sets up the function pointers*/
Initialize(Vec Y,void * s)1026 PetscErrorCode Initialize(Vec Y, void* s)
1027 {
1028 PetscErrorCode ierr;
1029 char *p = (char*) s;
1030 PetscScalar *y;
1031 PetscReal t0;
1032 PetscInt N = GetSize((const char *)s);
1033 PetscBool flg;
1034
1035 PetscFunctionBegin;
1036 VecZeroEntries(Y);
1037 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
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 ierr = PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);CHKERRQ(ierr);
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 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1132 PetscFunctionReturn(0);
1133 }
1134
1135 /* Calculates the exact solution to problems that have one */
ExactSolution(Vec Y,void * s,PetscReal t,PetscBool * flag)1136 PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag)
1137 {
1138 PetscErrorCode ierr;
1139 char *p = (char*) s;
1140 PetscScalar *y;
1141
1142 PetscFunctionBegin;
1143 if (!strcmp(p,"hull1972a1")) {
1144 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
1145 y[0] = PetscExpReal(-t);
1146 *flag = PETSC_TRUE;
1147 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1148 } else if (!strcmp(p,"hull1972a2")) {
1149 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
1150 y[0] = 1.0/PetscSqrtReal(t+1);
1151 *flag = PETSC_TRUE;
1152 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1153 } else if (!strcmp(p,"hull1972a3")) {
1154 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
1155 y[0] = PetscExpReal(PetscSinReal(t));
1156 *flag = PETSC_TRUE;
1157 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1158 } else if (!strcmp(p,"hull1972a4")) {
1159 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
1160 y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0));
1161 *flag = PETSC_TRUE;
1162 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1163 } else if (!strcmp(p,"kulik2013i")) {
1164 ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
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 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
1171 } else {
1172 ierr = VecSet(Y,0);CHKERRQ(ierr);
1173 *flag = PETSC_FALSE;
1174 }
1175 PetscFunctionReturn(0);
1176 }
1177
1178 /* Solves the specified ODE and computes the error if exact solution is available */
SolveODE(char * ptype,PetscReal dt,PetscReal tfinal,PetscInt maxiter,PetscReal * error,PetscBool * exact_flag)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 */
1191
1192 PetscFunctionBegin;
1193 N = GetSize((const char *)&ptype[0]);
1194 if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n");
1195 ierr = VecCreate(PETSC_COMM_WORLD,&Y);CHKERRQ(ierr);
1196 ierr = VecSetSizes(Y,N,PETSC_DECIDE);CHKERRQ(ierr);
1197 ierr = VecSetUp(Y);CHKERRQ(ierr);
1198 ierr = VecSet(Y,0);CHKERRQ(ierr);
1199
1200 /* Initialize the problem */
1201 ierr = Initialize(Y,&ptype[0]);CHKERRQ(ierr);
1202
1203 /* Create and initialize the time-integrator */
1204 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
1205 /* Default time integration options */
1206 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr);
1207 ierr = TSSetMaxSteps(ts,maxiter);CHKERRQ(ierr);
1208 ierr = TSSetMaxTime(ts,tfinal);CHKERRQ(ierr);
1209 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
1210 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
1211 /* Read command line options for time integration */
1212 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
1213 /* Set solution vector */
1214 ierr = TSSetSolution(ts,Y);CHKERRQ(ierr);
1215 /* Specify left/right-hand side functions */
1216 ierr = TSGetType(ts,&time_scheme);CHKERRQ(ierr);
1217
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 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);CHKERRQ(ierr);
1221 ierr = MatCreate(PETSC_COMM_WORLD,&Jac);CHKERRQ(ierr);
1222 ierr = MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
1223 ierr = MatSetFromOptions(Jac);CHKERRQ(ierr);
1224 ierr = MatSetUp(Jac);CHKERRQ(ierr);
1225 ierr = TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);CHKERRQ(ierr);
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 ierr = TSSetIFunction(ts,NULL,IFunction,&ptype[0]);CHKERRQ(ierr);
1230 ierr = MatCreate(PETSC_COMM_WORLD,&Jac);CHKERRQ(ierr);
1231 ierr = MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
1232 ierr = MatSetFromOptions(Jac);CHKERRQ(ierr);
1233 ierr = MatSetUp(Jac);CHKERRQ(ierr);
1234 ierr = TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);CHKERRQ(ierr);
1235 }
1236
1237 /* Solve */
1238 ierr = TSSolve(ts,Y);CHKERRQ(ierr);
1239 ierr = TSGetTime(ts,&final_time);CHKERRQ(ierr);
1240
1241 /* Get the estimated error, if available */
1242 ierr = VecDuplicate(Y,&Yerr);CHKERRQ(ierr);
1243 ierr = VecZeroEntries(Yerr);CHKERRQ(ierr);
1244 ierr = TSGetTimeError(ts,0,&Yerr);CHKERRQ(ierr);
1245 ierr = VecNorm(Yerr,NORM_2,&err_norm);CHKERRQ(ierr);
1246 ierr = VecDestroy(&Yerr);CHKERRQ(ierr);
1247 ierr = PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);CHKERRQ(ierr);
1248
1249 /* Exact solution */
1250 ierr = VecDuplicate(Y,&Yex);CHKERRQ(ierr);
1251 if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) {
1252 ierr = 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);CHKERRQ(ierr);
1253 }
1254 ierr = ExactSolution(Yex,&ptype[0],final_time,exact_flag);CHKERRQ(ierr);
1255
1256 /* Calculate Error */
1257 ierr = VecAYPX(Yex,-1.0,Y);CHKERRQ(ierr);
1258 ierr = VecNorm(Yex,NORM_2,error);CHKERRQ(ierr);
1259 *error = PetscSqrtReal(((*error)*(*error))/N);
1260
1261 /* Clean up and finalize */
1262 ierr = MatDestroy(&Jac);CHKERRQ(ierr);
1263 ierr = TSDestroy(&ts);CHKERRQ(ierr);
1264 ierr = VecDestroy(&Yex);CHKERRQ(ierr);
1265 ierr = VecDestroy(&Y);CHKERRQ(ierr);
1266
1267 PetscFunctionReturn(0);
1268 }
1269
main(int argc,char ** argv)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;
1284
1285 /* Initialize program */
1286 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
1287
1288 /* Check if running with only 1 proc */
1289 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
1290 if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
1291
1292 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);CHKERRQ(ierr);
1293 ierr = PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);CHKERRQ(ierr);
1294 ierr = PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);CHKERRQ(ierr);
1295 ierr = PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);CHKERRQ(ierr);
1296 ierr = PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);CHKERRQ(ierr);
1297 ierr = PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);CHKERRQ(ierr);
1298 ierr = PetscOptionsEnd();CHKERRQ(ierr);
1299
1300 ierr = PetscMalloc1(n_refine,&error);CHKERRQ(ierr);
1301 for (r = 0,dt = dt_initial; r < n_refine; r++) {
1302 error[r] = 0;
1303 if (r > 0) dt /= refine_fac;
1304
1305 ierr = 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]));CHKERRQ(ierr);
1306 ierr = SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);CHKERRQ(ierr);
1312 } else {
1313 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);CHKERRQ(ierr);
1314 }
1315 }
1316 }
1317 ierr = PetscFree(error);CHKERRQ(ierr);
1318 ierr = PetscFinalize();
1319 return ierr;
1320 }
1321
1322 /*TEST
1323
1324 test:
1325 suffix: 2
1326 args: -ts_type glee -final_time 5 -ts_adapt_type none
1327 timeoutfactor: 3
1328 requires: !single
1329
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
1335
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
1341
1342 TEST*/
1343