Actual source code: ex40.c
petsc-3.7.5 2017-01-01
1: static char help[] = "Serial bouncing ball example to test TS event feature.\n";
3: /*
4: The dynamics of the bouncing ball is described by the ODE
5: u1_t = u2
6: u2_t = -9.8
8: There are two events set in this example. The first one checks for the ball hitting the
9: ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
10: a factor of 0.9. The second event sets a limit on the number of ball bounces.
11: */
13: #include <petscts.h>
15: typedef struct {
16: PetscInt maxbounces;
17: PetscInt nbounces;
18: } AppCtx;
22: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
23: {
24: AppCtx *app=(AppCtx*)ctx;
25: PetscErrorCode ierr;
26: const PetscScalar *u;
29: /* Event for ball height */
30: VecGetArrayRead(U,&u);
31: fvalue[0] = u[0];
32: /* Event for number of bounces */
33: fvalue[1] = app->maxbounces - app->nbounces;
34: VecRestoreArrayRead(U,&u);
35: return(0);
36: }
40: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
41: {
42: AppCtx *app=(AppCtx*)ctx;
44: PetscScalar *u;
47: if (event_list[0] == 0) {
48: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds\n",(double)t);
49: /* Set new initial conditions with .9 attenuation */
50: VecGetArray(U,&u);
51: u[0] = 0.0;
52: u[1] = -0.9*u[1];
53: VecRestoreArray(U,&u);
54: } else if (event_list[0] == 1) {
55: PetscPrintf(PETSC_COMM_SELF,"Ball bounced %D times\n",app->nbounces);
56: }
57: app->nbounces++;
58: return(0);
59: }
63: /*
64: Defines the ODE passed to the ODE solver
65: */
66: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
67: {
68: PetscErrorCode ierr;
69: PetscScalar *f;
70: const PetscScalar *u,*udot;
73: /* The next three lines allow us to access the entries of the vectors directly */
74: VecGetArrayRead(U,&u);
75: VecGetArrayRead(Udot,&udot);
76: VecGetArray(F,&f);
78: f[0] = udot[0] - u[1];
79: f[1] = udot[1] + 9.8;
81: VecRestoreArrayRead(U,&u);
82: VecRestoreArrayRead(Udot,&udot);
83: VecRestoreArray(F,&f);
84: return(0);
85: }
89: /*
90: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
91: */
92: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
93: {
94: PetscErrorCode ierr;
95: PetscInt rowcol[] = {0,1};
96: PetscScalar J[2][2];
97: const PetscScalar *u,*udot;
100: VecGetArrayRead(U,&u);
101: VecGetArrayRead(Udot,&udot);
103: J[0][0] = a; J[0][1] = -1;
104: J[1][0] = 0.0; J[1][1] = a;
105: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
107: VecRestoreArrayRead(U,&u);
108: VecRestoreArrayRead(Udot,&udot);
110: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
112: if (A != B) {
113: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
115: }
116: return(0);
117: }
121: int main(int argc,char **argv)
122: {
123: TS ts; /* ODE integrator */
124: Vec U; /* solution will be stored here */
125: Mat A; /* Jacobian matrix */
127: PetscMPIInt size;
128: PetscInt n = 2;
129: PetscScalar *u;
130: AppCtx app;
131: PetscInt direction[2];
132: PetscBool terminate[2];
133: TSAdapt adapt;
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Initialize program
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscInitialize(&argc,&argv,(char*)0,help);
139: MPI_Comm_size(PETSC_COMM_WORLD,&size);
140: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
142: app.nbounces = 0;
143: app.maxbounces = 10;
144: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");
145: PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
146: PetscOptionsEnd();
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Create necessary matrix and vectors
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: MatCreate(PETSC_COMM_WORLD,&A);
152: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
153: MatSetType(A,MATDENSE);
154: MatSetFromOptions(A);
155: MatSetUp(A);
157: MatCreateVecs(A,&U,NULL);
159: VecGetArray(U,&u);
160: u[0] = 0.0;
161: u[1] = 20.0;
162: VecRestoreArray(U,&u);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Create timestepping solver context
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: TSCreate(PETSC_COMM_WORLD,&ts);
168: TSSetSaveTrajectory(ts);
169: TSSetProblemType(ts,TS_NONLINEAR);
170: TSSetType(ts,TSROSW);
171: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
172: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Set initial conditions
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: TSSetSolution(ts,U);
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Set solver options
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: TSSetDuration(ts,1000,30.0);
183: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
184: TSSetInitialTimeStep(ts,0.0,0.1);
186: /* Set directions and terminate flags for the two events */
187: direction[0] = -1; direction[1] = -1;
188: terminate[0] = PETSC_FALSE; terminate[1] = PETSC_TRUE;
189: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
191: /* The adapative time step controller could take very large timesteps resulting in
192: the same event occuring multiple times in the same interval. A maximum step size
193: limit is enforced here to avoid this issue.
194: */
195: TSGetAdapt(ts,&adapt);
196: TSAdaptSetStepLimits(adapt,0.0,0.5);
198: TSSetFromOptions(ts);
200: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: Run timestepping solver
202: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203: TSSolve(ts,U);
205: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206: Free work space. All PETSc objects should be destroyed when they are no longer needed.
207: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208: MatDestroy(&A);
209: VecDestroy(&U);
210: TSDestroy(&ts);
212: PetscFinalize();
213: return(0);
214: }