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: PetscErrorCode ierr;
20: const PetscScalar *u;
23: /* Event for ball height */
24: VecGetArrayRead(U,&u);
25: fvalue[0] = u[0];
26: VecRestoreArrayRead(U,&u);
27: return(0);
28: }
30: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
31: {
33: PetscScalar *u;
34: PetscMPIInt rank;
37: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
38: if (nevents) {
39: PetscPrintf(PETSC_COMM_SELF,"Ball hit the ground at t = %5.2f seconds -> Processor[%d]\n",(double)t,rank);
40: /* Set new initial conditions with .9 attenuation */
41: VecGetArray(U,&u);
42: u[0] = 1.0*rank;
43: u[1] = -0.9*u[1];
44: VecRestoreArray(U,&u);
45: }
46: return(0);
47: }
49: /*
50: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
51: */
52: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
53: {
54: PetscErrorCode ierr;
55: PetscScalar *f;
56: const PetscScalar *u;
59: /* The next three lines allow us to access the entries of the vectors directly */
60: VecGetArrayRead(U,&u);
61: VecGetArray(F,&f);
63: f[0] = u[1];
64: f[1] = - 9.8;
66: VecRestoreArrayRead(U,&u);
67: VecRestoreArray(F,&f);
68: return(0);
69: }
71: /*
72: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning the Jacobian.
73: */
74: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
75: {
76: PetscErrorCode ierr;
77: PetscInt rowcol[2],rstart;
78: PetscScalar J[2][2];
79: const PetscScalar *u;
82: VecGetArrayRead(U,&u);
84: MatGetOwnershipRange(B,&rstart,NULL);
85: rowcol[0] = rstart; rowcol[1] = rstart+1;
87: J[0][0] = 0.0; J[0][1] = 1.0;
88: J[1][0] = 0.0; J[1][1] = 0.0;
89: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
91: VecRestoreArrayRead(U,&u);
92: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
94: if (A != B) {
95: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
97: }
98: return(0);
99: }
101: /*
102: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
103: */
104: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
105: {
106: PetscErrorCode ierr;
107: PetscScalar *f;
108: const PetscScalar *u,*udot;
111: /* The next three lines allow us to access the entries of the vectors directly */
112: VecGetArrayRead(U,&u);
113: VecGetArrayRead(Udot,&udot);
114: VecGetArray(F,&f);
116: f[0] = udot[0] - u[1];
117: f[1] = udot[1] + 9.8;
119: VecRestoreArrayRead(U,&u);
120: VecRestoreArrayRead(Udot,&udot);
121: VecRestoreArray(F,&f);
122: return(0);
123: }
125: /*
126: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
127: */
128: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
129: {
130: PetscErrorCode ierr;
131: PetscInt rowcol[2],rstart;
132: PetscScalar J[2][2];
133: const PetscScalar *u,*udot;
136: VecGetArrayRead(U,&u);
137: VecGetArrayRead(Udot,&udot);
139: MatGetOwnershipRange(B,&rstart,NULL);
140: rowcol[0] = rstart; rowcol[1] = rstart+1;
142: J[0][0] = a; J[0][1] = -1.0;
143: J[1][0] = 0.0; J[1][1] = a;
144: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
146: VecRestoreArrayRead(U,&u);
147: VecRestoreArrayRead(Udot,&udot);
149: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
150: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
151: if (A != B) {
152: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
154: }
155: return(0);
156: }
158: int main(int argc,char **argv)
159: {
160: TS ts; /* ODE integrator */
161: Vec U; /* solution will be stored here */
163: PetscMPIInt rank;
164: PetscInt n = 2;
165: PetscScalar *u;
166: PetscInt direction=-1;
167: PetscBool terminate=PETSC_FALSE;
168: PetscBool rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
169: TSAdapt adapt;
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Initialize program
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
175: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
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: TSSetIFunction(ts,NULL,IFunction,NULL);
198: TSSetIJacobian(ts,NULL,NULL,IJacobian,NULL);
199: }
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Set initial conditions
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: VecCreate(PETSC_COMM_WORLD,&U);
205: VecSetSizes(U,n,PETSC_DETERMINE);
206: VecSetUp(U);
207: VecGetArray(U,&u);
208: u[0] = 1.0*rank;
209: u[1] = 20.0;
210: VecRestoreArray(U,&u);
211: TSSetSolution(ts,U);
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Set solver options
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: TSSetSaveTrajectory(ts);
217: TSSetMaxTime(ts,30.0);
218: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
219: TSSetTimeStep(ts,0.1);
220: /* The adapative time step controller could take very large timesteps resulting in
221: the same event occuring multiple times in the same interval. A maximum step size
222: limit is enforced here to avoid this issue. */
223: TSGetAdapt(ts,&adapt);
224: TSAdaptSetType(adapt,TSADAPTBASIC);
225: TSAdaptSetStepLimits(adapt,0.0,0.5);
227: /* Set direction and terminate flag for the event */
228: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
230: TSSetFromOptions(ts);
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Run timestepping solver
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: TSSolve(ts,U);
237: if (hist) { /* replay following history */
238: TSTrajectory tj;
239: PetscReal tf,t0,dt;
241: TSGetTime(ts,&tf);
242: TSSetMaxTime(ts,tf);
243: TSSetStepNumber(ts,0);
244: TSRestartStep(ts);
245: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
246: TSSetFromOptions(ts);
247: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
248: TSGetAdapt(ts,&adapt);
249: TSAdaptSetType(adapt,TSADAPTHISTORY);
250: TSGetTrajectory(ts,&tj);
251: TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
252: TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
253: /* this example fails with single (or smaller) precision */
254: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
255: TSAdaptSetType(adapt,TSADAPTBASIC);
256: TSAdaptSetStepLimits(adapt,0.0,0.5);
257: TSSetFromOptions(ts);
258: #endif
259: TSSetTime(ts,t0);
260: TSSetTimeStep(ts,dt);
261: TSResetTrajectory(ts);
262: VecGetArray(U,&u);
263: u[0] = 1.0*rank;
264: u[1] = 20.0;
265: VecRestoreArray(U,&u);
266: TSSolve(ts,U);
267: }
268: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: Free work space. All PETSc objects should be destroyed when they are no longer needed.
270: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: VecDestroy(&U);
272: TSDestroy(&ts);
274: PetscFinalize();
275: return ierr;
276: }
278: /*TEST
280: test:
281: suffix: a
282: nsize: 2
283: args: -ts_trajectory_type memory -snes_stol 1e-4
284: filter: sort -b
286: test:
287: suffix: b
288: nsize: 2
289: args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
290: filter: sort -b
292: test:
293: suffix: c
294: nsize: 2
295: args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
296: filter: sort -b
298: test:
299: suffix: d
300: nsize: 2
301: args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
302: filter: sort -b
304: test:
305: suffix: e
306: nsize: 2
307: args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
308: filter: sort -b
310: test:
311: suffix: f
312: nsize: 2
313: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
314: filter: sort -b
316: test:
317: suffix: g
318: nsize: 2
319: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
320: filter: sort -b
322: test:
323: suffix: h
324: nsize: 2
325: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
326: filter: sort -b
327: output_file: output/ex41_g.out
329: test:
330: suffix: i
331: nsize: 2
332: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
333: filter: sort -b
334: output_file: output/ex41_g.out
336: test:
337: suffix: j
338: nsize: 2
339: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
340: filter: sort -b
341: output_file: output/ex41_g.out
343: TEST*/