Actual source code: ex1.c

      1: static char help[] = "Solves the motion of spring.\n\
      2: Input parameters include:\n";

      4: /* ------------------------------------------------------------------------

      6:   This program solves the motion of spring by Hooke's law
      7:   x' = f(t,v) = v
      8:   v' = g(t,x) = -omega^2*x
      9:   on the interval 0 <= t <= 0.1, with the initial conditions
     10:     x(0) = 0.2, x'(0) = v(0) = 0,
     11:   and
     12:     omega = 64.
     13:   The exact solution is
     14:     x(t) = A*sin(t*omega) + B*cos(t*omega)
     15:   where A and B are constants that can be determined from the initial conditions.
     16:   In this case, B=0.2, A=0.

     18:   Notes:
     19:   This code demonstrates the TS solver interface to solve a separable Hamiltonian
     20:   system, which can be split into two subsystems involving two coupling components,
     21:   named generailized momentum and generailized position respectively.
     22:   Using a symplectic intergrator can preserve energy
     23:   E = (v^2+omega^2*x^2-omega^2*h*v*x)/2
     24:   ------------------------------------------------------------------------- */

     26: #include <petscts.h>
     27: #include <petscvec.h>

     29: typedef struct _n_User *User;
     30: struct _n_User {
     31:   PetscReal omega;
     32:   PetscInt  nts; /* print the energy at each nts time steps */
     33: };

     35: /*
     36:   User-defined routines.
     37:   The first RHS function provides f(t,x), the residual for the generalized momentum,
     38:   and the second one provides g(t,v), the residual for the generalized position.
     39: */
     40: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
     41: {
     42:   User               user = (User)ctx;
     43:   const PetscScalar *x;
     44:   PetscScalar       *vres;

     46:   PetscFunctionBeginUser;
     47:   PetscCall(VecGetArrayRead(X, &x));
     48:   PetscCall(VecGetArray(Vres, &vres));
     49:   vres[0] = -user->omega * user->omega * x[0];
     50:   PetscCall(VecRestoreArray(Vres, &vres));
     51:   PetscCall(VecRestoreArrayRead(X, &x));
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
     56: {
     57:   const PetscScalar *v;
     58:   PetscScalar       *xres;

     60:   PetscFunctionBeginUser;
     61:   PetscCall(VecGetArray(Xres, &xres));
     62:   PetscCall(VecGetArrayRead(V, &v));
     63:   xres[0] = v[0];
     64:   PetscCall(VecRestoreArrayRead(V, &v));
     65:   PetscCall(VecRestoreArray(Xres, &xres));
     66:   PetscFunctionReturn(PETSC_SUCCESS);
     67: }

     69: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
     70: {
     71:   User               user = (User)ctx;
     72:   const PetscScalar *u;
     73:   PetscScalar       *r;

     75:   PetscFunctionBeginUser;
     76:   PetscCall(VecGetArrayRead(U, &u));
     77:   PetscCall(VecGetArray(R, &r));
     78:   r[0] = u[1];
     79:   r[1] = -user->omega * user->omega * u[0];
     80:   PetscCall(VecRestoreArrayRead(U, &u));
     81:   PetscCall(VecRestoreArray(R, &r));
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
     86: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
     87: {
     88:   const PetscScalar *u;
     89:   PetscReal          dt;
     90:   PetscScalar        energy, menergy;
     91:   User               user = (User)ctx;

     93:   PetscFunctionBeginUser;
     94:   if (step % user->nts == 0) {
     95:     PetscCall(TSGetTimeStep(ts, &dt));
     96:     PetscCall(VecGetArrayRead(U, &u));
     97:     menergy = (u[1] * u[1] + user->omega * user->omega * u[0] * u[0] - user->omega * user->omega * dt * u[0] * u[1]) / 2.;
     98:     energy  = (u[1] * u[1] + user->omega * user->omega * u[0] * u[0]) / 2.;
     99:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At time %.6lf, Energy = %8g, Modified Energy = %8g\n", (double)t, (double)energy, (double)menergy));
    100:     PetscCall(VecRestoreArrayRead(U, &u));
    101:   }
    102:   PetscFunctionReturn(PETSC_SUCCESS);
    103: }

    105: int main(int argc, char **argv)
    106: {
    107:   TS             ts; /* nonlinear solver */
    108:   Vec            U;  /* solution, residual vectors */
    109:   IS             is1, is2;
    110:   PetscInt       nindices[1];
    111:   PetscReal      ftime   = 0.1;
    112:   PetscBool      monitor = PETSC_FALSE;
    113:   PetscScalar   *u_ptr;
    114:   PetscMPIInt    size;
    115:   struct _n_User user;

    117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118:      Initialize program
    119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    120:   PetscFunctionBeginUser;
    121:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    122:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    123:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

    125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    126:     Set runtime options
    127:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    128:   user.omega = 64.;
    129:   user.nts   = 100;
    130:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
    131:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL);
    132:   PetscCall(PetscOptionsReal("-omega", "parameter", "<64>", user.omega, &user.omega, NULL));
    133:   PetscCall(PetscOptionsInt("-next_output", "time steps for next output point", "<100>", user.nts, &user.nts, NULL));
    134:   PetscOptionsEnd();

    136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    137:     Create necessary matrix and vectors, solve same ODE on every process
    138:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    139:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &U));
    140:   nindices[0] = 0;
    141:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, nindices, PETSC_COPY_VALUES, &is1));
    142:   nindices[0] = 1;
    143:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, nindices, PETSC_COPY_VALUES, &is2));

    145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    146:      Create timestepping solver context
    147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    148:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    149:   PetscCall(TSSetType(ts, TSBASICSYMPLECTIC));
    150:   PetscCall(TSRHSSplitSetIS(ts, "position", is1));
    151:   PetscCall(TSRHSSplitSetIS(ts, "momentum", is2));
    152:   PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user));
    153:   PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user));
    154:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));

    156:   PetscCall(TSSetMaxTime(ts, ftime));
    157:   PetscCall(TSSetTimeStep(ts, 0.0001));
    158:   PetscCall(TSSetMaxSteps(ts, 1000));
    159:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    160:   if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));

    162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    163:      Set initial conditions
    164:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    165:   PetscCall(VecGetArray(U, &u_ptr));
    166:   u_ptr[0] = 0.2;
    167:   u_ptr[1] = 0.0;
    168:   PetscCall(VecRestoreArray(U, &u_ptr));

    170:   PetscCall(TSSetTime(ts, 0.0));

    172:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    173:      Set runtime options
    174:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    175:   PetscCall(TSSetFromOptions(ts));

    177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    178:      Solve nonlinear system
    179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    180:   PetscCall(TSSolve(ts, U));
    181:   PetscCall(TSGetSolveTime(ts, &ftime));
    182:   PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));

    184:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "The exact solution at time %.6lf is [%g %g]\n", (double)ftime, (double)(0.2 * PetscCosReal(user.omega * ftime)), (double)(-0.2 * user.omega * PetscSinReal(user.omega * ftime))));

    186:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    187:      Free work space.  All PETSc objects should be destroyed when they
    188:      are no longer needed.
    189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    190:   PetscCall(VecDestroy(&U));
    191:   PetscCall(TSDestroy(&ts));
    192:   PetscCall(ISDestroy(&is1));
    193:   PetscCall(ISDestroy(&is2));
    194:   PetscCall(PetscFinalize());
    195:   return 0;
    196: }

    198: /*TEST
    199:    build:
    200:      requires: !single !complex

    202:    test:
    203:      args: -ts_basicsymplectic_type 1 -monitor

    205:    test:
    206:      suffix: 2
    207:      args: -ts_basicsymplectic_type 2 -monitor

    209: TEST*/