1 static char help[] = "Single-DOF oscillator formulated as a second-order system.\n";
2
3 #include <petscts.h>
4
5 typedef struct {
6 PetscReal Omega; /* natural frequency */
7 PetscReal Xi; /* damping coefficient */
8 PetscReal u0, v0; /* initial conditions */
9 PetscBool use_pred; /* whether to use a predictor callback */
10 } UserParams;
11
Exact(PetscReal t,PetscReal omega,PetscReal xi,PetscReal u0,PetscReal v0,PetscReal * ut,PetscReal * vt)12 static void Exact(PetscReal t, PetscReal omega, PetscReal xi, PetscReal u0, PetscReal v0, PetscReal *ut, PetscReal *vt)
13 {
14 PetscReal u, v;
15 if (xi < 1) {
16 PetscReal a = xi * omega;
17 PetscReal w = PetscSqrtReal(1 - xi * xi) * omega;
18 PetscReal C1 = (v0 + a * u0) / w;
19 PetscReal C2 = u0;
20 u = PetscExpReal(-a * t) * (C1 * PetscSinReal(w * t) + C2 * PetscCosReal(w * t));
21 v = (-a * PetscExpReal(-a * t) * (C1 * PetscSinReal(w * t) + C2 * PetscCosReal(w * t)) + w * PetscExpReal(-a * t) * (C1 * PetscCosReal(w * t) - C2 * PetscSinReal(w * t)));
22 } else if (xi > 1) {
23 PetscReal w = PetscSqrtReal(xi * xi - 1) * omega;
24 PetscReal C1 = (w * u0 + xi * u0 + v0) / (2 * w);
25 PetscReal C2 = (w * u0 - xi * u0 - v0) / (2 * w);
26 u = C1 * PetscExpReal((-xi + w) * t) + C2 * PetscExpReal((-xi - w) * t);
27 v = C1 * (-xi + w) * PetscExpReal((-xi + w) * t) + C2 * (-xi - w) * PetscExpReal((-xi - w) * t);
28 } else {
29 PetscReal a = xi * omega;
30 PetscReal C1 = v0 + a * u0;
31 PetscReal C2 = u0;
32 u = (C1 * t + C2) * PetscExpReal(-a * t);
33 v = (C1 - a * (C1 * t + C2)) * PetscExpReal(-a * t);
34 }
35 if (ut) *ut = u;
36 if (vt) *vt = v;
37 }
38
Solution(TS ts,PetscReal t,Vec X,PetscCtx ctx)39 PetscErrorCode Solution(TS ts, PetscReal t, Vec X, PetscCtx ctx)
40 {
41 UserParams *user = (UserParams *)ctx;
42 PetscReal u, v;
43 PetscScalar *x;
44
45 PetscFunctionBeginUser;
46 Exact(t, user->Omega, user->Xi, user->u0, user->v0, &u, &v);
47 PetscCall(VecGetArray(X, &x));
48 x[0] = (PetscScalar)u;
49 PetscCall(VecRestoreArray(X, &x));
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
Predictor(TS ts,Vec X0,Vec V0,Vec A0,Vec X1,PetscCtx ctx)53 PetscErrorCode Predictor(TS ts, Vec X0, Vec V0, Vec A0, Vec X1, PetscCtx ctx)
54 {
55 PetscReal dt, accel_fac;
56
57 PetscFunctionBeginUser;
58 PetscCall(TSGetTimeStep(ts, &dt));
59 accel_fac = 0.5 * dt * dt;
60 /* X1 = X0 + dt*V0 + accel_fac*A0 */
61 PetscCall(VecCopy(X0, X1));
62 PetscCall(VecAXPY(X1, dt, V0));
63 PetscCall(VecAXPY(X1, accel_fac, A0));
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
Residual1(TS ts,PetscReal t,Vec U,Vec A,Vec R,PetscCtx ctx)67 PetscErrorCode Residual1(TS ts, PetscReal t, Vec U, Vec A, Vec R, PetscCtx ctx)
68 {
69 UserParams *user = (UserParams *)ctx;
70 PetscReal Omega = user->Omega;
71 const PetscScalar *u, *a;
72 PetscScalar *r;
73
74 PetscFunctionBeginUser;
75 PetscCall(VecGetArrayRead(U, &u));
76 PetscCall(VecGetArrayRead(A, &a));
77 PetscCall(VecGetArrayWrite(R, &r));
78
79 r[0] = a[0] + (Omega * Omega) * u[0];
80
81 PetscCall(VecRestoreArrayRead(U, &u));
82 PetscCall(VecRestoreArrayRead(A, &a));
83 PetscCall(VecRestoreArrayWrite(R, &r));
84 PetscCall(VecAssemblyBegin(R));
85 PetscCall(VecAssemblyEnd(R));
86 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
Tangent1(TS ts,PetscReal t,Vec U,Vec A,PetscReal shiftA,Mat J,Mat P,PetscCtx ctx)89 PetscErrorCode Tangent1(TS ts, PetscReal t, Vec U, Vec A, PetscReal shiftA, Mat J, Mat P, PetscCtx ctx)
90 {
91 UserParams *user = (UserParams *)ctx;
92 PetscReal Omega = user->Omega;
93 PetscReal T = 0;
94
95 PetscFunctionBeginUser;
96 T = shiftA + (Omega * Omega);
97
98 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES));
99 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
100 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
101 if (J != P) {
102 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
103 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
104 }
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
Residual2(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec R,PetscCtx ctx)108 PetscErrorCode Residual2(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, PetscCtx ctx)
109 {
110 UserParams *user = (UserParams *)ctx;
111 PetscReal Omega = user->Omega, Xi = user->Xi;
112 const PetscScalar *u, *v, *a;
113 PetscScalar *r;
114
115 PetscFunctionBeginUser;
116 PetscCall(VecGetArrayRead(U, &u));
117 PetscCall(VecGetArrayRead(V, &v));
118 PetscCall(VecGetArrayRead(A, &a));
119 PetscCall(VecGetArrayWrite(R, &r));
120
121 r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0];
122
123 PetscCall(VecRestoreArrayRead(U, &u));
124 PetscCall(VecRestoreArrayRead(V, &v));
125 PetscCall(VecRestoreArrayRead(A, &a));
126 PetscCall(VecRestoreArrayWrite(R, &r));
127 PetscCall(VecAssemblyBegin(R));
128 PetscCall(VecAssemblyEnd(R));
129 PetscFunctionReturn(PETSC_SUCCESS);
130 }
131
Tangent2(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P,PetscCtx ctx)132 PetscErrorCode Tangent2(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, PetscCtx ctx)
133 {
134 UserParams *user = (UserParams *)ctx;
135 PetscReal Omega = user->Omega, Xi = user->Xi;
136 PetscReal T = 0;
137
138 PetscFunctionBeginUser;
139 T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega);
140
141 PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES));
142 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
143 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
144 if (J != P) {
145 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
146 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
147 }
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
main(int argc,char * argv[])151 int main(int argc, char *argv[])
152 {
153 PetscMPIInt size;
154 TS ts;
155 Vec R;
156 Mat J;
157 Vec U, V;
158 PetscScalar *u, *v;
159 UserParams user = {/*Omega=*/1, /*Xi=*/0, /*u0=*/1, /*,v0=*/0, /*,use_pred=*/PETSC_FALSE};
160
161 PetscFunctionBeginUser;
162 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
163 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
164 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
165
166 PetscOptionsBegin(PETSC_COMM_SELF, "", "ex43 options", "");
167 PetscCall(PetscOptionsReal("-frequency", "Natual frequency", __FILE__, user.Omega, &user.Omega, NULL));
168 PetscCall(PetscOptionsReal("-damping", "Damping coefficient", __FILE__, user.Xi, &user.Xi, NULL));
169 PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL));
170 PetscCall(PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL));
171 PetscCall(PetscOptionsBool("-use_pred", "Use a custom predictor", __FILE__, user.use_pred, &user.use_pred, NULL));
172 PetscOptionsEnd();
173
174 PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
175 PetscCall(TSSetType(ts, TSALPHA2));
176 PetscCall(TSSetMaxTime(ts, 5 * (2 * PETSC_PI)));
177 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
178 PetscCall(TSSetTimeStep(ts, 0.01));
179
180 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R));
181 PetscCall(VecSetUp(R));
182 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J));
183 PetscCall(MatSetUp(J));
184 if (user.Xi) {
185 PetscCall(TSSetI2Function(ts, R, Residual2, &user));
186 PetscCall(TSSetI2Jacobian(ts, J, J, Tangent2, &user));
187 } else {
188 PetscCall(TSSetIFunction(ts, R, Residual1, &user));
189 PetscCall(TSSetIJacobian(ts, J, J, Tangent1, &user));
190 }
191 PetscCall(VecDestroy(&R));
192 PetscCall(MatDestroy(&J));
193 PetscCall(TSSetSolutionFunction(ts, Solution, &user));
194
195 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U));
196 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &V));
197 PetscCall(VecGetArrayWrite(U, &u));
198 PetscCall(VecGetArrayWrite(V, &v));
199 u[0] = user.u0;
200 v[0] = user.v0;
201 PetscCall(VecRestoreArrayWrite(U, &u));
202 PetscCall(VecRestoreArrayWrite(V, &v));
203
204 if (user.use_pred) PetscCall(TSAlpha2SetPredictor(ts, Predictor, NULL));
205
206 PetscCall(TS2SetSolution(ts, U, V));
207 PetscCall(TSSetFromOptions(ts));
208 PetscCall(TSSolve(ts, NULL));
209
210 PetscCall(VecDestroy(&U));
211 PetscCall(VecDestroy(&V));
212 PetscCall(TSDestroy(&ts));
213 PetscCall(PetscFinalize());
214 return 0;
215 }
216
217 /*TEST
218
219 test:
220 suffix: a
221 args: -ts_max_steps 10 -ts_view -use_pred
222 requires: !single
223
224 test:
225 suffix: b
226 args: -ts_max_steps 10 -ts_rtol 0 -ts_atol 1e-5 -ts_adapt_type basic -ts_adapt_monitor
227 requires: !single
228
229 TEST*/
230