Actual source code: ex31.c

  1: static char help[] = "Solves the ordinary differential equations (IVPs) using explicit and implicit time-integration methods.\n";

  3: /*

  5:   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)

 14: */

 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"

 36:  From Constantinescu, E. "Estimating Global Errors in Time Stepping" ArXiv e-prints,
 37:        https://arxiv.org/abs/1503.05166, 2016

 39:     Kulikov2013I -> "kulik2013i"

 41: */

 43: #include <petscts.h>

 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 *);

 51: /* Returns the size of the system of equations depending on problem specification */
 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: }

 65: /****************************************************************/

 67: /* Problem specific functions */

 69: /* Hull, 1972, Problem A1 */

 71: PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
 72: {
 73:   PetscScalar       *f;
 74:   const PetscScalar *y;

 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: }

 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;

 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: }

100: PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
101: {
102:   const PetscScalar *y;
103:   PetscScalar       *f;

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: }

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;

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: }

131: /* Hull, 1972, Problem A2 */

133: PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
134: {
135:   const PetscScalar *y;
136:   PetscScalar       *f;

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: }

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;

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: }

163: PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
164: {
165:   PetscScalar       *f;
166:   const PetscScalar *y;

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: }

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;

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: }

195: /* Hull, 1972, Problem A3 */

197: PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
198: {
199:   const PetscScalar *y;
200:   PetscScalar       *f;

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: }

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);

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: }

226: PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
227: {
228:   PetscScalar       *f;
229:   const PetscScalar *y;

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: }

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);

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: }

257: /* Hull, 1972, Problem A4 */

259: PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
260: {
261:   PetscScalar       *f;
262:   const PetscScalar *y;

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: }

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;

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: }

289: PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
290: {
291:   PetscScalar       *f;
292:   const PetscScalar *y;

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: }

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;

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: }

321: /* Hull, 1972, Problem A5 */

323: PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
324: {
325:   PetscScalar       *f;
326:   const PetscScalar *y;

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: }

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;

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: }

353: PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
354: {
355:   PetscScalar       *f;
356:   const PetscScalar *y;

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: }

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;

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: }

385: /* Hull, 1972, Problem B1 */

387: PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s)
388: {
389:   PetscScalar       *f;
390:   const PetscScalar *y;

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: }

402: PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
403: {
404:   PetscScalar       *f;
405:   const PetscScalar *y;

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: }

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];

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: }

438: /* Hull, 1972, Problem B2 */

440: PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s)
441: {
442:   PetscScalar       *f;
443:   const PetscScalar *y;

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: }

456: PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
457: {
458:   PetscScalar       *f;
459:   const PetscScalar *y;

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: }

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];

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: }

498: /* Hull, 1972, Problem B3 */

500: PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s)
501: {
502:   PetscScalar       *f;
503:   const PetscScalar *y;

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: }

516: PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
517: {
518:   PetscScalar       *f;
519:   const PetscScalar *y;

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: }

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];

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: }

558: /* Hull, 1972, Problem B4 */

560: PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s)
561: {
562:   PetscScalar       *f;
563:   const PetscScalar *y;

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: }

576: PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
577: {
578:   PetscScalar       *f;
579:   const PetscScalar *y;

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: }

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;

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: }

620: /* Hull, 1972, Problem B5 */

622: PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s)
623: {
624:   PetscScalar       *f;
625:   const PetscScalar *y;

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: }

638: PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
639: {
640:   PetscScalar       *f;
641:   const PetscScalar *y;

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: }

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];

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: }

680: /* Kulikov, 2013, Problem I */

682: PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s)
683: {
684:   PetscScalar       *f;
685:   const PetscScalar *y;

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: }

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;

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: }

735: PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s)
736: {
737:   PetscScalar       *f;
738:   const PetscScalar *y;

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: }

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;

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: }

790: /* Hull, 1972, Problem C1 */

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;

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: }

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;

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: }

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];

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: }

864: /* Hull, 1972, Problem C2 */

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;

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: }

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;

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: }

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];

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: }

938: /* Hull, 1972, Problem C3 and C4 */

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;

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: }

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;

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: }

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];

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: }

1018: /***************************************************************************/

1020: /* Sets the initial solution for the IVP and sets up the function pointers*/
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;

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: }

1128: /* Calculates the exact solution to problems that have one */
1129: PetscErrorCode ExactSolution(Vec Y, void *s, PetscReal t, PetscBool *flag)
1130: {
1131:   char        *p = (char *)s;
1132:   PetscScalar *y;

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: }

1170: /* Solves the specified ODE and computes the error if exact solution is available */
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  */

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));

1191:   /* Initialize the problem */
1192:   PetscCall(Initialize(Y, &ptype[0]));

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));

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:   }

1228:   /* Solve */
1229:   PetscCall(TSSolve(ts, Y));
1230:   PetscCall(TSGetTime(ts, &final_time));

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));

1240:   /* Exact solution */
1241:   PetscCall(VecDuplicate(Y, &Yex));
1242:   if (PetscAbsReal(final_time - tfinal) > 2. * PETSC_MACHINE_EPSILON) {
1243:     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));
1244:   }
1245:   PetscCall(ExactSolution(Yex, &ptype[0], final_time, exact_flag));

1247:   /* Calculate Error */
1248:   PetscCall(VecAYPX(Yex, -1.0, Y));
1249:   PetscCall(VecNorm(Yex, NORM_2, error));
1250:   *error = PetscSqrtReal(((*error) * (*error)) / N);

1252:   /* Clean up and finalize */
1253:   PetscCall(MatDestroy(&Jac));
1254:   PetscCall(TSDestroy(&ts));
1255:   PetscCall(VecDestroy(&Yex));
1256:   PetscCall(VecDestroy(&Y));
1257:   PetscFunctionReturn(PETSC_SUCCESS);
1258: }

1260: int main(int argc, char **argv)
1261: {
1262:   char        ptype[256] = "hull1972a1"; /* Problem specification                                */
1263:   PetscInt    n_refine   = 1;            /* Number of refinement levels for convergence analysis */
1264:   PetscReal   refine_fac = 2.0;          /* Refinement factor for dt                             */
1265:   PetscReal   dt_initial = 0.01;         /* Initial default value of dt                          */
1266:   PetscReal   dt;
1267:   PetscReal   tfinal  = 20.0;   /* Final time for the time-integration                  */
1268:   PetscInt    maxiter = 100000; /* Maximum number of time-integration iterations        */
1269:   PetscReal  *error;            /* Array to store the errors for convergence analysis   */
1270:   PetscMPIInt size;             /* No of processors                                     */
1271:   PetscBool   flag;             /* Flag denoting availability of exact solution         */
1272:   PetscInt    r, N;

1274:   /* Initialize program */
1275:   PetscFunctionBeginUser;
1276:   PetscCall(GetSize(&ptype[0], &N));
1277:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

1279:   /* Check if running with only 1 proc */
1280:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1281:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

1283:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL);
1284:   PetscCall(PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL));
1285:   PetscCall(PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL));
1286:   PetscCall(PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL));
1287:   PetscCall(PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL));
1288:   PetscCall(PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL));
1289:   PetscOptionsEnd();

1291:   PetscCall(PetscMalloc1(n_refine, &error));
1292:   for (r = 0, dt = dt_initial; r < n_refine; r++) {
1293:     error[r] = 0;
1294:     if (r > 0) dt /= refine_fac;

1296:     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));
1297:     PetscCall(SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag));
1298:     if (flag) {
1299:       /* If exact solution available for the specified ODE */
1300:       if (r > 0) {
1301:         PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac));
1302:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error           = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate));
1303:       } else {
1304:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error           = %E.\n", (double)error[r]));
1305:       }
1306:     }
1307:   }
1308:   PetscCall(PetscFree(error));
1309:   PetscCall(PetscFinalize());
1310:   return 0;
1311: }

1313: /*TEST

1315:     test:
1316:       suffix: 2
1317:       args: -ts_type glee -final_time 5 -ts_adapt_type none
1318:       timeoutfactor: 3
1319:       requires: !single

1321:     test:
1322:       suffix: 3
1323:       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
1324:       timeoutfactor: 3
1325:       requires: !single

1327:     test:
1328:       suffix: 4
1329:       # the test is so sensitive that I could not even replace VecMAXPY with a sequence of VECAXPY, so I just disable this GEMV optimization
1330:       args: -vec_maxpy_use_gemv 0
1331:       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
1332:       timeoutfactor: 3
1333:       requires: !single !__float128

1335: TEST*/