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