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