Actual source code: ex40.c
petsc-3.6.1 2015-08-06
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;
26: PetscErrorCode ierr;
27: const PetscScalar *u;
30: /* Event for ball height */
31: VecGetArrayRead(U,&u);
32: fvalue[0] = u[0];
33: /* Event for number of bounces */
34: fvalue[1] = app->maxbounces - app->nbounces;
35: VecRestoreArrayRead(U,&u);
36: return(0);
37: }
41: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,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",(double)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: {
70: PetscErrorCode ierr;
71: PetscScalar *f;
72: const PetscScalar *u,*udot;
75: /* The next three lines allow us to access the entries of the vectors directly */
76: VecGetArrayRead(U,&u);
77: VecGetArrayRead(Udot,&udot);
78: VecGetArray(F,&f);
80: f[0] = udot[0] - u[1];
81: f[1] = udot[1] + 9.8;
83: VecRestoreArrayRead(U,&u);
84: VecRestoreArrayRead(Udot,&udot);
85: VecRestoreArray(F,&f);
86: return(0);
87: }
91: /*
92: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
93: */
94: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
95: {
96: PetscErrorCode ierr;
97: PetscInt rowcol[] = {0,1};
98: PetscScalar J[2][2];
99: const PetscScalar *u,*udot;
102: VecGetArrayRead(U,&u);
103: VecGetArrayRead(Udot,&udot);
105: J[0][0] = a; J[0][1] = -1;
106: J[1][0] = 0.0; J[1][1] = a;
107: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
109: VecRestoreArrayRead(U,&u);
110: VecRestoreArrayRead(Udot,&udot);
112: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
114: if (A != B) {
115: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
117: }
118: return(0);
119: }
123: int main(int argc,char **argv)
124: {
125: TS ts; /* ODE integrator */
126: Vec U; /* solution will be stored here */
127: Mat A; /* Jacobian matrix */
129: PetscMPIInt size;
130: PetscInt n = 2;
131: PetscScalar *u;
132: AppCtx app;
133: PetscInt direction[2];
134: PetscBool terminate[2];
135: TSAdapt adapt;
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Initialize program
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: PetscInitialize(&argc,&argv,(char*)0,help);
141: MPI_Comm_size(PETSC_COMM_WORLD,&size);
142: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
144: app.nbounces = 0;
145: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21 options","");
146: {
147: app.maxbounces = 10;
148: PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
149: }
150: PetscOptionsEnd();
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Create necessary matrix and vectors
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: MatCreate(PETSC_COMM_WORLD,&A);
156: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
157: MatSetType(A,MATDENSE);
158: MatSetFromOptions(A);
159: MatSetUp(A);
161: MatCreateVecs(A,&U,NULL);
163: VecGetArray(U,&u);
164: u[0] = 0.0;
165: u[1] = 20.0;
166: VecRestoreArray(U,&u);
168: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: Create timestepping solver context
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: TSCreate(PETSC_COMM_WORLD,&ts);
172: TSSetSaveTrajectory(ts);
173: TSSetProblemType(ts,TS_NONLINEAR);
174: TSSetType(ts,TSROSW);
175: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
176: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Set initial conditions
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: TSSetSolution(ts,U);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Set solver options
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: TSSetDuration(ts,1000,30.0);
187: TSSetInitialTimeStep(ts,0.0,0.1);
188: TSSetFromOptions(ts);
189:
190: /* Set directions and terminate flags for the two events */
191: direction[0] = -1; direction[1] = -1;
192: terminate[0] = PETSC_FALSE; terminate[1] = PETSC_TRUE;
193: TSSetEventMonitor(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
195: TSGetAdapt(ts,&adapt);
196: /* The adapative time step controller could take very large timesteps resulting in
197: the same event occuring multiple times in the same interval. A max. step
198: limit is enforced here to avoid this issue.
199: */
200: TSAdaptSetStepLimits(adapt,0.0,0.5);
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Run timestepping solver
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: TSSolve(ts,U);
206: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207: Free work space. All PETSc objects should be destroyed when they are no longer needed.
208: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209:
210: MatDestroy(&A);
211: VecDestroy(&U);
212: TSDestroy(&ts);
213:
214: PetscFinalize();
215: return(0);
216: }