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