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;
45: Exact(t, user->Omega, user->Xi, user->u0, user->v0, &u, &v);
46: VecGetArray(X, &x);
47: x[0] = (PetscScalar)u;
48: VecRestoreArray(X, &x);
49: return 0;
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;
60: VecGetArrayRead(U, &u);
61: VecGetArrayRead(A, &a);
62: VecGetArrayWrite(R, &r);
64: r[0] = a[0] + (Omega * Omega) * u[0];
66: VecRestoreArrayRead(U, &u);
67: VecRestoreArrayRead(A, &a);
68: VecRestoreArrayWrite(R, &r);
69: VecAssemblyBegin(R);
70: VecAssemblyEnd(R);
71: return 0;
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;
82: T = shiftA + (Omega * Omega);
84: MatSetValue(P, 0, 0, T, INSERT_VALUES);
85: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
87: if (J != P) {
88: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
90: }
91: return 0;
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;
102: VecGetArrayRead(U, &u);
103: VecGetArrayRead(V, &v);
104: VecGetArrayRead(A, &a);
105: VecGetArrayWrite(R, &r);
107: r[0] = a[0] + (2 * Xi * Omega) * v[0] + (Omega * Omega) * u[0];
109: VecRestoreArrayRead(U, &u);
110: VecRestoreArrayRead(V, &v);
111: VecRestoreArrayRead(A, &a);
112: VecRestoreArrayWrite(R, &r);
113: VecAssemblyBegin(R);
114: VecAssemblyEnd(R);
115: return 0;
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;
126: T = shiftA + shiftV * (2 * Xi * Omega) + (Omega * Omega);
128: MatSetValue(P, 0, 0, T, INSERT_VALUES);
129: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
131: if (J != P) {
132: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
134: }
135: return 0;
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};
149: PetscInitialize(&argc, &argv, NULL, help);
150: MPI_Comm_size(PETSC_COMM_WORLD, &size);
153: PetscOptionsBegin(PETSC_COMM_SELF, "", "ex43 options", "");
154: PetscOptionsReal("-frequency", "Natual frequency", __FILE__, user.Omega, &user.Omega, NULL);
155: PetscOptionsReal("-damping", "Damping coefficient", __FILE__, user.Xi, &user.Xi, NULL);
156: PetscOptionsReal("-initial_u", "Initial displacement", __FILE__, user.u0, &user.u0, NULL);
157: PetscOptionsReal("-initial_v", "Initial velocity", __FILE__, user.v0, &user.v0, NULL);
158: PetscOptionsEnd();
160: TSCreate(PETSC_COMM_SELF, &ts);
161: TSSetType(ts, TSALPHA2);
162: TSSetMaxTime(ts, 5 * (2 * PETSC_PI));
163: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
164: TSSetTimeStep(ts, 0.01);
166: VecCreateSeq(PETSC_COMM_SELF, 1, &R);
167: VecSetUp(R);
168: MatCreateSeqDense(PETSC_COMM_SELF, 1, 1, NULL, &J);
169: MatSetUp(J);
170: if (user.Xi) {
171: TSSetI2Function(ts, R, Residual2, &user);
172: TSSetI2Jacobian(ts, J, J, Tangent2, &user);
173: } else {
174: TSSetIFunction(ts, R, Residual1, &user);
175: TSSetIJacobian(ts, J, J, Tangent1, &user);
176: }
177: VecDestroy(&R);
178: MatDestroy(&J);
179: TSSetSolutionFunction(ts, Solution, &user);
181: VecCreateSeq(PETSC_COMM_SELF, 1, &U);
182: VecCreateSeq(PETSC_COMM_SELF, 1, &V);
183: VecGetArrayWrite(U, &u);
184: VecGetArrayWrite(V, &v);
185: u[0] = user.u0;
186: v[0] = user.v0;
187: VecRestoreArrayWrite(U, &u);
188: VecRestoreArrayWrite(V, &v);
190: TS2SetSolution(ts, U, V);
191: TSSetFromOptions(ts);
192: TSSolve(ts, NULL);
194: VecDestroy(&U);
195: VecDestroy(&V);
196: TSDestroy(&ts);
197: 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*/