Actual source code: ex1.c


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

  4: /*F

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

 11: F*/

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

 23: #include <petscts.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

179: /*TEST

181:    build:
182:      requires: !complex !single

184:    test:
185:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10

187:    test:
188:       suffix: 2
189:       args: -ts_max_steps 10

191: TEST*/