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 = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr); 43 for (j = 0; j < numComps; j++) { 44 ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);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);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 = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr); 90 ierr = DMNetworkGetVariableOffset(networkdm,vto,&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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr); 195 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vtx[v],&goffset);CHKERRQ(ierr); 196 ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);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);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 = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr); 252 ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr); 253 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);CHKERRQ(ierr); 254 ierr = DMNetworkGetVariableGlobalOffset(networkdm,vto,&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);CHKERRQ(ierr); 266 ierr = DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust);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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetVariableOffset(networkdm,vtx[i],&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);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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetSubnetworkInfo(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; 432 DM networkdm; 433 PetscInt componentkey[4]; 434 UserCtx_Power User; 435 PetscLogStage stage1,stage2; 436 PetscMPIInt rank; 437 PetscInt nsubnet = 2; 438 PetscInt numVertices[2],numEdges[2]; 439 PetscInt *edgelist[2]; 440 PetscInt nv,ne; 441 const PetscInt *vtx; 442 const PetscInt *edges; 443 PetscInt i,j; 444 PetscInt genj,loadj; 445 Vec X,F; 446 Mat J; 447 SNES snes; 448 449 450 ierr = PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr; 451 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 452 { 453 /* 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 */ 454 /* this is an experiment to see how the analyzer reacts */ 455 const PetscMPIInt crank = rank; 456 457 /* Create an empty network object */ 458 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 459 460 /* Register the components in the network */ 461 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);CHKERRQ(ierr); 462 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);CHKERRQ(ierr); 463 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);CHKERRQ(ierr); 464 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);CHKERRQ(ierr); 465 466 ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr); 467 PetscLogStagePush(stage1); 468 /* READ THE DATA */ 469 if (!crank) { 470 /* Only rank 0 reads the data */ 471 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);CHKERRQ(ierr); 472 /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */ 473 474 /* READ DATA FOR THE FIRST SUBNETWORK */ 475 ierr = PetscNew(&pfdata1);CHKERRQ(ierr); 476 ierr = PFReadMatPowerData(pfdata1,pfdata_file);CHKERRQ(ierr); 477 User.Sbase = pfdata1->sbase; 478 479 numEdges1 = pfdata1->nbranch; 480 numVertices1 = pfdata1->nbus; 481 482 ierr = PetscMalloc1(2*numEdges1,&edgelist1);CHKERRQ(ierr); 483 ierr = GetListofEdges_Power(pfdata1,edgelist1);CHKERRQ(ierr); 484 485 /* READ DATA FOR THE SECOND SUBNETWORK */ 486 ierr = PetscNew(&pfdata2);CHKERRQ(ierr); 487 ierr = PFReadMatPowerData(pfdata2,pfdata_file);CHKERRQ(ierr); 488 User.Sbase = pfdata2->sbase; 489 490 numEdges2 = pfdata2->nbranch; 491 numVertices2 = pfdata2->nbus; 492 493 ierr = PetscMalloc1(2*numEdges2,&edgelist2);CHKERRQ(ierr); 494 ierr = GetListofEdges_Power(pfdata2,edgelist2);CHKERRQ(ierr); 495 } 496 497 PetscLogStagePop(); 498 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 499 ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr); 500 PetscLogStagePush(stage2); 501 502 /* Set number of nodes/edges */ 503 numVertices[0] = numVertices1; numVertices[1] = numVertices2; 504 numEdges[0] = numEdges1; numEdges[1] = numEdges2; 505 ierr = DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,0,NULL);CHKERRQ(ierr); 506 507 edgelist[0] = edgelist1; edgelist[1] = edgelist2; 508 509 /* Add edge connectivity */ 510 ierr = DMNetworkSetEdgeList(networkdm,edgelist,NULL);CHKERRQ(ierr); 511 512 /* Set up the network layout */ 513 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 514 515 /* Add network components only process 0 has any data to add*/ 516 if (!crank) { 517 genj=0; loadj=0; 518 519 /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */ 520 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 521 522 for (i = 0; i < ne; i++) { 523 ierr = DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i]);CHKERRQ(ierr); 524 } 525 526 for (i = 0; i < nv; i++) { 527 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i]);CHKERRQ(ierr); 528 if (pfdata1->bus[i].ngen) { 529 for (j = 0; j < pfdata1->bus[i].ngen; j++) { 530 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++]);CHKERRQ(ierr); 531 } 532 } 533 if (pfdata1->bus[i].nload) { 534 for (j=0; j < pfdata1->bus[i].nload; j++) { 535 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++]);CHKERRQ(ierr); 536 } 537 } 538 /* Add number of variables */ 539 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr); 540 } 541 542 genj=0; loadj=0; 543 544 /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */ 545 ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 546 547 for (i = 0; i < ne; i++) { 548 ierr = DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i]);CHKERRQ(ierr); 549 } 550 551 for (i = 0; i < nv; i++) { 552 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i]);CHKERRQ(ierr); 553 if (pfdata2->bus[i].ngen) { 554 for (j = 0; j < pfdata2->bus[i].ngen; j++) { 555 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++]);CHKERRQ(ierr); 556 } 557 } 558 if (pfdata2->bus[i].nload) { 559 for (j=0; j < pfdata2->bus[i].nload; j++) { 560 ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++]);CHKERRQ(ierr); 561 } 562 } 563 /* Add number of variables */ 564 ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr); 565 } 566 } 567 568 /* Set up DM for use */ 569 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 570 571 if (!crank) { 572 ierr = PetscFree(edgelist1);CHKERRQ(ierr); 573 ierr = PetscFree(edgelist2);CHKERRQ(ierr); 574 } 575 576 if (!crank) { 577 ierr = PetscFree(pfdata1->bus);CHKERRQ(ierr); 578 ierr = PetscFree(pfdata1->gen);CHKERRQ(ierr); 579 ierr = PetscFree(pfdata1->branch);CHKERRQ(ierr); 580 ierr = PetscFree(pfdata1->load);CHKERRQ(ierr); 581 ierr = PetscFree(pfdata1);CHKERRQ(ierr); 582 583 ierr = PetscFree(pfdata2->bus);CHKERRQ(ierr); 584 ierr = PetscFree(pfdata2->gen);CHKERRQ(ierr); 585 ierr = PetscFree(pfdata2->branch);CHKERRQ(ierr); 586 ierr = PetscFree(pfdata2->load);CHKERRQ(ierr); 587 ierr = PetscFree(pfdata2);CHKERRQ(ierr); 588 } 589 590 /* Distribute networkdm to multiple processes */ 591 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 592 593 PetscLogStagePop(); 594 595 /* Broadcast Sbase to all processors */ 596 ierr = MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 597 598 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 599 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 600 601 ierr = DMCreateMatrix(networkdm,&J);CHKERRQ(ierr); 602 ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 603 604 ierr = SetInitialValues(networkdm,X,&User);CHKERRQ(ierr); 605 606 /* HOOK UP SOLVER */ 607 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 608 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 609 ierr = SNESSetFunction(snes,F,FormFunction,&User);CHKERRQ(ierr); 610 ierr = SNESSetJacobian(snes,J,J,FormJacobian,&User);CHKERRQ(ierr); 611 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 612 613 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 614 /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 615 616 ierr = VecDestroy(&X);CHKERRQ(ierr); 617 ierr = VecDestroy(&F);CHKERRQ(ierr); 618 ierr = MatDestroy(&J);CHKERRQ(ierr); 619 620 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 621 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 622 } 623 ierr = PetscFinalize(); 624 return ierr; 625 } 626 627 /*TEST 628 629 build: 630 depends: PFReadData.c pffunctions.c 631 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 632 633 634 test: 635 args: -snes_rtol 1.e-3 636 localrunfiles: poweroptions case9.m 637 output_file: output/power2_1.out 638 requires: double !complex 639 640 test: 641 suffix: 2 642 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 643 nsize: 4 644 localrunfiles: poweroptions case9.m 645 output_file: output/power2_1.out 646 requires: double !complex 647 648 TEST*/ 649