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 /* Kulikov, 2013, Problem I */ 705 706 PetscErrorCode RHSFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec F, void *s) 707 { 708 PetscErrorCode ierr; 709 PetscScalar *f; 710 const PetscScalar *y; 711 712 PetscFunctionBegin; 713 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 714 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 715 f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3]; 716 f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.)); 717 f[2] = 2.*t*y[3]; 718 f[3] = -2.*t*PetscLogScalar(y[0]); 719 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 720 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 721 PetscFunctionReturn(0); 722 } 723 724 PetscErrorCode RHSJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Mat A, Mat B, void *s) 725 { 726 PetscErrorCode ierr; 727 const PetscScalar *y; 728 PetscInt row[4] = {0,1,2,3}; 729 PetscScalar value[4][4]; 730 PetscScalar m1,m2; 731 PetscFunctionBegin; 732 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 733 m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.)); 734 m2=2.*t*PetscPowScalar(y[1],1./5.); 735 value[0][0] = 0. ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2; 736 m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.)); 737 m2=10.*t*PetscExpScalar(5.0*(y[2]-1.)); 738 value[1][0] = 0.; value[1][1] = 0. ; value[1][2] = m1; value[1][3] = m2; 739 value[2][0] = 0.; value[2][1] = 0.; value[2][2] = 0.; value[2][3] = 2*t; 740 value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = 0.; 741 ierr = MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr); 742 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 743 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 744 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 745 PetscFunctionReturn(0); 746 } 747 748 PetscErrorCode IFunction_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 749 { 750 PetscErrorCode ierr; 751 PetscScalar *f; 752 const PetscScalar *y; 753 754 PetscFunctionBegin; 755 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 756 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 757 f[0] = 2.*t*PetscPowScalar(y[1],1./5.)*y[3]; 758 f[1] = 10.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.)); 759 f[2] = 2.*t*y[3]; 760 f[3] = -2.*t*PetscLogScalar(y[0]); 761 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 762 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 763 /* Left hand side = ydot - f(y) */ 764 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr); 765 PetscFunctionReturn(0); 766 } 767 768 PetscErrorCode IJacobian_Kulikov2013I(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 769 { 770 PetscErrorCode ierr; 771 const PetscScalar *y; 772 PetscInt row[4] = {0,1,2,3}; 773 PetscScalar value[4][4]; 774 PetscScalar m1,m2; 775 776 PetscFunctionBegin; 777 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 778 m1=(2.*t*y[3])/(5.*PetscPowScalar(y[1],4./5.)); 779 m2=2.*t*PetscPowScalar(y[1],1./5.); 780 value[0][0] = a ; value[0][1] = m1; value[0][2] = 0.; value[0][3] = m2; 781 m1=50.*t*y[3]*PetscExpScalar(5.0*(y[2]-1.)); 782 m2=10.*t*PetscExpScalar(5.0*(y[2]-1.)); 783 value[1][0] = 0.; value[1][1] = a ; value[1][2] = m1; value[1][3] = m2; 784 value[2][0] = 0.; value[2][1] = 0.; value[2][2] = a; value[2][3] = 2*t; 785 value[3][0] = -2.*t/y[0];value[3][1] = 0.; value[3][2] = 0.; value[3][3] = a; 786 ierr = MatSetValues(A,4,&row[0],4,&row[0],&value[0][0],INSERT_VALUES);CHKERRQ(ierr); 787 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 788 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 789 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 790 PetscFunctionReturn(0); 791 } 792 793 /* Hull, 1972, Problem C1 */ 794 795 PetscErrorCode RHSFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec F, void *s) 796 { 797 PetscErrorCode ierr; 798 PetscScalar *f; 799 const PetscScalar *y; 800 PetscInt N,i; 801 802 PetscFunctionBegin; 803 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 804 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 805 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 806 f[0] = -y[0]; 807 for (i = 1; i < N-1; i++) { 808 f[i] = y[i-1] - y[i]; 809 } 810 f[N-1] = y[N-2]; 811 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 812 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 813 PetscFunctionReturn(0); 814 } 815 816 PetscErrorCode IFunction_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 817 { 818 PetscErrorCode ierr; 819 PetscScalar *f; 820 const PetscScalar *y; 821 PetscInt N,i; 822 823 PetscFunctionBegin; 824 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 825 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 826 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 827 f[0] = -y[0]; 828 for (i = 1; i < N-1; i++) { 829 f[i] = y[i-1] - y[i]; 830 } 831 f[N-1] = y[N-2]; 832 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 833 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 834 /* Left hand side = ydot - f(y) */ 835 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 PetscErrorCode IJacobian_Hull1972C1(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 840 { 841 PetscErrorCode ierr; 842 const PetscScalar *y; 843 PetscInt N,i,col[2]; 844 PetscScalar value[2]; 845 846 PetscFunctionBegin; 847 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 848 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 849 i = 0; 850 value[0] = a+1; col[0] = 0; 851 value[1] = 0; col[1] = 1; 852 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 853 for (i = 0; i < N; i++) { 854 value[0] = -1; col[0] = i-1; 855 value[1] = a+1; col[1] = i; 856 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 857 } 858 i = N-1; 859 value[0] = -1; col[0] = N-2; 860 value[1] = a; col[1] = N-1; 861 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 862 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 863 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 864 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 868 /* Hull, 1972, Problem C2 */ 869 870 PetscErrorCode RHSFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec F, void *s) 871 { 872 PetscErrorCode ierr; 873 const PetscScalar *y; 874 PetscScalar *f; 875 PetscInt N,i; 876 877 PetscFunctionBegin; 878 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 879 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 880 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 881 f[0] = -y[0]; 882 for (i = 1; i < N-1; i++) { 883 f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i]; 884 } 885 f[N-1] = (PetscReal)(N-1)*y[N-2]; 886 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 887 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 888 PetscFunctionReturn(0); 889 } 890 891 PetscErrorCode IFunction_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 892 { 893 PetscErrorCode ierr; 894 PetscScalar *f; 895 const PetscScalar *y; 896 PetscInt N,i; 897 898 PetscFunctionBegin; 899 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 900 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 901 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 902 f[0] = -y[0]; 903 for (i = 1; i < N-1; i++) { 904 f[i] = (PetscReal)i*y[i-1] - (PetscReal)(i+1)*y[i]; 905 } 906 f[N-1] = (PetscReal)(N-1)*y[N-2]; 907 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 908 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 909 /* Left hand side = ydot - f(y) */ 910 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr); 911 PetscFunctionReturn(0); 912 } 913 914 PetscErrorCode IJacobian_Hull1972C2(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 915 { 916 PetscErrorCode ierr; 917 const PetscScalar *y; 918 PetscInt N,i,col[2]; 919 PetscScalar value[2]; 920 921 PetscFunctionBegin; 922 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 923 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 924 i = 0; 925 value[0] = a+1; col[0] = 0; 926 value[1] = 0; col[1] = 1; 927 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 928 for (i = 0; i < N; i++) { 929 value[0] = -(PetscReal) i; col[0] = i-1; 930 value[1] = a+(PetscReal)(i+1); col[1] = i; 931 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 932 } 933 i = N-1; 934 value[0] = -(PetscReal) (N-1); col[0] = N-2; 935 value[1] = a; col[1] = N-1; 936 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 937 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 938 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 939 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 940 PetscFunctionReturn(0); 941 } 942 943 /* Hull, 1972, Problem C3 and C4 */ 944 945 PetscErrorCode RHSFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec F, void *s) 946 { 947 PetscErrorCode ierr; 948 PetscScalar *f; 949 const PetscScalar *y; 950 PetscInt N,i; 951 952 PetscFunctionBegin; 953 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 954 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 955 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 956 f[0] = -2.0*y[0] + y[1]; 957 for (i = 1; i < N-1; i++) { 958 f[i] = y[i-1] - 2.0*y[i] + y[i+1]; 959 } 960 f[N-1] = y[N-2] - 2.0*y[N-1]; 961 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 962 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 963 PetscFunctionReturn(0); 964 } 965 966 PetscErrorCode IFunction_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *s) 967 { 968 PetscErrorCode ierr; 969 PetscScalar *f; 970 const PetscScalar *y; 971 PetscInt N,i; 972 973 PetscFunctionBegin; 974 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 975 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 976 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 977 f[0] = -2.0*y[0] + y[1]; 978 for (i = 1; i < N-1; i++) { 979 f[i] = y[i-1] - 2.0*y[i] + y[i+1]; 980 } 981 f[N-1] = y[N-2] - 2.0*y[N-1]; 982 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 983 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 984 /* Left hand side = ydot - f(y) */ 985 ierr = VecAYPX(F,-1.0,Ydot);CHKERRQ(ierr); 986 PetscFunctionReturn(0); 987 } 988 989 PetscErrorCode IJacobian_Hull1972C34(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *s) 990 { 991 PetscErrorCode ierr; 992 const PetscScalar *y; 993 PetscScalar value[3]; 994 PetscInt N,i,col[3]; 995 996 PetscFunctionBegin; 997 ierr = VecGetSize (Y,&N);CHKERRQ(ierr); 998 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 999 for (i = 0; i < N; i++) { 1000 if (i == 0) { 1001 value[0] = a+2; col[0] = i; 1002 value[1] = -1; col[1] = i+1; 1003 value[2] = 0; col[2] = i+2; 1004 } else if (i == N-1) { 1005 value[0] = 0; col[0] = i-2; 1006 value[1] = -1; col[1] = i-1; 1007 value[2] = a+2; col[2] = i; 1008 } else { 1009 value[0] = -1; col[0] = i-1; 1010 value[1] = a+2; col[1] = i; 1011 value[2] = -1; col[2] = i+1; 1012 } 1013 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 1014 } 1015 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1016 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1017 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 1018 PetscFunctionReturn(0); 1019 } 1020 1021 /***************************************************************************/ 1022 1023 /* Sets the initial solution for the IVP and sets up the function pointers*/ 1024 PetscErrorCode Initialize(Vec Y, void* s) 1025 { 1026 PetscErrorCode ierr; 1027 char *p = (char*) s; 1028 PetscScalar *y; 1029 PetscReal t0; 1030 PetscInt N = GetSize((const char *)s); 1031 PetscBool flg; 1032 1033 PetscFunctionBegin; 1034 VecZeroEntries(Y); 1035 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1036 if (!strcmp(p,"hull1972a1")) { 1037 y[0] = 1.0; 1038 RHSFunction = RHSFunction_Hull1972A1; 1039 RHSJacobian = RHSJacobian_Hull1972A1; 1040 IFunction = IFunction_Hull1972A1; 1041 IJacobian = IJacobian_Hull1972A1; 1042 } else if (!strcmp(p,"hull1972a2")) { 1043 y[0] = 1.0; 1044 RHSFunction = RHSFunction_Hull1972A2; 1045 RHSJacobian = RHSJacobian_Hull1972A2; 1046 IFunction = IFunction_Hull1972A2; 1047 IJacobian = IJacobian_Hull1972A2; 1048 } else if (!strcmp(p,"hull1972a3")) { 1049 y[0] = 1.0; 1050 RHSFunction = RHSFunction_Hull1972A3; 1051 RHSJacobian = RHSJacobian_Hull1972A3; 1052 IFunction = IFunction_Hull1972A3; 1053 IJacobian = IJacobian_Hull1972A3; 1054 } else if (!strcmp(p,"hull1972a4")) { 1055 y[0] = 1.0; 1056 RHSFunction = RHSFunction_Hull1972A4; 1057 RHSJacobian = RHSJacobian_Hull1972A4; 1058 IFunction = IFunction_Hull1972A4; 1059 IJacobian = IJacobian_Hull1972A4; 1060 } else if (!strcmp(p,"hull1972a5")) { 1061 y[0] = 4.0; 1062 RHSFunction = RHSFunction_Hull1972A5; 1063 RHSJacobian = RHSJacobian_Hull1972A5; 1064 IFunction = IFunction_Hull1972A5; 1065 IJacobian = IJacobian_Hull1972A5; 1066 } else if (!strcmp(p,"hull1972b1")) { 1067 y[0] = 1.0; 1068 y[1] = 3.0; 1069 RHSFunction = RHSFunction_Hull1972B1; 1070 IFunction = IFunction_Hull1972B1; 1071 IJacobian = IJacobian_Hull1972B1; 1072 } else if (!strcmp(p,"hull1972b2")) { 1073 y[0] = 2.0; 1074 y[1] = 0.0; 1075 y[2] = 1.0; 1076 RHSFunction = RHSFunction_Hull1972B2; 1077 IFunction = IFunction_Hull1972B2; 1078 IJacobian = IJacobian_Hull1972B2; 1079 } else if (!strcmp(p,"hull1972b3")) { 1080 y[0] = 1.0; 1081 y[1] = 0.0; 1082 y[2] = 0.0; 1083 RHSFunction = RHSFunction_Hull1972B3; 1084 IFunction = IFunction_Hull1972B3; 1085 IJacobian = IJacobian_Hull1972B3; 1086 } else if (!strcmp(p,"hull1972b4")) { 1087 y[0] = 3.0; 1088 y[1] = 0.0; 1089 y[2] = 0.0; 1090 RHSFunction = RHSFunction_Hull1972B4; 1091 IFunction = IFunction_Hull1972B4; 1092 IJacobian = IJacobian_Hull1972B4; 1093 } else if (!strcmp(p,"hull1972b5")) { 1094 y[0] = 0.0; 1095 y[1] = 1.0; 1096 y[2] = 1.0; 1097 RHSFunction = RHSFunction_Hull1972B5; 1098 IFunction = IFunction_Hull1972B5; 1099 IJacobian = IJacobian_Hull1972B5; 1100 } else if (!strcmp(p,"kulik2013i")) { 1101 t0=0.; 1102 y[0] = PetscExpReal(PetscSinReal(t0*t0)); 1103 y[1] = PetscExpReal(5.*PetscSinReal(t0*t0)); 1104 y[2] = PetscSinReal(t0*t0)+1.0; 1105 y[3] = PetscCosReal(t0*t0); 1106 RHSFunction = RHSFunction_Kulikov2013I; 1107 RHSJacobian = RHSJacobian_Kulikov2013I; 1108 IFunction = IFunction_Kulikov2013I; 1109 IJacobian = IJacobian_Kulikov2013I; 1110 } else if (!strcmp(p,"hull1972c1")) { 1111 y[0] = 1.0; 1112 RHSFunction = RHSFunction_Hull1972C1; 1113 IFunction = IFunction_Hull1972C1; 1114 IJacobian = IJacobian_Hull1972C1; 1115 } else if (!strcmp(p,"hull1972c2")) { 1116 y[0] = 1.0; 1117 RHSFunction = RHSFunction_Hull1972C2; 1118 IFunction = IFunction_Hull1972C2; 1119 IJacobian = IJacobian_Hull1972C2; 1120 } else if ((!strcmp(p,"hull1972c3")) 1121 ||(!strcmp(p,"hull1972c4"))) { 1122 y[0] = 1.0; 1123 RHSFunction = RHSFunction_Hull1972C34; 1124 IFunction = IFunction_Hull1972C34; 1125 IJacobian = IJacobian_Hull1972C34; 1126 } 1127 ierr = PetscOptionsGetScalarArray(NULL,NULL,"-yinit",y,&N,&flg);CHKERRQ(ierr); 1128 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)); 1129 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1130 PetscFunctionReturn(0); 1131 } 1132 1133 /* Calculates the exact solution to problems that have one */ 1134 PetscErrorCode ExactSolution(Vec Y, void* s, PetscReal t, PetscBool *flag) 1135 { 1136 PetscErrorCode ierr; 1137 char *p = (char*) s; 1138 PetscScalar *y; 1139 1140 PetscFunctionBegin; 1141 if (!strcmp(p,"hull1972a1")) { 1142 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1143 y[0] = PetscExpReal(-t); 1144 *flag = PETSC_TRUE; 1145 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1146 } else if (!strcmp(p,"hull1972a2")) { 1147 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1148 y[0] = 1.0/PetscSqrtReal(t+1); 1149 *flag = PETSC_TRUE; 1150 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1151 } else if (!strcmp(p,"hull1972a3")) { 1152 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1153 y[0] = PetscExpReal(PetscSinReal(t)); 1154 *flag = PETSC_TRUE; 1155 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1156 } else if (!strcmp(p,"hull1972a4")) { 1157 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1158 y[0] = 20.0/(1+19.0*PetscExpReal(-t/4.0)); 1159 *flag = PETSC_TRUE; 1160 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1161 } else if (!strcmp(p,"kulik2013i")) { 1162 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 1163 y[0] = PetscExpReal(PetscSinReal(t*t)); 1164 y[1] = PetscExpReal(5.*PetscSinReal(t*t)); 1165 y[2] = PetscSinReal(t*t)+1.0; 1166 y[3] = PetscCosReal(t*t); 1167 *flag = PETSC_TRUE; 1168 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 1169 } else { 1170 ierr = VecSet(Y,0);CHKERRQ(ierr); 1171 *flag = PETSC_FALSE; 1172 } 1173 PetscFunctionReturn(0); 1174 } 1175 1176 /* Solves the specified ODE and computes the error if exact solution is available */ 1177 PetscErrorCode SolveODE(char* ptype, PetscReal dt, PetscReal tfinal, PetscInt maxiter, PetscReal *error, PetscBool *exact_flag) 1178 { 1179 PetscErrorCode ierr; /* Error code */ 1180 TS ts; /* time-integrator */ 1181 Vec Y; /* Solution vector */ 1182 Vec Yex; /* Exact solution */ 1183 PetscInt N; /* Size of the system of equations */ 1184 TSType time_scheme; /* Type of time-integration scheme */ 1185 Mat Jac = NULL; /* Jacobian matrix */ 1186 Vec Yerr; /* Auxiliary solution vector */ 1187 PetscReal err_norm; /* Estimated error norm */ 1188 PetscReal final_time; /* Actual final time from the integrator */ 1189 1190 PetscFunctionBegin; 1191 N = GetSize((const char *)&ptype[0]); 1192 if (N < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"Illegal problem specification.\n"); 1193 ierr = VecCreate(PETSC_COMM_WORLD,&Y);CHKERRQ(ierr); 1194 ierr = VecSetSizes(Y,N,PETSC_DECIDE);CHKERRQ(ierr); 1195 ierr = VecSetUp(Y);CHKERRQ(ierr); 1196 ierr = VecSet(Y,0);CHKERRQ(ierr); 1197 1198 /* Initialize the problem */ 1199 ierr = Initialize(Y,&ptype[0]);CHKERRQ(ierr); 1200 1201 /* Create and initialize the time-integrator */ 1202 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 1203 /* Default time integration options */ 1204 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 1205 ierr = TSSetMaxSteps(ts,maxiter);CHKERRQ(ierr); 1206 ierr = TSSetMaxTime(ts,tfinal);CHKERRQ(ierr); 1207 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 1208 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 1209 /* Read command line options for time integration */ 1210 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 1211 /* Set solution vector */ 1212 ierr = TSSetSolution(ts,Y);CHKERRQ(ierr); 1213 /* Specify left/right-hand side functions */ 1214 ierr = TSGetType(ts,&time_scheme);CHKERRQ(ierr); 1215 1216 if ((!strcmp(time_scheme,TSEULER)) || (!strcmp(time_scheme,TSRK)) || (!strcmp(time_scheme,TSSSP) || (!strcmp(time_scheme,TSGLEE)))) { 1217 /* Explicit time-integration -> specify right-hand side function ydot = f(y) */ 1218 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&ptype[0]);CHKERRQ(ierr); 1219 ierr = MatCreate(PETSC_COMM_WORLD,&Jac);CHKERRQ(ierr); 1220 ierr = MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); 1221 ierr = MatSetFromOptions(Jac);CHKERRQ(ierr); 1222 ierr = MatSetUp(Jac);CHKERRQ(ierr); 1223 ierr = TSSetRHSJacobian(ts,Jac,Jac,RHSJacobian,&ptype[0]);CHKERRQ(ierr); 1224 } else if ((!strcmp(time_scheme,TSTHETA)) || (!strcmp(time_scheme,TSBEULER)) || (!strcmp(time_scheme,TSCN)) || (!strcmp(time_scheme,TSALPHA)) || (!strcmp(time_scheme,TSARKIMEX))) { 1225 /* Implicit time-integration -> specify left-hand side function ydot-f(y) = 0 */ 1226 /* and its Jacobian function */ 1227 ierr = TSSetIFunction(ts,NULL,IFunction,&ptype[0]);CHKERRQ(ierr); 1228 ierr = MatCreate(PETSC_COMM_WORLD,&Jac);CHKERRQ(ierr); 1229 ierr = MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); 1230 ierr = MatSetFromOptions(Jac);CHKERRQ(ierr); 1231 ierr = MatSetUp(Jac);CHKERRQ(ierr); 1232 ierr = TSSetIJacobian(ts,Jac,Jac,IJacobian,&ptype[0]);CHKERRQ(ierr); 1233 } 1234 1235 /* Solve */ 1236 ierr = TSSolve(ts,Y);CHKERRQ(ierr); 1237 ierr = TSGetTime(ts,&final_time);CHKERRQ(ierr); 1238 1239 /* Get the estimated error, if available */ 1240 ierr = VecDuplicate(Y,&Yerr);CHKERRQ(ierr); 1241 ierr = VecZeroEntries(Yerr);CHKERRQ(ierr); 1242 ierr = TSGetTimeError(ts,0,&Yerr);CHKERRQ(ierr); 1243 ierr = VecNorm(Yerr,NORM_2,&err_norm);CHKERRQ(ierr); 1244 ierr = VecDestroy(&Yerr);CHKERRQ(ierr); 1245 ierr = PetscPrintf(PETSC_COMM_WORLD,"Estimated Error = %E.\n",err_norm);CHKERRQ(ierr); 1246 1247 /* Exact solution */ 1248 ierr = VecDuplicate(Y,&Yex);CHKERRQ(ierr); 1249 if (PetscAbsScalar(final_time-tfinal)>2.*PETSC_MACHINE_EPSILON) { 1250 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); 1251 } 1252 ierr = ExactSolution(Yex,&ptype[0],final_time,exact_flag);CHKERRQ(ierr); 1253 1254 /* Calculate Error */ 1255 ierr = VecAYPX(Yex,-1.0,Y);CHKERRQ(ierr); 1256 ierr = VecNorm(Yex,NORM_2,error);CHKERRQ(ierr); 1257 *error = PetscSqrtReal(((*error)*(*error))/N); 1258 1259 /* Clean up and finalize */ 1260 ierr = MatDestroy(&Jac);CHKERRQ(ierr); 1261 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1262 ierr = VecDestroy(&Yex);CHKERRQ(ierr); 1263 ierr = VecDestroy(&Y);CHKERRQ(ierr); 1264 1265 PetscFunctionReturn(0); 1266 } 1267 1268 int main(int argc, char **argv) 1269 { 1270 PetscErrorCode ierr; /* Error code */ 1271 char ptype[256] = "hull1972a1"; /* Problem specification */ 1272 PetscInt n_refine = 1; /* Number of refinement levels for convergence analysis */ 1273 PetscReal refine_fac = 2.0; /* Refinement factor for dt */ 1274 PetscReal dt_initial = 0.01; /* Initial default value of dt */ 1275 PetscReal dt; 1276 PetscReal tfinal = 20.0; /* Final time for the time-integration */ 1277 PetscInt maxiter = 100000; /* Maximum number of time-integration iterations */ 1278 PetscReal *error; /* Array to store the errors for convergence analysis */ 1279 PetscMPIInt size; /* No of processors */ 1280 PetscBool flag; /* Flag denoting availability of exact solution */ 1281 PetscInt r; 1282 1283 /* Initialize program */ 1284 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 1285 1286 /* Check if running with only 1 proc */ 1287 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 1288 if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 1289 1290 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex31",NULL);CHKERRQ(ierr); 1291 ierr = PetscOptionsString("-problem","Problem specification","<hull1972a1>",ptype,ptype,sizeof(ptype),NULL);CHKERRQ(ierr); 1292 ierr = PetscOptionsInt("-refinement_levels","Number of refinement levels for convergence analysis","<1>",n_refine,&n_refine,NULL);CHKERRQ(ierr); 1293 ierr = PetscOptionsReal("-refinement_factor","Refinement factor for dt","<2.0>",refine_fac,&refine_fac,NULL);CHKERRQ(ierr); 1294 ierr = PetscOptionsReal("-dt","Time step size (for convergence analysis, initial time step)","<0.01>",dt_initial,&dt_initial,NULL);CHKERRQ(ierr); 1295 ierr = PetscOptionsReal("-final_time","Final time for the time-integration","<20.0>",tfinal,&tfinal,NULL);CHKERRQ(ierr); 1296 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1297 1298 ierr = PetscMalloc1(n_refine,&error);CHKERRQ(ierr); 1299 for (r = 0,dt = dt_initial; r < n_refine; r++) { 1300 error[r] = 0; 1301 if (r > 0) dt /= refine_fac; 1302 1303 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); 1304 ierr = SolveODE(&ptype[0],dt,tfinal,maxiter,&error[r],&flag);CHKERRQ(ierr); 1305 if (flag) { 1306 /* If exact solution available for the specified ODE */ 1307 if (r > 0) { 1308 PetscReal conv_rate = (PetscLogReal(error[r]) - PetscLogReal(error[r-1])) / (-PetscLogReal(refine_fac)); 1309 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error = %E,\tConvergence rate = %f.\n",(double)error[r],(double)conv_rate);CHKERRQ(ierr); 1310 } else { 1311 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error = %E.\n",error[r]);CHKERRQ(ierr); 1312 } 1313 } 1314 } 1315 ierr = PetscFree(error);CHKERRQ(ierr); 1316 ierr = PetscFinalize(); 1317 return ierr; 1318 } 1319 1320 /*TEST 1321 1322 test: 1323 suffix: 2 1324 args: -ts_type glee -final_time 5 -ts_adapt_type none 1325 timeoutfactor: 3 1326 requires: !single 1327 1328 test: 1329 suffix: 3 1330 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 1331 timeoutfactor: 3 1332 requires: !single 1333 1334 test: 1335 suffix: 4 1336 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 1337 timeoutfactor: 3 1338 requires: !single !__float128 1339 1340 TEST*/ 1341