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:   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:   VecGetArrayRead(U, &u);
 45:   VecGetArrayRead(Udot, &udot);
 46:   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:   VecRestoreArrayRead(U, &u);
 51:   VecRestoreArrayRead(Udot, &udot);
 52:   VecRestoreArray(F, &f);
 53:   return 0;
 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:   VecGetArrayRead(U, &u);
 66:   VecGetArrayRead(Udot, &udot);
 67:   J[0][0] = 2.0 * ctx->H * a / ctx->omega_s;
 68:   J[0][1] = -ctx->E * ctx->V * PetscCosScalar(u[1]) / ctx->X;
 69:   J[1][0] = -1.0;
 70:   J[1][1] = a;
 71:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
 72:   VecRestoreArrayRead(U, &u);
 73:   VecRestoreArrayRead(Udot, &udot);

 75:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 77:   if (A != B) {
 78:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 79:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 80:   }
 81:   return 0;
 82: }

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

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Initialize program
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   PetscInitialize(&argc, &argv, (char *)0, help);
 99:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

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

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

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

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

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

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

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

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

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

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

177: /*TEST

179:    build:
180:      requires: !complex !single

182:    test:
183:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10

185:    test:
186:       suffix: 2
187:       args: -ts_max_steps 10

189: TEST*/