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