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