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