Actual source code: ex9.c


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

  4: /*F

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

 11:   Ensemble of initial conditions
 12:    ./ex2 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3      -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly

 14:   Fault at .1 seconds
 15:    ./ex2           -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly

 17:   Initial conditions same as when fault is ended
 18:    ./ex2 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05  -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly

 20: F*/

 22: /*
 23:    Include "petscts.h" so that we can use TS solvers.  Note that this
 24:    file automatically includes:
 25:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 26:      petscmat.h - matrices
 27:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 28:      petscviewer.h - viewers               petscpc.h  - preconditioners
 29:      petscksp.h   - linear solvers
 30: */

 32: #include <petscts.h>

 34: typedef struct {
 35:   PetscScalar H, D, omega_b, omega_s, Pmax, Pm, E, V, X;
 36:   PetscReal   tf, tcl;
 37: } AppCtx;

 39: /*
 40:      Defines the ODE passed to the ODE solver
 41: */
 42: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
 43: {
 44:   const PetscScalar *u;
 45:   PetscScalar       *f, Pmax;

 47:   /*  The next three lines allow us to access the entries of the vectors directly */
 48:   VecGetArrayRead(U, &u);
 49:   VecGetArray(F, &f);
 50:   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
 51:   else Pmax = ctx->Pmax;

 53:   f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
 54:   f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);

 56:   VecRestoreArrayRead(U, &u);
 57:   VecRestoreArray(F, &f);
 58:   return 0;
 59: }

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

 70:   VecGetArrayRead(U, &u);
 71:   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
 72:   else Pmax = ctx->Pmax;

 74:   J[0][0] = 0;
 75:   J[0][1] = ctx->omega_b;
 76:   J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
 77:   J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);

 79:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
 80:   VecRestoreArrayRead(U, &u);

 82:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 84:   if (A != B) {
 85:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 86:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 87:   }
 88:   return 0;
 89: }

 91: int main(int argc, char **argv)
 92: {
 93:   TS           ts; /* ODE integrator */
 94:   Vec          U;  /* solution will be stored here */
 95:   Mat          A;  /* Jacobian matrix */
 96:   PetscMPIInt  size;
 97:   PetscInt     n = 2;
 98:   AppCtx       ctx;
 99:   PetscScalar *u;
100:   PetscReal    du[2]    = {0.0, 0.0};
101:   PetscBool    ensemble = PETSC_FALSE, flg1, flg2;

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:      Initialize program
105:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107:   PetscInitialize(&argc, &argv, (char *)0, help);
108:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:     Create necessary matrix and vectors
113:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   MatCreate(PETSC_COMM_WORLD, &A);
115:   MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
116:   MatSetType(A, MATDENSE);
117:   MatSetFromOptions(A);
118:   MatSetUp(A);

120:   MatCreateVecs(A, &U, NULL);

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:     Set runtime options
124:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
126:   {
127:     ctx.omega_b = 1.0;
128:     ctx.omega_s = 2.0 * PETSC_PI * 60.0;
129:     ctx.H       = 5.0;
130:     PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL);
131:     ctx.D = 5.0;
132:     PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL);
133:     ctx.E    = 1.1378;
134:     ctx.V    = 1.0;
135:     ctx.X    = 0.545;
136:     ctx.Pmax = ctx.E * ctx.V / ctx.X;
137:     PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL);
138:     ctx.Pm = 0.9;
139:     PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL);
140:     ctx.tf  = 1.0;
141:     ctx.tcl = 1.05;
142:     PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL);
143:     PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL);
144:     PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL);
145:     if (ensemble) {
146:       ctx.tf  = -1;
147:       ctx.tcl = -1;
148:     }

150:     VecGetArray(U, &u);
151:     u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
152:     u[1] = 1.0;
153:     PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1);
154:     n = 2;
155:     PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2);
156:     u[0] += du[0];
157:     u[1] += du[1];
158:     VecRestoreArray(U, &u);
159:     if (flg1 || flg2) {
160:       ctx.tf  = -1;
161:       ctx.tcl = -1;
162:     }
163:   }
164:   PetscOptionsEnd();

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:      Create timestepping solver context
168:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:   TSCreate(PETSC_COMM_WORLD, &ts);
170:   TSSetProblemType(ts, TS_NONLINEAR);
171:   TSSetType(ts, TSTHETA);
172:   TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx);
173:   TSSetRHSJacobian(ts, A, A, (TSRHSJacobian)RHSJacobian, &ctx);

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Set initial conditions
177:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   TSSetSolution(ts, U);

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Set solver options
182:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183:   TSSetMaxTime(ts, 35.0);
184:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
185:   TSSetTimeStep(ts, .01);
186:   TSSetFromOptions(ts);

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:      Solve nonlinear system
190:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191:   if (ensemble) {
192:     for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
193:       VecGetArray(U, &u);
194:       u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
195:       u[1] = ctx.omega_s;
196:       u[0] += du[0];
197:       u[1] += du[1];
198:       VecRestoreArray(U, &u);
199:       TSSetTimeStep(ts, .01);
200:       TSSolve(ts, U);
201:     }
202:   } else {
203:     TSSolve(ts, U);
204:   }
205:   VecView(U, PETSC_VIEWER_STDOUT_WORLD);
206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
208:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209:   MatDestroy(&A);
210:   VecDestroy(&U);
211:   TSDestroy(&ts);
212:   PetscFinalize();
213:   return 0;
214: }

216: /*TEST

218:    build:
219:      requires: !complex

221:    test:

223: TEST*/