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 = DMNetworkGetSubnetworkInfo(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 = DMNetworkGetSubnetworkInfo(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);CHKERRQ(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);CHKERRQ(ierr); 131 ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr); 132 PetscLogStagePush(stage2); 133 /* Set number of nodes/edges */ 134 ierr = DMNetworkSetSizes(networkdm,1,&numVertices,&numEdges,0,NULL);CHKERRQ(ierr); 135 /* Add edge connectivity */ 136 ierr = DMNetworkSetEdgeList(networkdm,&edges,NULL);CHKERRQ(ierr); 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]);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]);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++]);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++]);CHKERRQ(ierr); 162 } 163 } 164 /* Add number of variables */ 165 ierr = DMNetworkAddNumVariables(networkdm,i,2);CHKERRQ(ierr); 166 } 167 } 168 169 /* Set up DM for use */ 170 ierr = DMSetUp(networkdm);CHKERRQ(ierr); 171 172 if (!crank) { 173 ierr = PetscFree(pfdata->bus);CHKERRQ(ierr); 174 ierr = PetscFree(pfdata->gen);CHKERRQ(ierr); 175 ierr = PetscFree(pfdata->branch);CHKERRQ(ierr); 176 ierr = PetscFree(pfdata->load);CHKERRQ(ierr); 177 ierr = PetscFree(pfdata);CHKERRQ(ierr); 178 } 179 180 /* Distribute networkdm to multiple processes */ 181 ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr); 182 183 PetscLogStagePop(); 184 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 185 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 186 187 #if 0 188 EDGE_Power edge; 189 PetscInt key,kk,numComponents; 190 VERTEX_Power bus; 191 GEN gen; 192 LOAD load; 193 194 for (i = eStart; i < eEnd; i++) { 195 ierr = DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);CHKERRQ(ierr); 196 ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr); 197 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);CHKERRQ(ierr); 198 } 199 200 for (i = vStart; i < vEnd; i++) { 201 ierr = DMNetworkGetNumComponents(networkdm,i,&numComponents);CHKERRQ(ierr); 202 for (kk=0; kk < numComponents; kk++) { 203 ierr = DMNetworkGetComponent(networkdm,i,kk,&key,&component);CHKERRQ(ierr); 204 if (key == 1) { 205 bus = (VERTEX_Power)(component); 206 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);CHKERRQ(ierr); 207 } else if (key == 2) { 208 gen = (GEN)(component); 209 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);CHKERRQ(ierr); 210 } else if (key == 3) { 211 load = (LOAD)(component); 212 ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);CHKERRQ(ierr); 213 } 214 } 215 } 216 #endif 217 /* Broadcast Sbase to all processors */ 218 ierr = MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 219 220 ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr); 221 ierr = VecDuplicate(X,&F);CHKERRQ(ierr); 222 223 ierr = DMCreateMatrix(networkdm,&J);CHKERRQ(ierr); 224 ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 225 226 ierr = SetInitialValues(networkdm,X,&User);CHKERRQ(ierr); 227 228 /* HOOK UP SOLVER */ 229 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 230 ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr); 231 ierr = SNESSetFunction(snes,F,FormFunction,&User);CHKERRQ(ierr); 232 ierr = SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);CHKERRQ(ierr); 233 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 234 235 ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); 236 /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 237 238 ierr = VecDestroy(&X);CHKERRQ(ierr); 239 ierr = VecDestroy(&F);CHKERRQ(ierr); 240 ierr = MatDestroy(&J);CHKERRQ(ierr); 241 242 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 243 ierr = DMDestroy(&networkdm);CHKERRQ(ierr); 244 } 245 ierr = PetscFinalize(); 246 return ierr; 247 } 248 249 /*TEST 250 251 build: 252 depends: PFReadData.c pffunctions.c 253 requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED) 254 255 256 test: 257 args: -snes_rtol 1.e-3 258 localrunfiles: poweroptions case9.m 259 output_file: output/power_1.out 260 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 261 262 test: 263 suffix: 2 264 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 265 nsize: 4 266 localrunfiles: poweroptions case9.m 267 output_file: output/power_1.out 268 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 269 270 TEST*/ 271