Actual source code: ex41.c
petsc-3.9.4 2018-09-11
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,"Processor [%d]: Ball hit the ground at t = %5.2f seconds\n",rank,(double)t);
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;
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: 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: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: Free work space. All PETSc objects should be destroyed when they are no longer needed.
239: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240: VecDestroy(&U);
241: TSDestroy(&ts);
243: PetscFinalize();
244: return ierr;
245: }
247: /*TEST
249: test:
250: suffix: a
251: nsize: 2
252: args: -snes_stol 1e-4
253: filter: sort -b
254: filter_output: sort -b
256: test:
257: suffix: b
258: nsize: 2
259: args: -ts_type arkimex -snes_stol 1e-4
260: filter: sort -b
261: filter_output: sort -b
263: test:
264: suffix: c
265: nsize: 2
266: args: -ts_type theta -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
267: filter: sort -b
268: filter_output: sort -b
270: test:
271: suffix: d
272: nsize: 2
273: args: -ts_type alpha -ts_adapt_type basic -ts_atol 1e-1 -snes_stol 1e-4
274: filter: sort -b
275: filter_output: sort -b
277: test:
278: suffix: e
279: nsize: 2
280: args: -ts_type bdf -ts_adapt_dt_max 0.015 -ts_max_steps 3000
281: filter: sort -b
282: filter_output: sort -b
284: test:
285: suffix: f
286: nsize: 2
287: args: -rhs-form -ts_type rk -ts_rk_type 3bs
288: filter: sort -b
289: filter_output: sort -b
291: test:
292: suffix: g
293: nsize: 2
294: args: -rhs-form -ts_type rk -ts_rk_type 5bs
295: filter: sort -b
296: filter_output: sort -b
298: TEST*/