1 static char help[] = "Solves the motion of spring.\n\
2 Input parameters include:\n";
3
4 /* ------------------------------------------------------------------------
5
6 This program solves the motion of spring by Hooke's law
7 x' = f(t,v) = v
8 v' = g(t,x) = -omega^2*x
9 on the interval 0 <= t <= 0.1, with the initial conditions
10 x(0) = 0.2, x'(0) = v(0) = 0,
11 and
12 omega = 64.
13 The exact solution is
14 x(t) = A*sin(t*omega) + B*cos(t*omega)
15 where A and B are constants that can be determined from the initial conditions.
16 In this case, B=0.2, A=0.
17
18 Notes:
19 This code demonstrates the TS solver interface to solve a separable Hamiltonian
20 system, which can be split into two subsystems involving two coupling components,
21 named generailized momentum and generailized position respectively.
22 Using a symplectic intergrator can preserve energy
23 E = (v^2+omega^2*x^2-omega^2*h*v*x)/2
24 ------------------------------------------------------------------------- */
25
26 #include <petscts.h>
27 #include <petscvec.h>
28
29 typedef struct _n_User *User;
30 struct _n_User {
31 PetscReal omega;
32 PetscInt nts; /* print the energy at each nts time steps */
33 };
34
35 /*
36 User-defined routines.
37 The first RHS function provides f(t,x), the residual for the generalized momentum,
38 and the second one provides g(t,v), the residual for the generalized position.
39 */
RHSFunction2(TS ts,PetscReal t,Vec X,Vec Vres,PetscCtx ctx)40 static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, PetscCtx ctx)
41 {
42 User user = (User)ctx;
43 const PetscScalar *x;
44 PetscScalar *vres;
45
46 PetscFunctionBeginUser;
47 PetscCall(VecGetArrayRead(X, &x));
48 PetscCall(VecGetArray(Vres, &vres));
49 vres[0] = -user->omega * user->omega * x[0];
50 PetscCall(VecRestoreArray(Vres, &vres));
51 PetscCall(VecRestoreArrayRead(X, &x));
52 PetscFunctionReturn(PETSC_SUCCESS);
53 }
54
RHSFunction1(TS ts,PetscReal t,Vec V,Vec Xres,PetscCtx ctx)55 static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, PetscCtx ctx)
56 {
57 const PetscScalar *v;
58 PetscScalar *xres;
59
60 PetscFunctionBeginUser;
61 PetscCall(VecGetArray(Xres, &xres));
62 PetscCall(VecGetArrayRead(V, &v));
63 xres[0] = v[0];
64 PetscCall(VecRestoreArrayRead(V, &v));
65 PetscCall(VecRestoreArray(Xres, &xres));
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68
RHSFunction(TS ts,PetscReal t,Vec U,Vec R,PetscCtx ctx)69 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec R, PetscCtx ctx)
70 {
71 User user = (User)ctx;
72 const PetscScalar *u;
73 PetscScalar *r;
74
75 PetscFunctionBeginUser;
76 PetscCall(VecGetArrayRead(U, &u));
77 PetscCall(VecGetArray(R, &r));
78 r[0] = u[1];
79 r[1] = -user->omega * user->omega * u[0];
80 PetscCall(VecRestoreArrayRead(U, &u));
81 PetscCall(VecRestoreArray(R, &r));
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
85 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
Monitor(TS ts,PetscInt step,PetscReal t,Vec U,PetscCtx ctx)86 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
87 {
88 const PetscScalar *u;
89 PetscReal dt;
90 PetscScalar energy, menergy;
91 User user = (User)ctx;
92
93 PetscFunctionBeginUser;
94 if (step % user->nts == 0) {
95 PetscCall(TSGetTimeStep(ts, &dt));
96 PetscCall(VecGetArrayRead(U, &u));
97 menergy = (u[1] * u[1] + user->omega * user->omega * u[0] * u[0] - user->omega * user->omega * dt * u[0] * u[1]) / 2.;
98 energy = (u[1] * u[1] + user->omega * user->omega * u[0] * u[0]) / 2.;
99 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At time %.6lf, Energy = %8g, Modified Energy = %8g\n", (double)t, (double)energy, (double)menergy));
100 PetscCall(VecRestoreArrayRead(U, &u));
101 }
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
main(int argc,char ** argv)105 int main(int argc, char **argv)
106 {
107 TS ts; /* nonlinear solver */
108 Vec U; /* solution, residual vectors */
109 IS is1, is2;
110 PetscInt nindices[1];
111 PetscReal ftime = 0.1;
112 PetscBool monitor = PETSC_FALSE;
113 PetscScalar *u_ptr;
114 PetscMPIInt size;
115 struct _n_User user;
116
117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 Initialize program
119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120 PetscFunctionBeginUser;
121 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
122 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
123 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
124
125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 Set runtime options
127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128 user.omega = 64.;
129 user.nts = 100;
130 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
131 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL);
132 PetscCall(PetscOptionsReal("-omega", "parameter", "<64>", user.omega, &user.omega, NULL));
133 PetscCall(PetscOptionsInt("-next_output", "time steps for next output point", "<100>", user.nts, &user.nts, NULL));
134 PetscOptionsEnd();
135
136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 Create necessary matrix and vectors, solve same ODE on every process
138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &U));
140 nindices[0] = 0;
141 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, nindices, PETSC_COPY_VALUES, &is1));
142 nindices[0] = 1;
143 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, nindices, PETSC_COPY_VALUES, &is2));
144
145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 Create timestepping solver context
147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
149 PetscCall(TSSetType(ts, TSBASICSYMPLECTIC));
150 PetscCall(TSRHSSplitSetIS(ts, "position", is1));
151 PetscCall(TSRHSSplitSetIS(ts, "momentum", is2));
152 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user));
153 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user));
154 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
155
156 PetscCall(TSSetMaxTime(ts, ftime));
157 PetscCall(TSSetTimeStep(ts, 0.0001));
158 PetscCall(TSSetMaxSteps(ts, 1000));
159 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
160 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
161
162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163 Set initial conditions
164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165 PetscCall(VecGetArray(U, &u_ptr));
166 u_ptr[0] = 0.2;
167 u_ptr[1] = 0.0;
168 PetscCall(VecRestoreArray(U, &u_ptr));
169
170 PetscCall(TSSetTime(ts, 0.0));
171
172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173 Set runtime options
174 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175 PetscCall(TSSetFromOptions(ts));
176
177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178 Solve nonlinear system
179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180 PetscCall(TSSolve(ts, U));
181 PetscCall(TSGetSolveTime(ts, &ftime));
182 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
183
184 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "The exact solution at time %.6lf is [%g %g]\n", (double)ftime, (double)(0.2 * PetscCosReal(user.omega * ftime)), (double)(-0.2 * user.omega * PetscSinReal(user.omega * ftime))));
185
186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187 Free work space. All PETSc objects should be destroyed when they
188 are no longer needed.
189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190 PetscCall(VecDestroy(&U));
191 PetscCall(TSDestroy(&ts));
192 PetscCall(ISDestroy(&is1));
193 PetscCall(ISDestroy(&is2));
194 PetscCall(PetscFinalize());
195 return 0;
196 }
197
198 /*TEST
199 build:
200 requires: !single !complex
201
202 test:
203 args: -ts_basicsymplectic_type 1 -monitor
204
205 test:
206 suffix: 2
207 args: -ts_basicsymplectic_type 2 -monitor
208
209 TEST*/
210