1 2 static char help[] = "Nonlinear DAE benchmark problems.\n"; 3 4 /* 5 Include "petscts.h" so that we can use TS solvers. Note that this 6 file automatically includes: 7 petscsys.h - base PETSc routines petscvec.h - vectors 8 petscmat.h - matrices 9 petscis.h - index sets petscksp.h - Krylov subspace methods 10 petscviewer.h - viewers petscpc.h - preconditioners 11 petscksp.h - linear solvers 12 */ 13 #include <petscts.h> 14 15 typedef struct _Problem* Problem; 16 struct _Problem { 17 PetscErrorCode (*destroy)(Problem); 18 TSIFunction function; 19 TSIJacobian jacobian; 20 PetscErrorCode (*solution)(PetscReal,Vec,void*); 21 MPI_Comm comm; 22 PetscReal final_time; 23 PetscInt n; 24 PetscBool hasexact; 25 void *data; 26 }; 27 28 /* 29 Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996 30 */ 31 static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 32 { 33 PetscErrorCode ierr; 34 PetscScalar *f; 35 const PetscScalar *x,*xdot; 36 37 PetscFunctionBeginUser; 38 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 39 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 40 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 41 f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2]; 42 f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]); 43 f[2] = xdot[2] - 3e7*PetscSqr(x[1]); 44 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 45 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 46 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 51 { 52 PetscErrorCode ierr; 53 PetscInt rowcol[] = {0,1,2}; 54 PetscScalar J[3][3]; 55 const PetscScalar *x,*xdot; 56 57 PetscFunctionBeginUser; 58 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 59 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 60 J[0][0] = a + 0.04; J[0][1] = -1e4*x[2]; J[0][2] = -1e4*x[1]; 61 J[1][0] = -0.04; J[1][1] = a + 1e4*x[2] + 3e7*2*x[1]; J[1][2] = 1e4*x[1]; 62 J[2][0] = 0; J[2][1] = -3e7*2*x[1]; J[2][2] = a; 63 ierr = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 64 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 65 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 66 67 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 68 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 69 if (A != B) { 70 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 71 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 72 } 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx) 77 { 78 PetscErrorCode ierr; 79 PetscScalar *x; 80 81 PetscFunctionBeginUser; 82 if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented"); 83 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 84 x[0] = 1; 85 x[1] = 0; 86 x[2] = 0; 87 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 91 static PetscErrorCode RoberCreate(Problem p) 92 { 93 94 PetscFunctionBeginUser; 95 p->destroy = 0; 96 p->function = &RoberFunction; 97 p->jacobian = &RoberJacobian; 98 p->solution = &RoberSolution; 99 p->final_time = 1e11; 100 p->n = 3; 101 PetscFunctionReturn(0); 102 } 103 104 /* 105 Stiff scalar valued problem 106 */ 107 108 typedef struct { 109 PetscReal lambda; 110 } CECtx; 111 112 static PetscErrorCode CEDestroy(Problem p) 113 { 114 PetscErrorCode ierr; 115 116 PetscFunctionBeginUser; 117 ierr = PetscFree(p->data);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 122 { 123 PetscErrorCode ierr; 124 PetscReal l = ((CECtx*)ctx)->lambda; 125 PetscScalar *f; 126 const PetscScalar *x,*xdot; 127 128 PetscFunctionBeginUser; 129 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 130 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 131 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 132 f[0] = xdot[0] + l*(x[0] - PetscCosReal(t)); 133 #if 0 134 ierr = PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);CHKERRQ(ierr); 135 #endif 136 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 137 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 138 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 143 { 144 PetscReal l = ((CECtx*)ctx)->lambda; 145 PetscErrorCode ierr; 146 PetscInt rowcol[] = {0}; 147 PetscScalar J[1][1]; 148 const PetscScalar *x,*xdot; 149 150 PetscFunctionBeginUser; 151 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 152 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 153 J[0][0] = a + l; 154 ierr = MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 155 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 156 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 157 158 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 159 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 160 if (A != B) { 161 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 162 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 163 } 164 PetscFunctionReturn(0); 165 } 166 167 static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx) 168 { 169 PetscReal l = ((CECtx*)ctx)->lambda; 170 PetscErrorCode ierr; 171 PetscScalar *x; 172 173 PetscFunctionBeginUser; 174 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 175 x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t); 176 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode CECreate(Problem p) 181 { 182 PetscErrorCode ierr; 183 CECtx *ce; 184 185 PetscFunctionBeginUser; 186 ierr = PetscMalloc(sizeof(CECtx),&ce);CHKERRQ(ierr); 187 p->data = (void*)ce; 188 189 p->destroy = &CEDestroy; 190 p->function = &CEFunction; 191 p->jacobian = &CEJacobian; 192 p->solution = &CESolution; 193 p->final_time = 10; 194 p->n = 1; 195 p->hasexact = PETSC_TRUE; 196 197 ce->lambda = 10; 198 ierr = PetscOptionsBegin(p->comm,NULL,"CE options","");CHKERRQ(ierr); 199 { 200 ierr = PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);CHKERRQ(ierr); 201 } 202 ierr = PetscOptionsEnd();CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 /* 207 * Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner 208 */ 209 static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 210 { 211 PetscErrorCode ierr; 212 PetscScalar *f; 213 const PetscScalar *x,*xdot; 214 215 PetscFunctionBeginUser; 216 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 217 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 218 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 219 f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1])); 220 f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]); 221 f[2] = xdot[2] - 0.161*(x[0] - x[2]); 222 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 223 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 224 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 228 static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 229 { 230 PetscErrorCode ierr; 231 PetscInt rowcol[] = {0,1,2}; 232 PetscScalar J[3][3]; 233 const PetscScalar *x,*xdot; 234 235 PetscFunctionBeginUser; 236 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 237 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 238 J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]); 239 J[0][1] = -77.27*(1. - x[0]); 240 J[0][2] = 0; 241 J[1][0] = 1./77.27*x[1]; 242 J[1][1] = a + 1./77.27*(1. + x[0]); 243 J[1][2] = -1./77.27; 244 J[2][0] = -0.161; 245 J[2][1] = 0; 246 J[2][2] = a + 0.161; 247 ierr = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 248 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 249 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 250 251 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 253 if (A != B) { 254 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 255 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 256 } 257 PetscFunctionReturn(0); 258 } 259 260 static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx) 261 { 262 PetscErrorCode ierr; 263 PetscScalar *x; 264 265 PetscFunctionBeginUser; 266 if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented"); 267 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 268 x[0] = 1; 269 x[1] = 2; 270 x[2] = 3; 271 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 static PetscErrorCode OregoCreate(Problem p) 276 { 277 278 PetscFunctionBeginUser; 279 p->destroy = 0; 280 p->function = &OregoFunction; 281 p->jacobian = &OregoJacobian; 282 p->solution = &OregoSolution; 283 p->final_time = 360; 284 p->n = 3; 285 PetscFunctionReturn(0); 286 } 287 288 289 /* 290 * User-defined monitor for comparing to exact solutions when possible 291 */ 292 typedef struct { 293 MPI_Comm comm; 294 Problem problem; 295 Vec x; 296 } MonitorCtx; 297 298 static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx) 299 { 300 PetscErrorCode ierr; 301 MonitorCtx *mon = (MonitorCtx*)ctx; 302 PetscReal h,nrm_x,nrm_exact,nrm_diff; 303 304 PetscFunctionBeginUser; 305 if (!mon->problem->solution) PetscFunctionReturn(0); 306 ierr = (*mon->problem->solution)(t,mon->x,mon->problem->data);CHKERRQ(ierr); 307 ierr = VecNorm(x,NORM_2,&nrm_x);CHKERRQ(ierr); 308 ierr = VecNorm(mon->x,NORM_2,&nrm_exact);CHKERRQ(ierr); 309 ierr = VecAYPX(mon->x,-1,x);CHKERRQ(ierr); 310 ierr = VecNorm(mon->x,NORM_2,&nrm_diff);CHKERRQ(ierr); 311 ierr = TSGetTimeStep(ts,&h);CHKERRQ(ierr); 312 if (step < 0) { 313 ierr = PetscPrintf(mon->comm,"Interpolated final solution ");CHKERRQ(ierr); 314 } 315 ierr = PetscPrintf(mon->comm,"step %4D t=%12.8e h=% 8.2e |x|=%9.2e |x_e|=%9.2e |x-x_e|=%9.2e\n",step,(double)t,(double)h,(double)nrm_x,(double)nrm_exact,(double)nrm_diff);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318 319 320 int main(int argc,char **argv) 321 { 322 PetscFunctionList plist = NULL; 323 char pname[256]; 324 TS ts; /* nonlinear solver */ 325 Vec x,r; /* solution, residual vectors */ 326 Mat A; /* Jacobian matrix */ 327 Problem problem; 328 PetscBool use_monitor; 329 PetscInt steps,nonlinits,linits,snesfails,rejects; 330 PetscReal ftime; 331 MonitorCtx mon; 332 PetscErrorCode ierr; 333 PetscMPIInt size; 334 335 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 336 Initialize program 337 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 338 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 339 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 340 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 341 342 /* Register the available problems */ 343 ierr = PetscFunctionListAdd(&plist,"rober",&RoberCreate);CHKERRQ(ierr); 344 ierr = PetscFunctionListAdd(&plist,"ce",&CECreate);CHKERRQ(ierr); 345 ierr = PetscFunctionListAdd(&plist,"orego",&OregoCreate);CHKERRQ(ierr); 346 ierr = PetscStrcpy(pname,"ce");CHKERRQ(ierr); 347 348 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 349 Set runtime options 350 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 351 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");CHKERRQ(ierr); 352 { 353 ierr = PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);CHKERRQ(ierr); 354 use_monitor = PETSC_FALSE; 355 ierr = PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);CHKERRQ(ierr); 356 } 357 ierr = PetscOptionsEnd();CHKERRQ(ierr); 358 359 /* Create the new problem */ 360 ierr = PetscNew(&problem);CHKERRQ(ierr); 361 problem->comm = MPI_COMM_WORLD; 362 { 363 PetscErrorCode (*pcreate)(Problem); 364 365 ierr = PetscFunctionListFind(plist,pname,&pcreate);CHKERRQ(ierr); 366 if (!pcreate) SETERRQ1(PETSC_COMM_SELF,1,"No problem '%s'",pname); 367 ierr = (*pcreate)(problem);CHKERRQ(ierr); 368 } 369 370 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 371 Create necessary matrix and vectors 372 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 373 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 374 ierr = MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 375 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 376 ierr = MatSetUp(A);CHKERRQ(ierr); 377 378 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 379 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 380 381 mon.comm = PETSC_COMM_WORLD; 382 mon.problem = problem; 383 ierr = VecDuplicate(x,&mon.x);CHKERRQ(ierr); 384 385 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 386 Create timestepping solver context 387 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 388 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 389 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 390 ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); /* Rosenbrock-W */ 391 ierr = TSSetIFunction(ts,NULL,problem->function,problem->data);CHKERRQ(ierr); 392 ierr = TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);CHKERRQ(ierr); 393 ierr = TSSetMaxTime(ts,problem->final_time);CHKERRQ(ierr); 394 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 395 ierr = TSSetMaxStepRejections(ts,10);CHKERRQ(ierr); 396 ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr); /* unlimited */ 397 if (use_monitor) { 398 ierr = TSMonitorSet(ts,&MonitorError,&mon,NULL);CHKERRQ(ierr); 399 } 400 401 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 402 Set initial conditions 403 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 404 ierr = (*problem->solution)(0,x,problem->data);CHKERRQ(ierr); 405 ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 406 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 407 408 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 409 Set runtime options 410 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 411 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 412 413 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 414 Solve nonlinear system 415 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 416 ierr = TSSolve(ts,x);CHKERRQ(ierr); 417 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 418 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 419 ierr = TSGetSNESFailures(ts,&snesfails);CHKERRQ(ierr); 420 ierr = TSGetStepRejections(ts,&rejects);CHKERRQ(ierr); 421 ierr = TSGetSNESIterations(ts,&nonlinits);CHKERRQ(ierr); 422 ierr = TSGetKSPIterations(ts,&linits);CHKERRQ(ierr); 423 ierr = PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",steps,rejects,snesfails,(double)ftime,nonlinits,linits);CHKERRQ(ierr); 424 425 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 426 Free work space. All PETSc objects should be destroyed when they 427 are no longer needed. 428 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 429 ierr = MatDestroy(&A);CHKERRQ(ierr); 430 ierr = VecDestroy(&x);CHKERRQ(ierr); 431 ierr = VecDestroy(&r);CHKERRQ(ierr); 432 ierr = VecDestroy(&mon.x);CHKERRQ(ierr); 433 ierr = TSDestroy(&ts);CHKERRQ(ierr); 434 if (problem->destroy) { 435 ierr = (*problem->destroy)(problem);CHKERRQ(ierr); 436 } 437 ierr = PetscFree(problem);CHKERRQ(ierr); 438 ierr = PetscFunctionListDestroy(&plist);CHKERRQ(ierr); 439 440 ierr = PetscFinalize(); 441 return ierr; 442 } 443 444 /*TEST 445 446 test: 447 requires: !complex 448 args: -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex 449 450 test: 451 suffix: 2 452 requires: !single !complex 453 args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 0 -ts_adapt_time_step_increase_delay 4 454 455 test: 456 suffix: 3 457 requires: !single !complex 458 args: -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 1 459 460 TEST*/ 461