Actual source code: power.c

petsc-3.14.3 2021-01-09
Report Typos and Errors
  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: /* T
 10:    Concepts: DMNetwork
 11:    Concepts: PETSc SNES solver
 12: */

 14: #include "power.h"
 15: #include <petscdmnetwork.h>

 17: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
 18: {
 20:   DM             networkdm;
 21:   UserCtx_Power  *User=(UserCtx_Power*)appctx;
 22:   Vec            localX,localF;
 23:   PetscInt       nv,ne;
 24:   const PetscInt *vtx,*edges;

 27:   SNESGetDM(snes,&networkdm);
 28:   DMGetLocalVector(networkdm,&localX);
 29:   DMGetLocalVector(networkdm,&localF);
 30:   VecSet(F,0.0);
 31:   VecSet(localF,0.0);

 33:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 34:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 36:   DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
 37:   FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,User);

 39:   DMRestoreLocalVector(networkdm,&localX);

 41:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
 42:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
 43:   DMRestoreLocalVector(networkdm,&localF);
 44:   return(0);
 45: }

 47: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
 48: {
 50:   PetscInt       vStart,vEnd,nv,ne;
 51:   const PetscInt *vtx,*edges;
 52:   Vec            localX;
 53:   UserCtx_Power  *user_power=(UserCtx_Power*)appctx;

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

 58:   DMGetLocalVector(networkdm,&localX);

 60:   VecSet(X,0.0);
 61:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 62:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 64:   DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
 65:   SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,user_power);

 67:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
 68:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
 69:   DMRestoreLocalVector(networkdm,&localX);
 70:   return(0);
 71: }

 73: int main(int argc,char ** argv)
 74: {
 75:   PetscErrorCode   ierr;
 76:   char             pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
 77:   PFDATA           *pfdata;
 78:   PetscInt         numEdges=0,numVertices=0;
 79:   PetscInt         *edges = NULL;
 80:   PetscInt         i;
 81:   DM               networkdm;
 82:   UserCtx_Power    User;
 83: #if defined(PETSC_USE_LOG)
 84:   PetscLogStage    stage1,stage2;
 85: #endif
 86:   PetscMPIInt      rank;
 87:   PetscInt         eStart, eEnd, vStart, vEnd,j;
 88:   PetscInt         genj,loadj;
 89:   Vec              X,F;
 90:   Mat              J;
 91:   SNES             snes;

 93:   PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr;
 94:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 95:   {
 96:     /* 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 */
 97:     /* this is an experiment to see how the analyzer reacts */
 98:     const PetscMPIInt crank = rank;

100:     /* Create an empty network object */
101:     DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
102:     /* Register the components in the network */
103:     DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);
104:     DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);
105:     DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);
106:     DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);

108:     PetscLogStageRegister("Read Data",&stage1);
109:     PetscLogStagePush(stage1);
110:     /* READ THE DATA */
111:     if (!crank) {
112:       /*    READ DATA */
113:       /* Only rank 0 reads the data */
114:       PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
115:       PetscNew(&pfdata);
116:       PFReadMatPowerData(pfdata,pfdata_file);
117:       User.Sbase = pfdata->sbase;

119:       numEdges = pfdata->nbranch;
120:       numVertices = pfdata->nbus;

122:       PetscMalloc1(2*numEdges,&edges);
123:       GetListofEdges_Power(pfdata,edges);
124:     }

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

129:     PetscLogStagePop();
130:     MPI_Barrier(PETSC_COMM_WORLD);
131:     PetscLogStageRegister("Create network",&stage2);
132:     PetscLogStagePush(stage2);
133:     /* Set number of nodes/edges */
134:     DMNetworkSetSizes(networkdm,1,&numVertices,&numEdges,0,NULL);
135:     /* Add edge connectivity */
136:     DMNetworkSetEdgeList(networkdm,&edges,NULL);
137:     /* Set up the network layout */
138:     DMNetworkLayoutSetUp(networkdm);

140:     if (!crank) {
141:       PetscFree(edges);
142:     }

144:     /* Add network components only process 0 has any data to add*/
145:     if (!crank) {
146:       genj=0; loadj=0;
147:       DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
148:       for (i = eStart; i < eEnd; i++) {
149:         DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart]);
150:       }
151:       DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
152:       for (i = vStart; i < vEnd; i++) {
153:         DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart]);
154:         if (pfdata->bus[i-vStart].ngen) {
155:           for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) {
156:             DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++]);
157:           }
158:         }
159:         if (pfdata->bus[i-vStart].nload) {
160:           for (j=0; j < pfdata->bus[i-vStart].nload; j++) {
161:             DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++]);
162:           }
163:         }
164:         /* Add number of variables */
165:         DMNetworkAddNumVariables(networkdm,i,2);
166:       }
167:     }

169:     /* Set up DM for use */
170:     DMSetUp(networkdm);

172:     if (!crank) {
173:       PetscFree(pfdata->bus);
174:       PetscFree(pfdata->gen);
175:       PetscFree(pfdata->branch);
176:       PetscFree(pfdata->load);
177:       PetscFree(pfdata);
178:     }

180:     /* Distribute networkdm to multiple processes */
181:     DMNetworkDistribute(&networkdm,0);

183:     PetscLogStagePop();
184:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
185:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

187: #if 0
188:     EDGE_Power     edge;
189:     PetscInt       key,kk,numComponents;
190:     VERTEX_Power   bus;
191:     GEN            gen;
192:     LOAD           load;

194:     for (i = eStart; i < eEnd; i++) {
195:       DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
196:       DMNetworkGetNumComponents(networkdm,i,&numComponents);
197:       PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
198:     }

200:     for (i = vStart; i < vEnd; i++) {
201:       DMNetworkGetNumComponents(networkdm,i,&numComponents);
202:       for (kk=0; kk < numComponents; kk++) {
203:         DMNetworkGetComponent(networkdm,i,kk,&key,&component);
204:         if (key == 1) {
205:           bus = (VERTEX_Power)(component);
206:           PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
207:         } else if (key == 2) {
208:           gen = (GEN)(component);
209:           PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
210:         } else if (key == 3) {
211:           load = (LOAD)(component);
212:           PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
213:         }
214:       }
215:     }
216: #endif
217:     /* Broadcast Sbase to all processors */
218:     MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);

220:     DMCreateGlobalVector(networkdm,&X);
221:     VecDuplicate(X,&F);

223:     DMCreateMatrix(networkdm,&J);
224:     MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

228:     /* HOOK UP SOLVER */
229:     SNESCreate(PETSC_COMM_WORLD,&snes);
230:     SNESSetDM(snes,networkdm);
231:     SNESSetFunction(snes,F,FormFunction,&User);
232:     SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);
233:     SNESSetFromOptions(snes);

235:     SNESSolve(snes,NULL,X);
236:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

238:     VecDestroy(&X);
239:     VecDestroy(&F);
240:     MatDestroy(&J);

242:     SNESDestroy(&snes);
243:     DMDestroy(&networkdm);
244:   }
245:   PetscFinalize();
246:   return ierr;
247: }

249: /*TEST

251:    build:
252:      depends: PFReadData.c pffunctions.c
253:      requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)


256:    test:
257:      args: -snes_rtol 1.e-3
258:      localrunfiles: poweroptions case9.m
259:      output_file: output/power_1.out
260:      requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)

262:    test:
263:      suffix: 2
264:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
265:      nsize: 4
266:      localrunfiles: poweroptions case9.m
267:      output_file: output/power_1.out
268:      requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)

270: TEST*/