Actual source code: ex41.c
1: static char help[] = "Parallel 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: Each processor is assigned one ball.
10: The event function routine checks for the ball hitting the
11: ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
12: a factor of 0.9 and its height set to 1.0*rank.
13: */
15: #include <petscts.h>
17: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
18: {
19: const PetscScalar *u;
21: /* Event for ball height */
22: VecGetArrayRead(U,&u);
23: fvalue[0] = u[0];
24: VecRestoreArrayRead(U,&u);
25: return 0;
26: }
28: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
29: {
30: PetscScalar *u;
31: PetscMPIInt rank;
33: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
34: if (nevents) {
35: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n",(double)t,rank);
36: /* Set new initial conditions with .9 attenuation */
37: VecGetArray(U,&u);
38: u[0] = 1.0*rank;
39: u[1] = -0.9*u[1];
40: VecRestoreArray(U,&u);
41: }
42: return 0;
43: }
45: /*
46: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
47: */
48: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
49: {
50: PetscScalar *f;
51: const PetscScalar *u;
53: /* The next three lines allow us to access the entries of the vectors directly */
54: VecGetArrayRead(U,&u);
55: VecGetArray(F,&f);
57: f[0] = u[1];
58: f[1] = - 9.8;
60: VecRestoreArrayRead(U,&u);
61: VecRestoreArray(F,&f);
62: return 0;
63: }
65: /*
66: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
67: */
68: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
69: {
70: PetscInt rowcol[2],rstart;
71: PetscScalar J[2][2];
72: const PetscScalar *u;
74: VecGetArrayRead(U,&u);
76: MatGetOwnershipRange(B,&rstart,NULL);
77: rowcol[0] = rstart; rowcol[1] = rstart+1;
79: J[0][0] = 0.0; J[0][1] = 1.0;
80: J[1][0] = 0.0; J[1][1] = 0.0;
81: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
83: VecRestoreArrayRead(U,&u);
84: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
86: if (A != B) {
87: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
89: }
90: return 0;
91: }
93: /*
94: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
95: */
96: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
97: {
98: PetscScalar *f;
99: const PetscScalar *u,*udot;
101: /* The next three lines allow us to access the entries of the vectors directly */
102: VecGetArrayRead(U,&u);
103: VecGetArrayRead(Udot,&udot);
104: VecGetArray(F,&f);
106: f[0] = udot[0] - u[1];
107: f[1] = udot[1] + 9.8;
109: VecRestoreArrayRead(U,&u);
110: VecRestoreArrayRead(Udot,&udot);
111: VecRestoreArray(F,&f);
112: return 0;
113: }
115: /*
116: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
117: */
118: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
119: {
120: PetscInt rowcol[2],rstart;
121: PetscScalar J[2][2];
122: const PetscScalar *u,*udot;
124: VecGetArrayRead(U,&u);
125: VecGetArrayRead(Udot,&udot);
127: MatGetOwnershipRange(B,&rstart,NULL);
128: rowcol[0] = rstart; rowcol[1] = rstart+1;
130: J[0][0] = a; J[0][1] = -1.0;
131: J[1][0] = 0.0; J[1][1] = a;
132: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
134: VecRestoreArrayRead(U,&u);
135: VecRestoreArrayRead(Udot,&udot);
137: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
139: if (A != B) {
140: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142: }
143: return 0;
144: }
146: int main(int argc,char **argv)
147: {
148: TS ts; /* ODE integrator */
149: Vec U; /* solution will be stored here */
150: PetscMPIInt rank;
151: PetscInt n = 2;
152: PetscScalar *u;
153: PetscInt direction=-1;
154: PetscBool terminate=PETSC_FALSE;
155: PetscBool rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
156: TSAdapt adapt;
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Initialize program
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: PetscInitialize(&argc,&argv,(char*)0,help);
162: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Create timestepping solver context
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: TSCreate(PETSC_COMM_WORLD,&ts);
168: TSSetType(ts,TSROSW);
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Set ODE routines
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: TSSetProblemType(ts,TS_NONLINEAR);
174: /* Users are advised against the following branching and code duplication.
175: For problems without a mass matrix like the one at hand, the RHSFunction
176: (and companion RHSJacobian) interface is enough to support both explicit
177: and implicit timesteppers. This tutorial example also deals with the
178: IFunction/IJacobian interface for demonstration and testing purposes. */
179: PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
180: if (rhs_form) {
181: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
182: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
183: } else {
184: TSSetIFunction(ts,NULL,IFunction,NULL);
185: TSSetIJacobian(ts,NULL,NULL,IJacobian,NULL);
186: }
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Set initial conditions
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: VecCreate(PETSC_COMM_WORLD,&U);
192: VecSetSizes(U,n,PETSC_DETERMINE);
193: VecSetUp(U);
194: VecGetArray(U,&u);
195: u[0] = 1.0*rank;
196: u[1] = 20.0;
197: VecRestoreArray(U,&u);
198: TSSetSolution(ts,U);
200: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: Set solver options
202: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203: TSSetSaveTrajectory(ts);
204: TSSetMaxTime(ts,30.0);
205: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
206: TSSetTimeStep(ts,0.1);
207: /* The adaptive time step controller could take very large timesteps resulting in
208: the same event occurring multiple times in the same interval. A maximum step size
209: limit is enforced here to avoid this issue. */
210: TSGetAdapt(ts,&adapt);
211: TSAdaptSetType(adapt,TSADAPTBASIC);
212: TSAdaptSetStepLimits(adapt,0.0,0.5);
214: /* Set direction and terminate flag for the event */
215: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
217: TSSetFromOptions(ts);
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Run timestepping solver
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222: TSSolve(ts,U);
224: if (hist) { /* replay following history */
225: TSTrajectory tj;
226: PetscReal tf,t0,dt;
228: TSGetTime(ts,&tf);
229: TSSetMaxTime(ts,tf);
230: TSSetStepNumber(ts,0);
231: TSRestartStep(ts);
232: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
233: TSSetFromOptions(ts);
234: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
235: TSGetAdapt(ts,&adapt);
236: TSAdaptSetType(adapt,TSADAPTHISTORY);
237: TSGetTrajectory(ts,&tj);
238: TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
239: TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
240: /* this example fails with single (or smaller) precision */
241: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
242: TSAdaptSetType(adapt,TSADAPTBASIC);
243: TSAdaptSetStepLimits(adapt,0.0,0.5);
244: TSSetFromOptions(ts);
245: #endif
246: TSSetTime(ts,t0);
247: TSSetTimeStep(ts,dt);
248: TSResetTrajectory(ts);
249: VecGetArray(U,&u);
250: u[0] = 1.0*rank;
251: u[1] = 20.0;
252: VecRestoreArray(U,&u);
253: TSSolve(ts,U);
254: }
255: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: Free work space. All PETSc objects should be destroyed when they are no longer needed.
257: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
258: VecDestroy(&U);
259: TSDestroy(&ts);
261: PetscFinalize();
262: return 0;
263: }
265: /*TEST
267: test:
268: suffix: a
269: nsize: 2
270: args: -ts_trajectory_type memory -snes_stol 1e-4
271: filter: sort -b
273: test:
274: suffix: b
275: nsize: 2
276: args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
277: filter: sort -b
279: test:
280: suffix: c
281: nsize: 2
282: args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
283: filter: sort -b
285: test:
286: suffix: d
287: nsize: 2
288: args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
289: filter: sort -b
291: test:
292: suffix: e
293: nsize: 2
294: args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
295: filter: sort -b
297: test:
298: suffix: f
299: nsize: 2
300: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
301: filter: sort -b
303: test:
304: suffix: g
305: nsize: 2
306: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
307: filter: sort -b
309: test:
310: suffix: h
311: nsize: 2
312: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
313: filter: sort -b
314: output_file: output/ex41_g.out
316: test:
317: suffix: i
318: nsize: 2
319: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
320: filter: sort -b
321: output_file: output/ex41_g.out
323: test:
324: suffix: j
325: nsize: 2
326: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
327: filter: sort -b
328: output_file: output/ex41_g.out
330: TEST*/