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