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