Actual source code: ex40.c
petsc-3.12.5 2020-03-29
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,hist=PETSC_TRUE;
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: PetscOptionsBool("-test_adapthistory","","",hist,&hist,NULL);
187: PetscOptionsEnd();
189: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: Create timestepping solver context
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192: TSCreate(PETSC_COMM_WORLD,&ts);
193: TSSetType(ts,TSROSW);
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Set ODE routines
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: TSSetProblemType(ts,TS_NONLINEAR);
199: /* Users are advised against the following branching and code duplication.
200: For problems without a mass matrix like the one at hand, the RHSFunction
201: (and companion RHSJacobian) interface is enough to support both explicit
202: and implicit timesteppers. This tutorial example also deals with the
203: IFunction/IJacobian interface for demonstration and testing purposes. */
204: PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
205: if (rhs_form) {
206: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
207: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
208: } else {
209: Mat A; /* Jacobian matrix */
210: MatCreate(PETSC_COMM_WORLD,&A);
211: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
212: MatSetType(A,MATDENSE);
213: MatSetFromOptions(A);
214: MatSetUp(A);
215: TSSetIFunction(ts,NULL,IFunction,NULL);
216: TSSetIJacobian(ts,A,A,IJacobian,NULL);
217: MatDestroy(&A);
218: }
220: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: Set initial conditions
222: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223: VecCreate(PETSC_COMM_WORLD,&U);
224: VecSetSizes(U,n,PETSC_DETERMINE);
225: VecSetUp(U);
226: VecGetArray(U,&u);
227: u[0] = 0.0;
228: u[1] = 20.0;
229: VecRestoreArray(U,&u);
230: TSSetSolution(ts,U);
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Set solver options
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: TSSetSaveTrajectory(ts);
236: TSSetMaxTime(ts,30.0);
237: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
238: TSSetTimeStep(ts,0.1);
239: /* The adapative time step controller could take very large timesteps resulting in
240: the same event occuring multiple times in the same interval. A maximum step size
241: limit is enforced here to avoid this issue. */
242: TSGetAdapt(ts,&adapt);
243: TSAdaptSetStepLimits(adapt,0.0,0.5);
245: /* Set directions and terminate flags for the two events */
246: direction[0] = -1; direction[1] = -1;
247: terminate[0] = PETSC_FALSE; terminate[1] = PETSC_TRUE;
248: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&app);
250: TSSetFromOptions(ts);
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Run timestepping solver
254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255: TSSolve(ts,U);
257: if (hist) { /* replay following history */
258: TSTrajectory tj;
259: PetscReal tf,t0,dt;
261: app.nbounces = 0;
262: TSGetTime(ts,&tf);
263: TSSetMaxTime(ts,tf);
264: TSSetStepNumber(ts,0);
265: TSRestartStep(ts);
266: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
267: TSSetFromOptions(ts);
268: TSGetAdapt(ts,&adapt);
269: TSAdaptSetType(adapt,TSADAPTHISTORY);
270: TSGetTrajectory(ts,&tj);
271: TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
272: TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
273: /* this example fails with single (or smaller) precision */
274: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
275: TSAdaptSetType(adapt,TSADAPTBASIC);
276: TSAdaptSetStepLimits(adapt,0.0,0.5);
277: TSSetFromOptions(ts);
278: #endif
279: TSSetTime(ts,t0);
280: TSSetTimeStep(ts,dt);
281: TSResetTrajectory(ts);
282: VecGetArray(U,&u);
283: u[0] = 0.0;
284: u[1] = 20.0;
285: VecRestoreArray(U,&u);
286: TSSolve(ts,U);
287: }
288: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: Free work space. All PETSc objects should be destroyed when they are no longer needed.
290: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291: VecDestroy(&U);
292: TSDestroy(&ts);
294: PetscFinalize();
295: return ierr;
296: }
298: /*TEST
300: test:
301: suffix: a
302: args: -snes_stol 1e-4 -ts_trajectory_dirname ex40_a_dir
303: output_file: output/ex40.out
305: test:
306: suffix: b
307: args: -ts_type arkimex -snes_stol 1e-4 -ts_trajectory_dirname ex40_b_dir
308: output_file: output/ex40.out
310: test:
311: suffix: c
312: args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_c_dir
313: output_file: output/ex40.out
315: test:
316: suffix: d
317: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4 -ts_trajectory_dirname ex40_d_dir
318: output_file: output/ex40.out
320: test:
321: suffix: e
322: args: -ts_type bdf -ts_adapt_dt_max 0.025 -ts_max_steps 1500 -ts_trajectory_dirname ex40_e_dir
323: output_file: output/ex40.out
325: test:
326: suffix: f
327: args: -rhs-form -ts_type rk -ts_rk_type 3bs -ts_trajectory_dirname ex40_f_dir
328: output_file: output/ex40.out
330: test:
331: suffix: g
332: args: -rhs-form -ts_type rk -ts_rk_type 5bs -ts_trajectory_dirname ex40_g_dir
333: output_file: output/ex40.out
335: test:
336: suffix: h
337: args: -rhs-form -ts_type rk -ts_rk_type 6vr -ts_trajectory_dirname ex40_h_dir
338: output_file: output/ex40.out
340: test:
341: suffix: i
342: args: -rhs-form -ts_type rk -ts_rk_type 7vr -ts_trajectory_dirname ex40_i_dir
343: output_file: output/ex40.out
345: test:
346: suffix: j
347: args: -rhs-form -ts_type rk -ts_rk_type 8vr -ts_trajectory_dirname ex40_j_dir
348: output_file: output/ex40.out
350: TEST*/