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