Actual source code: ex40.c
petsc-3.10.5 2019-03-28
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;
20: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
21: {
22: AppCtx *app=(AppCtx*)ctx;
23: PetscErrorCode ierr;
24: const PetscScalar *u;
27: /* Event for ball height */
28: VecGetArrayRead(U,&u);
29: fvalue[0] = u[0];
30: /* Event for number of bounces */
31: fvalue[1] = app->maxbounces - app->nbounces;
32: VecRestoreArrayRead(U,&u);
33: return(0);
34: }
36: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
37: {
38: AppCtx *app=(AppCtx*)ctx;
40: PetscScalar *u;
43: if (event_list[0] == 0) {
44: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds\n",(double)t);
45: /* Set new initial conditions with .9 attenuation */
46: VecGetArray(U,&u);
47: u[0] = 0.0;
48: u[1] = -0.9*u[1];
49: VecRestoreArray(U,&u);
50: } else if (event_list[0] == 1) {
51: PetscPrintf(PETSC_COMM_SELF,"Ball bounced %D times\n",app->nbounces);
52: }
53: app->nbounces++;
54: return(0);
55: }
57: /*
58: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
59: */
60: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
61: {
62: PetscErrorCode ierr;
63: PetscScalar *f;
64: const PetscScalar *u;
67: /* The next three lines allow us to access the entries of the vectors directly */
68: VecGetArrayRead(U,&u);
69: VecGetArray(F,&f);
71: f[0] = u[1];
72: f[1] = - 9.8;
74: VecRestoreArrayRead(U,&u);
75: VecRestoreArray(F,&f);
76: return(0);
77: }
79: /*
80: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
81: */
82: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
83: {
84: PetscErrorCode ierr;
85: PetscInt rowcol[] = {0,1};
86: PetscScalar J[2][2];
87: const PetscScalar *u;
90: VecGetArrayRead(U,&u);
92: J[0][0] = 0.0; J[0][1] = 1.0;
93: J[1][0] = 0.0; J[1][1] = 0.0;
94: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
96: VecRestoreArrayRead(U,&u);
98: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
100: if (A != B) {
101: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
102: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
103: }
104: return(0);
105: }
107: /*
108: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
109: */
110: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
111: {
112: PetscErrorCode ierr;
113: PetscScalar *f;
114: const PetscScalar *u,*udot;
117: /* The next three lines allow us to access the entries of the vectors directly */
118: VecGetArrayRead(U,&u);
119: VecGetArrayRead(Udot,&udot);
120: VecGetArray(F,&f);
122: f[0] = udot[0] - u[1];
123: f[1] = udot[1] + 9.8;
125: VecRestoreArrayRead(U,&u);
126: VecRestoreArrayRead(Udot,&udot);
127: VecRestoreArray(F,&f);
128: return(0);
129: }
131: /*
132: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
133: */
134: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
135: {
136: PetscErrorCode ierr;
137: PetscInt rowcol[] = {0,1};
138: PetscScalar J[2][2];
139: const PetscScalar *u,*udot;
142: VecGetArrayRead(U,&u);
143: VecGetArrayRead(Udot,&udot);
145: J[0][0] = a; J[0][1] = -1.0;
146: J[1][0] = 0.0; J[1][1] = a;
147: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
149: VecRestoreArrayRead(U,&u);
150: VecRestoreArrayRead(Udot,&udot);
152: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
154: if (A != B) {
155: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
157: }
158: return(0);
159: }
161: int main(int argc,char **argv)
162: {
163: TS ts; /* ODE integrator */
164: Vec U; /* solution will be stored here */
166: PetscMPIInt size;
167: PetscInt n = 2;
168: PetscScalar *u;
169: AppCtx app;
170: PetscInt direction[2];
171: PetscBool terminate[2];
172: PetscBool rhs_form=PETSC_FALSE;
173: TSAdapt adapt;
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Initialize program
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
179: MPI_Comm_size(PETSC_COMM_WORLD,&size);
180: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
182: app.nbounces = 0;
183: app.maxbounces = 10;
184: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex40 options","");
185: PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
186: PetscOptionsEnd();
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Create timestepping solver context
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: TSCreate(PETSC_COMM_WORLD,&ts);
192: TSSetType(ts,TSROSW);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Set ODE routines
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: TSSetProblemType(ts,TS_NONLINEAR);
198: /* Users are advised against the following branching and code duplication.
199: For problems without a mass matrix like the one at hand, the RHSFunction
200: (and companion RHSJacobian) interface is enough to support both explicit
201: and implicit timesteppers. This tutorial example also deals with the
202: IFunction/IJacobian interface for demonstration and testing purposes. */
203: PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
204: if (rhs_form) {
205: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
206: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
207: } else {
208: Mat A; /* Jacobian matrix */
209: MatCreate(PETSC_COMM_WORLD,&A);
210: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
211: MatSetType(A,MATDENSE);
212: MatSetFromOptions(A);
213: MatSetUp(A);
214: TSSetIFunction(ts,NULL,IFunction,NULL);
215: TSSetIJacobian(ts,A,A,IJacobian,NULL);
216: MatDestroy(&A);
217: }
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Set initial conditions
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222: VecCreate(PETSC_COMM_WORLD,&U);
223: VecSetSizes(U,n,PETSC_DETERMINE);
224: VecSetUp(U);
225: VecGetArray(U,&u);
226: u[0] = 0.0;
227: u[1] = 20.0;
228: VecRestoreArray(U,&u);
229: TSSetSolution(ts,U);
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Set solver options
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: TSSetSaveTrajectory(ts);
235: TSSetMaxTime(ts,30.0);
236: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
237: TSSetTimeStep(ts,0.1);
238: /* The adapative time step controller could take very large timesteps resulting in
239: the same event occuring multiple times in the same interval. A maximum step size
240: limit is enforced here to avoid this issue. */
241: TSGetAdapt(ts,&adapt);
242: TSAdaptSetStepLimits(adapt,0.0,0.5);
244: /* Set directions and terminate flags for the two events */
245: direction[0] = -1; direction[1] = -1;
246: terminate[0] = PETSC_FALSE; terminate[1] = PETSC_TRUE;
247: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
249: TSSetFromOptions(ts);
251: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: Run timestepping solver
253: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: TSSolve(ts,U);
256: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257: Free work space. All PETSc objects should be destroyed when they are no longer needed.
258: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
259: VecDestroy(&U);
260: TSDestroy(&ts);
262: PetscFinalize();
263: return ierr;
264: }
266: /*TEST
268: test:
269: suffix: a
270: args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
271: output_file: output/ex40.out
273: test:
274: suffix: b
275: args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
276: output_file: output/ex40.out
278: test:
279: suffix: c
280: args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
281: output_file: output/ex40.out
283: test:
284: suffix: d
285: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
286: output_file: output/ex40.out
288: test:
289: suffix: e
290: args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
291: output_file: output/ex40.out
293: test:
294: suffix: f
295: args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
296: output_file: output/ex40.out
298: test:
299: suffix: g
300: args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
301: output_file: output/ex40.out
303: TEST*/