1 static char help[] = "This example demonstrates the use of DMNetwork interface 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 See 'Evaluation of overlapping restricted additive schwarz preconditioning for parallel solution \n\ 4 of very large power flow problems' https://dl.acm.org/citation.cfm?id=2536784).\n\ 5 The data file format used is from the MatPower package (http://www.pserc.cornell.edu//matpower/).\n\ 6 Run this program: mpiexec -n <n> ./pf\n\ 7 mpiexec -n <n> ./pfc \n"; 8 9 /* T 10 Concepts: DMNetwork 11 Concepts: PETSc SNES solver 12 */ 13 14 #include "power.h" 15 #include <petscdmnetwork.h> 16 17 PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx) 18 { 19 PetscErrorCode ierr; 20 DM networkdm; 21 UserCtx_Power *User=(UserCtx_Power*)appctx; 22 Vec localX,localF; 23 PetscInt nv,ne; 24 const PetscInt *vtx,*edges; 25 26 PetscFunctionBegin; 27 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 28 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 29 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 30 ierr = VecSet(F,0.0);CHKERRQ(ierr); 31 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 32 33 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 34 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 35 36 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 37 ierr = FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,User);CHKERRQ(ierr); 38 39 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 40 41 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 42 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 43 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx) 48 { 49 PetscErrorCode ierr; 50 PetscInt vStart,vEnd,nv,ne; 51 const PetscInt *vtx,*edges; 52 Vec localX; 53 UserCtx_Power *user_power=(UserCtx_Power*)appctx; 54 55 PetscFunctionBegin; 56 ierr = DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);CHKERRQ(ierr); 57 58 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 59 60 ierr = VecSet(X,0.0);CHKERRQ(ierr); 61 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 62 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 63 64 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 65 ierr = SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,user_power);CHKERRQ(ierr); 66 67 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 68 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 69 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 int main(int argc,char ** argv) 74 { 75 PetscErrorCode ierr; 76 char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m"; 77 PFDATA *pfdata; 78 PetscInt numEdges=0,numVertices=0; 79 PetscInt *edges = NULL; 80 PetscInt i; 81 DM networkdm; 82 UserCtx_Power User; 83 #if defined(PETSC_USE_LOG) 84 PetscLogStage stage1,stage2; 85 #endif 86 PetscMPIInt rank; 87 PetscInt eStart, eEnd, vStart, vEnd,j; 88 PetscInt genj,loadj; 89 Vec X,F; 90 Mat J; 91 SNES snes; 92 93 ierr = PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr; 94 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 95 { 96 /* 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 */ 97 /* this is an experiment to see how the analyzer reacts */ 98 const PetscMPIInt crank = rank; 99 100 /* Create an empty network object */ 101 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 102 /* Register the components in the network */ 103 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);CHKERRQ(ierr); 104 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);CHKERRQ(ierr); 105 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);CHKERRQ(ierr); 106 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);CHKERRQ(ierr); 107 108 ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr); 109 PetscLogStagePush(stage1); 110 /* READ THE DATA */ 111 if (!crank) { 112 /* READ DATA */ 113 /* Only rank 0 reads the data */ 114 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);CHKERRQ(ierr); 115 ierr = PetscNew(&pfdata);CHKERRQ(ierr); 116 ierr = PFReadMatPowerData(pfdata,pfdata_file);CHKERRQ(ierr); 117 User.Sbase = pfdata->sbase; 118 119 numEdges = pfdata->nbranch; 120 numVertices = pfdata->nbus; 121 122 ierr = PetscMalloc1(2*numEdges,&edges);CHKERRQ(ierr); 123 ierr = GetListofEdges_Power(pfdata,edges);CHKERRQ(ierr); 124 } 125 126 /* If external option activated. Introduce error in jacobian */ 127 ierr = PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);CHKERRQ(ierr); 128 129 PetscLogStagePop(); 130 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 131 ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr); 132 PetscLogStagePush(stage2); 133 /* Set number of nodes/edges */ 134 ierr = DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);CHKERRQ(ierr); 135 ierr = DMNetworkAddSubnetwork(networkdm,"",numVertices,numEdges,edges,NULL);CHKERRQ(ierr); 136 137 /* Set up the network layout */ 138 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 139 140 if (!crank) { 141 ierr = PetscFree(edges);CHKERRQ(ierr); 142 } 143 144 /* Add network components only process 0 has any data to add */ 145 if (!crank) { 146 genj=0; loadj=0; 147 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 148 for (i = eStart; i < eEnd; i++) { 149 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart],0);CHKERRQ(ierr); 150 } 151 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 152 for (i = vStart; i < vEnd; i++) { 153 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart],2);CHKERRQ(ierr); 154 if (pfdata->bus[i-vStart].ngen) { 155 for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) { 156 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++],0);CHKERRQ(ierr); 157 } 158 } 159 if (pfdata->bus[i-vStart].nload) { 160 for (j=0; j < pfdata->bus[i-vStart].nload; j++) { 161 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++],0);CHKERRQ(ierr); 162 } 163 } 164 } 165 } 166 167 /* Set up DM for use */ 168 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 169 170 if (!crank) { 171 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 172 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 173 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 174 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 175 ierr = PetscFree(pfdata);CHKERRQ(ierr); 176 } 177 178 /* Distribute networkdm to multiple processes */ 179 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 180 181 PetscLogStagePop(); 182 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 183 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 184 185 #if 0 186 EDGE_Power edge; 187 PetscInt key,kk,numComponents; 188 VERTEX_Power bus; 189 GEN gen; 190 LOAD load; 191 192 for (i = eStart; i < eEnd; i++) { 193 ierr = DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);CHKERRQ(ierr); 194 ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr); 195 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);CHKERRQ(ierr); 196 } 197 198 for (i = vStart; i < vEnd; i++) { 199 ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr); 200 for (kk=0; kk < numComponents; kk++) { 201 ierr = DMNetworkGetComponent(networkdm,i,kk,&key,&component);CHKERRQ(ierr); 202 if (key == 1) { 203 bus = (VERTEX_Power)(component); 204 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);CHKERRQ(ierr); 205 } else if (key == 2) { 206 gen = (GEN)(component); 207 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);CHKERRQ(ierr); 208 } else if (key == 3) { 209 load = (LOAD)(component); 210 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);CHKERRQ(ierr); 211 } 212 } 213 } 214 #endif 215 /* Broadcast Sbase to all processors */ 216 ierr = MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRMPI(ierr); 217 218 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 219 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 220 221 ierr = DMCreateMatrix(networkdm,&J);CHKERRQ(ierr); 222 ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 223 224 ierr = SetInitialValues(networkdm,X,&User);CHKERRQ(ierr); 225 226 /* HOOK UP SOLVER */ 227 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 228 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 229 ierr = SNESSetFunction(snes,F,FormFunction,&User);CHKERRQ(ierr); 230 ierr = SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);CHKERRQ(ierr); 231 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 232 233 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 234 /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 235 236 ierr = VecDestroy(&X);CHKERRQ(ierr); 237 ierr = VecDestroy(&F);CHKERRQ(ierr); 238 ierr = MatDestroy(&J);CHKERRQ(ierr); 239 240 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 241 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 242 } 243 ierr = PetscFinalize(); 244 return ierr; 245 } 246 247 /*TEST 248 249 build: 250 depends: PFReadData.c pffunctions.c 251 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 252 253 test: 254 args: -snes_rtol 1.e-3 255 localrunfiles: poweroptions case9.m 256 output_file: output/power_1.out 257 258 test: 259 suffix: 2 260 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 261 nsize: 4 262 localrunfiles: poweroptions case9.m 263 output_file: output/power_1.out 264 265 TEST*/ 266