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 56 if (!strcmp(p, "hull1972a1") || !strcmp(p, "hull1972a2") || !strcmp(p, "hull1972a3") || !strcmp(p, "hull1972a4") || !strcmp(p, "hull1972a5")) *sz = 1; 57 else if (!strcmp(p, "hull1972b1")) *sz = 2; 58 else if (!strcmp(p, "hull1972b2") || !strcmp(p, "hull1972b3") || !strcmp(p, "hull1972b4") || !strcmp(p, "hull1972b5")) *sz = 3; 59 else if (!strcmp(p, "kulik2013i")) *sz = 4; 60 else if (!strcmp(p, "hull1972c1") || !strcmp(p, "hull1972c2") || !strcmp(p, "hull1972c3")) *sz = 10; 61 else if (!strcmp(p, "hull1972c4")) *sz = 52; 62 else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Incorrect problem name %s", p); 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 /****************************************************************/ 67 68 /* Problem specific functions */ 69 70 /* Hull, 1972, Problem A1 */ 71 72 PetscErrorCode RHSFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec F, void *s) 73 { 74 PetscScalar *f; 75 const PetscScalar *y; 76 77 PetscFunctionBeginUser; 78 PetscCall(VecGetArrayRead(Y, &y)); 79 PetscCall(VecGetArray(F, &f)); 80 f[0] = -y[0]; 81 PetscCall(VecRestoreArrayRead(Y, &y)); 82 PetscCall(VecRestoreArray(F, &f)); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 PetscErrorCode RHSJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 87 { 88 const PetscScalar *y; 89 PetscInt row = 0, col = 0; 90 PetscScalar value = -1.0; 91 92 PetscFunctionBeginUser; 93 PetscCall(VecGetArrayRead(Y, &y)); 94 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 95 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 96 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 97 PetscCall(VecRestoreArrayRead(Y, &y)); 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 PetscErrorCode IFunction_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 102 { 103 const PetscScalar *y; 104 PetscScalar *f; 105 106 PetscFunctionBeginUser; 107 PetscCall(VecGetArrayRead(Y, &y)); 108 PetscCall(VecGetArray(F, &f)); 109 f[0] = -y[0]; 110 PetscCall(VecRestoreArrayRead(Y, &y)); 111 PetscCall(VecRestoreArray(F, &f)); 112 /* Left hand side = ydot - f(y) */ 113 PetscCall(VecAYPX(F, -1.0, Ydot)); 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 PetscErrorCode IJacobian_Hull1972A1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 118 { 119 const PetscScalar *y; 120 PetscInt row = 0, col = 0; 121 PetscScalar value = a - 1.0; 122 123 PetscFunctionBeginUser; 124 PetscCall(VecGetArrayRead(Y, &y)); 125 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 126 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 128 PetscCall(VecRestoreArrayRead(Y, &y)); 129 PetscFunctionReturn(PETSC_SUCCESS); 130 } 131 132 /* Hull, 1972, Problem A2 */ 133 134 PetscErrorCode RHSFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec F, void *s) 135 { 136 const PetscScalar *y; 137 PetscScalar *f; 138 139 PetscFunctionBeginUser; 140 PetscCall(VecGetArrayRead(Y, &y)); 141 PetscCall(VecGetArray(F, &f)); 142 f[0] = -0.5 * y[0] * y[0] * y[0]; 143 PetscCall(VecRestoreArrayRead(Y, &y)); 144 PetscCall(VecRestoreArray(F, &f)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 PetscErrorCode RHSJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 149 { 150 const PetscScalar *y; 151 PetscInt row = 0, col = 0; 152 PetscScalar value; 153 154 PetscFunctionBeginUser; 155 PetscCall(VecGetArrayRead(Y, &y)); 156 value = -0.5 * 3.0 * y[0] * y[0]; 157 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 158 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 159 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 160 PetscCall(VecRestoreArrayRead(Y, &y)); 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 PetscErrorCode IFunction_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 165 { 166 PetscScalar *f; 167 const PetscScalar *y; 168 169 PetscFunctionBeginUser; 170 PetscCall(VecGetArrayRead(Y, &y)); 171 PetscCall(VecGetArray(F, &f)); 172 f[0] = -0.5 * y[0] * y[0] * y[0]; 173 PetscCall(VecRestoreArrayRead(Y, &y)); 174 PetscCall(VecRestoreArray(F, &f)); 175 /* Left hand side = ydot - f(y) */ 176 PetscCall(VecAYPX(F, -1.0, Ydot)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 PetscErrorCode IJacobian_Hull1972A2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 181 { 182 const PetscScalar *y; 183 PetscInt row = 0, col = 0; 184 PetscScalar value; 185 186 PetscFunctionBeginUser; 187 PetscCall(VecGetArrayRead(Y, &y)); 188 value = a + 0.5 * 3.0 * y[0] * y[0]; 189 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 190 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 191 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 192 PetscCall(VecRestoreArrayRead(Y, &y)); 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 /* Hull, 1972, Problem A3 */ 197 198 PetscErrorCode RHSFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec F, void *s) 199 { 200 const PetscScalar *y; 201 PetscScalar *f; 202 203 PetscFunctionBeginUser; 204 PetscCall(VecGetArrayRead(Y, &y)); 205 PetscCall(VecGetArray(F, &f)); 206 f[0] = y[0] * PetscCosReal(t); 207 PetscCall(VecRestoreArrayRead(Y, &y)); 208 PetscCall(VecRestoreArray(F, &f)); 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 PetscErrorCode RHSJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 213 { 214 const PetscScalar *y; 215 PetscInt row = 0, col = 0; 216 PetscScalar value = PetscCosReal(t); 217 218 PetscFunctionBeginUser; 219 PetscCall(VecGetArrayRead(Y, &y)); 220 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 221 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 222 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 223 PetscCall(VecRestoreArrayRead(Y, &y)); 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 PetscErrorCode IFunction_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 228 { 229 PetscScalar *f; 230 const PetscScalar *y; 231 232 PetscFunctionBeginUser; 233 PetscCall(VecGetArrayRead(Y, &y)); 234 PetscCall(VecGetArray(F, &f)); 235 f[0] = y[0] * PetscCosReal(t); 236 PetscCall(VecRestoreArrayRead(Y, &y)); 237 PetscCall(VecRestoreArray(F, &f)); 238 /* Left hand side = ydot - f(y) */ 239 PetscCall(VecAYPX(F, -1.0, Ydot)); 240 PetscFunctionReturn(PETSC_SUCCESS); 241 } 242 243 PetscErrorCode IJacobian_Hull1972A3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 244 { 245 const PetscScalar *y; 246 PetscInt row = 0, col = 0; 247 PetscScalar value = a - PetscCosReal(t); 248 249 PetscFunctionBeginUser; 250 PetscCall(VecGetArrayRead(Y, &y)); 251 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 252 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 253 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 254 PetscCall(VecRestoreArrayRead(Y, &y)); 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 /* Hull, 1972, Problem A4 */ 259 260 PetscErrorCode RHSFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec F, void *s) 261 { 262 PetscScalar *f; 263 const PetscScalar *y; 264 265 PetscFunctionBeginUser; 266 PetscCall(VecGetArrayRead(Y, &y)); 267 PetscCall(VecGetArray(F, &f)); 268 f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]); 269 PetscCall(VecRestoreArrayRead(Y, &y)); 270 PetscCall(VecRestoreArray(F, &f)); 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 PetscErrorCode RHSJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 275 { 276 const PetscScalar *y; 277 PetscInt row = 0, col = 0; 278 PetscScalar value; 279 280 PetscFunctionBeginUser; 281 PetscCall(VecGetArrayRead(Y, &y)); 282 value = 0.25 * (1.0 - 0.05 * y[0]) - (0.25 * y[0]) * 0.05; 283 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 284 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 285 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 286 PetscCall(VecRestoreArrayRead(Y, &y)); 287 PetscFunctionReturn(PETSC_SUCCESS); 288 } 289 290 PetscErrorCode IFunction_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 291 { 292 PetscScalar *f; 293 const PetscScalar *y; 294 295 PetscFunctionBeginUser; 296 PetscCall(VecGetArrayRead(Y, &y)); 297 PetscCall(VecGetArray(F, &f)); 298 f[0] = (0.25 * y[0]) * (1.0 - 0.05 * y[0]); 299 PetscCall(VecRestoreArrayRead(Y, &y)); 300 PetscCall(VecRestoreArray(F, &f)); 301 /* Left hand side = ydot - f(y) */ 302 PetscCall(VecAYPX(F, -1.0, Ydot)); 303 PetscFunctionReturn(PETSC_SUCCESS); 304 } 305 306 PetscErrorCode IJacobian_Hull1972A4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 307 { 308 const PetscScalar *y; 309 PetscInt row = 0, col = 0; 310 PetscScalar value; 311 312 PetscFunctionBeginUser; 313 PetscCall(VecGetArrayRead(Y, &y)); 314 value = a - 0.25 * (1.0 - 0.05 * y[0]) + (0.25 * y[0]) * 0.05; 315 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 316 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 317 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 318 PetscCall(VecRestoreArrayRead(Y, &y)); 319 PetscFunctionReturn(PETSC_SUCCESS); 320 } 321 322 /* Hull, 1972, Problem A5 */ 323 324 PetscErrorCode RHSFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec F, void *s) 325 { 326 PetscScalar *f; 327 const PetscScalar *y; 328 329 PetscFunctionBeginUser; 330 PetscCall(VecGetArrayRead(Y, &y)); 331 PetscCall(VecGetArray(F, &f)); 332 f[0] = (y[0] - t) / (y[0] + t); 333 PetscCall(VecRestoreArrayRead(Y, &y)); 334 PetscCall(VecRestoreArray(F, &f)); 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337 338 PetscErrorCode RHSJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 339 { 340 const PetscScalar *y; 341 PetscInt row = 0, col = 0; 342 PetscScalar value; 343 344 PetscFunctionBeginUser; 345 PetscCall(VecGetArrayRead(Y, &y)); 346 value = 2 * t / ((t + y[0]) * (t + y[0])); 347 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 348 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 349 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 350 PetscCall(VecRestoreArrayRead(Y, &y)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 PetscErrorCode IFunction_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 355 { 356 PetscScalar *f; 357 const PetscScalar *y; 358 359 PetscFunctionBeginUser; 360 PetscCall(VecGetArrayRead(Y, &y)); 361 PetscCall(VecGetArray(F, &f)); 362 f[0] = (y[0] - t) / (y[0] + t); 363 PetscCall(VecRestoreArrayRead(Y, &y)); 364 PetscCall(VecRestoreArray(F, &f)); 365 /* Left hand side = ydot - f(y) */ 366 PetscCall(VecAYPX(F, -1.0, Ydot)); 367 PetscFunctionReturn(PETSC_SUCCESS); 368 } 369 370 PetscErrorCode IJacobian_Hull1972A5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 371 { 372 const PetscScalar *y; 373 PetscInt row = 0, col = 0; 374 PetscScalar value; 375 376 PetscFunctionBeginUser; 377 PetscCall(VecGetArrayRead(Y, &y)); 378 value = a - 2 * t / ((t + y[0]) * (t + y[0])); 379 PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES)); 380 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 381 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 382 PetscCall(VecRestoreArrayRead(Y, &y)); 383 PetscFunctionReturn(PETSC_SUCCESS); 384 } 385 386 /* Hull, 1972, Problem B1 */ 387 388 PetscErrorCode RHSFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec F, void *s) 389 { 390 PetscScalar *f; 391 const PetscScalar *y; 392 393 PetscFunctionBeginUser; 394 PetscCall(VecGetArrayRead(Y, &y)); 395 PetscCall(VecGetArray(F, &f)); 396 f[0] = 2.0 * (y[0] - y[0] * y[1]); 397 f[1] = -(y[1] - y[0] * y[1]); 398 PetscCall(VecRestoreArrayRead(Y, &y)); 399 PetscCall(VecRestoreArray(F, &f)); 400 PetscFunctionReturn(PETSC_SUCCESS); 401 } 402 403 PetscErrorCode IFunction_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 404 { 405 PetscScalar *f; 406 const PetscScalar *y; 407 408 PetscFunctionBeginUser; 409 PetscCall(VecGetArrayRead(Y, &y)); 410 PetscCall(VecGetArray(F, &f)); 411 f[0] = 2.0 * (y[0] - y[0] * y[1]); 412 f[1] = -(y[1] - y[0] * y[1]); 413 PetscCall(VecRestoreArrayRead(Y, &y)); 414 PetscCall(VecRestoreArray(F, &f)); 415 /* Left hand side = ydot - f(y) */ 416 PetscCall(VecAYPX(F, -1.0, Ydot)); 417 PetscFunctionReturn(PETSC_SUCCESS); 418 } 419 420 PetscErrorCode IJacobian_Hull1972B1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 421 { 422 const PetscScalar *y; 423 PetscInt row[2] = {0, 1}; 424 PetscScalar value[2][2]; 425 426 PetscFunctionBeginUser; 427 PetscCall(VecGetArrayRead(Y, &y)); 428 value[0][0] = a - 2.0 * (1.0 - y[1]); 429 value[0][1] = 2.0 * y[0]; 430 value[1][0] = -y[1]; 431 value[1][1] = a + 1.0 - y[0]; 432 PetscCall(MatSetValues(A, 2, &row[0], 2, &row[0], &value[0][0], INSERT_VALUES)); 433 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 434 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 435 PetscCall(VecRestoreArrayRead(Y, &y)); 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 /* Hull, 1972, Problem B2 */ 440 441 PetscErrorCode RHSFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec F, void *s) 442 { 443 PetscScalar *f; 444 const PetscScalar *y; 445 446 PetscFunctionBeginUser; 447 PetscCall(VecGetArrayRead(Y, &y)); 448 PetscCall(VecGetArray(F, &f)); 449 f[0] = -y[0] + y[1]; 450 f[1] = y[0] - 2.0 * y[1] + y[2]; 451 f[2] = y[1] - y[2]; 452 PetscCall(VecRestoreArrayRead(Y, &y)); 453 PetscCall(VecRestoreArray(F, &f)); 454 PetscFunctionReturn(PETSC_SUCCESS); 455 } 456 457 PetscErrorCode IFunction_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 458 { 459 PetscScalar *f; 460 const PetscScalar *y; 461 462 PetscFunctionBeginUser; 463 PetscCall(VecGetArrayRead(Y, &y)); 464 PetscCall(VecGetArray(F, &f)); 465 f[0] = -y[0] + y[1]; 466 f[1] = y[0] - 2.0 * y[1] + y[2]; 467 f[2] = y[1] - y[2]; 468 PetscCall(VecRestoreArrayRead(Y, &y)); 469 PetscCall(VecRestoreArray(F, &f)); 470 /* Left hand side = ydot - f(y) */ 471 PetscCall(VecAYPX(F, -1.0, Ydot)); 472 PetscFunctionReturn(PETSC_SUCCESS); 473 } 474 475 PetscErrorCode IJacobian_Hull1972B2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 476 { 477 const PetscScalar *y; 478 PetscInt row[3] = {0, 1, 2}; 479 PetscScalar value[3][3]; 480 481 PetscFunctionBeginUser; 482 PetscCall(VecGetArrayRead(Y, &y)); 483 value[0][0] = a + 1.0; 484 value[0][1] = -1.0; 485 value[0][2] = 0; 486 value[1][0] = -1.0; 487 value[1][1] = a + 2.0; 488 value[1][2] = -1.0; 489 value[2][0] = 0; 490 value[2][1] = -1.0; 491 value[2][2] = a + 1.0; 492 PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES)); 493 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 494 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 495 PetscCall(VecRestoreArrayRead(Y, &y)); 496 PetscFunctionReturn(PETSC_SUCCESS); 497 } 498 499 /* Hull, 1972, Problem B3 */ 500 501 PetscErrorCode RHSFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec F, void *s) 502 { 503 PetscScalar *f; 504 const PetscScalar *y; 505 506 PetscFunctionBeginUser; 507 PetscCall(VecGetArrayRead(Y, &y)); 508 PetscCall(VecGetArray(F, &f)); 509 f[0] = -y[0]; 510 f[1] = y[0] - y[1] * y[1]; 511 f[2] = y[1] * y[1]; 512 PetscCall(VecRestoreArrayRead(Y, &y)); 513 PetscCall(VecRestoreArray(F, &f)); 514 PetscFunctionReturn(PETSC_SUCCESS); 515 } 516 517 PetscErrorCode IFunction_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 518 { 519 PetscScalar *f; 520 const PetscScalar *y; 521 522 PetscFunctionBeginUser; 523 PetscCall(VecGetArrayRead(Y, &y)); 524 PetscCall(VecGetArray(F, &f)); 525 f[0] = -y[0]; 526 f[1] = y[0] - y[1] * y[1]; 527 f[2] = y[1] * y[1]; 528 PetscCall(VecRestoreArrayRead(Y, &y)); 529 PetscCall(VecRestoreArray(F, &f)); 530 /* Left hand side = ydot - f(y) */ 531 PetscCall(VecAYPX(F, -1.0, Ydot)); 532 PetscFunctionReturn(PETSC_SUCCESS); 533 } 534 535 PetscErrorCode IJacobian_Hull1972B3(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 536 { 537 const PetscScalar *y; 538 PetscInt row[3] = {0, 1, 2}; 539 PetscScalar value[3][3]; 540 541 PetscFunctionBeginUser; 542 PetscCall(VecGetArrayRead(Y, &y)); 543 value[0][0] = a + 1.0; 544 value[0][1] = 0; 545 value[0][2] = 0; 546 value[1][0] = -1.0; 547 value[1][1] = a + 2.0 * y[1]; 548 value[1][2] = 0; 549 value[2][0] = 0; 550 value[2][1] = -2.0 * y[1]; 551 value[2][2] = a; 552 PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES)); 553 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 554 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 555 PetscCall(VecRestoreArrayRead(Y, &y)); 556 PetscFunctionReturn(PETSC_SUCCESS); 557 } 558 559 /* Hull, 1972, Problem B4 */ 560 561 PetscErrorCode RHSFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec F, void *s) 562 { 563 PetscScalar *f; 564 const PetscScalar *y; 565 566 PetscFunctionBeginUser; 567 PetscCall(VecGetArrayRead(Y, &y)); 568 PetscCall(VecGetArray(F, &f)); 569 f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 570 f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 571 f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 572 PetscCall(VecRestoreArrayRead(Y, &y)); 573 PetscCall(VecRestoreArray(F, &f)); 574 PetscFunctionReturn(PETSC_SUCCESS); 575 } 576 577 PetscErrorCode IFunction_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 578 { 579 PetscScalar *f; 580 const PetscScalar *y; 581 582 PetscFunctionBeginUser; 583 PetscCall(VecGetArrayRead(Y, &y)); 584 PetscCall(VecGetArray(F, &f)); 585 f[0] = -y[1] - y[0] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 586 f[1] = y[0] - y[1] * y[2] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 587 f[2] = y[0] / PetscSqrtScalar(y[0] * y[0] + y[1] * y[1]); 588 PetscCall(VecRestoreArrayRead(Y, &y)); 589 PetscCall(VecRestoreArray(F, &f)); 590 /* Left hand side = ydot - f(y) */ 591 PetscCall(VecAYPX(F, -1.0, Ydot)); 592 PetscFunctionReturn(PETSC_SUCCESS); 593 } 594 595 PetscErrorCode IJacobian_Hull1972B4(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 596 { 597 const PetscScalar *y; 598 PetscInt row[3] = {0, 1, 2}; 599 PetscScalar value[3][3], fac, fac2; 600 601 PetscFunctionBeginUser; 602 PetscCall(VecGetArrayRead(Y, &y)); 603 fac = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -1.5); 604 fac2 = PetscPowScalar(y[0] * y[0] + y[1] * y[1], -0.5); 605 value[0][0] = a + (y[1] * y[1] * y[2]) * fac; 606 value[0][1] = 1.0 - (y[0] * y[1] * y[2]) * fac; 607 value[0][2] = y[0] * fac2; 608 value[1][0] = -1.0 - y[0] * y[1] * y[2] * fac; 609 value[1][1] = a + y[0] * y[0] * y[2] * fac; 610 value[1][2] = y[1] * fac2; 611 value[2][0] = -y[1] * y[1] * fac; 612 value[2][1] = y[0] * y[1] * fac; 613 value[2][2] = a; 614 PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES)); 615 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 616 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 617 PetscCall(VecRestoreArrayRead(Y, &y)); 618 PetscFunctionReturn(PETSC_SUCCESS); 619 } 620 621 /* Hull, 1972, Problem B5 */ 622 623 PetscErrorCode RHSFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec F, void *s) 624 { 625 PetscScalar *f; 626 const PetscScalar *y; 627 628 PetscFunctionBeginUser; 629 PetscCall(VecGetArrayRead(Y, &y)); 630 PetscCall(VecGetArray(F, &f)); 631 f[0] = y[1] * y[2]; 632 f[1] = -y[0] * y[2]; 633 f[2] = -0.51 * y[0] * y[1]; 634 PetscCall(VecRestoreArrayRead(Y, &y)); 635 PetscCall(VecRestoreArray(F, &f)); 636 PetscFunctionReturn(PETSC_SUCCESS); 637 } 638 639 PetscErrorCode IFunction_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 640 { 641 PetscScalar *f; 642 const PetscScalar *y; 643 644 PetscFunctionBeginUser; 645 PetscCall(VecGetArrayRead(Y, &y)); 646 PetscCall(VecGetArray(F, &f)); 647 f[0] = y[1] * y[2]; 648 f[1] = -y[0] * y[2]; 649 f[2] = -0.51 * y[0] * y[1]; 650 PetscCall(VecRestoreArrayRead(Y, &y)); 651 PetscCall(VecRestoreArray(F, &f)); 652 /* Left hand side = ydot - f(y) */ 653 PetscCall(VecAYPX(F, -1.0, Ydot)); 654 PetscFunctionReturn(PETSC_SUCCESS); 655 } 656 657 PetscErrorCode IJacobian_Hull1972B5(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 658 { 659 const PetscScalar *y; 660 PetscInt row[3] = {0, 1, 2}; 661 PetscScalar value[3][3]; 662 663 PetscFunctionBeginUser; 664 PetscCall(VecGetArrayRead(Y, &y)); 665 value[0][0] = a; 666 value[0][1] = -y[2]; 667 value[0][2] = -y[1]; 668 value[1][0] = y[2]; 669 value[1][1] = a; 670 value[1][2] = y[0]; 671 value[2][0] = 0.51 * y[1]; 672 value[2][1] = 0.51 * y[0]; 673 value[2][2] = a; 674 PetscCall(MatSetValues(A, 3, &row[0], 3, &row[0], &value[0][0], INSERT_VALUES)); 675 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 676 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 677 PetscCall(VecRestoreArrayRead(Y, &y)); 678 PetscFunctionReturn(PETSC_SUCCESS); 679 } 680 681 /* Kulikov, 2013, Problem I */ 682 683 PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s) 684 { 685 PetscScalar *f; 686 const PetscScalar *y; 687 688 PetscFunctionBeginUser; 689 PetscCall(VecGetArrayRead(Y, &y)); 690 PetscCall(VecGetArray(F, &f)); 691 f[0] = 2. * t * PetscPowScalar(y[1], 1. / 5.) * y[3]; 692 f[1] = 10. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.)); 693 f[2] = 2. * t * y[3]; 694 f[3] = -2. * t * PetscLogScalar(y[0]); 695 PetscCall(VecRestoreArrayRead(Y, &y)); 696 PetscCall(VecRestoreArray(F, &f)); 697 PetscFunctionReturn(PETSC_SUCCESS); 698 } 699 700 PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 701 { 702 const PetscScalar *y; 703 PetscInt row[4] = {0, 1, 2, 3}; 704 PetscScalar value[4][4]; 705 PetscScalar m1, m2; 706 PetscFunctionBeginUser; 707 PetscCall(VecGetArrayRead(Y, &y)); 708 m1 = (2. * t * y[3]) / (5. * PetscPowScalar(y[1], 4. / 5.)); 709 m2 = 2. * t * PetscPowScalar(y[1], 1. / 5.); 710 value[0][0] = 0.; 711 value[0][1] = m1; 712 value[0][2] = 0.; 713 value[0][3] = m2; 714 m1 = 50. * t * y[3] * PetscExpScalar(5.0 * (y[2] - 1.)); 715 m2 = 10. * t * PetscExpScalar(5.0 * (y[2] - 1.)); 716 value[1][0] = 0.; 717 value[1][1] = 0.; 718 value[1][2] = m1; 719 value[1][3] = m2; 720 value[2][0] = 0.; 721 value[2][1] = 0.; 722 value[2][2] = 0.; 723 value[2][3] = 2 * t; 724 value[3][0] = -2. * t / y[0]; 725 value[3][1] = 0.; 726 value[3][2] = 0.; 727 value[3][3] = 0.; 728 PetscCall(MatSetValues(A, 4, &row[0], 4, &row[0], &value[0][0], INSERT_VALUES)); 729 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 730 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 731 PetscCall(VecRestoreArrayRead(Y, &y)); 732 PetscFunctionReturn(PETSC_SUCCESS); 733 } 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) { 1243 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference between the prescribed final time %g and the actual final time, %g.\n", (double)tfinal, (double)final_time)); 1244 } 1245 PetscCall(ExactSolution(Yex, &ptype[0], final_time, exact_flag)); 1246 1247 /* Calculate Error */ 1248 PetscCall(VecAYPX(Yex, -1.0, Y)); 1249 PetscCall(VecNorm(Yex, NORM_2, error)); 1250 *error = PetscSqrtReal(((*error) * (*error)) / N); 1251 1252 /* Clean up and finalize */ 1253 PetscCall(MatDestroy(&Jac)); 1254 PetscCall(TSDestroy(&ts)); 1255 PetscCall(VecDestroy(&Yex)); 1256 PetscCall(VecDestroy(&Y)); 1257 1258 PetscFunctionReturn(PETSC_SUCCESS); 1259 } 1260 1261 int main(int argc, char **argv) 1262 { 1263 char ptype[256] = "hull1972a1"; /* Problem specification */ 1264 PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */ 1265 PetscReal refine_fac = 2.0; /* Refinement factor for dt */ 1266 PetscReal dt_initial = 0.01; /* Initial default value of dt */ 1267 PetscReal dt; 1268 PetscReal tfinal = 20.0; /* Final time for the time-integration */ 1269 PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */ 1270 PetscReal *error; /* Array to store the errors for convergence analysis */ 1271 PetscMPIInt size; /* No of processors */ 1272 PetscBool flag; /* Flag denoting availability of exact solution */ 1273 PetscInt r, N; 1274 1275 /* Initialize program */ 1276 PetscFunctionBeginUser; 1277 PetscCall(GetSize(&ptype[0], &N)); 1278 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 1279 1280 /* Check if running with only 1 proc */ 1281 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 1282 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 1283 1284 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex31", NULL); 1285 PetscCall(PetscOptionsString("-problem", "Problem specification", "<hull1972a1>", ptype, ptype, sizeof(ptype), NULL)); 1286 PetscCall(PetscOptionsInt("-refinement_levels", "Number of refinement levels for convergence analysis", "<1>", n_refine, &n_refine, NULL)); 1287 PetscCall(PetscOptionsReal("-refinement_factor", "Refinement factor for dt", "<2.0>", refine_fac, &refine_fac, NULL)); 1288 PetscCall(PetscOptionsReal("-dt", "Time step size (for convergence analysis, initial time step)", "<0.01>", dt_initial, &dt_initial, NULL)); 1289 PetscCall(PetscOptionsReal("-final_time", "Final time for the time-integration", "<20.0>", tfinal, &tfinal, NULL)); 1290 PetscOptionsEnd(); 1291 1292 PetscCall(PetscMalloc1(n_refine, &error)); 1293 for (r = 0, dt = dt_initial; r < n_refine; r++) { 1294 error[r] = 0; 1295 if (r > 0) dt /= refine_fac; 1296 1297 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)); 1298 PetscCall(SolveODE(&ptype[0], dt, tfinal, maxiter, &error[r], &flag)); 1299 if (flag) { 1300 /* If exact solution available for the specified ODE */ 1301 if (r > 0) { 1302 PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r - 1])) / (-PetscLogReal(refine_fac)); 1303 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E,\tConvergence rate = %f.\n", (double)error[r], (double)conv_rate)); 1304 } else { 1305 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error = %E.\n", (double)error[r])); 1306 } 1307 } 1308 } 1309 PetscCall(PetscFree(error)); 1310 PetscCall(PetscFinalize()); 1311 return 0; 1312 } 1313 1314 /*TEST 1315 1316 test: 1317 suffix: 2 1318 args: -ts_type glee -final_time 5 -ts_adapt_type none 1319 timeoutfactor: 3 1320 requires: !single 1321 1322 test: 1323 suffix: 3 1324 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 1325 timeoutfactor: 3 1326 requires: !single 1327 1328 test: 1329 suffix: 4 1330 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 1331 timeoutfactor: 3 1332 requires: !single !__float128 1333 1334 TEST*/ 1335