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