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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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 #if defined(PETSC_USE_LOG) 442 PetscLogStage stage1, stage2; 443 #endif 444 PetscMPIInt rank; 445 PetscInt nsubnet = 2, nv, ne, i, j, genj, loadj; 446 const PetscInt *vtx, *edges; 447 Vec X, F; 448 Mat J; 449 SNES snes; 450 451 PetscFunctionBeginUser; 452 PetscCall(PetscInitialize(&argc, &argv, "poweroptions", help)); 453 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 454 { 455 /* 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 */ 456 /* this is an experiment to see how the analyzer reacts */ 457 const PetscMPIInt crank = rank; 458 459 /* Create an empty network object */ 460 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm)); 461 462 /* Register the components in the network */ 463 PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &componentkey[0])); 464 PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &componentkey[1])); 465 PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &componentkey[2])); 466 PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &componentkey[3])); 467 468 PetscCall(PetscLogStageRegister("Read Data", &stage1)); 469 PetscLogStagePush(stage1); 470 /* READ THE DATA */ 471 if (!crank) { 472 /* Only rank 0 reads the data */ 473 PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL)); 474 /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */ 475 476 /* READ DATA FOR THE FIRST SUBNETWORK */ 477 PetscCall(PetscNew(&pfdata1)); 478 PetscCall(PFReadMatPowerData(pfdata1, pfdata_file)); 479 User.Sbase = pfdata1->sbase; 480 481 numEdges1 = pfdata1->nbranch; 482 PetscCall(PetscMalloc1(2 * numEdges1, &edgelist1)); 483 PetscCall(GetListofEdges_Power(pfdata1, edgelist1)); 484 485 /* READ DATA FOR THE SECOND SUBNETWORK */ 486 PetscCall(PetscNew(&pfdata2)); 487 PetscCall(PFReadMatPowerData(pfdata2, pfdata_file)); 488 User.Sbase = pfdata2->sbase; 489 490 numEdges2 = pfdata2->nbranch; 491 PetscCall(PetscMalloc1(2 * numEdges2, &edgelist2)); 492 PetscCall(GetListofEdges_Power(pfdata2, edgelist2)); 493 } 494 495 PetscLogStagePop(); 496 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 497 PetscCall(PetscLogStageRegister("Create network", &stage2)); 498 PetscLogStagePush(stage2); 499 500 /* Set number of nodes/edges and edge connectivity */ 501 PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, nsubnet)); 502 PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges1, edgelist1, NULL)); 503 PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges2, edgelist2, NULL)); 504 505 /* Set up the network layout */ 506 PetscCall(DMNetworkLayoutSetUp(networkdm)); 507 508 /* Add network components only process 0 has any data to add*/ 509 if (!crank) { 510 genj = 0; 511 loadj = 0; 512 513 /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */ 514 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 515 516 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata1->branch[i], 0)); 517 518 for (i = 0; i < nv; i++) { 519 PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata1->bus[i], 2)); 520 if (pfdata1->bus[i].ngen) { 521 for (j = 0; j < pfdata1->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata1->gen[genj++], 0)); 522 } 523 if (pfdata1->bus[i].nload) { 524 for (j = 0; j < pfdata1->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata1->load[loadj++], 0)); 525 } 526 } 527 528 genj = 0; 529 loadj = 0; 530 531 /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */ 532 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges)); 533 534 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata2->branch[i], 0)); 535 536 for (i = 0; i < nv; i++) { 537 PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata2->bus[i], 2)); 538 if (pfdata2->bus[i].ngen) { 539 for (j = 0; j < pfdata2->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata2->gen[genj++], 0)); 540 } 541 if (pfdata2->bus[i].nload) { 542 for (j = 0; j < pfdata2->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata2->load[loadj++], 0)); 543 } 544 } 545 } 546 547 /* Set up DM for use */ 548 PetscCall(DMSetUp(networkdm)); 549 550 if (!crank) { 551 PetscCall(PetscFree(edgelist1)); 552 PetscCall(PetscFree(edgelist2)); 553 } 554 555 if (!crank) { 556 PetscCall(PetscFree(pfdata1->bus)); 557 PetscCall(PetscFree(pfdata1->gen)); 558 PetscCall(PetscFree(pfdata1->branch)); 559 PetscCall(PetscFree(pfdata1->load)); 560 PetscCall(PetscFree(pfdata1)); 561 562 PetscCall(PetscFree(pfdata2->bus)); 563 PetscCall(PetscFree(pfdata2->gen)); 564 PetscCall(PetscFree(pfdata2->branch)); 565 PetscCall(PetscFree(pfdata2->load)); 566 PetscCall(PetscFree(pfdata2)); 567 } 568 569 /* Distribute networkdm to multiple processes */ 570 PetscCall(DMNetworkDistribute(&networkdm, 0)); 571 572 PetscLogStagePop(); 573 574 /* Broadcast Sbase to all processors */ 575 PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD)); 576 577 PetscCall(DMCreateGlobalVector(networkdm, &X)); 578 PetscCall(VecDuplicate(X, &F)); 579 580 PetscCall(DMCreateMatrix(networkdm, &J)); 581 PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 582 583 PetscCall(SetInitialValues(networkdm, X, &User)); 584 585 /* HOOK UP SOLVER */ 586 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 587 PetscCall(SNESSetDM(snes, networkdm)); 588 PetscCall(SNESSetFunction(snes, F, FormFunction, &User)); 589 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &User)); 590 PetscCall(SNESSetFromOptions(snes)); 591 592 PetscCall(SNESSolve(snes, NULL, X)); 593 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 594 595 PetscCall(VecDestroy(&X)); 596 PetscCall(VecDestroy(&F)); 597 PetscCall(MatDestroy(&J)); 598 599 PetscCall(SNESDestroy(&snes)); 600 PetscCall(DMDestroy(&networkdm)); 601 } 602 PetscCall(PetscFinalize()); 603 return 0; 604 } 605 606 /*TEST 607 608 build: 609 depends: PFReadData.c pffunctions.c 610 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 611 612 test: 613 args: -snes_rtol 1.e-3 614 localrunfiles: poweroptions case9.m 615 output_file: output/power_1.out 616 617 test: 618 suffix: 2 619 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 620 nsize: 4 621 localrunfiles: poweroptions case9.m 622 output_file: output/power_1.out 623 624 TEST*/ 625