Actual source code: ex41.c
petsc-3.12.5 2020-03-29
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(A,&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);
93: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
95: if (A != B) {
96: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
98: }
99: return(0);
100: }
102: /*
103: Defines the ODE passed to the ODE solver in implicit form: F(U_t,U) = 0
104: */
105: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
106: {
107: PetscErrorCode ierr;
108: PetscScalar *f;
109: const PetscScalar *u,*udot;
112: /* The next three lines allow us to access the entries of the vectors directly */
113: VecGetArrayRead(U,&u);
114: VecGetArrayRead(Udot,&udot);
115: VecGetArray(F,&f);
117: f[0] = udot[0] - u[1];
118: f[1] = udot[1] + 9.8;
120: VecRestoreArrayRead(U,&u);
121: VecRestoreArrayRead(Udot,&udot);
122: VecRestoreArray(F,&f);
123: return(0);
124: }
126: /*
127: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
128: */
129: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
130: {
131: PetscErrorCode ierr;
132: PetscInt rowcol[2],rstart;
133: PetscScalar J[2][2];
134: const PetscScalar *u,*udot;
137: VecGetArrayRead(U,&u);
138: VecGetArrayRead(Udot,&udot);
140: MatGetOwnershipRange(A,&rstart,NULL);
141: rowcol[0] = rstart; rowcol[1] = rstart+1;
143: J[0][0] = a; J[0][1] = -1.0;
144: J[1][0] = 0.0; J[1][1] = a;
145: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
147: VecRestoreArrayRead(U,&u);
148: VecRestoreArrayRead(Udot,&udot);
150: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
152: if (A != B) {
153: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
155: }
156: return(0);
157: }
159: int main(int argc,char **argv)
160: {
161: TS ts; /* ODE integrator */
162: Vec U; /* solution will be stored here */
164: PetscMPIInt rank;
165: PetscInt n = 2;
166: PetscScalar *u;
167: PetscInt direction=-1;
168: PetscBool terminate=PETSC_FALSE;
169: PetscBool rhs_form=PETSC_FALSE,hist=PETSC_TRUE;
170: TSAdapt adapt;
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Initialize program
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
176: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Create timestepping solver context
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: TSCreate(PETSC_COMM_WORLD,&ts);
182: TSSetType(ts,TSROSW);
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Set ODE routines
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: TSSetProblemType(ts,TS_NONLINEAR);
188: /* Users are advised against the following branching and code duplication.
189: For problems without a mass matrix like the one at hand, the RHSFunction
190: (and companion RHSJacobian) interface is enough to support both explicit
191: and implicit timesteppers. This tutorial example also deals with the
192: IFunction/IJacobian interface for demonstration and testing purposes. */
193: PetscOptionsGetBool(NULL,NULL,"-rhs-form",&rhs_form,NULL);
194: if (rhs_form) {
195: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
196: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,NULL);
197: } else {
198: TSSetIFunction(ts,NULL,IFunction,NULL);
199: TSSetIJacobian(ts,NULL,NULL,IJacobian,NULL);
200: }
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Set initial conditions
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: VecCreate(PETSC_COMM_WORLD,&U);
206: VecSetSizes(U,n,PETSC_DETERMINE);
207: VecSetUp(U);
208: VecGetArray(U,&u);
209: u[0] = 1.0*rank;
210: u[1] = 20.0;
211: VecRestoreArray(U,&u);
212: TSSetSolution(ts,U);
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Set solver options
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: TSSetSaveTrajectory(ts);
218: TSSetMaxTime(ts,30.0);
219: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
220: TSSetTimeStep(ts,0.1);
221: /* The adapative time step controller could take very large timesteps resulting in
222: the same event occuring multiple times in the same interval. A maximum step size
223: limit is enforced here to avoid this issue. */
224: TSGetAdapt(ts,&adapt);
225: TSAdaptSetType(adapt,TSADAPTBASIC);
226: TSAdaptSetStepLimits(adapt,0.0,0.5);
228: /* Set direction and terminate flag for the event */
229: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
231: TSSetFromOptions(ts);
233: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234: Run timestepping solver
235: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236: TSSolve(ts,U);
238: if (hist) { /* replay following history */
239: TSTrajectory tj;
240: PetscReal tf,t0,dt;
242: TSGetTime(ts,&tf);
243: TSSetMaxTime(ts,tf);
244: TSSetStepNumber(ts,0);
245: TSRestartStep(ts);
246: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
247: TSSetFromOptions(ts);
248: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
249: TSGetAdapt(ts,&adapt);
250: TSAdaptSetType(adapt,TSADAPTHISTORY);
251: TSGetTrajectory(ts,&tj);
252: TSAdaptHistorySetTrajectory(adapt,tj,PETSC_FALSE);
253: TSAdaptHistoryGetStep(adapt,0,&t0,&dt);
254: /* this example fails with single (or smaller) precision */
255: #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL__FP16)
256: TSAdaptSetType(adapt,TSADAPTBASIC);
257: TSAdaptSetStepLimits(adapt,0.0,0.5);
258: TSSetFromOptions(ts);
259: #endif
260: TSSetTime(ts,t0);
261: TSSetTimeStep(ts,dt);
262: TSResetTrajectory(ts);
263: VecGetArray(U,&u);
264: u[0] = 1.0*rank;
265: u[1] = 20.0;
266: VecRestoreArray(U,&u);
267: TSSolve(ts,U);
268: }
269: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: Free work space. All PETSc objects should be destroyed when they are no longer needed.
271: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
272: VecDestroy(&U);
273: TSDestroy(&ts);
275: PetscFinalize();
276: return ierr;
277: }
279: /*TEST
281: test:
282: suffix: a
283: nsize: 2
284: args: -ts_trajectory_type memory -snes_stol 1e-4
285: filter: sort -b
287: test:
288: suffix: b
289: nsize: 2
290: args: -ts_trajectory_type memory -ts_type arkimex -snes_stol 1e-4
291: filter: sort -b
293: test:
294: suffix: c
295: nsize: 2
296: args: -ts_trajectory_type memory -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
297: filter: sort -b
299: test:
300: suffix: d
301: nsize: 2
302: args: -ts_trajectory_type memory -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
303: filter: sort -b
305: test:
306: suffix: e
307: nsize: 2
308: args: -ts_trajectory_type memory -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
309: filter: sort -b
311: test:
312: suffix: f
313: nsize: 2
314: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 3bs
315: filter: sort -b
317: test:
318: suffix: g
319: nsize: 2
320: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 5bs
321: filter: sort -b
323: test:
324: suffix: h
325: nsize: 2
326: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 6vr
327: filter: sort -b
328: output_file: output/ex41_g.out
330: test:
331: suffix: i
332: nsize: 2
333: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 7vr
334: filter: sort -b
335: output_file: output/ex41_g.out
337: test:
338: suffix: j
339: nsize: 2
340: args: -ts_trajectory_type memory -rhs-form -ts_type rk -ts_rk_type 8vr
341: filter: sort -b
342: output_file: output/ex41_g.out
344: TEST*/