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