Actual source code: power.c

  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";

  9: #include "power.h"
 10: #include <petscdmnetwork.h>

 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;

 20:   SNESGetDM(snes,&networkdm);
 21:   DMGetLocalVector(networkdm,&localX);
 22:   DMGetLocalVector(networkdm,&localF);
 23:   VecSet(F,0.0);
 24:   VecSet(localF,0.0);

 26:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 27:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 29:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
 30:   FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,User);

 32:   DMRestoreLocalVector(networkdm,&localX);

 34:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
 35:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
 36:   DMRestoreLocalVector(networkdm,&localF);
 37:   return 0;
 38: }

 40: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
 41: {
 42:   PetscInt       vStart,vEnd,nv,ne;
 43:   const PetscInt *vtx,*edges;
 44:   Vec            localX;
 45:   UserCtx_Power  *user_power=(UserCtx_Power*)appctx;

 47:   DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);

 49:   DMGetLocalVector(networkdm,&localX);

 51:   VecSet(X,0.0);
 52:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 53:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 55:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
 56:   SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,user_power);

 58:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
 59:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
 60:   DMRestoreLocalVector(networkdm,&localX);
 61:   return 0;
 62: }

 64: int main(int argc,char ** argv)
 65: {
 66:   char             pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
 67:   PFDATA           *pfdata;
 68:   PetscInt         numEdges=0;
 69:   PetscInt         *edges = NULL;
 70:   PetscInt         i;
 71:   DM               networkdm;
 72:   UserCtx_Power    User;
 73: #if defined(PETSC_USE_LOG)
 74:   PetscLogStage    stage1,stage2;
 75: #endif
 76:   PetscMPIInt      rank;
 77:   PetscInt         eStart, eEnd, vStart, vEnd,j;
 78:   PetscInt         genj,loadj;
 79:   Vec              X,F;
 80:   Mat              J;
 81:   SNES             snes;

 83:   PetscInitialize(&argc,&argv,"poweroptions",help);
 84:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 85:   {
 86:     /* 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 */
 87:     /* this is an experiment to see how the analyzer reacts */
 88:     const PetscMPIInt crank = rank;

 90:     /* Create an empty network object */
 91:     DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
 92:     /* Register the components in the network */
 93:     DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);
 94:     DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);
 95:     DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);
 96:     DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);

 98:     PetscLogStageRegister("Read Data",&stage1);
 99:     PetscLogStagePush(stage1);
100:     /* READ THE DATA */
101:     if (!crank) {
102:       /*    READ DATA */
103:       /* Only rank 0 reads the data */
104:       PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
105:       PetscNew(&pfdata);
106:       PFReadMatPowerData(pfdata,pfdata_file);
107:       User.Sbase = pfdata->sbase;

109:       numEdges = pfdata->nbranch;
110:       PetscMalloc1(2*numEdges,&edges);
111:       GetListofEdges_Power(pfdata,edges);
112:     }

114:     /* If external option activated. Introduce error in jacobian */
115:     PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);

117:     PetscLogStagePop();
118:     MPI_Barrier(PETSC_COMM_WORLD);
119:     PetscLogStageRegister("Create network",&stage2);
120:     PetscLogStagePush(stage2);
121:     /* Set number of nodes/edges */
122:     DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
123:     DMNetworkAddSubnetwork(networkdm,"",numEdges,edges,NULL);

125:     /* Set up the network layout */
126:     DMNetworkLayoutSetUp(networkdm);

128:     if (!crank) {
129:       PetscFree(edges);
130:     }

132:     /* Add network components only process 0 has any data to add */
133:     if (!crank) {
134:       genj=0; loadj=0;
135:       DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
136:       for (i = eStart; i < eEnd; i++) {
137:         DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart],0);
138:       }
139:       DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
140:       for (i = vStart; i < vEnd; i++) {
141:         DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart],2);
142:         if (pfdata->bus[i-vStart].ngen) {
143:           for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) {
144:             DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++],0);
145:           }
146:         }
147:         if (pfdata->bus[i-vStart].nload) {
148:           for (j=0; j < pfdata->bus[i-vStart].nload; j++) {
149:             DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++],0);
150:           }
151:         }
152:       }
153:     }

155:     /* Set up DM for use */
156:     DMSetUp(networkdm);

158:     if (!crank) {
159:       PetscFree(pfdata->bus);
160:       PetscFree(pfdata->gen);
161:       PetscFree(pfdata->branch);
162:       PetscFree(pfdata->load);
163:       PetscFree(pfdata);
164:     }

166:     /* Distribute networkdm to multiple processes */
167:     DMNetworkDistribute(&networkdm,0);

169:     PetscLogStagePop();
170:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
171:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

173: #if 0
174:     EDGE_Power     edge;
175:     PetscInt       key,kk,numComponents;
176:     VERTEX_Power   bus;
177:     GEN            gen;
178:     LOAD           load;

180:     for (i = eStart; i < eEnd; i++) {
181:       DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
182:       DMNetworkGetNumComponents(networkdm,i,&numComponents);
183:       PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
184:     }

186:     for (i = vStart; i < vEnd; i++) {
187:       DMNetworkGetNumComponents(networkdm,i,&numComponents);
188:       for (kk=0; kk < numComponents; kk++) {
189:         DMNetworkGetComponent(networkdm,i,kk,&key,&component);
190:         if (key == 1) {
191:           bus = (VERTEX_Power)(component);
192:           PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
193:         } else if (key == 2) {
194:           gen = (GEN)(component);
195:           PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
196:         } else if (key == 3) {
197:           load = (LOAD)(component);
198:           PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
199:         }
200:       }
201:     }
202: #endif
203:     /* Broadcast Sbase to all processors */
204:     MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);

206:     DMCreateGlobalVector(networkdm,&X);
207:     VecDuplicate(X,&F);

209:     DMCreateMatrix(networkdm,&J);
210:     MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

212:     SetInitialValues(networkdm,X,&User);

214:     /* HOOK UP SOLVER */
215:     SNESCreate(PETSC_COMM_WORLD,&snes);
216:     SNESSetDM(snes,networkdm);
217:     SNESSetFunction(snes,F,FormFunction,&User);
218:     SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);
219:     SNESSetFromOptions(snes);

221:     SNESSolve(snes,NULL,X);
222:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

224:     VecDestroy(&X);
225:     VecDestroy(&F);
226:     MatDestroy(&J);

228:     SNESDestroy(&snes);
229:     DMDestroy(&networkdm);
230:   }
231:   PetscFinalize();
232:   return 0;
233: }

235: /*TEST

237:    build:
238:      depends: PFReadData.c pffunctions.c
239:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

241:    test:
242:      args: -snes_rtol 1.e-3
243:      localrunfiles: poweroptions case9.m
244:      output_file: output/power_1.out

246:    test:
247:      suffix: 2
248:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
249:      nsize: 4
250:      localrunfiles: poweroptions case9.m
251:      output_file: output/power_1.out

253: TEST*/