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