Actual source code: ex3sa.c
2: static char help[] = "Adjoint and tangent linear sensitivity analysis of the 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: F*/
13: /*
14: This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities.
15: It computes the sensitivities of an integral cost function
16: \int c*max(0,\theta(t)-u_s)^beta dt
17: w.r.t. initial conditions and the parameter P_m.
18: Backward Euler method is used for time integration.
19: The discontinuities are detected with TSEvent.
20: */
22: #include <petscts.h>
23: #include "ex3.h"
25: int main(int argc, char **argv)
26: {
27: TS ts, quadts; /* ODE integrator */
28: Vec U; /* solution will be stored here */
29: PetscMPIInt size;
30: PetscInt n = 2;
31: AppCtx ctx;
32: PetscScalar *u;
33: PetscReal du[2] = {0.0, 0.0};
34: PetscBool ensemble = PETSC_FALSE, flg1, flg2;
35: PetscReal ftime;
36: PetscInt steps;
37: PetscScalar *x_ptr, *y_ptr, *s_ptr;
38: Vec lambda[1], q, mu[1];
39: PetscInt direction[2];
40: PetscBool terminate[2];
41: Mat qgrad;
42: Mat sp; /* Forward sensitivity matrix */
43: SAMethod sa;
45: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: Initialize program
47: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
48: PetscFunctionBeginUser;
49: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
50: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
51: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create necessary matrix and vectors
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jac));
57: PetscCall(MatSetSizes(ctx.Jac, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
58: PetscCall(MatSetType(ctx.Jac, MATDENSE));
59: PetscCall(MatSetFromOptions(ctx.Jac));
60: PetscCall(MatSetUp(ctx.Jac));
61: PetscCall(MatCreateVecs(ctx.Jac, &U, NULL));
62: PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jacp));
63: PetscCall(MatSetSizes(ctx.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
64: PetscCall(MatSetFromOptions(ctx.Jacp));
65: PetscCall(MatSetUp(ctx.Jacp));
66: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &ctx.DRDP));
67: PetscCall(MatSetUp(ctx.DRDP));
68: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &ctx.DRDU));
69: PetscCall(MatSetUp(ctx.DRDU));
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Set runtime options
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
75: {
76: ctx.beta = 2;
77: ctx.c = 10000.0;
78: ctx.u_s = 1.0;
79: ctx.omega_s = 1.0;
80: ctx.omega_b = 120.0 * PETSC_PI;
81: ctx.H = 5.0;
82: PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
83: ctx.D = 5.0;
84: PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
85: ctx.E = 1.1378;
86: ctx.V = 1.0;
87: ctx.X = 0.545;
88: ctx.Pmax = ctx.E * ctx.V / ctx.X;
89: ctx.Pmax_ini = ctx.Pmax;
90: PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
91: ctx.Pm = 1.1;
92: PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
93: ctx.tf = 0.1;
94: ctx.tcl = 0.2;
95: PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
96: PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
97: PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
98: if (ensemble) {
99: ctx.tf = -1;
100: ctx.tcl = -1;
101: }
103: PetscCall(VecGetArray(U, &u));
104: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
105: u[1] = 1.0;
106: PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
107: n = 2;
108: PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
109: u[0] += du[0];
110: u[1] += du[1];
111: PetscCall(VecRestoreArray(U, &u));
112: if (flg1 || flg2) {
113: ctx.tf = -1;
114: ctx.tcl = -1;
115: }
116: sa = SA_ADJ;
117: PetscCall(PetscOptionsEnum("-sa_method", "Sensitivity analysis method (adj or tlm)", "", SAMethods, (PetscEnum)sa, (PetscEnum *)&sa, NULL));
118: }
119: PetscOptionsEnd();
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create timestepping solver context
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
125: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
126: PetscCall(TSSetType(ts, TSBEULER));
127: PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx));
128: PetscCall(TSSetRHSJacobian(ts, ctx.Jac, ctx.Jac, (TSRHSJacobian)RHSJacobian, &ctx));
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Set initial conditions
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: PetscCall(TSSetSolution(ts, U));
135: /* Set RHS JacobianP */
136: PetscCall(TSSetRHSJacobianP(ts, ctx.Jacp, RHSJacobianP, &ctx));
138: PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
139: PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, &ctx));
140: PetscCall(TSSetRHSJacobian(quadts, ctx.DRDU, ctx.DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &ctx));
141: PetscCall(TSSetRHSJacobianP(quadts, ctx.DRDP, DRDPJacobianTranspose, &ctx));
142: if (sa == SA_ADJ) {
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Save trajectory of solution so that TSAdjointSolve() may be used
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: PetscCall(TSSetSaveTrajectory(ts));
147: PetscCall(MatCreateVecs(ctx.Jac, &lambda[0], NULL));
148: PetscCall(MatCreateVecs(ctx.Jacp, &mu[0], NULL));
149: PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
150: }
152: if (sa == SA_TLM) {
153: PetscScalar val[2];
154: PetscInt row[] = {0, 1}, col[] = {0};
156: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &qgrad));
157: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &sp));
158: PetscCall(TSForwardSetSensitivities(ts, 1, sp));
159: PetscCall(TSForwardSetSensitivities(quadts, 1, qgrad));
160: val[0] = 1. / PetscSqrtScalar(1. - (ctx.Pm / ctx.Pmax) * (ctx.Pm / ctx.Pmax)) / ctx.Pmax;
161: val[1] = 0.0;
162: PetscCall(MatSetValues(sp, 2, row, 1, col, val, INSERT_VALUES));
163: PetscCall(MatAssemblyBegin(sp, MAT_FINAL_ASSEMBLY));
164: PetscCall(MatAssemblyEnd(sp, MAT_FINAL_ASSEMBLY));
165: }
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Set solver options
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(TSSetMaxTime(ts, 1.0));
171: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
172: PetscCall(TSSetTimeStep(ts, 0.03125));
173: PetscCall(TSSetFromOptions(ts));
175: direction[0] = direction[1] = 1;
176: terminate[0] = terminate[1] = PETSC_FALSE;
178: PetscCall(TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&ctx));
180: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: Solve nonlinear system
182: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: if (ensemble) {
184: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
185: PetscCall(VecGetArray(U, &u));
186: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
187: u[1] = ctx.omega_s;
188: u[0] += du[0];
189: u[1] += du[1];
190: PetscCall(VecRestoreArray(U, &u));
191: PetscCall(TSSetTimeStep(ts, 0.03125));
192: PetscCall(TSSolve(ts, U));
193: }
194: } else {
195: PetscCall(TSSolve(ts, U));
196: }
197: PetscCall(TSGetSolveTime(ts, &ftime));
198: PetscCall(TSGetStepNumber(ts, &steps));
200: if (sa == SA_ADJ) {
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Adjoint model starts here
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: /* Set initial conditions for the adjoint integration */
205: PetscCall(VecGetArray(lambda[0], &y_ptr));
206: y_ptr[0] = 0.0;
207: y_ptr[1] = 0.0;
208: PetscCall(VecRestoreArray(lambda[0], &y_ptr));
210: PetscCall(VecGetArray(mu[0], &x_ptr));
211: x_ptr[0] = 0.0;
212: PetscCall(VecRestoreArray(mu[0], &x_ptr));
214: PetscCall(TSAdjointSolve(ts));
216: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n"));
217: PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
218: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n mu: d[Psi(tf)]/d[pm]\n"));
219: PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
220: PetscCall(TSGetCostIntegral(ts, &q));
221: PetscCall(VecGetArray(q, &x_ptr));
222: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
223: PetscCall(VecRestoreArray(q, &x_ptr));
224: PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
225: PetscCall(VecGetArray(mu[0], &x_ptr));
226: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)x_ptr[0]));
227: PetscCall(VecRestoreArray(mu[0], &x_ptr));
228: PetscCall(VecDestroy(&lambda[0]));
229: PetscCall(VecDestroy(&mu[0]));
230: }
231: if (sa == SA_TLM) {
232: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n trajectory sensitivity: d[phi(tf)]/d[pm] d[omega(tf)]/d[pm]\n"));
233: PetscCall(MatView(sp, PETSC_VIEWER_STDOUT_WORLD));
234: PetscCall(TSGetCostIntegral(ts, &q));
235: PetscCall(VecGetArray(q, &s_ptr));
236: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(s_ptr[0] - ctx.Pm)));
237: PetscCall(VecRestoreArray(q, &s_ptr));
238: PetscCall(MatDenseGetArray(qgrad, &s_ptr));
239: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)s_ptr[0]));
240: PetscCall(MatDenseRestoreArray(qgrad, &s_ptr));
241: PetscCall(MatDestroy(&qgrad));
242: PetscCall(MatDestroy(&sp));
243: }
244: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: Free work space. All PETSc objects should be destroyed when they are no longer needed.
246: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
247: PetscCall(MatDestroy(&ctx.Jac));
248: PetscCall(MatDestroy(&ctx.Jacp));
249: PetscCall(MatDestroy(&ctx.DRDU));
250: PetscCall(MatDestroy(&ctx.DRDP));
251: PetscCall(VecDestroy(&U));
252: PetscCall(TSDestroy(&ts));
253: PetscCall(PetscFinalize());
254: return 0;
255: }
257: /*TEST
259: build:
260: requires: !complex !single
262: test:
263: args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu
265: test:
266: suffix: 2
267: args: -sa_method tlm -ts_type cn -pc_type lu
269: test:
270: suffix: 3
271: args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp
273: TEST*/