Actual source code: ex43.c
1: static char help[] = "Single-DOF oscillator formulated as a second-order system.\n";
3: #include <petscts.h>
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;
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: }
39: PetscErrorCode Solution(TS ts, PetscReal t, Vec X, void *ctx)
40: {
41: UserParams *user = (UserParams *)ctx;
42: PetscReal u, v;
43: PetscScalar *x;
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: }
53: PetscErrorCode Predictor(TS ts, Vec X0, Vec V0, Vec A0, Vec X1, void *ctx)
54: {
55: PetscReal dt, accel_fac;
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: }
67: PetscErrorCode Residual1(TS ts, PetscReal t, Vec U, Vec A, Vec R, void *ctx)
68: {
69: UserParams *user = (UserParams *)ctx;
70: PetscReal Omega = user->Omega;
71: const PetscScalar *u, *a;
72: PetscScalar *r;
74: PetscFunctionBeginUser;
75: PetscCall(VecGetArrayRead(U, &u));
76: PetscCall(VecGetArrayRead(A, &a));
77: PetscCall(VecGetArrayWrite(R, &r));
79: r[0] = a[0] + (Omega * Omega) * u[0];
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: }
89: PetscErrorCode Tangent1(TS ts, PetscReal t, Vec U, Vec A, PetscReal shiftA, Mat J, Mat P, void *ctx)
90: {
91: UserParams *user = (UserParams *)ctx;
92: PetscReal Omega = user->Omega;
93: PetscReal T = 0;
95: PetscFunctionBeginUser;
96: T = shiftA + (Omega * Omega);
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: }
108: PetscErrorCode Residual2(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, void *ctx)
109: {
110: UserParams *user = (UserParams *)ctx;
111: PetscReal Omega = user->Omega, Xi = user->Xi;
112: const PetscScalar *u, *v, *a;
113: PetscScalar *r;
115: PetscFunctionBeginUser;
116: PetscCall(VecGetArrayRead(U, &u));
117: PetscCall(VecGetArrayRead(V, &v));
118: PetscCall(VecGetArrayRead(A, &a));
119: PetscCall(VecGetArrayWrite(R, &r));
121: r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0];
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: }
132: PetscErrorCode Tangent2(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx)
133: {
134: UserParams *user = (UserParams *)ctx;
135: PetscReal Omega = user->Omega, Xi = user->Xi;
136: PetscReal T = 0;
138: PetscFunctionBeginUser;
139: T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega);
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: }
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};
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");
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();
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));
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));
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));
204: if (user.use_pred) PetscCall(TSAlpha2SetPredictor(ts, Predictor, NULL));
206: PetscCall(TS2SetSolution(ts, U, V));
207: PetscCall(TSSetFromOptions(ts));
208: PetscCall(TSSolve(ts, NULL));
210: PetscCall(VecDestroy(&U));
211: PetscCall(VecDestroy(&V));
212: PetscCall(TSDestroy(&ts));
213: PetscCall(PetscFinalize());
214: return 0;
215: }
217: /*TEST
219: test:
220: suffix: a
221: args: -ts_max_steps 10 -ts_view -use_pred
222: requires: !single
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
229: TEST*/