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