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