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 PetscCall(SNESGetDM(snes,&networkdm)); 27 PetscCall(DMGetLocalVector(networkdm,&localX)); 28 PetscCall(DMGetLocalVector(networkdm,&localF)); 29 PetscCall(VecSet(F,0.0)); 30 PetscCall(VecSet(localF,0.0)); 31 32 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 33 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 34 35 PetscCall(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 36 PetscCall(FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,User)); 37 38 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 39 40 PetscCall(DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F)); 41 PetscCall(DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F)); 42 PetscCall(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 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart, &vEnd)); 55 56 PetscCall(DMGetLocalVector(networkdm,&localX)); 57 58 PetscCall(VecSet(X,0.0)); 59 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 60 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 61 62 PetscCall(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 63 PetscCall(SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,user_power)); 64 65 PetscCall(DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X)); 66 PetscCall(DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X)); 67 PetscCall(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 PetscCall(PetscInitialize(&argc,&argv,"poweroptions",help)); 91 PetscCallMPI(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 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD,&networkdm)); 99 /* Register the components in the network */ 100 PetscCall(DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch)); 101 PetscCall(DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus)); 102 PetscCall(DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen)); 103 PetscCall(DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load)); 104 105 PetscCall(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 PetscCall(PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL)); 112 PetscCall(PetscNew(&pfdata)); 113 PetscCall(PFReadMatPowerData(pfdata,pfdata_file)); 114 User.Sbase = pfdata->sbase; 115 116 numEdges = pfdata->nbranch; 117 PetscCall(PetscMalloc1(2*numEdges,&edges)); 118 PetscCall(GetListofEdges_Power(pfdata,edges)); 119 } 120 121 /* If external option activated. Introduce error in jacobian */ 122 PetscCall(PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error)); 123 124 PetscLogStagePop(); 125 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 126 PetscCall(PetscLogStageRegister("Create network",&stage2)); 127 PetscLogStagePush(stage2); 128 /* Set number of nodes/edges */ 129 PetscCall(DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1)); 130 PetscCall(DMNetworkAddSubnetwork(networkdm,"",numEdges,edges,NULL)); 131 132 /* Set up the network layout */ 133 PetscCall(DMNetworkLayoutSetUp(networkdm)); 134 135 if (!crank) { 136 PetscCall(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 PetscCall(DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd)); 143 for (i = eStart; i < eEnd; i++) { 144 PetscCall(DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart],0)); 145 } 146 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 147 for (i = vStart; i < vEnd; i++) { 148 PetscCall(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 PetscCall(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 PetscCall(DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++],0)); 157 } 158 } 159 } 160 } 161 162 /* Set up DM for use */ 163 PetscCall(DMSetUp(networkdm)); 164 165 if (!crank) { 166 PetscCall(PetscFree(pfdata->bus)); 167 PetscCall(PetscFree(pfdata->gen)); 168 PetscCall(PetscFree(pfdata->branch)); 169 PetscCall(PetscFree(pfdata->load)); 170 PetscCall(PetscFree(pfdata)); 171 } 172 173 /* Distribute networkdm to multiple processes */ 174 PetscCall(DMNetworkDistribute(&networkdm,0)); 175 176 PetscLogStagePop(); 177 PetscCall(DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd)); 178 PetscCall(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 PetscCall(DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge)); 189 PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents)); 190 PetscCall(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 PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents)); 195 for (kk=0; kk < numComponents; kk++) { 196 PetscCall(DMNetworkGetComponent(networkdm,i,kk,&key,&component)); 197 if (key == 1) { 198 bus = (VERTEX_Power)(component); 199 PetscCall(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 PetscCall(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 PetscCall(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 PetscCallMPI(MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD)); 212 213 PetscCall(DMCreateGlobalVector(networkdm,&X)); 214 PetscCall(VecDuplicate(X,&F)); 215 216 PetscCall(DMCreateMatrix(networkdm,&J)); 217 PetscCall(MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 218 219 PetscCall(SetInitialValues(networkdm,X,&User)); 220 221 /* HOOK UP SOLVER */ 222 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 223 PetscCall(SNESSetDM(snes,networkdm)); 224 PetscCall(SNESSetFunction(snes,F,FormFunction,&User)); 225 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian_Power,&User)); 226 PetscCall(SNESSetFromOptions(snes)); 227 228 PetscCall(SNESSolve(snes,NULL,X)); 229 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 230 231 PetscCall(VecDestroy(&X)); 232 PetscCall(VecDestroy(&F)); 233 PetscCall(MatDestroy(&J)); 234 235 PetscCall(SNESDestroy(&snes)); 236 PetscCall(DMDestroy(&networkdm)); 237 } 238 PetscCall(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