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: } UserParams;
11: static void Exact(PetscReal t, PetscReal omega, PetscReal xi, PetscReal u0, PetscReal v0, PetscReal *ut, PetscReal *vt)
12: {
13: PetscReal u, v;
14: if (xi < 1) {
15: PetscReal a = xi * omega;
16: PetscReal w = PetscSqrtReal(1 - xi * xi) * omega;
17: PetscReal C1 = (v0 + a * u0) / w;
18: PetscReal C2 = u0;
19: u = PetscExpReal(-a * t) * (C1 * PetscSinReal(w * t) + C2 * PetscCosReal(w * t));
20: 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)));
21: } else if (xi > 1) {
22: PetscReal w = PetscSqrtReal(xi * xi - 1) * omega;
23: PetscReal C1 = (w * u0 + xi * u0 + v0) / (2 * w);
24: PetscReal C2 = (w * u0 - xi * u0 - v0) / (2 * w);
25: u = C1 * PetscExpReal((-xi + w) * t) + C2 * PetscExpReal((-xi - w) * t);
26: v = C1 * (-xi + w) * PetscExpReal((-xi + w) * t) + C2 * (-xi - w) * PetscExpReal((-xi - w) * t);
27: } else {
28: PetscReal a = xi * omega;
29: PetscReal C1 = v0 + a * u0;
30: PetscReal C2 = u0;
31: u = (C1 * t + C2) * PetscExpReal(-a * t);
32: v = (C1 - a * (C1 * t + C2)) * PetscExpReal(-a * t);
33: }
34: if (ut) *ut = u;
35: if (vt) *vt = v;
36: }
38: PetscErrorCode Solution(TS ts, PetscReal t, Vec X, void *ctx)
39: {
40: UserParams *user = (UserParams *)ctx;
41: PetscReal u, v;
42: PetscScalar *x;
44: PetscFunctionBeginUser;
45: Exact(t, user->Omega, user->Xi, user->u0, user->v0, &u, &v);
46: PetscCall(VecGetArray(X, &x));
47: x[0] = (PetscScalar)u;
48: PetscCall(VecRestoreArray(X, &x));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: PetscErrorCode Residual1(TS ts, PetscReal t, Vec U, Vec A, Vec R, void *ctx)
53: {
54: UserParams *user = (UserParams *)ctx;
55: PetscReal Omega = user->Omega;
56: const PetscScalar *u, *a;
57: PetscScalar *r;
59: PetscFunctionBeginUser;
60: PetscCall(VecGetArrayRead(U, &u));
61: PetscCall(VecGetArrayRead(A, &a));
62: PetscCall(VecGetArrayWrite(R, &r));
64: r[0] = a[0] + (Omega * Omega) * u[0];
66: PetscCall(VecRestoreArrayRead(U, &u));
67: PetscCall(VecRestoreArrayRead(A, &a));
68: PetscCall(VecRestoreArrayWrite(R, &r));
69: PetscCall(VecAssemblyBegin(R));
70: PetscCall(VecAssemblyEnd(R));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: PetscErrorCode Tangent1(TS ts, PetscReal t, Vec U, Vec A, PetscReal shiftA, Mat J, Mat P, void *ctx)
75: {
76: UserParams *user = (UserParams *)ctx;
77: PetscReal Omega = user->Omega;
78: PetscReal T = 0;
80: PetscFunctionBeginUser;
82: T = shiftA + (Omega * Omega);
84: PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES));
85: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
87: if (J != P) {
88: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
89: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
90: }
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode Residual2(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec R, void *ctx)
95: {
96: UserParams *user = (UserParams *)ctx;
97: PetscReal Omega = user->Omega, Xi = user->Xi;
98: const PetscScalar *u, *v, *a;
99: PetscScalar *r;
101: PetscFunctionBeginUser;
102: PetscCall(VecGetArrayRead(U, &u));
103: PetscCall(VecGetArrayRead(V, &v));
104: PetscCall(VecGetArrayRead(A, &a));
105: PetscCall(VecGetArrayWrite(R, &r));
107: r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0];
109: PetscCall(VecRestoreArrayRead(U, &u));
110: PetscCall(VecRestoreArrayRead(V, &v));
111: PetscCall(VecRestoreArrayRead(A, &a));
112: PetscCall(VecRestoreArrayWrite(R, &r));
113: PetscCall(VecAssemblyBegin(R));
114: PetscCall(VecAssemblyEnd(R));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: PetscErrorCode Tangent2(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P, void *ctx)
119: {
120: UserParams *user = (UserParams *)ctx;
121: PetscReal Omega = user->Omega, Xi = user->Xi;
122: PetscReal T = 0;
124: PetscFunctionBeginUser;
126: T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega);
128: PetscCall(MatSetValue(P, 0, 0, T, INSERT_VALUES));
129: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
130: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
131: if (J != P) {
132: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
133: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
134: }
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: int main(int argc, char *argv[])
139: {
140: PetscMPIInt size;
141: TS ts;
142: Vec R;
143: Mat J;
144: Vec U, V;
145: PetscScalar *u, *v;
146: UserParams user = {/*Omega=*/1, /*Xi=*/0, /*u0=*/1, /*,v0=*/0};
148: PetscFunctionBeginUser;
149: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
150: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
151: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
153: PetscOptionsBegin(PETSC_COMM_SELF, "", "ex43 options", "");
154: PetscCall(PetscOptionsReal("-frequency", "Natual frequency", __FILE__, user.Omega, &user.Omega, NULL));
155: PetscCall(PetscOptionsReal("-damping", "Damping coefficient", __FILE__, user.Xi, &user.Xi, NULL));
156: PetscCall(PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL));
157: PetscCall(PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL));
158: PetscOptionsEnd();
160: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
161: PetscCall(TSSetType(ts, TSALPHA2));
162: PetscCall(TSSetMaxTime(ts, 5 * (2 * PETSC_PI)));
163: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
164: PetscCall(TSSetTimeStep(ts, 0.01));
166: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &R));
167: PetscCall(VecSetUp(R));
168: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J));
169: PetscCall(MatSetUp(J));
170: if (user.Xi) {
171: PetscCall(TSSetI2Function(ts, R, Residual2, &user));
172: PetscCall(TSSetI2Jacobian(ts, J, J, Tangent2, &user));
173: } else {
174: PetscCall(TSSetIFunction(ts, R, Residual1, &user));
175: PetscCall(TSSetIJacobian(ts, J, J, Tangent1, &user));
176: }
177: PetscCall(VecDestroy(&R));
178: PetscCall(MatDestroy(&J));
179: PetscCall(TSSetSolutionFunction(ts, Solution, &user));
181: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &U));
182: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &V));
183: PetscCall(VecGetArrayWrite(U, &u));
184: PetscCall(VecGetArrayWrite(V, &v));
185: u[0] = user.u0;
186: v[0] = user.v0;
187: PetscCall(VecRestoreArrayWrite(U, &u));
188: PetscCall(VecRestoreArrayWrite(V, &v));
190: PetscCall(TS2SetSolution(ts, U, V));
191: PetscCall(TSSetFromOptions(ts));
192: PetscCall(TSSolve(ts, NULL));
194: PetscCall(VecDestroy(&U));
195: PetscCall(VecDestroy(&V));
196: PetscCall(TSDestroy(&ts));
197: PetscCall(PetscFinalize());
198: return 0;
199: }
201: /*TEST
203: test:
204: suffix: a
205: args: -ts_max_steps 10 -ts_view
206: requires: !single
208: test:
209: suffix: b
210: args: -ts_max_steps 10 -ts_rtol 0 -ts_atol 1e-5 -ts_adapt_type basic -ts_adapt_monitor
211: requires: !single
213: TEST*/