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 is one event set in this example, which 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. On reaching the limit on the number of ball bounces, 11 the TS run is requested to terminate from the PostEvent() callback. 12 */ 13 14 #include <petscts.h> 15 16 typedef struct { 17 PetscReal Cd; /* drag coefficient */ 18 PetscReal Cr; /* restitution coefficient */ 19 PetscInt bounces; 20 PetscInt maxbounces; 21 } AppCtx; 22 23 static PetscErrorCode Event(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx) 24 { 25 Vec V; 26 const PetscScalar *u; 27 28 PetscFunctionBeginUser; 29 /* Event for ball height */ 30 PetscCall(TS2GetSolution(ts, &U, &V)); 31 PetscCall(VecGetArrayRead(U, &u)); 32 fvalue[0] = PetscRealPart(u[0]); 33 PetscCall(VecRestoreArrayRead(U, &u)); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx) 38 { 39 AppCtx *app = (AppCtx *)ctx; 40 Vec V; 41 PetscScalar *u, *v; 42 PetscMPIInt rank; 43 PetscBool inflag = PETSC_FALSE, outflag; 44 45 PetscFunctionBeginUser; 46 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 47 if (nevents > 0) { 48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball hit the ground at t = %5.2f seconds\n", rank, (double)t)); 49 /* Set new initial conditions with attenuation Cr */ 50 PetscCall(TS2GetSolution(ts, &U, &V)); 51 PetscCall(VecGetArray(U, &u)); 52 PetscCall(VecGetArray(V, &v)); 53 u[0] = 0.0; 54 v[0] = -app->Cr * v[0]; 55 PetscCall(VecRestoreArray(U, &u)); 56 PetscCall(VecRestoreArray(V, &v)); 57 app->bounces++; 58 } 59 if (app->bounces >= app->maxbounces) { // 'app->bounces' may be different on different processes 60 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Processor [%d]: Ball bounced %" PetscInt_FMT " times\n", rank, app->bounces)); 61 inflag = PETSC_TRUE; // current process requested to terminate 62 } 63 PetscCallMPI(MPIU_Allreduce(&inflag, &outflag, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)ts))); 64 if (outflag) PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER)); // request TS to terminate, sync on all ranks 65 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(P, MAT_FINAL_ASSEMBLY)); 108 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 109 if (J != P) { 110 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 111 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 112 } 113 PetscFunctionReturn(PETSC_SUCCESS); 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[1]; 126 PetscBool terminate[1]; 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 PetscCall(TSSetEventHandler(ts, 1, direction, terminate, Event, PostEvent, &app)); // each process has one event-function defined 157 158 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 0, NULL, &J)); 159 PetscCall(MatSetFromOptions(J)); 160 PetscCall(MatSetUp(J)); 161 PetscCall(MatCreateVecs(J, NULL, &F)); 162 PetscCall(TSSetI2Function(ts, F, I2Function, &app)); 163 PetscCall(TSSetI2Jacobian(ts, J, J, I2Jacobian, &app)); 164 PetscCall(VecDestroy(&F)); 165 PetscCall(MatDestroy(&J)); 166 167 PetscCall(TSGetI2Jacobian(ts, &J, NULL, NULL, NULL)); 168 PetscCall(MatCreateVecs(J, &U, NULL)); 169 PetscCall(MatCreateVecs(J, &V, NULL)); 170 PetscCall(VecGetArray(U, &u)); 171 PetscCall(VecGetArray(V, &v)); 172 u[0] = 5.0 * rank; 173 v[0] = 20.0; 174 PetscCall(VecRestoreArray(U, &u)); 175 PetscCall(VecRestoreArray(V, &v)); 176 177 PetscCall(TS2SetSolution(ts, U, V)); 178 PetscCall(TSSetFromOptions(ts)); 179 PetscCall(TSSolve(ts, NULL)); 180 181 PetscCall(VecDestroy(&U)); 182 PetscCall(VecDestroy(&V)); 183 PetscCall(TSDestroy(&ts)); 184 185 PetscCall(PetscFinalize()); 186 return 0; 187 } 188 189 /*TEST 190 191 test: 192 suffix: a 193 args: -ts_alpha_radius {{1.0 0.5}} -ts_max_time 50 194 output_file: output/ex44.out 195 196 test: 197 suffix: b 198 args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50 199 output_file: output/ex44.out 200 201 test: 202 suffix: bmf 203 args: -snes_mf_operator -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50 204 output_file: output/ex44.out 205 206 test: 207 suffix: 2 208 nsize: 2 209 args: -ts_rtol 0 -ts_atol 1e-1 -ts_adapt_type basic -ts_max_time 50 210 output_file: output/ex44_2.out 211 filter: sort -b 212 filter_output: sort -b 213 214 test: 215 requires: !single 216 args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor 217 args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false 218 219 TEST*/ 220