Actual source code: ex41.c
petsc-3.6.1 2015-08-06
2: static char help[] = "Parallel bouncing ball example to test TS event feature.\n";
4: /*
5: The dynamics of the bouncing ball is described by the ODE
6: u1_t = u2
7: u2_t = -9.8
8:
9: Each processor is assigned one ball.
10:
11: The event function routine checks for the ball hitting the
12: ground (u1 = 0). Every time the ball hits the ground, its velocity u2 is attenuated by
13: a factor of 0.9 and its height set to 1.0*rank.
14: */
16: #include <petscts.h>
20: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,void *ctx)
21: {
22: PetscErrorCode ierr;
23: const PetscScalar *u;
26: /* Event for ball height */
27: VecGetArrayRead(U,&u);
28: fvalue[0] = u[0];
29: VecRestoreArrayRead(U,&u);
30: return(0);
31: }
35: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)
36: {
38: PetscScalar *u;
39: PetscMPIInt rank;
42: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
43: if (nevents) {
44: VecGetArray(U,&u);
45: PetscPrintf(PETSC_COMM_SELF,"Processor [%d]: Ball hit the ground at t = %g seconds\n",rank,(double)t);
46: /* Set new initial conditions with .9 attenuation */
47: u[0] = 1.0*rank;
48: u[1] = -0.9*u[1];
49: VecRestoreArray(U,&u);
50: }
51: TSSetSolution(ts,U);
52: return(0);
53: }
57: /*
58: Defines the ODE passed to the ODE solver
59: */
60: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
61: {
62: PetscErrorCode ierr;
63: PetscScalar *f;
64: const PetscScalar *u,*udot;
67: /* The next three lines allow us to access the entries of the vectors directly */
68: VecGetArrayRead(U,&u);
69: VecGetArrayRead(Udot,&udot);
70: VecGetArray(F,&f);
72: f[0] = udot[0] - u[1];
73: f[1] = udot[1] + 9.8;
75: VecRestoreArrayRead(U,&u);
76: VecRestoreArrayRead(Udot,&udot);
77: VecRestoreArray(F,&f);
78: return(0);
79: }
83: /*
84: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
85: */
86: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
87: {
88: PetscErrorCode ierr;
89: PetscInt rowcol[2];
90: const PetscScalar *u,*udot;
91: PetscScalar J[2][2];
92: PetscInt rstart;
95: VecGetArrayRead(U,&u);
96: VecGetArrayRead(Udot,&udot);
98: MatGetOwnershipRange(A,&rstart,NULL);
99: rowcol[0] = rstart; rowcol[1] = rstart+1;
101: J[0][0] = a; J[0][1] = -1;
102: J[1][0] = 0.0; J[1][1] = a;
103: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
105: VecRestoreArrayRead(U,&u);
106: VecRestoreArrayRead(Udot,&udot);
108: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
110: if (A != B) {
111: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
112: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
113: }
114: return(0);
115: }
116:
119: int main(int argc,char **argv)
120: {
121: TS ts; /* ODE integrator */
122: Vec U; /* solution will be stored here */
123: Mat A; /* Jacobian matrix */
125: PetscMPIInt rank;
126: PetscInt n = 2,direction=-1;
127: PetscScalar *u;
128: PetscBool terminate=PETSC_FALSE;
129: TSAdapt adapt;
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Initialize program
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: PetscInitialize(&argc,&argv,(char*)0,help);
135: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create necessary matrix and vectors
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: MatCreate(PETSC_COMM_WORLD,&A);
141: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
142: MatSetType(A,MATDENSE);
143: MatSetFromOptions(A);
144: MatSetUp(A);
146: MatCreateVecs(A,&U,NULL);
148: VecGetArray(U,&u);
149: u[0] = 1.0*rank;
150: u[1] = 20.0;
151: VecRestoreArray(U,&u);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Create timestepping solver context
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: TSCreate(PETSC_COMM_WORLD,&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: TSSetInitialTimeStep(ts,0.0,0.1);
172: TSSetFromOptions(ts);
173:
174: TSSetEventMonitor(ts,1,&direction,&terminate,EventFunction,PostEventFunction,NULL);
176: TSGetAdapt(ts,&adapt);
177: /* The adapative time step controller could take very large timesteps resulting in
178: the same event occuring multiple times in the same interval. A max. step
179: limit is enforced here to avoid this issue.
180: */
181: TSAdaptSetStepLimits(adapt,0.0,0.5);
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Run timestepping solver
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: TSSolve(ts,U);
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Free work space. All PETSc objects should be destroyed when they are no longer needed.
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190:
191: MatDestroy(&A);
192: VecDestroy(&U);
193: TSDestroy(&ts);
194:
195: PetscFinalize();
196: return(0);
197: }