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