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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: PetscInitialize(&argc, &argv, (char *)0, help);
50: MPI_Comm_size(PETSC_COMM_WORLD, &size);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create necessary matrix and vectors
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: MatCreate(PETSC_COMM_WORLD, &ctx.Jac);
57: MatSetSizes(ctx.Jac, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
58: MatSetType(ctx.Jac, MATDENSE);
59: MatSetFromOptions(ctx.Jac);
60: MatSetUp(ctx.Jac);
61: MatCreateVecs(ctx.Jac, &U, NULL);
62: MatCreate(PETSC_COMM_WORLD, &ctx.Jacp);
63: MatSetSizes(ctx.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1);
64: MatSetFromOptions(ctx.Jacp);
65: MatSetUp(ctx.Jacp);
66: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &ctx.DRDP);
67: MatSetUp(ctx.DRDP);
68: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &ctx.DRDU);
69: 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: PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL);
83: ctx.D = 5.0;
84: 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: PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL);
91: ctx.Pm = 1.1;
92: PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL);
93: ctx.tf = 0.1;
94: ctx.tcl = 0.2;
95: PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL);
96: PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL);
97: PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL);
98: if (ensemble) {
99: ctx.tf = -1;
100: ctx.tcl = -1;
101: }
103: VecGetArray(U, &u);
104: u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
105: u[1] = 1.0;
106: PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1);
107: n = 2;
108: PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2);
109: u[0] += du[0];
110: u[1] += du[1];
111: VecRestoreArray(U, &u);
112: if (flg1 || flg2) {
113: ctx.tf = -1;
114: ctx.tcl = -1;
115: }
116: sa = SA_ADJ;
117: 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: TSCreate(PETSC_COMM_WORLD, &ts);
125: TSSetProblemType(ts, TS_NONLINEAR);
126: TSSetType(ts, TSBEULER);
127: TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx);
128: TSSetRHSJacobian(ts, ctx.Jac, ctx.Jac, (TSRHSJacobian)RHSJacobian, &ctx);
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Set initial conditions
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: TSSetSolution(ts, U);
135: /* Set RHS JacobianP */
136: TSSetRHSJacobianP(ts, ctx.Jacp, RHSJacobianP, &ctx);
138: TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts);
139: TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, &ctx);
140: TSSetRHSJacobian(quadts, ctx.DRDU, ctx.DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &ctx);
141: 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: TSSetSaveTrajectory(ts);
147: MatCreateVecs(ctx.Jac, &lambda[0], NULL);
148: MatCreateVecs(ctx.Jacp, &mu[0], NULL);
149: TSSetCostGradients(ts, 1, lambda, mu);
150: }
152: if (sa == SA_TLM) {
153: PetscScalar val[2];
154: PetscInt row[] = {0, 1}, col[] = {0};
156: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &qgrad);
157: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &sp);
158: TSForwardSetSensitivities(ts, 1, sp);
159: 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: MatSetValues(sp, 2, row, 1, col, val, INSERT_VALUES);
163: MatAssemblyBegin(sp, MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(sp, MAT_FINAL_ASSEMBLY);
165: }
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Set solver options
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: TSSetMaxTime(ts, 1.0);
171: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
172: TSSetTimeStep(ts, 0.03125);
173: TSSetFromOptions(ts);
175: direction[0] = direction[1] = 1;
176: terminate[0] = terminate[1] = PETSC_FALSE;
178: 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: 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: VecRestoreArray(U, &u);
191: TSSetTimeStep(ts, 0.03125);
192: TSSolve(ts, U);
193: }
194: } else {
195: TSSolve(ts, U);
196: }
197: TSGetSolveTime(ts, &ftime);
198: TSGetStepNumber(ts, &steps);
200: if (sa == SA_ADJ) {
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Adjoint model starts here
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: /* Set initial conditions for the adjoint integration */
205: VecGetArray(lambda[0], &y_ptr);
206: y_ptr[0] = 0.0;
207: y_ptr[1] = 0.0;
208: VecRestoreArray(lambda[0], &y_ptr);
210: VecGetArray(mu[0], &x_ptr);
211: x_ptr[0] = 0.0;
212: VecRestoreArray(mu[0], &x_ptr);
214: TSAdjointSolve(ts);
216: PetscPrintf(PETSC_COMM_WORLD, "\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n");
217: VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD);
218: PetscPrintf(PETSC_COMM_WORLD, "\n mu: d[Psi(tf)]/d[pm]\n");
219: VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD);
220: TSGetCostIntegral(ts, &q);
221: VecGetArray(q, &x_ptr);
222: PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm));
223: VecRestoreArray(q, &x_ptr);
224: ComputeSensiP(lambda[0], mu[0], &ctx);
225: VecGetArray(mu[0], &x_ptr);
226: PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)x_ptr[0]);
227: VecRestoreArray(mu[0], &x_ptr);
228: VecDestroy(&lambda[0]);
229: VecDestroy(&mu[0]);
230: }
231: if (sa == SA_TLM) {
232: PetscPrintf(PETSC_COMM_WORLD, "\n trajectory sensitivity: d[phi(tf)]/d[pm] d[omega(tf)]/d[pm]\n");
233: MatView(sp, PETSC_VIEWER_STDOUT_WORLD);
234: TSGetCostIntegral(ts, &q);
235: VecGetArray(q, &s_ptr);
236: PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(s_ptr[0] - ctx.Pm));
237: VecRestoreArray(q, &s_ptr);
238: MatDenseGetArray(qgrad, &s_ptr);
239: PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)s_ptr[0]);
240: MatDenseRestoreArray(qgrad, &s_ptr);
241: MatDestroy(&qgrad);
242: MatDestroy(&sp);
243: }
244: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: Free work space. All PETSc objects should be destroyed when they are no longer needed.
246: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
247: MatDestroy(&ctx.Jac);
248: MatDestroy(&ctx.Jacp);
249: MatDestroy(&ctx.DRDU);
250: MatDestroy(&ctx.DRDP);
251: VecDestroy(&U);
252: TSDestroy(&ts);
253: 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*/