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 CHKERRQ(VecGetArrayRead(localX,&xarr)); 26 CHKERRQ(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 CHKERRQ(DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex)); 40 CHKERRQ(DMNetworkGetNumComponents(networkdm,vtx[v],&numComps)); 41 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset)); 42 for (j = 0; j < numComps; j++) { 43 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMNetworkGetConnectedVertices(networkdm,e,&cone)); 85 vfrom = cone[0]; 86 vto = cone[1]; 87 88 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom)); 89 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(localX,&xarr)); 126 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&networkdm)); 139 CHKERRQ(DMGetLocalVector(networkdm,&localX)); 140 CHKERRQ(DMGetLocalVector(networkdm,&localF)); 141 CHKERRQ(VecSet(F,0.0)); 142 143 CHKERRQ(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 144 CHKERRQ(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 145 146 CHKERRQ(DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF)); 147 CHKERRQ(DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF)); 148 149 /* Form Function for first subnetwork */ 150 CHKERRQ(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 151 CHKERRQ(FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx)); 152 153 /* Form Function for second subnetwork */ 154 CHKERRQ(DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges)); 155 CHKERRQ(FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx)); 156 157 CHKERRQ(DMRestoreLocalVector(networkdm,&localX)); 158 159 CHKERRQ(DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F)); 160 CHKERRQ(DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F)); 161 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex)); 188 CHKERRQ(DMNetworkGetNumComponents(networkdm,vtx[v],&numComps)); 189 for (j = 0; j < numComps; j++) { 190 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset)); 191 CHKERRQ(DMNetworkGetGlobalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&goffset)); 192 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(J,2,row,2,col,values,ADD_VALUES)); 219 } 220 221 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMNetworkGetConnectedVertices(networkdm,e,&cone)); 244 vfrom = cone[0]; 245 vto = cone[1]; 246 247 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom)); 248 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto)); 249 CHKERRQ(DMNetworkGetGlobalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&goffsetfrom)); 250 CHKERRQ(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 CHKERRQ(DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf,NULL)); 262 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(J,1,row,1,col,values,ADD_VALUES)); 316 } 317 } 318 } 319 } 320 CHKERRQ(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 CHKERRQ(MatZeroEntries(J)); 333 334 CHKERRQ(SNESGetDM(snes,&networkdm)); 335 CHKERRQ(DMGetLocalVector(networkdm,&localX)); 336 337 CHKERRQ(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 338 CHKERRQ(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 339 340 /* Form Jacobian for first subnetwork */ 341 CHKERRQ(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 342 CHKERRQ(FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx)); 343 344 /* Form Jacobian for second subnetwork */ 345 CHKERRQ(DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges)); 346 CHKERRQ(FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx)); 347 348 CHKERRQ(DMRestoreLocalVector(networkdm,&localX)); 349 350 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 351 CHKERRQ(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 CHKERRQ(VecGetArray(localX,&xarr)); 367 for (i = 0; i < nv; i++) { 368 CHKERRQ(DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex)); 369 if (ghostvtex) continue; 370 371 CHKERRQ(DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset)); 372 CHKERRQ(DMNetworkGetNumComponents(networkdm,vtx[i],&numComps)); 373 for (j=0; j < numComps; j++) { 374 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMGetLocalVector(networkdm,&localX)); 399 400 CHKERRQ(VecSet(X,0.0)); 401 CHKERRQ(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 402 CHKERRQ(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 403 404 /* Set initial guess for first subnetwork */ 405 CHKERRQ(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 406 CHKERRQ(SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx)); 407 408 /* Set initial guess for second subnetwork */ 409 CHKERRQ(DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges)); 410 CHKERRQ(SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx)); 411 412 CHKERRQ(DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X)); 413 CHKERRQ(DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X)); 414 CHKERRQ(DMRestoreLocalVector(networkdm,&localX)); 415 PetscFunctionReturn(0); 416 } 417 418 int main(int argc,char ** argv) 419 { 420 PetscErrorCode ierr; 421 char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m"; 422 PFDATA *pfdata1,*pfdata2; 423 PetscInt numEdges1=0,numEdges2=0; 424 PetscInt *edgelist1 = NULL,*edgelist2 = NULL,componentkey[4]; 425 DM networkdm; 426 UserCtx_Power User; 427 #if defined(PETSC_USE_LOG) 428 PetscLogStage stage1,stage2; 429 #endif 430 PetscMPIInt rank; 431 PetscInt nsubnet = 2,nv,ne,i,j,genj,loadj; 432 const PetscInt *vtx,*edges; 433 Vec X,F; 434 Mat J; 435 SNES snes; 436 437 ierr = PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr; 438 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 439 { 440 /* 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 */ 441 /* this is an experiment to see how the analyzer reacts */ 442 const PetscMPIInt crank = rank; 443 444 /* Create an empty network object */ 445 CHKERRQ(DMNetworkCreate(PETSC_COMM_WORLD,&networkdm)); 446 447 /* Register the components in the network */ 448 CHKERRQ(DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0])); 449 CHKERRQ(DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1])); 450 CHKERRQ(DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2])); 451 CHKERRQ(DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3])); 452 453 CHKERRQ(PetscLogStageRegister("Read Data",&stage1)); 454 PetscLogStagePush(stage1); 455 /* READ THE DATA */ 456 if (!crank) { 457 /* Only rank 0 reads the data */ 458 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL)); 459 /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */ 460 461 /* READ DATA FOR THE FIRST SUBNETWORK */ 462 CHKERRQ(PetscNew(&pfdata1)); 463 CHKERRQ(PFReadMatPowerData(pfdata1,pfdata_file)); 464 User.Sbase = pfdata1->sbase; 465 466 numEdges1 = pfdata1->nbranch; 467 CHKERRQ(PetscMalloc1(2*numEdges1,&edgelist1)); 468 CHKERRQ(GetListofEdges_Power(pfdata1,edgelist1)); 469 470 /* READ DATA FOR THE SECOND SUBNETWORK */ 471 CHKERRQ(PetscNew(&pfdata2)); 472 CHKERRQ(PFReadMatPowerData(pfdata2,pfdata_file)); 473 User.Sbase = pfdata2->sbase; 474 475 numEdges2 = pfdata2->nbranch; 476 CHKERRQ(PetscMalloc1(2*numEdges2,&edgelist2)); 477 CHKERRQ(GetListofEdges_Power(pfdata2,edgelist2)); 478 } 479 480 PetscLogStagePop(); 481 CHKERRMPI(MPI_Barrier(PETSC_COMM_WORLD)); 482 CHKERRQ(PetscLogStageRegister("Create network",&stage2)); 483 PetscLogStagePush(stage2); 484 485 /* Set number of nodes/edges and edge connectivity */ 486 CHKERRQ(DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,nsubnet)); 487 CHKERRQ(DMNetworkAddSubnetwork(networkdm,"",numEdges1,edgelist1,NULL)); 488 CHKERRQ(DMNetworkAddSubnetwork(networkdm,"",numEdges2,edgelist2,NULL)); 489 490 /* Set up the network layout */ 491 CHKERRQ(DMNetworkLayoutSetUp(networkdm)); 492 493 /* Add network components only process 0 has any data to add*/ 494 if (!crank) { 495 genj=0; loadj=0; 496 497 /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */ 498 CHKERRQ(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 499 500 for (i = 0; i < ne; i++) { 501 CHKERRQ(DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i],0)); 502 } 503 504 for (i = 0; i < nv; i++) { 505 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i],2)); 506 if (pfdata1->bus[i].ngen) { 507 for (j = 0; j < pfdata1->bus[i].ngen; j++) { 508 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++],0)); 509 } 510 } 511 if (pfdata1->bus[i].nload) { 512 for (j=0; j < pfdata1->bus[i].nload; j++) { 513 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++],0)); 514 } 515 } 516 } 517 518 genj=0; loadj=0; 519 520 /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */ 521 CHKERRQ(DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges)); 522 523 for (i = 0; i < ne; i++) { 524 CHKERRQ(DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i],0)); 525 } 526 527 for (i = 0; i < nv; i++) { 528 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i],2)); 529 if (pfdata2->bus[i].ngen) { 530 for (j = 0; j < pfdata2->bus[i].ngen; j++) { 531 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++],0)); 532 } 533 } 534 if (pfdata2->bus[i].nload) { 535 for (j=0; j < pfdata2->bus[i].nload; j++) { 536 CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++],0)); 537 } 538 } 539 } 540 } 541 542 /* Set up DM for use */ 543 CHKERRQ(DMSetUp(networkdm)); 544 545 if (!crank) { 546 CHKERRQ(PetscFree(edgelist1)); 547 CHKERRQ(PetscFree(edgelist2)); 548 } 549 550 if (!crank) { 551 CHKERRQ(PetscFree(pfdata1->bus)); 552 CHKERRQ(PetscFree(pfdata1->gen)); 553 CHKERRQ(PetscFree(pfdata1->branch)); 554 CHKERRQ(PetscFree(pfdata1->load)); 555 CHKERRQ(PetscFree(pfdata1)); 556 557 CHKERRQ(PetscFree(pfdata2->bus)); 558 CHKERRQ(PetscFree(pfdata2->gen)); 559 CHKERRQ(PetscFree(pfdata2->branch)); 560 CHKERRQ(PetscFree(pfdata2->load)); 561 CHKERRQ(PetscFree(pfdata2)); 562 } 563 564 /* Distribute networkdm to multiple processes */ 565 CHKERRQ(DMNetworkDistribute(&networkdm,0)); 566 567 PetscLogStagePop(); 568 569 /* Broadcast Sbase to all processors */ 570 CHKERRMPI(MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD)); 571 572 CHKERRQ(DMCreateGlobalVector(networkdm,&X)); 573 CHKERRQ(VecDuplicate(X,&F)); 574 575 CHKERRQ(DMCreateMatrix(networkdm,&J)); 576 CHKERRQ(MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 577 578 CHKERRQ(SetInitialValues(networkdm,X,&User)); 579 580 /* HOOK UP SOLVER */ 581 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 582 CHKERRQ(SNESSetDM(snes,networkdm)); 583 CHKERRQ(SNESSetFunction(snes,F,FormFunction,&User)); 584 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian,&User)); 585 CHKERRQ(SNESSetFromOptions(snes)); 586 587 CHKERRQ(SNESSolve(snes,NULL,X)); 588 /* CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 589 590 CHKERRQ(VecDestroy(&X)); 591 CHKERRQ(VecDestroy(&F)); 592 CHKERRQ(MatDestroy(&J)); 593 594 CHKERRQ(SNESDestroy(&snes)); 595 CHKERRQ(DMDestroy(&networkdm)); 596 } 597 ierr = PetscFinalize(); 598 return ierr; 599 } 600 601 /*TEST 602 603 build: 604 depends: PFReadData.c pffunctions.c 605 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 606 607 test: 608 args: -snes_rtol 1.e-3 609 localrunfiles: poweroptions case9.m 610 output_file: output/power_1.out 611 612 test: 613 suffix: 2 614 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 615 nsize: 4 616 localrunfiles: poweroptions case9.m 617 output_file: output/power_1.out 618 619 TEST*/ 620