Actual source code: ex9adj.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: ./ex9 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
14: Fault at .1 seconds
15: ./ex9 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
17: Initial conditions same as when fault is ended
18: ./ex9 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -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, u_s, c;
36: PetscInt beta;
37: PetscReal tf, tcl;
38: } AppCtx;
40: PetscErrorCode PostStepFunction(TS ts)
41: {
42: Vec U;
43: PetscReal t;
44: const PetscScalar *u;
46: PetscFunctionBegin;
47: PetscCall(TSGetTime(ts, &t));
48: PetscCall(TSGetSolution(ts, &U));
49: PetscCall(VecGetArrayRead(U, &u));
50: PetscCall(PetscPrintf(PETSC_COMM_SELF, "delta(%3.2f) = %8.7f\n", (double)t, (double)u[0]));
51: PetscCall(VecRestoreArrayRead(U, &u));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /*
56: Defines the ODE passed to the ODE solver
57: */
58: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
59: {
60: PetscScalar *f, Pmax;
61: const PetscScalar *u;
63: PetscFunctionBegin;
64: /* The next three lines allow us to access the entries of the vectors directly */
65: PetscCall(VecGetArrayRead(U, &u));
66: PetscCall(VecGetArray(F, &f));
67: 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 */
68: else Pmax = ctx->Pmax;
70: f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
71: f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);
73: PetscCall(VecRestoreArrayRead(U, &u));
74: PetscCall(VecRestoreArray(F, &f));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*
79: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
80: */
81: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
82: {
83: PetscInt rowcol[] = {0, 1};
84: PetscScalar J[2][2], Pmax;
85: const PetscScalar *u;
87: PetscFunctionBegin;
88: PetscCall(VecGetArrayRead(U, &u));
89: 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 */
90: else Pmax = ctx->Pmax;
92: J[0][0] = 0;
93: J[0][1] = ctx->omega_b;
94: J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
95: J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);
97: PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
98: PetscCall(VecRestoreArrayRead(U, &u));
100: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
101: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
102: if (A != B) {
103: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
104: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
110: {
111: PetscInt row[] = {0, 1}, col[] = {0};
112: PetscScalar J[2][1];
113: AppCtx *ctx = (AppCtx *)ctx0;
115: PetscFunctionBeginUser;
116: J[0][0] = 0;
117: J[1][0] = ctx->omega_s / (2.0 * ctx->H);
118: PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
119: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
120: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
125: {
126: PetscScalar *r;
127: const PetscScalar *u;
129: PetscFunctionBegin;
130: PetscCall(VecGetArrayRead(U, &u));
131: PetscCall(VecGetArray(R, &r));
132: r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
133: PetscCall(VecRestoreArray(R, &r));
134: PetscCall(VecRestoreArrayRead(U, &u));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
139: {
140: PetscScalar ru[1];
141: const PetscScalar *u;
142: PetscInt row[] = {0}, col[] = {0};
144: PetscFunctionBegin;
145: PetscCall(VecGetArrayRead(U, &u));
146: ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
147: PetscCall(VecRestoreArrayRead(U, &u));
148: PetscCall(MatSetValues(DRDU, 1, row, 1, col, ru, INSERT_VALUES));
149: PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
150: PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, AppCtx *ctx)
155: {
156: PetscFunctionBegin;
157: PetscCall(MatZeroEntries(DRDP));
158: PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
159: PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
164: {
165: PetscScalar sensip;
166: const PetscScalar *x, *y;
168: PetscFunctionBegin;
169: PetscCall(VecGetArrayRead(lambda, &x));
170: PetscCall(VecGetArrayRead(mu, &y));
171: sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
172: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt parameter pm: %.7f \n", (double)sensip));
173: PetscCall(VecRestoreArrayRead(lambda, &x));
174: PetscCall(VecRestoreArrayRead(mu, &y));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: int main(int argc, char **argv)
179: {
180: TS ts, quadts; /* ODE integrator */
181: Vec U; /* solution will be stored here */
182: Mat A; /* Jacobian matrix */
183: Mat Jacp; /* Jacobian matrix */
184: Mat DRDU, DRDP;
185: PetscMPIInt size;
186: PetscInt n = 2;
187: AppCtx ctx;
188: PetscScalar *u;
189: PetscReal du[2] = {0.0, 0.0};
190: PetscBool ensemble = PETSC_FALSE, flg1, flg2;
191: PetscReal ftime;
192: PetscInt steps;
193: PetscScalar *x_ptr, *y_ptr;
194: Vec lambda[1], q, mu[1];
196: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: Initialize program
198: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199: PetscFunctionBeginUser;
200: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
201: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
202: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Create necessary matrix and vectors
206: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
208: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
209: PetscCall(MatSetType(A, MATDENSE));
210: PetscCall(MatSetFromOptions(A));
211: PetscCall(MatSetUp(A));
213: PetscCall(MatCreateVecs(A, &U, NULL));
215: PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
216: PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
217: PetscCall(MatSetFromOptions(Jacp));
218: PetscCall(MatSetUp(Jacp));
220: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &DRDP));
221: PetscCall(MatSetUp(DRDP));
222: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, &DRDU));
223: PetscCall(MatSetUp(DRDU));
225: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226: Set runtime options
227: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
228: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
229: {
230: ctx.beta = 2;
231: ctx.c = 10000.0;
232: ctx.u_s = 1.0;
233: ctx.omega_s = 1.0;
234: ctx.omega_b = 120.0 * PETSC_PI;
235: ctx.H = 5.0;
236: PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
237: ctx.D = 5.0;
238: PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
239: ctx.E = 1.1378;
240: ctx.V = 1.0;
241: ctx.X = 0.545;
242: ctx.Pmax = ctx.E * ctx.V / ctx.X;
243: PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
244: ctx.Pm = 1.1;
245: PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
246: ctx.tf = 0.1;
247: ctx.tcl = 0.2;
248: PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
249: PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
250: PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
251: if (ensemble) {
252: ctx.tf = -1;
253: ctx.tcl = -1;
254: }
256: PetscCall(VecGetArray(U, &u));
257: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
258: u[1] = 1.0;
259: PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
260: n = 2;
261: PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
262: u[0] += du[0];
263: u[1] += du[1];
264: PetscCall(VecRestoreArray(U, &u));
265: if (flg1 || flg2) {
266: ctx.tf = -1;
267: ctx.tcl = -1;
268: }
269: }
270: PetscOptionsEnd();
272: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273: Create timestepping solver context
274: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
276: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
277: PetscCall(TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
278: PetscCall(TSSetType(ts, TSRK));
279: PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx));
280: PetscCall(TSSetRHSJacobian(ts, A, A, (TSRHSJacobian)RHSJacobian, &ctx));
281: PetscCall(TSCreateQuadratureTS(ts, PETSC_TRUE, &quadts));
282: PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, &ctx));
283: PetscCall(TSSetRHSJacobian(quadts, DRDU, DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &ctx));
284: PetscCall(TSSetRHSJacobianP(quadts, DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, &ctx));
285: PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
286: PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &ctx));
288: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: Set initial conditions
290: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291: PetscCall(TSSetSolution(ts, U));
293: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
294: Save trajectory of solution so that TSAdjointSolve() may be used
295: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
296: PetscCall(TSSetSaveTrajectory(ts));
298: PetscCall(MatCreateVecs(A, &lambda[0], NULL));
299: /* Set initial conditions for the adjoint integration */
300: PetscCall(VecGetArray(lambda[0], &y_ptr));
301: y_ptr[0] = 0.0;
302: y_ptr[1] = 0.0;
303: PetscCall(VecRestoreArray(lambda[0], &y_ptr));
305: PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
306: PetscCall(VecGetArray(mu[0], &x_ptr));
307: x_ptr[0] = -1.0;
308: PetscCall(VecRestoreArray(mu[0], &x_ptr));
310: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
311: Set solver options
312: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
313: PetscCall(TSSetMaxTime(ts, 10.0));
314: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
315: PetscCall(TSSetTimeStep(ts, .01));
316: PetscCall(TSSetFromOptions(ts));
318: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
319: Solve nonlinear system
320: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
321: if (ensemble) {
322: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
323: PetscCall(VecGetArray(U, &u));
324: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
325: u[1] = ctx.omega_s;
326: u[0] += du[0];
327: u[1] += du[1];
328: PetscCall(VecRestoreArray(U, &u));
329: PetscCall(TSSetTimeStep(ts, .01));
330: PetscCall(TSSolve(ts, U));
331: }
332: } else {
333: PetscCall(TSSolve(ts, U));
334: }
335: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
336: PetscCall(TSGetSolveTime(ts, &ftime));
337: PetscCall(TSGetStepNumber(ts, &steps));
339: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340: Adjoint model starts here
341: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
342: /* Set initial conditions for the adjoint integration */
343: PetscCall(VecGetArray(lambda[0], &y_ptr));
344: y_ptr[0] = 0.0;
345: y_ptr[1] = 0.0;
346: PetscCall(VecRestoreArray(lambda[0], &y_ptr));
348: PetscCall(VecGetArray(mu[0], &x_ptr));
349: x_ptr[0] = -1.0;
350: PetscCall(VecRestoreArray(mu[0], &x_ptr));
352: PetscCall(TSAdjointSolve(ts));
354: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n"));
355: PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
356: PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
357: PetscCall(TSGetCostIntegral(ts, &q));
358: PetscCall(VecView(q, PETSC_VIEWER_STDOUT_WORLD));
359: PetscCall(VecGetArray(q, &x_ptr));
360: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
361: PetscCall(VecRestoreArray(q, &x_ptr));
363: PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
365: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
366: Free work space. All PETSc objects should be destroyed when they are no longer needed.
367: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
368: PetscCall(MatDestroy(&A));
369: PetscCall(MatDestroy(&Jacp));
370: PetscCall(MatDestroy(&DRDU));
371: PetscCall(MatDestroy(&DRDP));
372: PetscCall(VecDestroy(&U));
373: PetscCall(VecDestroy(&lambda[0]));
374: PetscCall(VecDestroy(&mu[0]));
375: PetscCall(TSDestroy(&ts));
376: PetscCall(PetscFinalize());
377: return 0;
378: }
380: /*TEST
382: build:
383: requires: !complex
385: test:
386: args: -viewer_binary_skip_info -ts_adapt_type none
388: TEST*/