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 CHKERRQ(VecGetArrayRead(U,&u)); 24 fvalue[0] = u[0]; 25 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 36 if (nevents) { 37 CHKERRQ(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 CHKERRQ(VecGetArray(U,&u)); 40 u[0] = 1.0*rank; 41 u[1] = -0.9*u[1]; 42 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 58 CHKERRQ(VecGetArray(F,&f)); 59 60 f[0] = u[1]; 61 f[1] = - 9.8; 62 63 CHKERRQ(VecRestoreArrayRead(U,&u)); 64 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 79 80 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 86 87 CHKERRQ(VecRestoreArrayRead(U,&u)); 88 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 89 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 90 if (A != B) { 91 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 92 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 108 CHKERRQ(VecGetArrayRead(Udot,&udot)); 109 CHKERRQ(VecGetArray(F,&f)); 110 111 f[0] = udot[0] - u[1]; 112 f[1] = udot[1] + 9.8; 113 114 CHKERRQ(VecRestoreArrayRead(U,&u)); 115 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 116 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 131 CHKERRQ(VecGetArrayRead(Udot,&udot)); 132 133 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 139 140 CHKERRQ(VecRestoreArrayRead(U,&u)); 141 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 142 143 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 144 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 145 if (A != B) { 146 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 147 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 168 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 169 170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171 Create timestepping solver context 172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 173 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 174 CHKERRQ(TSSetType(ts,TSROSW)); 175 176 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177 Set ODE routines 178 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179 CHKERRQ(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 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL)); 186 if (rhs_form) { 187 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 188 CHKERRQ(TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL)); 189 } else { 190 CHKERRQ(TSSetIFunction(ts,NULL,IFunction,NULL)); 191 CHKERRQ(TSSetIJacobian(ts,NULL,NULL,IJacobian,NULL)); 192 } 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Set initial conditions 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&U)); 198 CHKERRQ(VecSetSizes(U,n,PETSC_DETERMINE)); 199 CHKERRQ(VecSetUp(U)); 200 CHKERRQ(VecGetArray(U,&u)); 201 u[0] = 1.0*rank; 202 u[1] = 20.0; 203 CHKERRQ(VecRestoreArray(U,&u)); 204 CHKERRQ(TSSetSolution(ts,U)); 205 206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207 Set solver options 208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 209 CHKERRQ(TSSetSaveTrajectory(ts)); 210 CHKERRQ(TSSetMaxTime(ts,30.0)); 211 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 212 CHKERRQ(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 CHKERRQ(TSGetAdapt(ts,&adapt)); 217 CHKERRQ(TSAdaptSetType(adapt,TSADAPTBASIC)); 218 CHKERRQ(TSAdaptSetStepLimits(adapt,0.0,0.5)); 219 220 /* Set direction and terminate flag for the event */ 221 CHKERRQ(TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL)); 222 223 CHKERRQ(TSSetFromOptions(ts)); 224 225 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 226 Run timestepping solver 227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 228 CHKERRQ(TSSolve(ts,U)); 229 230 if (hist) { /* replay following history */ 231 TSTrajectory tj; 232 PetscReal tf,t0,dt; 233 234 CHKERRQ(TSGetTime(ts,&tf)); 235 CHKERRQ(TSSetMaxTime(ts,tf)); 236 CHKERRQ(TSSetStepNumber(ts,0)); 237 CHKERRQ(TSRestartStep(ts)); 238 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 239 CHKERRQ(TSSetFromOptions(ts)); 240 CHKERRQ(TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL)); 241 CHKERRQ(TSGetAdapt(ts,&adapt)); 242 CHKERRQ(TSAdaptSetType(adapt,TSADAPTHISTORY)); 243 CHKERRQ(TSGetTrajectory(ts,&tj)); 244 CHKERRQ(TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE)); 245 CHKERRQ(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 CHKERRQ(TSAdaptSetType(adapt,TSADAPTBASIC)); 249 CHKERRQ(TSAdaptSetStepLimits(adapt,0.0,0.5)); 250 CHKERRQ(TSSetFromOptions(ts)); 251 #endif 252 CHKERRQ(TSSetTime(ts,t0)); 253 CHKERRQ(TSSetTimeStep(ts,dt)); 254 CHKERRQ(TSResetTrajectory(ts)); 255 CHKERRQ(VecGetArray(U,&u)); 256 u[0] = 1.0*rank; 257 u[1] = 20.0; 258 CHKERRQ(VecRestoreArray(U,&u)); 259 CHKERRQ(TSSolve(ts,U)); 260 } 261 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 262 Free work space. All PETSc objects should be destroyed when they are no longer needed. 263 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 264 CHKERRQ(VecDestroy(&U)); 265 CHKERRQ(TSDestroy(&ts)); 266 267 CHKERRQ(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