1 static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks for solving a nonlinear electric power grid problem.\n\ 2 The available solver options are in the poweroptions file and the data files are in the datafiles directory.\n\ 3 The data file format used is from the MatPower package (http://www.pserc.cornell.edu//matpower/).\n\ 4 This example shows the use of subnetwork feature in DMNetwork. It creates duplicates of the same network which are treated as subnetworks.\n\ 5 Run this program: mpiexec -n <n> ./pf2\n\ 6 mpiexec -n <n> ./pf2 \n"; 7 8 #include "power.h" 9 #include <petscdmnetwork.h> 10 11 PetscErrorCode FormFunction_Subnet(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 12 { 13 UserCtx_Power *User = (UserCtx_Power *)appctx; 14 PetscInt e, v, vfrom, vto; 15 const PetscScalar *xarr; 16 PetscScalar *farr; 17 PetscInt offsetfrom, offsetto, offset; 18 19 PetscFunctionBegin; 20 PetscCall(VecGetArrayRead(localX, &xarr)); 21 PetscCall(VecGetArray(localF, &farr)); 22 23 for (v = 0; v < nv; v++) { 24 PetscInt i, j, key; 25 PetscScalar Vm; 26 PetscScalar Sbase = User->Sbase; 27 VERTEX_Power bus = NULL; 28 GEN gen; 29 LOAD load; 30 PetscBool ghostvtex; 31 PetscInt numComps; 32 void *component; 33 34 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex)); 35 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps)); 36 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset)); 37 for (j = 0; j < numComps; j++) { 38 PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL)); 39 if (key == 1) { 40 PetscInt nconnedges; 41 const PetscInt *connedges; 42 43 bus = (VERTEX_Power)(component); 44 /* Handle reference bus constrained dofs */ 45 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 46 farr[offset] = xarr[offset] - bus->va * PETSC_PI / 180.0; 47 farr[offset + 1] = xarr[offset + 1] - bus->vm; 48 break; 49 } 50 51 if (!ghostvtex) { 52 Vm = xarr[offset + 1]; 53 54 /* Shunt injections */ 55 farr[offset] += Vm * Vm * bus->gl / Sbase; 56 if (bus->ide != PV_BUS) farr[offset + 1] += -Vm * Vm * bus->bl / Sbase; 57 } 58 59 PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges)); 60 for (i = 0; i < nconnedges; i++) { 61 EDGE_Power branch; 62 PetscInt keye; 63 PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; 64 const PetscInt *cone; 65 PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf; 66 67 e = connedges[i]; 68 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL)); 69 if (!branch->status) continue; 70 Gff = branch->yff[0]; 71 Bff = branch->yff[1]; 72 Gft = branch->yft[0]; 73 Bft = branch->yft[1]; 74 Gtf = branch->ytf[0]; 75 Btf = branch->ytf[1]; 76 Gtt = branch->ytt[0]; 77 Btt = branch->ytt[1]; 78 79 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone)); 80 vfrom = cone[0]; 81 vto = cone[1]; 82 83 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom)); 84 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto)); 85 86 thetaf = xarr[offsetfrom]; 87 Vmf = xarr[offsetfrom + 1]; 88 thetat = xarr[offsetto]; 89 Vmt = xarr[offsetto + 1]; 90 thetaft = thetaf - thetat; 91 thetatf = thetat - thetaf; 92 93 if (vfrom == vtx[v]) { 94 farr[offsetfrom] += Gff * Vmf * Vmf + Vmf * Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); 95 farr[offsetfrom + 1] += -Bff * Vmf * Vmf + Vmf * Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 96 } else { 97 farr[offsetto] += Gtt * Vmt * Vmt + Vmt * Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); 98 farr[offsetto + 1] += -Btt * Vmt * Vmt + Vmt * Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 99 } 100 } 101 } else if (key == 2) { 102 if (!ghostvtex) { 103 gen = (GEN)(component); 104 if (!gen->status) continue; 105 farr[offset] += -gen->pg / Sbase; 106 farr[offset + 1] += -gen->qg / Sbase; 107 } 108 } else if (key == 3) { 109 if (!ghostvtex) { 110 load = (LOAD)(component); 111 farr[offset] += load->pl / Sbase; 112 farr[offset + 1] += load->ql / Sbase; 113 } 114 } 115 } 116 if (bus && bus->ide == PV_BUS) farr[offset + 1] = xarr[offset + 1] - bus->vm; 117 } 118 PetscCall(VecRestoreArrayRead(localX, &xarr)); 119 PetscCall(VecRestoreArray(localF, &farr)); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx) 124 { 125 DM networkdm; 126 Vec localX, localF; 127 PetscInt nv, ne; 128 const PetscInt *vtx, *edges; 129 130 PetscFunctionBegin; 131 PetscCall(SNESGetDM(snes, &networkdm)); 132 PetscCall(DMGetLocalVector(networkdm, &localX)); 133 PetscCall(DMGetLocalVector(networkdm, &localF)); 134 PetscCall(VecSet(F, 0.0)); 135 136 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 137 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 138 139 PetscCall(DMGlobalToLocalBegin(networkdm, F, INSERT_VALUES, localF)); 140 PetscCall(DMGlobalToLocalEnd(networkdm, F, INSERT_VALUES, localF)); 141 142 /* Form Function for first subnetwork */ 143 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 144 PetscCall(FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx)); 145 146 /* Form Function for second subnetwork */ 147 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges)); 148 PetscCall(FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx)); 149 150 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 151 152 PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F)); 153 PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F)); 154 PetscCall(DMRestoreLocalVector(networkdm, &localF)); 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 PetscErrorCode FormJacobian_Subnet(DM networkdm, Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 159 { 160 UserCtx_Power *User = (UserCtx_Power *)appctx; 161 PetscInt e, v, vfrom, vto; 162 const PetscScalar *xarr; 163 PetscInt offsetfrom, offsetto, goffsetfrom, goffsetto; 164 PetscInt row[2], col[8]; 165 PetscScalar values[8]; 166 167 PetscFunctionBegin; 168 PetscCall(VecGetArrayRead(localX, &xarr)); 169 170 for (v = 0; v < nv; v++) { 171 PetscInt i, j, key; 172 PetscInt offset, goffset; 173 PetscScalar Vm; 174 PetscScalar Sbase = User->Sbase; 175 VERTEX_Power bus; 176 PetscBool ghostvtex; 177 PetscInt numComps; 178 void *component; 179 180 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex)); 181 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps)); 182 for (j = 0; j < numComps; j++) { 183 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset)); 184 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &goffset)); 185 PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL)); 186 if (key == 1) { 187 PetscInt nconnedges; 188 const PetscInt *connedges; 189 190 bus = (VERTEX_Power)(component); 191 if (!ghostvtex) { 192 /* Handle reference bus constrained dofs */ 193 if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) { 194 row[0] = goffset; 195 row[1] = goffset + 1; 196 col[0] = goffset; 197 col[1] = goffset + 1; 198 col[2] = goffset; 199 col[3] = goffset + 1; 200 values[0] = 1.0; 201 values[1] = 0.0; 202 values[2] = 0.0; 203 values[3] = 1.0; 204 PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES)); 205 break; 206 } 207 208 Vm = xarr[offset + 1]; 209 210 /* Shunt injections */ 211 row[0] = goffset; 212 row[1] = goffset + 1; 213 col[0] = goffset; 214 col[1] = goffset + 1; 215 values[0] = values[1] = values[2] = values[3] = 0.0; 216 if (bus->ide != PV_BUS) { 217 values[1] = 2.0 * Vm * bus->gl / Sbase; 218 values[3] = -2.0 * Vm * bus->bl / Sbase; 219 } 220 PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES)); 221 } 222 223 PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges)); 224 for (i = 0; i < nconnedges; i++) { 225 EDGE_Power branch; 226 VERTEX_Power busf, bust; 227 PetscInt keyf, keyt; 228 PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; 229 const PetscInt *cone; 230 PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf; 231 232 e = connedges[i]; 233 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL)); 234 if (!branch->status) continue; 235 236 Gff = branch->yff[0]; 237 Bff = branch->yff[1]; 238 Gft = branch->yft[0]; 239 Bft = branch->yft[1]; 240 Gtf = branch->ytf[0]; 241 Btf = branch->ytf[1]; 242 Gtt = branch->ytt[0]; 243 Btt = branch->ytt[1]; 244 245 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone)); 246 vfrom = cone[0]; 247 vto = cone[1]; 248 249 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom)); 250 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto)); 251 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &goffsetfrom)); 252 PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vto, ALL_COMPONENTS, &goffsetto)); 253 254 if (goffsetto < 0) goffsetto = -goffsetto - 1; 255 256 thetaf = xarr[offsetfrom]; 257 Vmf = xarr[offsetfrom + 1]; 258 thetat = xarr[offsetto]; 259 Vmt = xarr[offsetto + 1]; 260 thetaft = thetaf - thetat; 261 thetatf = thetat - thetaf; 262 263 PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL)); 264 PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL)); 265 266 if (vfrom == vtx[v]) { 267 if (busf->ide != REF_BUS) { 268 /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */ 269 row[0] = goffsetfrom; 270 col[0] = goffsetfrom; 271 col[1] = goffsetfrom + 1; 272 col[2] = goffsetto; 273 col[3] = goffsetto + 1; 274 values[0] = Vmf * Vmt * (Gft * -PetscSinScalar(thetaft) + Bft * PetscCosScalar(thetaft)); /* df_dthetaf */ 275 values[1] = 2.0 * Gff * Vmf + Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmf */ 276 values[2] = Vmf * Vmt * (Gft * PetscSinScalar(thetaft) + Bft * -PetscCosScalar(thetaft)); /* df_dthetat */ 277 values[3] = Vmf * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmt */ 278 279 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 280 } 281 if (busf->ide != PV_BUS && busf->ide != REF_BUS) { 282 row[0] = goffsetfrom + 1; 283 col[0] = goffsetfrom; 284 col[1] = goffsetfrom + 1; 285 col[2] = goffsetto; 286 col[3] = goffsetto + 1; 287 /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */ 288 values[0] = Vmf * Vmt * (Bft * PetscSinScalar(thetaft) + Gft * PetscCosScalar(thetaft)); 289 values[1] = -2.0 * Bff * Vmf + Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 290 values[2] = Vmf * Vmt * (-Bft * PetscSinScalar(thetaft) + Gft * -PetscCosScalar(thetaft)); 291 values[3] = Vmf * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft)); 292 293 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 294 } 295 } else { 296 if (bust->ide != REF_BUS) { 297 row[0] = goffsetto; 298 col[0] = goffsetto; 299 col[1] = goffsetto + 1; 300 col[2] = goffsetfrom; 301 col[3] = goffsetfrom + 1; 302 /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */ 303 values[0] = Vmt * Vmf * (Gtf * -PetscSinScalar(thetatf) + Btf * PetscCosScalar(thetaft)); /* df_dthetat */ 304 values[1] = 2.0 * Gtt * Vmt + Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmt */ 305 values[2] = Vmt * Vmf * (Gtf * PetscSinScalar(thetatf) + Btf * -PetscCosScalar(thetatf)); /* df_dthetaf */ 306 values[3] = Vmt * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmf */ 307 308 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 309 } 310 if (bust->ide != PV_BUS && bust->ide != REF_BUS) { 311 row[0] = goffsetto + 1; 312 col[0] = goffsetto; 313 col[1] = goffsetto + 1; 314 col[2] = goffsetfrom; 315 col[3] = goffsetfrom + 1; 316 /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */ 317 values[0] = Vmt * Vmf * (Btf * PetscSinScalar(thetatf) + Gtf * PetscCosScalar(thetatf)); 318 values[1] = -2.0 * Btt * Vmt + Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 319 values[2] = Vmt * Vmf * (-Btf * PetscSinScalar(thetatf) + Gtf * -PetscCosScalar(thetatf)); 320 values[3] = Vmt * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf)); 321 322 PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES)); 323 } 324 } 325 } 326 if (!ghostvtex && bus->ide == PV_BUS) { 327 row[0] = goffset + 1; 328 col[0] = goffset + 1; 329 values[0] = 1.0; 330 PetscCall(MatSetValues(J, 1, row, 1, col, values, ADD_VALUES)); 331 } 332 } 333 } 334 } 335 PetscCall(VecRestoreArrayRead(localX, &xarr)); 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx) 340 { 341 DM networkdm; 342 Vec localX; 343 PetscInt ne, nv; 344 const PetscInt *vtx, *edges; 345 346 PetscFunctionBegin; 347 PetscCall(MatZeroEntries(J)); 348 349 PetscCall(SNESGetDM(snes, &networkdm)); 350 PetscCall(DMGetLocalVector(networkdm, &localX)); 351 352 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 353 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 354 355 /* Form Jacobian for first subnetwork */ 356 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 357 PetscCall(FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx)); 358 359 /* Form Jacobian for second subnetwork */ 360 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges)); 361 PetscCall(FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx)); 362 363 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 364 365 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 366 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 367 PetscFunctionReturn(PETSC_SUCCESS); 368 } 369 370 PetscErrorCode SetInitialValues_Subnet(DM networkdm, Vec localX, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx) 371 { 372 VERTEX_Power bus; 373 PetscInt i; 374 GEN gen; 375 PetscBool ghostvtex; 376 PetscScalar *xarr; 377 PetscInt key, numComps, j, offset; 378 void *component; 379 380 PetscFunctionBegin; 381 PetscCall(VecGetArray(localX, &xarr)); 382 for (i = 0; i < nv; i++) { 383 PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex)); 384 if (ghostvtex) continue; 385 386 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset)); 387 PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &numComps)); 388 for (j = 0; j < numComps; j++) { 389 PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, &component, NULL)); 390 if (key == 1) { 391 bus = (VERTEX_Power)(component); 392 xarr[offset] = bus->va * PETSC_PI / 180.0; 393 xarr[offset + 1] = bus->vm; 394 } else if (key == 2) { 395 gen = (GEN)(component); 396 if (!gen->status) continue; 397 xarr[offset + 1] = gen->vs; 398 break; 399 } 400 } 401 } 402 PetscCall(VecRestoreArray(localX, &xarr)); 403 PetscFunctionReturn(PETSC_SUCCESS); 404 } 405 406 PetscErrorCode SetInitialValues(DM networkdm, Vec X, void *appctx) 407 { 408 PetscInt nv, ne; 409 const PetscInt *vtx, *edges; 410 Vec localX; 411 412 PetscFunctionBegin; 413 PetscCall(DMGetLocalVector(networkdm, &localX)); 414 415 PetscCall(VecSet(X, 0.0)); 416 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 417 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 418 419 /* Set initial guess for first subnetwork */ 420 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 421 PetscCall(SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx)); 422 423 /* Set initial guess for second subnetwork */ 424 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges)); 425 PetscCall(SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx)); 426 427 PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X)); 428 PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X)); 429 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 430 PetscFunctionReturn(PETSC_SUCCESS); 431 } 432 433 int main(int argc, char **argv) 434 { 435 char pfdata_file[PETSC_MAX_PATH_LEN] = "case9.m"; 436 PFDATA *pfdata1, *pfdata2; 437 PetscInt numEdges1 = 0, numEdges2 = 0; 438 PetscInt *edgelist1 = NULL, *edgelist2 = NULL, componentkey[4]; 439 DM networkdm; 440 UserCtx_Power User; 441 PetscLogStage stage1, stage2; 442 PetscMPIInt rank; 443 PetscInt nsubnet = 2, nv, ne, i, j, genj, loadj; 444 const PetscInt *vtx, *edges; 445 Vec X, F; 446 Mat J; 447 SNES snes; 448 449 PetscFunctionBeginUser; 450 PetscCall(PetscInitialize(&argc, &argv, "poweroptions", help)); 451 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 452 { 453 /* introduce the const crank so the clang static analyzer realizes that if it enters any of the if (crank) then it must have entered the first */ 454 /* this is an experiment to see how the analyzer reacts */ 455 const PetscMPIInt crank = rank; 456 457 /* Create an empty network object */ 458 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm)); 459 460 /* Register the components in the network */ 461 PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &componentkey[0])); 462 PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &componentkey[1])); 463 PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &componentkey[2])); 464 PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &componentkey[3])); 465 466 PetscCall(PetscLogStageRegister("Read Data", &stage1)); 467 PetscCall(PetscLogStagePush(stage1)); 468 /* READ THE DATA */ 469 if (!crank) { 470 /* Only rank 0 reads the data */ 471 PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL)); 472 /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */ 473 474 /* READ DATA FOR THE FIRST SUBNETWORK */ 475 PetscCall(PetscNew(&pfdata1)); 476 PetscCall(PFReadMatPowerData(pfdata1, pfdata_file)); 477 User.Sbase = pfdata1->sbase; 478 479 numEdges1 = pfdata1->nbranch; 480 PetscCall(PetscMalloc1(2 * numEdges1, &edgelist1)); 481 PetscCall(GetListofEdges_Power(pfdata1, edgelist1)); 482 483 /* READ DATA FOR THE SECOND SUBNETWORK */ 484 PetscCall(PetscNew(&pfdata2)); 485 PetscCall(PFReadMatPowerData(pfdata2, pfdata_file)); 486 User.Sbase = pfdata2->sbase; 487 488 numEdges2 = pfdata2->nbranch; 489 PetscCall(PetscMalloc1(2 * numEdges2, &edgelist2)); 490 PetscCall(GetListofEdges_Power(pfdata2, edgelist2)); 491 } 492 493 PetscCall(PetscLogStagePop()); 494 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 495 PetscCall(PetscLogStageRegister("Create network", &stage2)); 496 PetscCall(PetscLogStagePush(stage2)); 497 498 /* Set number of nodes/edges and edge connectivity */ 499 PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, nsubnet)); 500 PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges1, edgelist1, NULL)); 501 PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges2, edgelist2, NULL)); 502 503 /* Set up the network layout */ 504 PetscCall(DMNetworkLayoutSetUp(networkdm)); 505 506 /* Add network components only process 0 has any data to add*/ 507 if (!crank) { 508 genj = 0; 509 loadj = 0; 510 511 /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */ 512 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 513 514 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata1->branch[i], 0)); 515 516 for (i = 0; i < nv; i++) { 517 PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata1->bus[i], 2)); 518 if (pfdata1->bus[i].ngen) { 519 for (j = 0; j < pfdata1->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata1->gen[genj++], 0)); 520 } 521 if (pfdata1->bus[i].nload) { 522 for (j = 0; j < pfdata1->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata1->load[loadj++], 0)); 523 } 524 } 525 526 genj = 0; 527 loadj = 0; 528 529 /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */ 530 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges)); 531 532 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata2->branch[i], 0)); 533 534 for (i = 0; i < nv; i++) { 535 PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata2->bus[i], 2)); 536 if (pfdata2->bus[i].ngen) { 537 for (j = 0; j < pfdata2->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata2->gen[genj++], 0)); 538 } 539 if (pfdata2->bus[i].nload) { 540 for (j = 0; j < pfdata2->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata2->load[loadj++], 0)); 541 } 542 } 543 } 544 545 /* Set up DM for use */ 546 PetscCall(DMSetUp(networkdm)); 547 548 if (!crank) { 549 PetscCall(PetscFree(edgelist1)); 550 PetscCall(PetscFree(edgelist2)); 551 } 552 553 if (!crank) { 554 PetscCall(PetscFree(pfdata1->bus)); 555 PetscCall(PetscFree(pfdata1->gen)); 556 PetscCall(PetscFree(pfdata1->branch)); 557 PetscCall(PetscFree(pfdata1->load)); 558 PetscCall(PetscFree(pfdata1)); 559 560 PetscCall(PetscFree(pfdata2->bus)); 561 PetscCall(PetscFree(pfdata2->gen)); 562 PetscCall(PetscFree(pfdata2->branch)); 563 PetscCall(PetscFree(pfdata2->load)); 564 PetscCall(PetscFree(pfdata2)); 565 } 566 567 /* Distribute networkdm to multiple processes */ 568 PetscCall(DMNetworkDistribute(&networkdm, 0)); 569 570 PetscCall(PetscLogStagePop()); 571 572 /* Broadcast Sbase to all processors */ 573 PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD)); 574 575 PetscCall(DMCreateGlobalVector(networkdm, &X)); 576 PetscCall(VecDuplicate(X, &F)); 577 578 PetscCall(DMCreateMatrix(networkdm, &J)); 579 PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 580 581 PetscCall(SetInitialValues(networkdm, X, &User)); 582 583 /* HOOK UP SOLVER */ 584 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 585 PetscCall(SNESSetDM(snes, networkdm)); 586 PetscCall(SNESSetFunction(snes, F, FormFunction, &User)); 587 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &User)); 588 PetscCall(SNESSetFromOptions(snes)); 589 590 PetscCall(SNESSolve(snes, NULL, X)); 591 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 592 593 PetscCall(VecDestroy(&X)); 594 PetscCall(VecDestroy(&F)); 595 PetscCall(MatDestroy(&J)); 596 597 PetscCall(SNESDestroy(&snes)); 598 PetscCall(DMDestroy(&networkdm)); 599 } 600 PetscCall(PetscFinalize()); 601 return 0; 602 } 603 604 /*TEST 605 606 build: 607 depends: PFReadData.c pffunctions.c 608 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 609 610 test: 611 args: -snes_rtol 1.e-3 612 localrunfiles: poweroptions case9.m 613 output_file: output/power_1.out 614 615 test: 616 suffix: 2 617 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 618 nsize: 4 619 localrunfiles: poweroptions case9.m 620 output_file: output/power_1.out 621 622 TEST*/ 623