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*/