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 const PetscScalar *u; 19 20 PetscFunctionBeginUser; 21 /* Event for ball height */ 22 PetscCall(VecGetArrayRead(U, &u)); 23 fvalue[0] = u[0]; 24 PetscCall(VecRestoreArrayRead(U, &u)); 25 PetscFunctionReturn(0); 26 } 27 28 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx) { 29 PetscScalar *u; 30 PetscMPIInt rank; 31 32 PetscFunctionBeginUser; 33 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 34 if (nevents) { 35 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n", (double)t, rank)); 36 /* Set new initial conditions with .9 attenuation */ 37 PetscCall(VecGetArray(U, &u)); 38 u[0] = 1.0 * rank; 39 u[1] = -0.9 * u[1]; 40 PetscCall(VecRestoreArray(U, &u)); 41 } 42 PetscFunctionReturn(0); 43 } 44 45 /* 46 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U) 47 */ 48 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx) { 49 PetscScalar *f; 50 const PetscScalar *u; 51 52 PetscFunctionBeginUser; 53 /* The next three lines allow us to access the entries of the vectors directly */ 54 PetscCall(VecGetArrayRead(U, &u)); 55 PetscCall(VecGetArray(F, &f)); 56 57 f[0] = u[1]; 58 f[1] = -9.8; 59 60 PetscCall(VecRestoreArrayRead(U, &u)); 61 PetscCall(VecRestoreArray(F, &f)); 62 PetscFunctionReturn(0); 63 } 64 65 /* 66 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian. 67 */ 68 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx) { 69 PetscInt rowcol[2], rstart; 70 PetscScalar J[2][2]; 71 const PetscScalar *u; 72 73 PetscFunctionBeginUser; 74 PetscCall(VecGetArrayRead(U, &u)); 75 76 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 77 rowcol[0] = rstart; 78 rowcol[1] = rstart + 1; 79 80 J[0][0] = 0.0; 81 J[0][1] = 1.0; 82 J[1][0] = 0.0; 83 J[1][1] = 0.0; 84 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 85 86 PetscCall(VecRestoreArrayRead(U, &u)); 87 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 88 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 89 if (A != B) { 90 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 91 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 92 } 93 PetscFunctionReturn(0); 94 } 95 96 /* 97 Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0 98 */ 99 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) { 100 PetscScalar *f; 101 const PetscScalar *u, *udot; 102 103 PetscFunctionBeginUser; 104 /* The next three lines allow us to access the entries of the vectors directly */ 105 PetscCall(VecGetArrayRead(U, &u)); 106 PetscCall(VecGetArrayRead(Udot, &udot)); 107 PetscCall(VecGetArray(F, &f)); 108 109 f[0] = udot[0] - u[1]; 110 f[1] = udot[1] + 9.8; 111 112 PetscCall(VecRestoreArrayRead(U, &u)); 113 PetscCall(VecRestoreArrayRead(Udot, &udot)); 114 PetscCall(VecRestoreArray(F, &f)); 115 PetscFunctionReturn(0); 116 } 117 118 /* 119 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 120 */ 121 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx) { 122 PetscInt rowcol[2], rstart; 123 PetscScalar J[2][2]; 124 const PetscScalar *u, *udot; 125 126 PetscFunctionBeginUser; 127 PetscCall(VecGetArrayRead(U, &u)); 128 PetscCall(VecGetArrayRead(Udot, &udot)); 129 130 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 131 rowcol[0] = rstart; 132 rowcol[1] = rstart + 1; 133 134 J[0][0] = a; 135 J[0][1] = -1.0; 136 J[1][0] = 0.0; 137 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 TS ts; /* ODE integrator */ 154 Vec U; /* solution will be stored here */ 155 PetscMPIInt rank; 156 PetscInt n = 2; 157 PetscScalar *u; 158 PetscInt direction = -1; 159 PetscBool terminate = PETSC_FALSE; 160 PetscBool rhs_form = PETSC_FALSE, hist = PETSC_TRUE; 161 TSAdapt adapt; 162 163 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164 Initialize program 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 166 PetscFunctionBeginUser; 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