1 static char help[] = "Parallel bouncing ball example formulated as a second-order system to test TS event feature.\n"; 2 3 /* 4 The dynamics of the bouncing ball with drag coefficient Cd is described by the ODE 5 6 u_tt = -9.8 - 1/2 Cd (u_t)^2 sign(u_t) 7 8 There are two events set in this example. The first one checks for the ball hitting the 9 ground (u = 0). Every time the ball hits the ground, its velocity u_t is attenuated by 10 a restitution coefficient Cr. The second event sets a limit on the number of ball bounces. 11 */ 12 13 #include <petscts.h> 14 15 typedef struct { 16 PetscReal Cd; /* drag coefficient */ 17 PetscReal Cr; /* restitution coefficient */ 18 PetscInt bounces; 19 PetscInt maxbounces; 20 } AppCtx; 21 22 static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscScalar *fvalue, void *ctx) 23 { 24 AppCtx *app = (AppCtx *)ctx; 25 Vec V; 26 const PetscScalar *u, *v; 27 28 PetscFunctionBeginUser; 29 /* Event for ball height */ 30 PetscCall(TS2GetSolution(ts, &U, &V)); 31 PetscCall(VecGetArrayRead(U, &u)); 32 PetscCall(VecGetArrayRead(V, &v)); 33 fvalue[0] = u[0]; 34 /* Event for number of bounces */ 35 fvalue[1] = app->maxbounces - app->bounces; 36 PetscCall(VecRestoreArrayRead(U, &u)); 37 PetscCall(VecRestoreArrayRead(V, &v)); 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx) 42 { 43 AppCtx *app = (AppCtx *)ctx; 44 Vec V; 45 PetscScalar *u, *v; 46 PetscMPIInt rank; 47 48 PetscFunctionBeginUser; 49 if (!nevents) PetscFunctionReturn(0); 50 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 51 if (event_list[0] == 0) { 52 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t)); 53 /* Set new initial conditions with .9 attenuation */ 54 PetscCall(TS2GetSolution(ts, &U, &V)); 55 PetscCall(VecGetArray(U, &u)); 56 PetscCall(VecGetArray(V, &v)); 57 u[0] = 0.0; 58 v[0] = -app->Cr * v[0]; 59 PetscCall(VecRestoreArray(U, &u)); 60 PetscCall(VecRestoreArray(V, &v)); 61 app->bounces++; 62 } else if (event_list[0] == 1) { 63 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces)); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 static PetscErrorCode I2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F, void *ctx) 69 { 70 AppCtx *app = (AppCtx *)ctx; 71 const PetscScalar *u, *v, *a; 72 PetscScalar Res, *f; 73 74 PetscFunctionBeginUser; 75 PetscCall(VecGetArrayRead(U, &u)); 76 PetscCall(VecGetArrayRead(V, &v)); 77 PetscCall(VecGetArrayRead(A, &a)); 78 Res = a[0] + 9.8 + 0.5 * app->Cd * v[0] * v[0] * PetscSignReal(PetscRealPart(v[0])); 79 PetscCall(VecRestoreArrayRead(U, &u)); 80 PetscCall(VecRestoreArrayRead(V, &v)); 81 PetscCall(VecRestoreArrayRead(A, &a)); 82 83 PetscCall(VecGetArray(F, &f)); 84 f[0] = Res; 85 PetscCall(VecRestoreArray(F, &f)); 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode I2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx) 90 { 91 AppCtx *app = (AppCtx *)ctx; 92 const PetscScalar *u, *v, *a; 93 PetscInt i; 94 PetscScalar Jac; 95 96 PetscFunctionBeginUser; 97 PetscCall(VecGetArrayRead(U, &u)); 98 PetscCall(VecGetArrayRead(V, &v)); 99 PetscCall(VecGetArrayRead(A, &a)); 100 Jac = shiftA + shiftV * app->Cd * v[0]; 101 PetscCall(VecRestoreArrayRead(U, &u)); 102 PetscCall(VecRestoreArrayRead(V, &v)); 103 PetscCall(VecRestoreArrayRead(A, &a)); 104 105 PetscCall(MatGetOwnershipRange(P, &i, NULL)); 106 PetscCall(MatSetValue(P, i, i, Jac, INSERT_VALUES)); 107 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 108 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 109 if (J != P) { 110 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 111 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 112 } 113 PetscFunctionReturn(0); 114 } 115 116 int main(int argc, char **argv) 117 { 118 TS ts; /* ODE integrator */ 119 Vec U, V; /* solution will be stored here */ 120 Vec F; /* residual vector */ 121 Mat J; /* Jacobian matrix */ 122 PetscMPIInt rank; 123 PetscScalar *u, *v; 124 AppCtx app; 125 PetscInt direction[2]; 126 PetscBool terminate[2]; 127 TSAdapt adapt; 128 129 PetscFunctionBeginUser; 130 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 131 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 132 133 app.Cd = 0.0; 134 app.Cr = 0.9; 135 app.bounces = 0; 136 app.maxbounces = 10; 137 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex44 options", ""); 138 PetscCall(PetscOptionsReal("-Cd", "Drag coefficient", "", app.Cd, &app.Cd, NULL)); 139 PetscCall(PetscOptionsReal("-Cr", "Restitution coefficient", "", app.Cr, &app.Cr, NULL)); 140 PetscCall(PetscOptionsInt("-maxbounces", "Maximum number of bounces", "", app.maxbounces, &app.maxbounces, NULL)); 141 PetscOptionsEnd(); 142 143 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 144 /*PetscCall(TSSetSaveTrajectory(ts));*/ 145 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 146 PetscCall(TSSetType(ts, TSALPHA2)); 147 148 PetscCall(TSSetMaxTime(ts, PETSC_INFINITY)); 149 PetscCall(TSSetTimeStep(ts, 0.1)); 150 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 151 PetscCall(TSGetAdapt(ts, &adapt)); 152 PetscCall(TSAdaptSetStepLimits(adapt, 0.0, 0.5)); 153 154 direction[0] = -1; 155 terminate[0] = PETSC_FALSE; 156 direction[1] = -1; 157 terminate[1] = PETSC_TRUE; 158 PetscCall(TSSetEventHandler(ts, 2, direction, terminate, Event, PostEvent, &app)); 159 160 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J)); 161 PetscCall(MatSetFromOptions(J)); 162 PetscCall(MatSetUp(J)); 163 PetscCall(MatCreateVecs(J, NULL, &F)); 164 PetscCall(TSSetI2Function(ts, F, I2Function, &app)); 165 PetscCall(TSSetI2Jacobian(ts, J, J, I2Jacobian, &app)); 166 PetscCall(VecDestroy(&F)); 167 PetscCall(MatDestroy(&J)); 168 169 PetscCall(TSGetI2Jacobian(ts, &J, NULL, NULL, NULL)); 170 PetscCall(MatCreateVecs(J, &U, NULL)); 171 PetscCall(MatCreateVecs(J, &V, NULL)); 172 PetscCall(VecGetArray(U, &u)); 173 PetscCall(VecGetArray(V, &v)); 174 u[0] = 5.0 * rank; 175 v[0] = 20.0; 176 PetscCall(VecRestoreArray(U, &u)); 177 PetscCall(VecRestoreArray(V, &v)); 178 179 PetscCall(TS2SetSolution(ts, U, V)); 180 PetscCall(TSSetFromOptions(ts)); 181 PetscCall(TSSolve(ts, NULL)); 182 183 PetscCall(VecDestroy(&U)); 184 PetscCall(VecDestroy(&V)); 185 PetscCall(TSDestroy(&ts)); 186 187 PetscCall(PetscFinalize()); 188 return 0; 189 } 190 191 /*TEST 192 193 test: 194 suffix: a 195 args: -ts_alpha_radius {{1.0 0.5}} 196 output_file: output/ex44.out 197 198 test: 199 suffix: b 200 args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic 201 output_file: output/ex44.out 202 203 test: 204 suffix: 2 205 nsize: 2 206 args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic 207 output_file: output/ex44_2.out 208 filter: sort -b 209 filter_output: sort -b 210 211 test: 212 requires: !single 213 args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor 214 args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false 215 216 TEST*/ 217