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