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*/