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 PetscLogStage stage1,stage2; 84 PetscMPIInt rank; 85 PetscInt eStart, eEnd, vStart, vEnd,j; 86 PetscInt genj,loadj; 87 Vec X,F; 88 Mat J; 89 SNES snes; 90 91 ierr = PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr; 92 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 93 { 94 /* 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 */ 95 /* this is an experiment to see how the analyzer reacts */ 96 const PetscMPIInt crank = rank; 97 98 /* Create an empty network object */ 99 ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr); 100 /* Register the components in the network */ 101 ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);CHKERRQ(ierr); 102 ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);CHKERRQ(ierr); 103 ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);CHKERRQ(ierr); 104 ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);CHKERRQ(ierr); 105 106 ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr); 107 PetscLogStagePush(stage1); 108 /* READ THE DATA */ 109 if (!crank) { 110 /* READ DATA */ 111 /* Only rank 0 reads the data */ 112 ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);CHKERRQ(ierr); 113 ierr = PetscNew(&pfdata);CHKERRQ(ierr); 114 ierr = PFReadMatPowerData(pfdata,pfdata_file);CHKERRQ(ierr); 115 User.Sbase = pfdata->sbase; 116 117 numEdges = pfdata->nbranch; 118 numVertices = pfdata->nbus; 119 120 ierr = PetscMalloc1(2*numEdges,&edges);CHKERRQ(ierr); 121 ierr = GetListofEdges_Power(pfdata,edges);CHKERRQ(ierr); 122 } 123 124 /* If external option activated. Introduce error in jacobian */ 125 ierr = PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);CHKERRQ(ierr); 126 127 PetscLogStagePop(); 128 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr); 129 ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr); 130 PetscLogStagePush(stage2); 131 /* Set number of nodes/edges */ 132 ierr = DMNetworkSetSizes(networkdm,1,&numVertices,&numEdges,0,NULL);CHKERRQ(ierr); 133 /* Add edge connectivity */ 134 ierr = DMNetworkSetEdgeList(networkdm,&edges,NULL);CHKERRQ(ierr); 135 /* Set up the network layout */ 136 ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr); 137 138 if (!crank) { 139 ierr = PetscFree(edges);CHKERRQ(ierr); 140 } 141 142 /* Add network components only process 0 has any data to add*/ 143 if (!crank) { 144 genj=0; loadj=0; 145 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 146 for (i = eStart; i < eEnd; i++) { 147 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart]);CHKERRQ(ierr); 148 } 149 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 150 for (i = vStart; i < vEnd; i++) { 151 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart]);CHKERRQ(ierr); 152 if (pfdata->bus[i-vStart].ngen) { 153 for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) { 154 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++]);CHKERRQ(ierr); 155 } 156 } 157 if (pfdata->bus[i-vStart].nload) { 158 for (j=0; j < pfdata->bus[i-vStart].nload; j++) { 159 ierr = DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++]);CHKERRQ(ierr); 160 } 161 } 162 /* Add number of variables */ 163 ierr = DMNetworkAddNumVariables(networkdm,i,2);CHKERRQ(ierr); 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);CHKERRQ(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 254 test: 255 args: -snes_rtol 1.e-3 256 localrunfiles: poweroptions case9.m 257 output_file: output/power_1.out 258 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 259 260 test: 261 suffix: 2 262 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 263 nsize: 4 264 localrunfiles: poweroptions case9.m 265 output_file: output/power_1.out 266 requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED) 267 268 TEST*/ 269