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