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