1 static char help[] = "Basic equation for generator stability analysis.\n";
2
3 /*F
4
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}
9
10 F*/
11
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 */
21
22 #include <petscts.h>
23
24 typedef struct {
25 PetscScalar H, omega_s, E, V, X;
26 PetscRandom rand;
27 } AppCtx;
28
29 /*
30 Defines the ODE passed to the ODE solver
31 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)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;
37
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;
49
50 PetscCall(VecRestoreArrayRead(U, &u));
51 PetscCall(VecRestoreArrayRead(Udot, &udot));
52 PetscCall(VecRestoreArray(F, &f));
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55
56 /*
57 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
58 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)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;
64
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));
75
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 }
84
main(int argc,char ** argv)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;
94
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");
102
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));
110
111 PetscCall(MatCreateVecs(A, &U, NULL));
112
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));
128
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();
136
137 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ctx.rand));
138 PetscCall(PetscRandomSetFromOptions(ctx.rand));
139
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));
148
149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150 Set initial conditions
151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152 PetscCall(TSSetSolution(ts, U));
153
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));
161
162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163 Solve nonlinear system
164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165 PetscCall(TSSolve(ts, U));
166
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 }
177
178 /*TEST
179
180 build:
181 requires: !complex !single
182
183 test:
184 args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
185 output_file: output/empty.out
186
187 test:
188 suffix: 2
189 args: -ts_max_steps 10
190 output_file: output/empty.out
191
192 TEST*/
193