Actual source code: ex1.c

      1: static char help[] = "Basic equation for generator stability analysis.\n";

      3: /*F

      5: \begin{eqnarray}
      6:                  \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - \frac{EV}{X} \sin(\theta) \\
      7:                  \frac{d \theta}{dt} = \omega - \omega_s
      8: \end{eqnarray}

     10: F*/

     12: /*
     13:    Include "petscts.h" so that we can use TS solvers.  Note that this
     14:    file automatically includes:
     15:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     16:      petscmat.h - matrices
     17:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     18:      petscviewer.h - viewers               petscpc.h  - preconditioners
     19:      petscksp.h   - linear solvers
     20: */

     22: #include <petscts.h>

     24: typedef struct {
     25:   PetscScalar H, omega_s, E, V, X;
     26:   PetscRandom rand;
     27: } AppCtx;

     29: /*
     30:      Defines the ODE passed to the ODE solver
     31: */
     32: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
     33: {
     34:   PetscScalar       *f, r;
     35:   const PetscScalar *u, *udot;
     36:   static PetscScalar R = .4;

     38:   PetscFunctionBegin;
     39:   PetscCall(PetscRandomGetValue(ctx->rand, &r));
     40:   if (r > .9) R = .5;
     41:   if (r < .1) R = .4;
     42:   R = .4;
     43:   /*  The next three lines allow us to access the entries of the vectors directly */
     44:   PetscCall(VecGetArrayRead(U, &u));
     45:   PetscCall(VecGetArrayRead(Udot, &udot));
     46:   PetscCall(VecGetArray(F, &f));
     47:   f[0] = 2.0 * ctx->H * udot[0] / ctx->omega_s + ctx->E * ctx->V * PetscSinScalar(u[1]) / ctx->X - R;
     48:   f[1] = udot[1] - u[0] + ctx->omega_s;

     50:   PetscCall(VecRestoreArrayRead(U, &u));
     51:   PetscCall(VecRestoreArrayRead(Udot, &udot));
     52:   PetscCall(VecRestoreArray(F, &f));
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: /*
     57:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
     58: */
     59: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
     60: {
     61:   PetscInt           rowcol[] = {0, 1};
     62:   PetscScalar        J[2][2];
     63:   const PetscScalar *u, *udot;

     65:   PetscFunctionBegin;
     66:   PetscCall(VecGetArrayRead(U, &u));
     67:   PetscCall(VecGetArrayRead(Udot, &udot));
     68:   J[0][0] = 2.0 * ctx->H * a / ctx->omega_s;
     69:   J[0][1] = -ctx->E * ctx->V * PetscCosScalar(u[1]) / ctx->X;
     70:   J[1][0] = -1.0;
     71:   J[1][1] = a;
     72:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
     73:   PetscCall(VecRestoreArrayRead(U, &u));
     74:   PetscCall(VecRestoreArrayRead(Udot, &udot));

     76:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     77:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     78:   if (A != B) {
     79:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     80:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     81:   }
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: int main(int argc, char **argv)
     86: {
     87:   TS           ts; /* ODE integrator */
     88:   Vec          U;  /* solution will be stored here */
     89:   Mat          A;  /* Jacobian matrix */
     90:   PetscMPIInt  size;
     91:   PetscInt     n = 2;
     92:   AppCtx       ctx;
     93:   PetscScalar *u;

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Initialize program
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   PetscFunctionBeginUser;
     99:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    100:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    101:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    104:     Create necessary matrix and vectors
    105:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    106:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    107:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
    108:   PetscCall(MatSetFromOptions(A));
    109:   PetscCall(MatSetUp(A));

    111:   PetscCall(MatCreateVecs(A, &U, NULL));

    113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114:     Set runtime options
    115:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    116:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Reaction options", "");
    117:   {
    118:     ctx.omega_s = 1.0;
    119:     PetscCall(PetscOptionsScalar("-omega_s", "", "", ctx.omega_s, &ctx.omega_s, NULL));
    120:     ctx.H = 1.0;
    121:     PetscCall(PetscOptionsScalar("-H", "", "", ctx.H, &ctx.H, NULL));
    122:     ctx.E = 1.0;
    123:     PetscCall(PetscOptionsScalar("-E", "", "", ctx.E, &ctx.E, NULL));
    124:     ctx.V = 1.0;
    125:     PetscCall(PetscOptionsScalar("-V", "", "", ctx.V, &ctx.V, NULL));
    126:     ctx.X = 1.0;
    127:     PetscCall(PetscOptionsScalar("-X", "", "", ctx.X, &ctx.X, NULL));

    129:     PetscCall(VecGetArray(U, &u));
    130:     u[0] = 1;
    131:     u[1] = .7;
    132:     PetscCall(VecRestoreArray(U, &u));
    133:     PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", U, NULL));
    134:   }
    135:   PetscOptionsEnd();

    137:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ctx.rand));
    138:   PetscCall(PetscRandomSetFromOptions(ctx.rand));

    140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141:      Create timestepping solver context
    142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    143:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    144:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    145:   PetscCall(TSSetType(ts, TSROSW));
    146:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx));
    147:   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));

    149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    150:      Set initial conditions
    151:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    152:   PetscCall(TSSetSolution(ts, U));

    154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    155:      Set solver options
    156:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    157:   PetscCall(TSSetMaxTime(ts, 2000.0));
    158:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    159:   PetscCall(TSSetTimeStep(ts, .001));
    160:   PetscCall(TSSetFromOptions(ts));

    162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    163:      Solve nonlinear system
    164:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    165:   PetscCall(TSSolve(ts, U));

    167:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    168:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    169:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    170:   PetscCall(MatDestroy(&A));
    171:   PetscCall(VecDestroy(&U));
    172:   PetscCall(TSDestroy(&ts));
    173:   PetscCall(PetscRandomDestroy(&ctx.rand));
    174:   PetscCall(PetscFinalize());
    175:   return 0;
    176: }

    178: /*TEST

    180:    build:
    181:      requires: !complex !single

    183:    test:
    184:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
    185:       output_file: output/empty.out

    187:    test:
    188:       suffix: 2
    189:       args: -ts_max_steps 10
    190:       output_file: output/empty.out

    192: TEST*/