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;

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

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

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

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

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

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

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

129:     if (!crank) PetscFree(edges);

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

149:     /* Set up DM for use */
150:     DMSetUp(networkdm);

152:     if (!crank) {
153:       PetscFree(pfdata->bus);
154:       PetscFree(pfdata->gen);
155:       PetscFree(pfdata->branch);
156:       PetscFree(pfdata->load);
157:       PetscFree(pfdata);
158:     }

160:     /* Distribute networkdm to multiple processes */
161:     DMNetworkDistribute(&networkdm, 0);

163:     PetscLogStagePop();
164:     DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
165:     DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);

167: #if 0
168:     EDGE_Power     edge;
169:     PetscInt       key,kk,numComponents;
170:     VERTEX_Power   bus;
171:     GEN            gen;
172:     LOAD           load;

174:     for (i = eStart; i < eEnd; i++) {
175:       DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
176:       DMNetworkGetNumComponents(networkdm,i,&numComponents);
177:       PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
178:     }

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

200:     DMCreateGlobalVector(networkdm, &X);
201:     VecDuplicate(X, &F);

203:     DMCreateMatrix(networkdm, &J);
204:     MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

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

208:     /* HOOK UP SOLVER */
209:     SNESCreate(PETSC_COMM_WORLD, &snes);
210:     SNESSetDM(snes, networkdm);
211:     SNESSetFunction(snes, F, FormFunction, &User);
212:     SNESSetJacobian(snes, J, J, FormJacobian_Power, &User);
213:     SNESSetFromOptions(snes);

215:     SNESSolve(snes, NULL, X);
216:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

218:     VecDestroy(&X);
219:     VecDestroy(&F);
220:     MatDestroy(&J);

222:     SNESDestroy(&snes);
223:     DMDestroy(&networkdm);
224:   }
225:   PetscFinalize();
226:   return 0;
227: }

229: /*TEST

231:    build:
232:      depends: PFReadData.c pffunctions.c
233:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

235:    test:
236:      args: -snes_rtol 1.e-3
237:      localrunfiles: poweroptions case9.m
238:      output_file: output/power_1.out

240:    test:
241:      suffix: 2
242:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
243:      nsize: 4
244:      localrunfiles: poweroptions case9.m
245:      output_file: output/power_1.out

247: TEST*/