Actual source code: ex41.c
petsc-3.7.5 2017-01-01
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>
19: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
20: {
21: PetscErrorCode ierr;
22: const PetscScalar *u;
25: /* Event for ball height */
26: VecGetArrayRead(U,&u);
27: fvalue[0] = u[0];
28: VecRestoreArrayRead(U,&u);
29: return(0);
30: }
34: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
35: {
37: PetscScalar *u;
38: PetscMPIInt rank;
41: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
42: if (nevents) {
43: PetscPrintf(PETSC_COMM_SELF,"Processor [%d]: Ball hit the ground at t = %5.2f seconds\n",rank,(double)t);
44: /* Set new initial conditions with .9 attenuation */
45: VecGetArray(U,&u);
46: u[0] = 1.0*rank;
47: u[1] = -0.9*u[1];
48: VecRestoreArray(U,&u);
49: }
50: return(0);
51: }
55: /*
56: Defines the ODE passed to the ODE solver
57: */
58: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
59: {
60: PetscErrorCode ierr;
61: PetscScalar *f;
62: const PetscScalar *u,*udot;
65: /* The next three lines allow us to access the entries of the vectors directly */
66: VecGetArrayRead(U,&u);
67: VecGetArrayRead(Udot,&udot);
68: VecGetArray(F,&f);
70: f[0] = udot[0] - u[1];
71: f[1] = udot[1] + 9.8;
73: VecRestoreArrayRead(U,&u);
74: VecRestoreArrayRead(Udot,&udot);
75: VecRestoreArray(F,&f);
76: return(0);
77: }
81: /*
82: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
83: */
84: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
85: {
86: PetscErrorCode ierr;
87: PetscInt rowcol[2];
88: const PetscScalar *u,*udot;
89: PetscScalar J[2][2];
90: PetscInt rstart;
93: VecGetArrayRead(U,&u);
94: VecGetArrayRead(Udot,&udot);
96: MatGetOwnershipRange(A,&rstart,NULL);
97: rowcol[0] = rstart; rowcol[1] = rstart+1;
99: J[0][0] = a; J[0][1] = -1;
100: J[1][0] = 0.0; J[1][1] = a;
101: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
103: VecRestoreArrayRead(U,&u);
104: VecRestoreArrayRead(Udot,&udot);
106: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108: if (A != B) {
109: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
111: }
112: return(0);
113: }
117: int main(int argc,char **argv)
118: {
119: TS ts; /* ODE integrator */
120: Vec U; /* solution will be stored here */
121: Mat A; /* Jacobian matrix */
123: PetscMPIInt rank;
124: PetscInt n = 2;
125: PetscScalar *u;
126: PetscInt direction=-1;
127: PetscBool terminate=PETSC_FALSE;
128: TSAdapt adapt;
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Initialize program
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: PetscInitialize(&argc,&argv,(char*)0,help);
134: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Create necessary matrix and vectors
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: MatCreate(PETSC_COMM_WORLD,&A);
140: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
141: MatSetType(A,MATDENSE);
142: MatSetFromOptions(A);
143: MatSetUp(A);
145: MatCreateVecs(A,&U,NULL);
147: VecGetArray(U,&u);
148: u[0] = 1.0*rank;
149: u[1] = 20.0;
150: VecRestoreArray(U,&u);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Create timestepping solver context
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: TSCreate(PETSC_COMM_WORLD,&ts);
156: TSSetSaveTrajectory(ts);
157: TSSetProblemType(ts,TS_NONLINEAR);
158: TSSetType(ts,TSROSW);
159: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,NULL);
160: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,NULL);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Set initial conditions
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: TSSetSolution(ts,U);
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Set solver options
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: TSSetDuration(ts,1000,30.0);
171: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
172: TSSetInitialTimeStep(ts,0.0,0.1);
174: TSSetEventHandler(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
176: /* The adapative time step controller could take very large timesteps resulting in
177: the same event occuring multiple times in the same interval. A maximum step size
178: limit is enforced here to avoid this issue.
179: */
180: TSGetAdapt(ts,&adapt);
181: TSAdaptSetStepLimits(adapt,0.0,0.5);
183: TSSetFromOptions(ts);
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Run timestepping solver
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188: TSSolve(ts,U);
190: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: Free work space. All PETSc objects should be destroyed when they are no longer needed.
192: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193: MatDestroy(&A);
194: VecDestroy(&U);
195: TSDestroy(&ts);
197: PetscFinalize();
198: return(0);
199: }