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 char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m"; 74 PFDATA *pfdata; 75 PetscInt numEdges=0; 76 PetscInt *edges = NULL; 77 PetscInt i; 78 DM networkdm; 79 UserCtx_Power User; 80 #if defined(PETSC_USE_LOG) 81 PetscLogStage stage1,stage2; 82 #endif 83 PetscMPIInt rank; 84 PetscInt eStart, eEnd, vStart, vEnd,j; 85 PetscInt genj,loadj; 86 Vec X,F; 87 Mat J; 88 SNES snes; 89 90 CHKERRQ(PetscInitialize(&argc,&argv,"poweroptions",help)); 91 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 92 { 93 /* 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 */ 94 /* this is an experiment to see how the analyzer reacts */ 95 const PetscMPIInt crank = rank; 96 97 /* Create an empty network object */ 98 CHKERRQ(DMNetworkCreate(PETSC_COMM_WORLD,&networkdm)); 99 /* Register the components in the network */ 100 CHKERRQ(DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch)); 101 CHKERRQ(DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus)); 102 CHKERRQ(DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen)); 103 CHKERRQ(DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load)); 104 105 CHKERRQ(PetscLogStageRegister("Read Data",&stage1)); 106 PetscLogStagePush(stage1); 107 /* READ THE DATA */ 108 if (!crank) { 109 /* READ DATA */ 110 /* Only rank 0 reads the data */ 111 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL)); 112 CHKERRQ(PetscNew(&pfdata)); 113 CHKERRQ(PFReadMatPowerData(pfdata,pfdata_file)); 114 User.Sbase = pfdata->sbase; 115 116 numEdges = pfdata->nbranch; 117 CHKERRQ(PetscMalloc1(2*numEdges,&edges)); 118 CHKERRQ(GetListofEdges_Power(pfdata,edges)); 119 } 120 121 /* If external option activated. Introduce error in jacobian */ 122 CHKERRQ(PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error)); 123 124 PetscLogStagePop(); 125 CHKERRMPI(MPI_Barrier(PETSC_COMM_WORLD)); 126 CHKERRQ(PetscLogStageRegister("Create network",&stage2)); 127 PetscLogStagePush(stage2); 128 /* Set number of nodes/edges */ 129 CHKERRQ(DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1)); 130 CHKERRQ(DMNetworkAddSubnetwork(networkdm,"",numEdges,edges,NULL)); 131 132 /* Set up the network layout */ 133 CHKERRQ(DMNetworkLayoutSetUp(networkdm)); 134 135 if (!crank) { 136 CHKERRQ(PetscFree(edges)); 137 } 138 139 /* Add network components only process 0 has any data to add */ 140 if (!crank) { 141 genj=0; loadj=0; 142 CHKERRQ(DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd)); 143 for (i = eStart; i < eEnd; i++) { 144 CHKERRQ(DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart],0)); 145 } 146 CHKERRQ(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 147 for (i = vStart; i < vEnd; i++) { 148 CHKERRQ(DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart],2)); 149 if (pfdata->bus[i-vStart].ngen) { 150 for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) { 151 CHKERRQ(DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++],0)); 152 } 153 } 154 if (pfdata->bus[i-vStart].nload) { 155 for (j=0; j < pfdata->bus[i-vStart].nload; j++) { 156 CHKERRQ(DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++],0)); 157 } 158 } 159 } 160 } 161 162 /* Set up DM for use */ 163 CHKERRQ(DMSetUp(networkdm)); 164 165 if (!crank) { 166 CHKERRQ(PetscFree(pfdata->bus)); 167 CHKERRQ(PetscFree(pfdata->gen)); 168 CHKERRQ(PetscFree(pfdata->branch)); 169 CHKERRQ(PetscFree(pfdata->load)); 170 CHKERRQ(PetscFree(pfdata)); 171 } 172 173 /* Distribute networkdm to multiple processes */ 174 CHKERRQ(DMNetworkDistribute(&networkdm,0)); 175 176 PetscLogStagePop(); 177 CHKERRQ(DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd)); 178 CHKERRQ(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 179 180 #if 0 181 EDGE_Power edge; 182 PetscInt key,kk,numComponents; 183 VERTEX_Power bus; 184 GEN gen; 185 LOAD load; 186 187 for (i = eStart; i < eEnd; i++) { 188 CHKERRQ(DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge)); 189 CHKERRQ(DMNetworkGetNumComponents(networkdm,i,&numComponents)); 190 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j)); 191 } 192 193 for (i = vStart; i < vEnd; i++) { 194 CHKERRQ(DMNetworkGetNumComponents(networkdm,i,&numComponents)); 195 for (kk=0; kk < numComponents; kk++) { 196 CHKERRQ(DMNetworkGetComponent(networkdm,i,kk,&key,&component)); 197 if (key == 1) { 198 bus = (VERTEX_Power)(component); 199 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i)); 200 } else if (key == 2) { 201 gen = (GEN)(component); 202 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg)); 203 } else if (key == 3) { 204 load = (LOAD)(component); 205 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql)); 206 } 207 } 208 } 209 #endif 210 /* Broadcast Sbase to all processors */ 211 CHKERRMPI(MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD)); 212 213 CHKERRQ(DMCreateGlobalVector(networkdm,&X)); 214 CHKERRQ(VecDuplicate(X,&F)); 215 216 CHKERRQ(DMCreateMatrix(networkdm,&J)); 217 CHKERRQ(MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 218 219 CHKERRQ(SetInitialValues(networkdm,X,&User)); 220 221 /* HOOK UP SOLVER */ 222 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 223 CHKERRQ(SNESSetDM(snes,networkdm)); 224 CHKERRQ(SNESSetFunction(snes,F,FormFunction,&User)); 225 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian_Power,&User)); 226 CHKERRQ(SNESSetFromOptions(snes)); 227 228 CHKERRQ(SNESSolve(snes,NULL,X)); 229 /* CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 230 231 CHKERRQ(VecDestroy(&X)); 232 CHKERRQ(VecDestroy(&F)); 233 CHKERRQ(MatDestroy(&J)); 234 235 CHKERRQ(SNESDestroy(&snes)); 236 CHKERRQ(DMDestroy(&networkdm)); 237 } 238 CHKERRQ(PetscFinalize()); 239 return 0; 240 } 241 242 /*TEST 243 244 build: 245 depends: PFReadData.c pffunctions.c 246 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 247 248 test: 249 args: -snes_rtol 1.e-3 250 localrunfiles: poweroptions case9.m 251 output_file: output/power_1.out 252 253 test: 254 suffix: 2 255 args: -snes_rtol 1.e-3 -petscpartitioner_type simple 256 nsize: 4 257 localrunfiles: poweroptions case9.m 258 output_file: output/power_1.out 259 260 TEST*/ 261