1 2 static char help[] = "Adjoint and tangent linear sensitivity analysis of the basic equation for generator stability analysis.\n"; 3 4 /*F 5 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} 10 11 F*/ 12 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 */ 21 22 #include <petscts.h> 23 #include "ex3.h" 24 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; 44 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"); 52 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)); 70 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 } 102 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(); 120 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)); 129 130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 Set initial conditions 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 PetscCall(TSSetSolution(ts, U)); 134 135 /* Set RHS JacobianP */ 136 PetscCall(TSSetRHSJacobianP(ts, ctx.Jacp, RHSJacobianP, &ctx)); 137 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 } 151 152 if (sa == SA_TLM) { 153 PetscScalar val[2]; 154 PetscInt row[] = {0, 1}, col[] = {0}; 155 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 } 166 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)); 174 175 direction[0] = direction[1] = 1; 176 terminate[0] = terminate[1] = PETSC_FALSE; 177 178 PetscCall(TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&ctx)); 179 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)); 199 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)); 209 210 PetscCall(VecGetArray(mu[0], &x_ptr)); 211 x_ptr[0] = 0.0; 212 PetscCall(VecRestoreArray(mu[0], &x_ptr)); 213 214 PetscCall(TSAdjointSolve(ts)); 215 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 } 256 257 /*TEST 258 259 build: 260 requires: !complex !single 261 262 test: 263 args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu 264 265 test: 266 suffix: 2 267 args: -sa_method tlm -ts_type cn -pc_type lu 268 269 test: 270 suffix: 3 271 args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp 272 273 TEST*/ 274