Actual source code: ex40.c
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: const PetscScalar *u;
25: /* Event for ball height */
26: VecGetArrayRead(U,&u);
27: fvalue[0] = u[0];
28: /* Event for number of bounces */
29: fvalue[1] = app->maxbounces - app->nbounces;
30: VecRestoreArrayRead(U,&u);
31: return 0;
32: }
34: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
35: {
36: AppCtx *app=(AppCtx*)ctx;
37: PetscScalar *u;
39: if (event_list[0] == 0) {
40: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds\n",(double)t);
41: /* Set new initial conditions with .9 attenuation */
42: VecGetArray(U,&u);
43: u[0] = 0.0;
44: u[1] = -0.9*u[1];
45: VecRestoreArray(U,&u);
46: } else if (event_list[0] == 1) {
47: PetscPrintf(PETSC_COMM_SELF,"Ball bounced %D times\n",app->nbounces);
48: }
49: app->nbounces++;
50: return 0;
51: }
53: /*
54: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
55: */
56: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
57: {
58: PetscScalar *f;
59: const PetscScalar *u;
61: /* The next three lines allow us to access the entries of the vectors directly */
62: VecGetArrayRead(U,&u);
63: VecGetArray(F,&f);
65: f[0] = u[1];
66: f[1] = - 9.8;
68: VecRestoreArrayRead(U,&u);
69: VecRestoreArray(F,&f);
70: return 0;
71: }
73: /*
74: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of the Jacobian.
75: */
76: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
77: {
78: PetscInt rowcol[] = {0,1};
79: PetscScalar J[2][2];
80: const PetscScalar *u;
82: VecGetArrayRead(U,&u);
84: J[0][0] = 0.0; J[0][1] = 1.0;
85: J[1][0] = 0.0; J[1][1] = 0.0;
86: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
88: VecRestoreArrayRead(U,&u);
90: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
92: if (A != B) {
93: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
95: }
96: return 0;
97: }
99: /*
100: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
101: */
102: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
103: {
104: PetscScalar *f;
105: const PetscScalar *u,*udot;
107: /* The next three lines allow us to access the entries of the vectors directly */
108: VecGetArrayRead(U,&u);
109: VecGetArrayRead(Udot,&udot);
110: VecGetArray(F,&f);
112: f[0] = udot[0] - u[1];
113: f[1] = udot[1] + 9.8;
115: VecRestoreArrayRead(U,&u);
116: VecRestoreArrayRead(Udot,&udot);
117: VecRestoreArray(F,&f);
118: return 0;
119: }
121: /*
122: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
123: */
124: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
125: {
126: PetscInt rowcol[] = {0,1};
127: PetscScalar J[2][2];
128: const PetscScalar *u,*udot;
130: VecGetArrayRead(U,&u);
131: VecGetArrayRead(Udot,&udot);
133: J[0][0] = a; J[0][1] = -1.0;
134: J[1][0] = 0.0; J[1][1] = a;
135: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
137: VecRestoreArrayRead(U,&u);
138: VecRestoreArrayRead(Udot,&udot);
140: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142: if (A != B) {
143: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
144: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
145: }
146: return 0;
147: }
149: int main(int argc,char **argv)
150: {
151: TS ts; /* ODE integrator */
152: Vec U; /* solution will be stored here */
154: PetscMPIInt size;
155: PetscInt n = 2;
156: PetscScalar *u;
157: AppCtx app;
158: PetscInt direction[2];
159: PetscBool terminate[2];
160: PetscBool rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
161: TSAdapt adapt;
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Initialize program
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: PetscInitialize(&argc,&argv,(char*)0,help);
167: MPI_Comm_size(PETSC_COMM_WORLD,&size);
170: app.nbounces = 0;
171: app.maxbounces = 10;
172: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex40 options","");
173: PetscOptionsInt("-maxbounces","","",app.maxbounces,&app.maxbounces,NULL);
174: PetscOptionsBool("-test_adapthistory","","",hist,&hist,NULL);
175: PetscOptionsEnd();
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Create timestepping solver context
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: TSCreate(PETSC_COMM_WORLD,&ts);
181: TSSetType(ts,TSROSW);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Set ODE routines
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: TSSetProblemType(ts,TS_NONLINEAR);
187: /* Users are advised against the following branching and code duplication.
188: For problems without a mass matrix like the one at hand, the RHSFunction
189: (and companion RHSJacobian) interface is enough to support both explicit
190: and implicit timesteppers. This tutorial example also deals with the
191: IFunction/IJacobian interface for demonstration and testing purposes. */
192: PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
193: if (rhs_form) {
194: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
195: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
196: } else {
197: Mat A; /* Jacobian matrix */
198: MatCreate(PETSC_COMM_WORLD,&A);
199: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
200: MatSetType(A,MATDENSE);
201: MatSetFromOptions(A);
202: MatSetUp(A);
203: TSSetIFunction(ts,NULL,IFunction,NULL);
204: TSSetIJacobian(ts,A,A,IJacobian,NULL);
205: MatDestroy(&A);
206: }
208: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: Set initial conditions
210: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
211: VecCreate(PETSC_COMM_WORLD,&U);
212: VecSetSizes(U,n,PETSC_DETERMINE);
213: VecSetUp(U);
214: VecGetArray(U,&u);
215: u[0] = 0.0;
216: u[1] = 20.0;
217: VecRestoreArray(U,&u);
218: TSSetSolution(ts,U);
220: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: Set solver options
222: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223: if (hist) TSSetSaveTrajectory(ts);
224: TSSetMaxTime(ts,30.0);
225: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
226: TSSetTimeStep(ts,0.1);
227: /* The adaptive time step controller could take very large timesteps resulting in
228: the same event occurring multiple times in the same interval. A maximum step size
229: limit is enforced here to avoid this issue. */
230: TSGetAdapt(ts,&adapt);
231: TSAdaptSetStepLimits(adapt,0.0,0.5);
233: /* Set directions and terminate flags for the two events */
234: direction[0] = -1; direction[1] = -1;
235: terminate[0] = PETSC_FALSE; terminate[1] = PETSC_TRUE;
236: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
238: TSSetFromOptions(ts);
240: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: Run timestepping solver
242: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243: TSSolve(ts,U);
245: if (hist) { /* replay following history */
246: TSTrajectory tj;
247: PetscReal tf,t0,dt;
249: app.nbounces = 0;
250: TSGetTime(ts,&tf);
251: TSSetMaxTime(ts,tf);
252: TSSetStepNumber(ts,0);
253: TSRestartStep(ts);
254: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
255: TSSetFromOptions(ts);
256: TSGetAdapt(ts,&adapt);
257: TSAdaptSetType(adapt,TSADAPTHISTORY);
258: TSGetTrajectory(ts,&tj);
259: TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
260: TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
261: /* this example fails with single (or smaller) precision */
262: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
263: TSAdaptSetType(adapt,TSADAPTBASIC);
264: TSAdaptSetStepLimits(adapt,0.0,0.5);
265: TSSetFromOptions(ts);
266: #endif
267: TSSetTime(ts,t0);
268: TSSetTimeStep(ts,dt);
269: TSResetTrajectory(ts);
270: VecGetArray(U,&u);
271: u[0] = 0.0;
272: u[1] = 20.0;
273: VecRestoreArray(U,&u);
274: TSSolve(ts,U);
275: }
276: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: Free work space. All PETSc objects should be destroyed when they are no longer needed.
278: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279: VecDestroy(&U);
280: TSDestroy(&ts);
282: PetscFinalize();
283: return 0;
284: }
286: /*TEST
288: test:
289: suffix: a
290: args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
291: output_file: output/ex40.out
293: test:
294: suffix: b
295: args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
296: output_file: output/ex40.out
298: test:
299: suffix: c
300: args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
301: output_file: output/ex40.out
303: test:
304: suffix: d
305: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
306: output_file: output/ex40.out
308: test:
309: suffix: e
310: args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
311: output_file: output/ex40.out
313: test:
314: suffix: f
315: args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
316: output_file: output/ex40.out
318: test:
319: suffix: g
320: args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
321: output_file: output/ex40.out
323: test:
324: suffix: h
325: args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
326: output_file: output/ex40.out
328: test:
329: suffix: i
330: args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
331: output_file: output/ex40.out
333: test:
334: suffix: j
335: args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
336: output_file: output/ex40.out
338: test:
339: suffix: k
340: args: -ts_type theta -ts_adapt_type dsp -ts_trajectory_dirname ex40_k_dir
341: output_file: output/ex40.out
343: test:
344: suffix: l
345: args: -rhs-form -ts_type rk -ts_rk_type 2a -ts_trajectory_dirname ex40_l_dir
346: args: -ts_adapt_type dsp -ts_adapt_always_accept {{false true}} -ts_adapt_dt_min 0.01
347: output_file: output/ex40.out
349: test:
350: suffix: m
351: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -test_adapthistory false
352: args: -ts_max_time 10 -ts_exact_final_time {{STEPOVER MATCHSTEP INTERPOLATE}}
354: test:
355: requires: !single
356: suffix: n
357: args: -test_adapthistory false
358: args: -ts_type alpha -ts_alpha_radius 1.0 -ts_view
359: args: -ts_dt 0.25 -ts_adapt_type basic -ts_adapt_wnormtype INFINITY -ts_adapt_monitor
360: args: -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false
362: test:
363: requires: !single
364: suffix: o
365: args: -rhs-form -ts_type rk -ts_rk_type 2b -ts_trajectory_dirname ex40_o_dir
366: output_file: output/ex40.out
367: TEST*/