1 static char help[] = "Serial bouncing ball example to test TS event feature.\n"; 2 3 /* 4 The dynamics of the bouncing ball is described by the ODE 5 u1_t = u2 6 u2_t = -9.8 7 8 There are two events set in this example. The first one checks for the ball hitting the 9 ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by 10 a factor of 0.9. The second event sets a limit on the number of ball bounces. 11 */ 12 13 #include <petscts.h> 14 15 typedef struct { 16 PetscInt maxbounces; 17 PetscInt nbounces; 18 } AppCtx; 19 20 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx) { 21 AppCtx *app = (AppCtx *)ctx; 22 const PetscScalar *u; 23 24 PetscFunctionBeginUser; 25 /* Event for ball height */ 26 PetscCall(VecGetArrayRead(U, &u)); 27 fvalue[0] = u[0]; 28 /* Event for number of bounces */ 29 fvalue[1] = app->maxbounces - app->nbounces; 30 PetscCall(VecRestoreArrayRead(U, &u)); 31 PetscFunctionReturn(0); 32 } 33 34 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx) { 35 AppCtx *app = (AppCtx *)ctx; 36 PetscScalar *u; 37 38 PetscFunctionBeginUser; 39 if (event_list[0] == 0) { 40 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds\n", (double)t)); 41 /* Set new initial conditions with .9 attenuation */ 42 PetscCall(VecGetArray(U, &u)); 43 u[0] = 0.0; 44 u[1] = -0.9 * u[1]; 45 PetscCall(VecRestoreArray(U, &u)); 46 } else if (event_list[0] == 1) { 47 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball bounced %" PetscInt_FMT " times\n", app->nbounces)); 48 } 49 app->nbounces++; 50 PetscFunctionReturn(0); 51 } 52 53 /* 54 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U) 55 */ 56 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx) { 57 PetscScalar *f; 58 const PetscScalar *u; 59 60 PetscFunctionBeginUser; 61 /* The next three lines allow us to access the entries of the vectors directly */ 62 PetscCall(VecGetArrayRead(U, &u)); 63 PetscCall(VecGetArray(F, &f)); 64 65 f[0] = u[1]; 66 f[1] = -9.8; 67 68 PetscCall(VecRestoreArrayRead(U, &u)); 69 PetscCall(VecRestoreArray(F, &f)); 70 PetscFunctionReturn(0); 71 } 72 73 /* 74 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian. 75 */ 76 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx) { 77 PetscInt rowcol[] = {0, 1}; 78 PetscScalar J[2][2]; 79 const PetscScalar *u; 80 81 PetscFunctionBeginUser; 82 PetscCall(VecGetArrayRead(U, &u)); 83 84 J[0][0] = 0.0; 85 J[0][1] = 1.0; 86 J[1][0] = 0.0; 87 J[1][1] = 0.0; 88 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 89 90 PetscCall(VecRestoreArrayRead(U, &u)); 91 92 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 93 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 94 if (A != B) { 95 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 96 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 97 } 98 PetscFunctionReturn(0); 99 } 100 101 /* 102 Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0 103 */ 104 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) { 105 PetscScalar *f; 106 const PetscScalar *u, *udot; 107 108 PetscFunctionBeginUser; 109 /* The next three lines allow us to access the entries of the vectors directly */ 110 PetscCall(VecGetArrayRead(U, &u)); 111 PetscCall(VecGetArrayRead(Udot, &udot)); 112 PetscCall(VecGetArray(F, &f)); 113 114 f[0] = udot[0] - u[1]; 115 f[1] = udot[1] + 9.8; 116 117 PetscCall(VecRestoreArrayRead(U, &u)); 118 PetscCall(VecRestoreArrayRead(Udot, &udot)); 119 PetscCall(VecRestoreArray(F, &f)); 120 PetscFunctionReturn(0); 121 } 122 123 /* 124 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 125 */ 126 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx) { 127 PetscInt rowcol[] = {0, 1}; 128 PetscScalar J[2][2]; 129 const PetscScalar *u, *udot; 130 131 PetscFunctionBeginUser; 132 PetscCall(VecGetArrayRead(U, &u)); 133 PetscCall(VecGetArrayRead(Udot, &udot)); 134 135 J[0][0] = a; 136 J[0][1] = -1.0; 137 J[1][0] = 0.0; 138 J[1][1] = a; 139 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 140 141 PetscCall(VecRestoreArrayRead(U, &u)); 142 PetscCall(VecRestoreArrayRead(Udot, &udot)); 143 144 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 145 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 146 if (A != B) { 147 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 int main(int argc, char **argv) { 154 TS ts; /* ODE integrator */ 155 Vec U; /* solution will be stored here */ 156 PetscMPIInt size; 157 PetscInt n = 2; 158 PetscScalar *u; 159 AppCtx app; 160 PetscInt direction[2]; 161 PetscBool terminate[2]; 162 PetscBool rhs_form = PETSC_FALSE, hist = PETSC_TRUE; 163 TSAdapt adapt; 164 165 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 166 Initialize program 167 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 168 PetscFunctionBeginUser; 169 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 170 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 171 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 172 173 app.nbounces = 0; 174 app.maxbounces = 10; 175 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex40 options", ""); 176 PetscCall(PetscOptionsInt("-maxbounces", "", "", app.maxbounces, &app.maxbounces, NULL)); 177 PetscCall(PetscOptionsBool("-test_adapthistory", "", "", hist, &hist, NULL)); 178 PetscOptionsEnd(); 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Create timestepping solver context 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 184 PetscCall(TSSetType(ts, TSROSW)); 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Set ODE routines 188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 189 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 190 /* Users are advised against the following branching and code duplication. 191 For problems without a mass matrix like the one at hand, the RHSFunction 192 (and companion RHSJacobian) interface is enough to support both explicit 193 and implicit timesteppers. This tutorial example also deals with the 194 IFunction/IJacobian interface for demonstration and testing purposes. */ 195 PetscCall(PetscOptionsGetBool(NULL, NULL, "-rhs-form", &rhs_form, NULL)); 196 if (rhs_form) { 197 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 198 PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, NULL)); 199 } else { 200 Mat A; /* Jacobian matrix */ 201 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 202 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 203 PetscCall(MatSetType(A, MATDENSE)); 204 PetscCall(MatSetFromOptions(A)); 205 PetscCall(MatSetUp(A)); 206 PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL)); 207 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL)); 208 PetscCall(MatDestroy(&A)); 209 } 210 211 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 212 Set initial conditions 213 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 214 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 215 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 216 PetscCall(VecSetUp(U)); 217 PetscCall(VecGetArray(U, &u)); 218 u[0] = 0.0; 219 u[1] = 20.0; 220 PetscCall(VecRestoreArray(U, &u)); 221 PetscCall(TSSetSolution(ts, U)); 222 223 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 224 Set solver options 225 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 226 if (hist) PetscCall(TSSetSaveTrajectory(ts)); 227 PetscCall(TSSetMaxTime(ts, 30.0)); 228 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 229 PetscCall(TSSetTimeStep(ts, 0.1)); 230 /* The adaptive time step controller could take very large timesteps resulting in 231 the same event occurring multiple times in the same interval. A maximum step size 232 limit is enforced here to avoid this issue. */ 233 PetscCall(TSGetAdapt(ts, &adapt)); 234 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5)); 235 236 /* Set directions and terminate flags for the two events */ 237 direction[0] = -1; 238 direction[1] = -1; 239 terminate[0] = PETSC_FALSE; 240 terminate[1] = PETSC_TRUE; 241 PetscCall(TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&app)); 242 243 PetscCall(TSSetFromOptions(ts)); 244 245 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 246 Run timestepping solver 247 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 248 PetscCall(TSSolve(ts, U)); 249 250 if (hist) { /* replay following history */ 251 TSTrajectory tj; 252 PetscReal tf, t0, dt; 253 254 app.nbounces = 0; 255 PetscCall(TSGetTime(ts, &tf)); 256 PetscCall(TSSetMaxTime(ts, tf)); 257 PetscCall(TSSetStepNumber(ts, 0)); 258 PetscCall(TSRestartStep(ts)); 259 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 260 PetscCall(TSSetFromOptions(ts)); 261 PetscCall(TSGetAdapt(ts, &adapt)); 262 PetscCall(TSAdaptSetType(adapt, TSADAPTHISTORY)); 263 PetscCall(TSGetTrajectory(ts, &tj)); 264 PetscCall(TSAdaptHistorySetTrajectory(adapt, tj, PETSC_FALSE)); 265 PetscCall(TSAdaptHistoryGetStep(adapt, 0, &t0, &dt)); 266 /* this example fails with single (or smaller) precision */ 267 #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16) 268 PetscCall(TSAdaptSetType(adapt, TSADAPTBASIC)); 269 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5)); 270 PetscCall(TSSetFromOptions(ts)); 271 #endif 272 PetscCall(TSSetTime(ts, t0)); 273 PetscCall(TSSetTimeStep(ts, dt)); 274 PetscCall(TSResetTrajectory(ts)); 275 PetscCall(VecGetArray(U, &u)); 276 u[0] = 0.0; 277 u[1] = 20.0; 278 PetscCall(VecRestoreArray(U, &u)); 279 PetscCall(TSSolve(ts, U)); 280 } 281 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 282 Free work space. All PETSc objects should be destroyed when they are no longer needed. 283 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 284 PetscCall(VecDestroy(&U)); 285 PetscCall(TSDestroy(&ts)); 286 287 PetscCall(PetscFinalize()); 288 return 0; 289 } 290 291 /*TEST 292 293 test: 294 suffix: a 295 args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir 296 output_file: output/ex40.out 297 298 test: 299 suffix: b 300 args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir 301 output_file: output/ex40.out 302 303 test: 304 suffix: c 305 args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir 306 output_file: output/ex40.out 307 308 test: 309 suffix: d 310 args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir 311 output_file: output/ex40.out 312 313 test: 314 suffix: e 315 args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir 316 output_file: output/ex40.out 317 318 test: 319 suffix: f 320 args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir 321 output_file: output/ex40.out 322 323 test: 324 suffix: g 325 args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir 326 output_file: output/ex40.out 327 328 test: 329 suffix: h 330 args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir 331 output_file: output/ex40.out 332 333 test: 334 suffix: i 335 args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir 336 output_file: output/ex40.out 337 338 test: 339 suffix: j 340 args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir 341 output_file: output/ex40.out 342 343 test: 344 suffix: k 345 args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir 346 output_file: output/ex40.out 347 348 test: 349 suffix: l 350 args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir 351 args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01 352 output_file: output/ex40.out 353 354 test: 355 suffix: m 356 args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false 357 args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}} 358 359 test: 360 requires: !single 361 suffix: n 362 args: -test_adapthistory false 363 args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view 364 args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor 365 args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false 366 367 test: 368 requires: !single 369 suffix: o 370 args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir 371 output_file: output/ex40.out 372 TEST*/ 373