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