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