Actual source code: ex40.c
petsc-3.5.4 2015-05-23
2: static char help[] = "Serial bouncing ball example to test TS event feature.\n";
4: /*
5: The dynamics of the bouncing ball is described by the ODE
6: u1_t = u2
7: u2_t = -9.8
9: There are two events set in this example. The first one checks for the ball hitting the
10: ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
11: a factor of 0.9. The second event sets a limit on the number of ball bounces.
12: */
14: #include <petscts.h>
16: typedef struct {
17: PetscInt maxbounces;
18: PetscInt nbounces;
19: } AppCtx;
23: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
24: {
25: AppCtx *app=(AppCtx*)ctx;
27: PetscScalar *u;
30: /* Event for ball height */
31: VecGetArray(U,&u);
32: fvalue[0] = u[0];
33: /* Event for number of bounces */
34: fvalue[1] = app->maxbounces - app->nbounces;
35: VecRestoreArray(U,&u);
36: return(0);
37: }
41: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,void* ctx)
42: {
43: AppCtx *app=(AppCtx*)ctx;
45: PetscScalar *u;
48: if (event_list[0] == 0) {
49: VecGetArray(U,&u);
50: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %f seconds\n",t);
51: /* Set new initial conditions with .9 attenuation */
52: u[0] = 0.0;
53: u[1] = -0.9*u[1];
54: VecRestoreArray(U,&u);
55: TSSetSolution(ts,U);
56: } else if (event_list[0] == 1) {
57: PetscPrintf(PETSC_COMM_SELF,"Ball bounced %d times\n",app->nbounces);
58: }
59: app->nbounces++;
60: return(0);
61: }
65: /*
66: Defines the ODE passed to the ODE solver
67: */
68: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
69: {
71: PetscScalar *u,*udot,*f;
74: /* The next three lines allow us to access the entries of the vectors directly */
75: VecGetArray(U,&u);
76: VecGetArray(Udot,&udot);
77: VecGetArray(F,&f);
79: f[0] = udot[0] - u[1];
80: f[1] = udot[1] + 9.8;
82: VecRestoreArray(U,&u);
83: VecRestoreArray(Udot,&udot);
84: VecRestoreArray(F,&f);
85: return(0);
86: }
90: /*
91: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
92: */
93: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
94: {
96: PetscInt rowcol[] = {0,1};
97: PetscScalar *u,*udot,J[2][2];
100: VecGetArray(U,&u);
101: VecGetArray(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: VecRestoreArray(U,&u);
108: VecRestoreArray(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: }
118:
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];
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Initialize program
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: PetscInitialize(&argc,&argv,(char*)0,help);
138: MPI_Comm_size(PETSC_COMM_WORLD,&size);
139: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
141: app.nbounces = 0;
142: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");
143: {
144: app.maxbounces = 10;
145: PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
146: }
147: PetscOptionsEnd();
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Create necessary matrix and vectors
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: MatCreate(PETSC_COMM_WORLD,&A);
153: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
154: MatSetType(A,MATDENSE);
155: MatSetFromOptions(A);
156: MatSetUp(A);
158: MatGetVecs(A,&U,NULL);
160: VecGetArray(U,&u);
161: u[0] = 0.0;
162: u[1] = 20.0;
163: VecRestoreArray(U,&u);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Create timestepping solver context
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: TSCreate(PETSC_COMM_WORLD,&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: TSSetInitialTimeStep(ts,0.0,0.1);
184: TSSetFromOptions(ts);
185:
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: TSSetEventMonitor(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
191: TSAdapt adapt;
192: TSGetAdapt(ts,&adapt);
193: /* The adapative time step controller could take very large timesteps resulting in
194: the same event occuring multiple times in the same interval. A max. step
195: limit is enforced here to avoid this issue.
196: */
197: TSAdaptSetStepLimits(adapt,0.0,0.5);
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Run timestepping solver
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: TSSolve(ts,U);
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Free work space. All PETSc objects should be destroyed when they are no longer needed.
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206:
207: MatDestroy(&A);
208: VecDestroy(&U);
209: TSDestroy(&ts);
210:
211: PetscFinalize();
212: return(0);
213: }