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 CECtx *ce; 174 175 PetscFunctionBeginUser; 176 PetscCall(PetscMalloc(sizeof(CECtx),&ce)); 177 p->data = (void*)ce; 178 179 p->destroy = &CEDestroy; 180 p->function = &CEFunction; 181 p->jacobian = &CEJacobian; 182 p->solution = &CESolution; 183 p->final_time = 10; 184 p->n = 1; 185 p->hasexact = PETSC_TRUE; 186 187 ce->lambda = 10; 188 PetscOptionsBegin(p->comm,NULL,"CE options",""); 189 { 190 PetscCall(PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL)); 191 } 192 PetscOptionsEnd(); 193 PetscFunctionReturn(0); 194 } 195 196 /* 197 Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner 198 */ 199 static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 200 { 201 PetscScalar *f; 202 const PetscScalar *x,*xdot; 203 204 PetscFunctionBeginUser; 205 PetscCall(VecGetArrayRead(X,&x)); 206 PetscCall(VecGetArrayRead(Xdot,&xdot)); 207 PetscCall(VecGetArray(F,&f)); 208 f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1])); 209 f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]); 210 f[2] = xdot[2] - 0.161*(x[0] - x[2]); 211 PetscCall(VecRestoreArrayRead(X,&x)); 212 PetscCall(VecRestoreArrayRead(Xdot,&xdot)); 213 PetscCall(VecRestoreArray(F,&f)); 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 218 { 219 PetscInt rowcol[] = {0,1,2}; 220 PetscScalar J[3][3]; 221 const PetscScalar *x,*xdot; 222 223 PetscFunctionBeginUser; 224 PetscCall(VecGetArrayRead(X,&x)); 225 PetscCall(VecGetArrayRead(Xdot,&xdot)); 226 J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]); 227 J[0][1] = -77.27*(1. - x[0]); 228 J[0][2] = 0; 229 J[1][0] = 1./77.27*x[1]; 230 J[1][1] = a + 1./77.27*(1. + x[0]); 231 J[1][2] = -1./77.27; 232 J[2][0] = -0.161; 233 J[2][1] = 0; 234 J[2][2] = a + 0.161; 235 PetscCall(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 236 PetscCall(VecRestoreArrayRead(X,&x)); 237 PetscCall(VecRestoreArrayRead(Xdot,&xdot)); 238 239 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 240 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 241 if (A != B) { 242 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 243 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 244 } 245 PetscFunctionReturn(0); 246 } 247 248 static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx) 249 { 250 PetscScalar *x; 251 252 PetscFunctionBeginUser; 253 PetscCheck(t == 0,PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented"); 254 PetscCall(VecGetArray(X,&x)); 255 x[0] = 1; 256 x[1] = 2; 257 x[2] = 3; 258 PetscCall(VecRestoreArray(X,&x)); 259 PetscFunctionReturn(0); 260 } 261 262 static PetscErrorCode OregoCreate(Problem p) 263 { 264 PetscFunctionBeginUser; 265 p->destroy = 0; 266 p->function = &OregoFunction; 267 p->jacobian = &OregoJacobian; 268 p->solution = &OregoSolution; 269 p->final_time = 360; 270 p->n = 3; 271 PetscFunctionReturn(0); 272 } 273 274 /* 275 User-defined monitor for comparing to exact solutions when possible 276 */ 277 typedef struct { 278 MPI_Comm comm; 279 Problem problem; 280 Vec x; 281 } MonitorCtx; 282 283 static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx) 284 { 285 MonitorCtx *mon = (MonitorCtx*)ctx; 286 PetscReal h,nrm_x,nrm_exact,nrm_diff; 287 288 PetscFunctionBeginUser; 289 if (!mon->problem->solution) PetscFunctionReturn(0); 290 PetscCall((*mon->problem->solution)(t,mon->x,mon->problem->data)); 291 PetscCall(VecNorm(x,NORM_2,&nrm_x)); 292 PetscCall(VecNorm(mon->x,NORM_2,&nrm_exact)); 293 PetscCall(VecAYPX(mon->x,-1,x)); 294 PetscCall(VecNorm(mon->x,NORM_2,&nrm_diff)); 295 PetscCall(TSGetTimeStep(ts,&h)); 296 if (step < 0) { 297 PetscCall(PetscPrintf(mon->comm,"Interpolated final solution ")); 298 } 299 PetscCall(PetscPrintf(mon->comm,"step %4" PetscInt_FMT " 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)); 300 PetscFunctionReturn(0); 301 } 302 303 int main(int argc,char **argv) 304 { 305 PetscFunctionList plist = NULL; 306 char pname[256]; 307 TS ts; /* nonlinear solver */ 308 Vec x,r; /* solution, residual vectors */ 309 Mat A; /* Jacobian matrix */ 310 Problem problem; 311 PetscBool use_monitor = PETSC_FALSE; 312 PetscBool use_result = PETSC_FALSE; 313 PetscInt steps,nonlinits,linits,snesfails,rejects; 314 PetscReal ftime; 315 MonitorCtx mon; 316 PetscMPIInt size; 317 318 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 319 Initialize program 320 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 321 PetscFunctionBeginUser; 322 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 323 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 324 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 325 326 /* Register the available problems */ 327 PetscCall(PetscFunctionListAdd(&plist,"rober",&RoberCreate)); 328 PetscCall(PetscFunctionListAdd(&plist,"ce",&CECreate)); 329 PetscCall(PetscFunctionListAdd(&plist,"orego",&OregoCreate)); 330 PetscCall(PetscStrcpy(pname,"ce")); 331 332 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 333 Set runtime options 334 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 335 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options",""); 336 { 337 PetscCall(PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL)); 338 use_monitor = PETSC_FALSE; 339 PetscCall(PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL)); 340 PetscCall(PetscOptionsBool("-monitor_result","Display result","",use_result,&use_result,NULL)); 341 } 342 PetscOptionsEnd(); 343 344 /* Create the new problem */ 345 PetscCall(PetscNew(&problem)); 346 problem->comm = MPI_COMM_WORLD; 347 { 348 PetscErrorCode (*pcreate)(Problem); 349 350 PetscCall(PetscFunctionListFind(plist,pname,&pcreate)); 351 PetscCheck(pcreate,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"No problem '%s'",pname); 352 PetscCall((*pcreate)(problem)); 353 } 354 355 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 356 Create necessary matrix and vectors 357 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 358 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 359 PetscCall(MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE)); 360 PetscCall(MatSetFromOptions(A)); 361 PetscCall(MatSetUp(A)); 362 363 PetscCall(MatCreateVecs(A,&x,NULL)); 364 PetscCall(VecDuplicate(x,&r)); 365 366 mon.comm = PETSC_COMM_WORLD; 367 mon.problem = problem; 368 PetscCall(VecDuplicate(x,&mon.x)); 369 370 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 371 Create timestepping solver context 372 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 373 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 374 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 375 PetscCall(TSSetType(ts,TSROSW)); /* Rosenbrock-W */ 376 PetscCall(TSSetIFunction(ts,NULL,problem->function,problem->data)); 377 PetscCall(TSSetIJacobian(ts,A,A,problem->jacobian,problem->data)); 378 PetscCall(TSSetMaxTime(ts,problem->final_time)); 379 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 380 PetscCall(TSSetMaxStepRejections(ts,10)); 381 PetscCall(TSSetMaxSNESFailures(ts,-1)); /* unlimited */ 382 if (use_monitor) { 383 PetscCall(TSMonitorSet(ts,&MonitorError,&mon,NULL)); 384 } 385 386 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 387 Set initial conditions 388 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 389 PetscCall((*problem->solution)(0,x,problem->data)); 390 PetscCall(TSSetTimeStep(ts,.001)); 391 PetscCall(TSSetSolution(ts,x)); 392 393 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 394 Set runtime options 395 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 396 PetscCall(TSSetFromOptions(ts)); 397 398 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 399 Solve nonlinear system 400 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 401 PetscCall(TSSolve(ts,x)); 402 PetscCall(TSGetSolveTime(ts,&ftime)); 403 PetscCall(TSGetStepNumber(ts,&steps)); 404 PetscCall(TSGetSNESFailures(ts,&snesfails)); 405 PetscCall(TSGetStepRejections(ts,&rejects)); 406 PetscCall(TSGetSNESIterations(ts,&nonlinits)); 407 PetscCall(TSGetKSPIterations(ts,&linits)); 408 if (use_result) { 409 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"steps %" PetscInt_FMT " (%" PetscInt_FMT " rejected, %" PetscInt_FMT " SNES fails), ftime %g, nonlinits %" PetscInt_FMT ", linits %" PetscInt_FMT "\n",steps,rejects,snesfails,(double)ftime,nonlinits,linits)); 410 } 411 412 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 413 Free work space. All PETSc objects should be destroyed when they 414 are no longer needed. 415 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 416 PetscCall(MatDestroy(&A)); 417 PetscCall(VecDestroy(&x)); 418 PetscCall(VecDestroy(&r)); 419 PetscCall(VecDestroy(&mon.x)); 420 PetscCall(TSDestroy(&ts)); 421 if (problem->destroy) PetscCall((*problem->destroy)(problem)); 422 PetscCall(PetscFree(problem)); 423 PetscCall(PetscFunctionListDestroy(&plist)); 424 425 PetscCall(PetscFinalize()); 426 return 0; 427 } 428 429 /*TEST 430 431 test: 432 requires: !complex 433 args: -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex 434 435 test: 436 suffix: 2 437 requires: !single !complex 438 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 439 440 test: 441 suffix: 3 442 requires: !single !complex 443 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 444 445 test: 446 suffix: 4 447 448 test: 449 suffix: 5 450 args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists 451 452 TEST*/ 453