1 /* function subroutines used by power.c */ 2 3 #include "power.h" 4 5 PetscErrorCode GetListofEdges_Power(PFDATA *pfdata, PetscInt *edgelist) 6 { 7 PetscInt i, fbus, tbus, nbranches = pfdata->nbranch; 8 EDGE_Power branch = pfdata->branch; 9 PetscBool netview = PETSC_FALSE; 10 11 PetscFunctionBegin; 12 PetscCall(PetscOptionsHasName(NULL, NULL, "-powernet_view", &netview)); 13 for (i = 0; i < nbranches; i++) { 14 fbus = branch[i].internal_i; 15 tbus = branch[i].internal_j; 16 edgelist[2 * i] = fbus; 17 edgelist[2 * i + 1] = tbus; 18 if (netview) PetscCall(PetscPrintf(PETSC_COMM_SELF, "branch %" PetscInt_FMT ", bus[%" PetscInt_FMT "] -> bus[%" PetscInt_FMT "]\n", i, fbus, tbus)); 19 } 20 if (netview) { 21 for (i = 0; i < pfdata->nbus; i++) { 22 if (pfdata->bus[i].ngen) { 23 PetscCall(PetscPrintf(PETSC_COMM_SELF, " bus %" PetscInt_FMT ": gen\n", i)); 24 } else if (pfdata->bus[i].nload) { 25 PetscCall(PetscPrintf(PETSC_COMM_SELF, " bus %" PetscInt_FMT ": load\n", i)); 26 } 27 } 28 } 29 PetscFunctionReturn(PETSC_SUCCESS); 30 } 31 32 PetscErrorCode FormJacobian_Power_private(DM networkdm, Vec localX, Mat J, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 33 { 34 const PetscScalar *xarr; 35 PetscInt i, v, row[2], col[8], e, vfrom, vto; 36 PetscInt offsetfrom, offsetto, goffsetfrom, goffsetto, numComps; 37 PetscScalar values[8]; 38 PetscInt j, key, offset, goffset; 39 PetscScalar Vm; 40 UserCtx_Power *user_power = (UserCtx_Power *)appctx; 41 PetscScalar Sbase = user_power->Sbase; 42 VERTEX_Power bus; 43 PetscBool ghostvtex; 44 void *component; 45 46 PetscFunctionBegin; 47 PetscCall(VecGetArrayRead(localX, &xarr)); 48 49 for (v = 0; v < nv; v++) { 50 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex)); 51 52 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps)); 53 for (j = 0; j < numComps; j++) { 54 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset)); 55 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &goffset)); 56 PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL)); 57 58 if (key == user_power->compkey_bus) { 59 PetscInt nconnedges; 60 const PetscInt *connedges; 61 62 bus = (VERTEX_Power)component; 63 if (!ghostvtex) { 64 /* Handle reference bus constrained dofs */ 65 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 66 row[0] = goffset; 67 row[1] = goffset + 1; 68 col[0] = goffset; 69 col[1] = goffset + 1; 70 col[2] = goffset; 71 col[3] = goffset + 1; 72 values[0] = 1.0; 73 values[1] = 0.0; 74 values[2] = 0.0; 75 values[3] = 1.0; 76 PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES)); 77 break; 78 } 79 80 Vm = xarr[offset + 1]; 81 82 /* Shunt injections */ 83 row[0] = goffset; 84 row[1] = goffset + 1; 85 col[0] = goffset; 86 col[1] = goffset + 1; 87 values[0] = values[1] = values[2] = values[3] = 0.0; 88 if (bus->ide != PV_BUS) { 89 values[1] = 2.0 * Vm * bus->gl / Sbase; 90 values[3] = -2.0 * Vm * bus->bl / Sbase; 91 } 92 PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES)); 93 } 94 95 PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges)); 96 97 for (i = 0; i < nconnedges; i++) { 98 EDGE_Power branch; 99 VERTEX_Power busf, bust; 100 PetscInt keyf, keyt; 101 PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; 102 const PetscInt *cone; 103 PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf; 104 105 e = connedges[i]; 106 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL)); 107 if (!branch->status) continue; 108 109 Gff = branch->yff[0]; 110 Bff = branch->yff[1]; 111 Gft = branch->yft[0]; 112 Bft = branch->yft[1]; 113 Gtf = branch->ytf[0]; 114 Btf = branch->ytf[1]; 115 Gtt = branch->ytt[0]; 116 Btt = branch->ytt[1]; 117 118 PetscCall(DMNetworkGetConnectedVertices(networkdm, edges[e], &cone)); 119 vfrom = cone[0]; 120 vto = cone[1]; 121 122 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom)); 123 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto)); 124 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &goffsetfrom)); 125 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vto, ALL_COMPONENTS, &goffsetto)); 126 127 if (goffsetto < 0) goffsetto = -goffsetto - 1; 128 129 thetaf = xarr[offsetfrom]; 130 Vmf = xarr[offsetfrom + 1]; 131 thetat = xarr[offsetto]; 132 Vmt = xarr[offsetto + 1]; 133 thetaft = thetaf - thetat; 134 thetatf = thetat - thetaf; 135 136 PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL)); 137 PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL)); 138 139 if (vfrom == vtx[v]) { 140 if (busf->ide != REF_BUS) { 141 /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */ 142 row[0] = goffsetfrom; 143 col[0] = goffsetfrom; 144 col[1] = goffsetfrom + 1; 145 col[2] = goffsetto; 146 col[3] = goffsetto + 1; 147 values[0] = Vmf * Vmt * (Gft * -PetscSinScalar(thetaft) + Bft * PetscCosScalar(thetaft)); /* df_dthetaf */ 148 values[1] = 2.0 * Gff * Vmf + Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmf */ 149 values[2] = Vmf * Vmt * (Gft * PetscSinScalar(thetaft) + Bft * -PetscCosScalar(thetaft)); /* df_dthetat */ 150 values[3] = Vmf * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmt */ 151 152 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 153 } 154 if (busf->ide != PV_BUS && busf->ide != REF_BUS) { 155 row[0] = goffsetfrom + 1; 156 col[0] = goffsetfrom; 157 col[1] = goffsetfrom + 1; 158 col[2] = goffsetto; 159 col[3] = goffsetto + 1; 160 /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */ 161 values[0] = Vmf * Vmt * (Bft * PetscSinScalar(thetaft) + Gft * PetscCosScalar(thetaft)); 162 values[1] = -2.0 * Bff * Vmf + Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 163 values[2] = Vmf * Vmt * (-Bft * PetscSinScalar(thetaft) + Gft * -PetscCosScalar(thetaft)); 164 values[3] = Vmf * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 165 166 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 167 } 168 } else { 169 if (bust->ide != REF_BUS) { 170 row[0] = goffsetto; 171 col[0] = goffsetto; 172 col[1] = goffsetto + 1; 173 col[2] = goffsetfrom; 174 col[3] = goffsetfrom + 1; 175 /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */ 176 values[0] = Vmt * Vmf * (Gtf * -PetscSinScalar(thetatf) + Btf * PetscCosScalar(thetaft)); /* df_dthetat */ 177 values[1] = 2.0 * Gtt * Vmt + Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmt */ 178 values[2] = Vmt * Vmf * (Gtf * PetscSinScalar(thetatf) + Btf * -PetscCosScalar(thetatf)); /* df_dthetaf */ 179 values[3] = Vmt * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmf */ 180 181 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 182 } 183 if (bust->ide != PV_BUS && bust->ide != REF_BUS) { 184 row[0] = goffsetto + 1; 185 col[0] = goffsetto; 186 col[1] = goffsetto + 1; 187 col[2] = goffsetfrom; 188 col[3] = goffsetfrom + 1; 189 /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */ 190 values[0] = Vmt * Vmf * (Btf * PetscSinScalar(thetatf) + Gtf * PetscCosScalar(thetatf)); 191 values[1] = -2.0 * Btt * Vmt + Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 192 values[2] = Vmt * Vmf * (-Btf * PetscSinScalar(thetatf) + Gtf * -PetscCosScalar(thetatf)); 193 values[3] = Vmt * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 194 195 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 196 } 197 } 198 } 199 if (!ghostvtex && bus->ide == PV_BUS) { 200 row[0] = goffset + 1; 201 col[0] = goffset + 1; 202 values[0] = 1.0; 203 if (user_power->jac_error) values[0] = 50.0; 204 PetscCall(MatSetValues(J, 1, row, 1, col, values, ADD_VALUES)); 205 } 206 } 207 } 208 } 209 210 PetscCall(VecRestoreArrayRead(localX, &xarr)); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 PetscErrorCode FormJacobian_Power(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx) 215 { 216 DM networkdm; 217 Vec localX; 218 PetscInt nv, ne; 219 const PetscInt *vtx, *edges; 220 221 PetscFunctionBegin; 222 PetscCall(MatZeroEntries(J)); 223 224 PetscCall(SNESGetDM(snes, &networkdm)); 225 PetscCall(DMGetLocalVector(networkdm, &localX)); 226 227 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 228 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 229 230 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 231 PetscCall(FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx)); 232 233 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 234 235 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 236 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 PetscErrorCode FormFunction_Power(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 241 { 242 UserCtx_Power *User = (UserCtx_Power *)appctx; 243 PetscInt e, v, vfrom, vto; 244 const PetscScalar *xarr; 245 PetscScalar *farr; 246 PetscInt offsetfrom, offsetto, offset, i, j, key, numComps; 247 PetscScalar Vm; 248 PetscScalar Sbase = User->Sbase; 249 VERTEX_Power bus = NULL; 250 GEN gen; 251 LOAD load; 252 PetscBool ghostvtex; 253 void *component; 254 255 PetscFunctionBegin; 256 PetscCall(VecGetArrayRead(localX, &xarr)); 257 PetscCall(VecGetArray(localF, &farr)); 258 259 for (v = 0; v < nv; v++) { 260 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex)); 261 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps)); 262 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset)); 263 264 for (j = 0; j < numComps; j++) { 265 PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL)); 266 if (key == User->compkey_bus) { 267 PetscInt nconnedges; 268 const PetscInt *connedges; 269 270 bus = (VERTEX_Power)component; 271 /* Handle reference bus constrained dofs */ 272 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 273 farr[offset] = xarr[offset] - bus->va * PETSC_PI / 180.0; 274 farr[offset + 1] = xarr[offset + 1] - bus->vm; 275 break; 276 } 277 278 if (!ghostvtex) { 279 Vm = xarr[offset + 1]; 280 281 /* Shunt injections */ 282 farr[offset] += Vm * Vm * bus->gl / Sbase; 283 if (bus->ide != PV_BUS) farr[offset + 1] += -Vm * Vm * bus->bl / Sbase; 284 } 285 286 PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges)); 287 for (i = 0; i < nconnedges; i++) { 288 EDGE_Power branch; 289 PetscInt keye; 290 PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; 291 const PetscInt *cone; 292 PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf; 293 294 e = connedges[i]; 295 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL)); 296 if (!branch->status) continue; 297 Gff = branch->yff[0]; 298 Bff = branch->yff[1]; 299 Gft = branch->yft[0]; 300 Bft = branch->yft[1]; 301 Gtf = branch->ytf[0]; 302 Btf = branch->ytf[1]; 303 Gtt = branch->ytt[0]; 304 Btt = branch->ytt[1]; 305 306 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone)); 307 vfrom = cone[0]; 308 vto = cone[1]; 309 310 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom)); 311 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto)); 312 313 thetaf = xarr[offsetfrom]; 314 Vmf = xarr[offsetfrom + 1]; 315 thetat = xarr[offsetto]; 316 Vmt = xarr[offsetto + 1]; 317 thetaft = thetaf - thetat; 318 thetatf = thetat - thetaf; 319 320 if (vfrom == vtx[v]) { 321 farr[offsetfrom] += Gff * Vmf * Vmf + Vmf * Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); 322 farr[offsetfrom + 1] += -Bff * Vmf * Vmf + Vmf * Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 323 } else { 324 farr[offsetto] += Gtt * Vmt * Vmt + Vmt * Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); 325 farr[offsetto + 1] += -Btt * Vmt * Vmt + Vmt * Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 326 } 327 } 328 } else if (key == User->compkey_gen) { 329 if (!ghostvtex) { 330 gen = (GEN)component; 331 if (!gen->status) continue; 332 farr[offset] += -gen->pg / Sbase; 333 farr[offset + 1] += -gen->qg / Sbase; 334 } 335 } else if (key == User->compkey_load) { 336 if (!ghostvtex) { 337 load = (LOAD)component; 338 farr[offset] += load->pl / Sbase; 339 farr[offset + 1] += load->ql / Sbase; 340 } 341 } 342 } 343 if (bus && bus->ide == PV_BUS) farr[offset + 1] = xarr[offset + 1] - bus->vm; 344 } 345 PetscCall(VecRestoreArrayRead(localX, &xarr)); 346 PetscCall(VecRestoreArray(localF, &farr)); 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349 350 PetscErrorCode SetInitialGuess_Power(DM networkdm, Vec localX, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 351 { 352 VERTEX_Power bus; 353 PetscInt i; 354 GEN gen; 355 PetscBool ghostvtex, sharedv; 356 PetscScalar *xarr; 357 PetscInt key, numComps, j, offset; 358 void *component; 359 PetscMPIInt rank; 360 MPI_Comm comm; 361 UserCtx_Power *User = (UserCtx_Power *)appctx; 362 363 PetscFunctionBegin; 364 PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm)); 365 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 366 PetscCall(VecGetArray(localX, &xarr)); 367 for (i = 0; i < nv; i++) { 368 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex)); 369 PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &sharedv)); 370 if (ghostvtex || sharedv) continue; 371 372 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset)); 373 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &numComps)); 374 for (j = 0; j < numComps; j++) { 375 PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, &component, NULL)); 376 if (key == User->compkey_bus) { 377 bus = (VERTEX_Power)component; 378 xarr[offset] = bus->va * PETSC_PI / 180.0; 379 xarr[offset + 1] = bus->vm; 380 } else if (key == User->compkey_gen) { 381 gen = (GEN)component; 382 if (!gen->status) continue; 383 xarr[offset + 1] = gen->vs; 384 break; 385 } 386 } 387 } 388 PetscCall(VecRestoreArray(localX, &xarr)); 389 PetscFunctionReturn(PETSC_SUCCESS); 390 } 391