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