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