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