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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 172 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 173 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 174 175 app.nbounces = 0; 176 app.maxbounces = 10; 177 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex40 options",""); 178 PetscCall(PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL)); 179 PetscCall(PetscOptionsBool("-test_adapthistory","","",hist,&hist,NULL)); 180 PetscOptionsEnd(); 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Create timestepping solver context 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 186 PetscCall(TSSetType(ts,TSROSW)); 187 188 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 189 Set ODE routines 190 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 191 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 192 /* Users are advised against the following branching and code duplication. 193 For problems without a mass matrix like the one at hand, the RHSFunction 194 (and companion RHSJacobian) interface is enough to support both explicit 195 and implicit timesteppers. This tutorial example also deals with the 196 IFunction/IJacobian interface for demonstration and testing purposes. */ 197 PetscCall(PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL)); 198 if (rhs_form) { 199 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 200 PetscCall(TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL)); 201 } else { 202 Mat A; /* Jacobian matrix */ 203 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 204 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 205 PetscCall(MatSetType(A,MATDENSE)); 206 PetscCall(MatSetFromOptions(A)); 207 PetscCall(MatSetUp(A)); 208 PetscCall(TSSetIFunction(ts,NULL,IFunction,NULL)); 209 PetscCall(TSSetIJacobian(ts,A,A,IJacobian,NULL)); 210 PetscCall(MatDestroy(&A)); 211 } 212 213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 214 Set initial conditions 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 PetscCall(VecCreate(PETSC_COMM_WORLD,&U)); 217 PetscCall(VecSetSizes(U,n,PETSC_DETERMINE)); 218 PetscCall(VecSetUp(U)); 219 PetscCall(VecGetArray(U,&u)); 220 u[0] = 0.0; 221 u[1] = 20.0; 222 PetscCall(VecRestoreArray(U,&u)); 223 PetscCall(TSSetSolution(ts,U)); 224 225 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 226 Set solver options 227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 228 if (hist) PetscCall(TSSetSaveTrajectory(ts)); 229 PetscCall(TSSetMaxTime(ts,30.0)); 230 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 231 PetscCall(TSSetTimeStep(ts,0.1)); 232 /* The adaptive time step controller could take very large timesteps resulting in 233 the same event occurring multiple times in the same interval. A maximum step size 234 limit is enforced here to avoid this issue. */ 235 PetscCall(TSGetAdapt(ts,&adapt)); 236 PetscCall(TSAdaptSetStepLimits(adapt,0.0,0.5)); 237 238 /* Set directions and terminate flags for the two events */ 239 direction[0] = -1; direction[1] = -1; 240 terminate[0] = PETSC_FALSE; 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