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