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