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