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