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