Actual source code: ex1.c
1: static char help[] = "Solves the motion of spring.\n\
2: Input parameters include:\n";
4: /* ------------------------------------------------------------------------
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.
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: ------------------------------------------------------------------------- */
26: #include <petscts.h>
27: #include <petscvec.h>
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: };
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: */
40: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
41: {
42: User user = (User)ctx;
43: const PetscScalar *x;
44: PetscScalar *vres;
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: }
55: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
56: {
57: const PetscScalar *v;
58: PetscScalar *xres;
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: }
69: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
70: {
71: User user = (User)ctx;
72: const PetscScalar *u;
73: PetscScalar *r;
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: }
85: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
86: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
87: {
88: const PetscScalar *u;
89: PetscReal dt;
90: PetscScalar energy, menergy;
91: User user = (User)ctx;
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: }
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;
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!");
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();
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));
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));
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));
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));
170: PetscCall(TSSetTime(ts, 0.0));
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Set runtime options
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: PetscCall(TSSetFromOptions(ts));
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Solve nonlinear system
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: PetscCall(TSSolve(ts, U));
181: PetscCall(TSGetSolveTime(ts, &ftime));
182: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
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))));
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: }
198: /*TEST
199: build:
200: requires: !single !complex
202: test:
203: args: -ts_basicsymplectic_type 1 -monitor
205: test:
206: suffix: 2
207: args: -ts_basicsymplectic_type 2 -monitor
209: TEST*/