Actual source code: ex9adj.c
1: static char help[] = "Basic equation for generator stability analysis.\n";
3: /*F
5: \begin{eqnarray}
6: \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
7: \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
8: \end{eqnarray}
10: Ensemble of initial conditions
11: ./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
13: Fault at .1 seconds
14: ./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
16: Initial conditions same as when fault is ended
17: ./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
19: F*/
21: /*
22: Include "petscts.h" so that we can use TS solvers. Note that this
23: file automatically includes:
24: petscsys.h - base PETSc routines petscvec.h - vectors
25: petscmat.h - matrices
26: petscis.h - index sets petscksp.h - Krylov subspace methods
27: petscviewer.h - viewers petscpc.h - preconditioners
28: petscksp.h - linear solvers
29: */
31: #include <petscts.h>
33: typedef struct {
34: PetscScalar H, D, omega_b, omega_s, Pmax, Pm, E, V, X, u_s, c;
35: PetscInt beta;
36: PetscReal tf, tcl;
37: } AppCtx;
39: PetscErrorCode PostStepFunction(TS ts)
40: {
41: Vec U;
42: PetscReal t;
43: const PetscScalar *u;
45: PetscFunctionBegin;
46: PetscCall(TSGetTime(ts, &t));
47: PetscCall(TSGetSolution(ts, &U));
48: PetscCall(VecGetArrayRead(U, &u));
49: PetscCall(PetscPrintf(PETSC_COMM_SELF, "delta(%3.2f) = %8.7f\n", (double)t, (double)u[0]));
50: PetscCall(VecRestoreArrayRead(U, &u));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*
55: Defines the ODE passed to the ODE solver
56: */
57: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
58: {
59: PetscScalar *f, Pmax;
60: const PetscScalar *u;
62: PetscFunctionBegin;
63: /* The next three lines allow us to access the entries of the vectors directly */
64: PetscCall(VecGetArrayRead(U, &u));
65: PetscCall(VecGetArray(F, &f));
66: 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 */
67: else Pmax = ctx->Pmax;
69: f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
70: f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);
72: PetscCall(VecRestoreArrayRead(U, &u));
73: PetscCall(VecRestoreArray(F, &f));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*
78: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
79: */
80: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
81: {
82: PetscInt rowcol[] = {0, 1};
83: PetscScalar J[2][2], Pmax;
84: const PetscScalar *u;
86: PetscFunctionBegin;
87: PetscCall(VecGetArrayRead(U, &u));
88: 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 */
89: else Pmax = ctx->Pmax;
91: J[0][0] = 0;
92: J[0][1] = ctx->omega_b;
93: J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
94: J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);
96: PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
97: PetscCall(VecRestoreArrayRead(U, &u));
99: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
100: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
101: if (A != B) {
102: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
103: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
104: }
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
109: {
110: PetscInt row[] = {0, 1}, col[] = {0};
111: PetscScalar J[2][1];
112: AppCtx *ctx = (AppCtx *)ctx0;
114: PetscFunctionBeginUser;
115: J[0][0] = 0;
116: J[1][0] = ctx->omega_s / (2.0 * ctx->H);
117: PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
118: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
119: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
124: {
125: PetscScalar *r;
126: const PetscScalar *u;
128: PetscFunctionBegin;
129: PetscCall(VecGetArrayRead(U, &u));
130: PetscCall(VecGetArray(R, &r));
131: r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
132: PetscCall(VecRestoreArray(R, &r));
133: PetscCall(VecRestoreArrayRead(U, &u));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
138: {
139: PetscScalar ru[1];
140: const PetscScalar *u;
141: PetscInt row[] = {0}, col[] = {0};
143: PetscFunctionBegin;
144: PetscCall(VecGetArrayRead(U, &u));
145: ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
146: PetscCall(VecRestoreArrayRead(U, &u));
147: PetscCall(MatSetValues(DRDU, 1, row, 1, col, ru, INSERT_VALUES));
148: PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
149: PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, AppCtx *ctx)
154: {
155: PetscFunctionBegin;
156: PetscCall(MatZeroEntries(DRDP));
157: PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
158: PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
163: {
164: PetscScalar sensip;
165: const PetscScalar *x, *y;
167: PetscFunctionBegin;
168: PetscCall(VecGetArrayRead(lambda, &x));
169: PetscCall(VecGetArrayRead(mu, &y));
170: sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
171: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt parameter pm: %.7f \n", (double)sensip));
172: PetscCall(VecRestoreArrayRead(lambda, &x));
173: PetscCall(VecRestoreArrayRead(mu, &y));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: int main(int argc, char **argv)
178: {
179: TS ts, quadts; /* ODE integrator */
180: Vec U; /* solution will be stored here */
181: Mat A; /* Jacobian matrix */
182: Mat Jacp; /* Jacobian matrix */
183: Mat DRDU, DRDP;
184: PetscMPIInt size;
185: PetscInt n = 2;
186: AppCtx ctx;
187: PetscScalar *u;
188: PetscReal du[2] = {0.0, 0.0};
189: PetscBool ensemble = PETSC_FALSE, flg1, flg2;
190: PetscReal ftime;
191: PetscInt steps;
192: PetscScalar *x_ptr, *y_ptr;
193: Vec lambda[1], q, mu[1];
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Initialize program
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: PetscFunctionBeginUser;
199: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
200: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
201: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Create necessary matrix and vectors
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
207: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
208: PetscCall(MatSetType(A, MATDENSE));
209: PetscCall(MatSetFromOptions(A));
210: PetscCall(MatSetUp(A));
212: PetscCall(MatCreateVecs(A, &U, NULL));
214: PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
215: PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
216: PetscCall(MatSetFromOptions(Jacp));
217: PetscCall(MatSetUp(Jacp));
219: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &DRDP));
220: PetscCall(MatSetUp(DRDP));
221: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, &DRDU));
222: PetscCall(MatSetUp(DRDU));
224: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225: Set runtime options
226: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
228: {
229: ctx.beta = 2;
230: ctx.c = 10000.0;
231: ctx.u_s = 1.0;
232: ctx.omega_s = 1.0;
233: ctx.omega_b = 120.0 * PETSC_PI;
234: ctx.H = 5.0;
235: PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
236: ctx.D = 5.0;
237: PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
238: ctx.E = 1.1378;
239: ctx.V = 1.0;
240: ctx.X = 0.545;
241: ctx.Pmax = ctx.E * ctx.V / ctx.X;
242: PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
243: ctx.Pm = 1.1;
244: PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
245: ctx.tf = 0.1;
246: ctx.tcl = 0.2;
247: PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
248: PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
249: PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
250: if (ensemble) {
251: ctx.tf = -1;
252: ctx.tcl = -1;
253: }
255: PetscCall(VecGetArray(U, &u));
256: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
257: u[1] = 1.0;
258: PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
259: n = 2;
260: PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
261: u[0] += du[0];
262: u[1] += du[1];
263: PetscCall(VecRestoreArray(U, &u));
264: if (flg1 || flg2) {
265: ctx.tf = -1;
266: ctx.tcl = -1;
267: }
268: }
269: PetscOptionsEnd();
271: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272: Create timestepping solver context
273: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
275: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
276: PetscCall(TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
277: PetscCall(TSSetType(ts, TSRK));
278: PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
279: PetscCall(TSSetRHSJacobian(ts, A, A, (TSRHSJacobianFn *)RHSJacobian, &ctx));
280: PetscCall(TSCreateQuadratureTS(ts, PETSC_TRUE, &quadts));
281: PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &ctx));
282: PetscCall(TSSetRHSJacobian(quadts, DRDU, DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, &ctx));
283: PetscCall(TSSetRHSJacobianP(quadts, DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, &ctx));
284: PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
285: PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &ctx));
287: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288: Set initial conditions
289: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290: PetscCall(TSSetSolution(ts, U));
292: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
293: Save trajectory of solution so that TSAdjointSolve() may be used
294: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
295: PetscCall(TSSetSaveTrajectory(ts));
297: PetscCall(MatCreateVecs(A, &lambda[0], NULL));
298: /* Set initial conditions for the adjoint integration */
299: PetscCall(VecGetArray(lambda[0], &y_ptr));
300: y_ptr[0] = 0.0;
301: y_ptr[1] = 0.0;
302: PetscCall(VecRestoreArray(lambda[0], &y_ptr));
304: PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
305: PetscCall(VecGetArray(mu[0], &x_ptr));
306: x_ptr[0] = -1.0;
307: PetscCall(VecRestoreArray(mu[0], &x_ptr));
309: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310: Set solver options
311: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312: PetscCall(TSSetMaxTime(ts, 10.0));
313: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
314: PetscCall(TSSetTimeStep(ts, .01));
315: PetscCall(TSSetFromOptions(ts));
317: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318: Solve nonlinear system
319: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320: if (ensemble) {
321: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
322: PetscCall(VecGetArray(U, &u));
323: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
324: u[1] = ctx.omega_s;
325: u[0] += du[0];
326: u[1] += du[1];
327: PetscCall(VecRestoreArray(U, &u));
328: PetscCall(TSSetTimeStep(ts, .01));
329: PetscCall(TSSolve(ts, U));
330: }
331: } else {
332: PetscCall(TSSolve(ts, U));
333: }
334: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
335: PetscCall(TSGetSolveTime(ts, &ftime));
336: PetscCall(TSGetStepNumber(ts, &steps));
338: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339: Adjoint model starts here
340: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
341: /* Set initial conditions for the adjoint integration */
342: PetscCall(VecGetArray(lambda[0], &y_ptr));
343: y_ptr[0] = 0.0;
344: y_ptr[1] = 0.0;
345: PetscCall(VecRestoreArray(lambda[0], &y_ptr));
347: PetscCall(VecGetArray(mu[0], &x_ptr));
348: x_ptr[0] = -1.0;
349: PetscCall(VecRestoreArray(mu[0], &x_ptr));
351: PetscCall(TSAdjointSolve(ts));
353: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n"));
354: PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
355: PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
356: PetscCall(TSGetCostIntegral(ts, &q));
357: PetscCall(VecView(q, PETSC_VIEWER_STDOUT_WORLD));
358: PetscCall(VecGetArray(q, &x_ptr));
359: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
360: PetscCall(VecRestoreArray(q, &x_ptr));
362: PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
364: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
365: Free work space. All PETSc objects should be destroyed when they are no longer needed.
366: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
367: PetscCall(MatDestroy(&A));
368: PetscCall(MatDestroy(&Jacp));
369: PetscCall(MatDestroy(&DRDU));
370: PetscCall(MatDestroy(&DRDP));
371: PetscCall(VecDestroy(&U));
372: PetscCall(VecDestroy(&lambda[0]));
373: PetscCall(VecDestroy(&mu[0]));
374: PetscCall(TSDestroy(&ts));
375: PetscCall(PetscFinalize());
376: return 0;
377: }
379: /*TEST
381: build:
382: requires: !complex
384: test:
385: args: -viewer_binary_skip_info -ts_adapt_type none
387: TEST*/