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