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 CHKERRQ(VecGetArrayRead(U,&u)); 28 fvalue[0] = u[0]; 29 /* Event for number of bounces */ 30 fvalue[1] = app->maxbounces - app->nbounces; 31 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecGetArray(U,&u)); 45 u[0] = 0.0; 46 u[1] = -0.9*u[1]; 47 CHKERRQ(VecRestoreArray(U,&u)); 48 } else if (event_list[0] == 1) { 49 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Ball bounced %D 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 CHKERRQ(VecGetArrayRead(U,&u)); 66 CHKERRQ(VecGetArray(F,&f)); 67 68 f[0] = u[1]; 69 f[1] = - 9.8; 70 71 CHKERRQ(VecRestoreArrayRead(U,&u)); 72 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 91 92 CHKERRQ(VecRestoreArrayRead(U,&u)); 93 94 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 95 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 96 if (A != B) { 97 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 98 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 114 CHKERRQ(VecGetArrayRead(Udot,&udot)); 115 CHKERRQ(VecGetArray(F,&f)); 116 117 f[0] = udot[0] - u[1]; 118 f[1] = udot[1] + 9.8; 119 120 CHKERRQ(VecRestoreArrayRead(U,&u)); 121 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 122 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 137 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 142 143 CHKERRQ(VecRestoreArrayRead(U,&u)); 144 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 145 146 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 147 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 148 if (A != B) { 149 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 150 CHKERRQ(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 PetscErrorCode ierr; 160 PetscMPIInt size; 161 PetscInt n = 2; 162 PetscScalar *u; 163 AppCtx app; 164 PetscInt direction[2]; 165 PetscBool terminate[2]; 166 PetscBool rhs_form=PETSC_FALSE,hist=PETSC_TRUE; 167 TSAdapt adapt; 168 169 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 Initialize program 171 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 172 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 173 CHKERRMPI(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 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex40 options","");CHKERRQ(ierr); 179 CHKERRQ(PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL)); 180 CHKERRQ(PetscOptionsBool("-test_adapthistory","","",hist,&hist,NULL)); 181 ierr = PetscOptionsEnd();CHKERRQ(ierr); 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Create timestepping solver context 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 187 CHKERRQ(TSSetType(ts,TSROSW)); 188 189 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 190 Set ODE routines 191 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 192 CHKERRQ(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 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL)); 199 if (rhs_form) { 200 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 201 CHKERRQ(TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL)); 202 } else { 203 Mat A; /* Jacobian matrix */ 204 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 205 CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 206 CHKERRQ(MatSetType(A,MATDENSE)); 207 CHKERRQ(MatSetFromOptions(A)); 208 CHKERRQ(MatSetUp(A)); 209 CHKERRQ(TSSetIFunction(ts,NULL,IFunction,NULL)); 210 CHKERRQ(TSSetIJacobian(ts,A,A,IJacobian,NULL)); 211 CHKERRQ(MatDestroy(&A)); 212 } 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Set initial conditions 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&U)); 218 CHKERRQ(VecSetSizes(U,n,PETSC_DETERMINE)); 219 CHKERRQ(VecSetUp(U)); 220 CHKERRQ(VecGetArray(U,&u)); 221 u[0] = 0.0; 222 u[1] = 20.0; 223 CHKERRQ(VecRestoreArray(U,&u)); 224 CHKERRQ(TSSetSolution(ts,U)); 225 226 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227 Set solver options 228 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 229 if (hist) CHKERRQ(TSSetSaveTrajectory(ts)); 230 CHKERRQ(TSSetMaxTime(ts,30.0)); 231 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 232 CHKERRQ(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 CHKERRQ(TSGetAdapt(ts,&adapt)); 237 CHKERRQ(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 CHKERRQ(TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app)); 243 244 CHKERRQ(TSSetFromOptions(ts)); 245 246 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 247 Run timestepping solver 248 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 249 CHKERRQ(TSSolve(ts,U)); 250 251 if (hist) { /* replay following history */ 252 TSTrajectory tj; 253 PetscReal tf,t0,dt; 254 255 app.nbounces = 0; 256 CHKERRQ(TSGetTime(ts,&tf)); 257 CHKERRQ(TSSetMaxTime(ts,tf)); 258 CHKERRQ(TSSetStepNumber(ts,0)); 259 CHKERRQ(TSRestartStep(ts)); 260 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 261 CHKERRQ(TSSetFromOptions(ts)); 262 CHKERRQ(TSGetAdapt(ts,&adapt)); 263 CHKERRQ(TSAdaptSetType(adapt,TSADAPTHISTORY)); 264 CHKERRQ(TSGetTrajectory(ts,&tj)); 265 CHKERRQ(TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE)); 266 CHKERRQ(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 CHKERRQ(TSAdaptSetType(adapt,TSADAPTBASIC)); 270 CHKERRQ(TSAdaptSetStepLimits(adapt,0.0,0.5)); 271 CHKERRQ(TSSetFromOptions(ts)); 272 #endif 273 CHKERRQ(TSSetTime(ts,t0)); 274 CHKERRQ(TSSetTimeStep(ts,dt)); 275 CHKERRQ(TSResetTrajectory(ts)); 276 CHKERRQ(VecGetArray(U,&u)); 277 u[0] = 0.0; 278 u[1] = 20.0; 279 CHKERRQ(VecRestoreArray(U,&u)); 280 CHKERRQ(TSSolve(ts,U)); 281 } 282 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 283 Free work space. All PETSc objects should be destroyed when they are no longer needed. 284 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 285 CHKERRQ(VecDestroy(&U)); 286 CHKERRQ(TSDestroy(&ts)); 287 288 CHKERRQ(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