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:   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));

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

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

 33:   PetscCall(DMRestoreLocalVector(networkdm, &localX));

 35:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
 36:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
 37:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 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;

 48:   PetscFunctionBegin;
 49:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

 51:   PetscCall(DMGetLocalVector(networkdm, &localX));

 53:   PetscCall(VecSet(X, 0.0));
 54:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
 55:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

 57:   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
 58:   PetscCall(SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, user_power));

 60:   PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
 61:   PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
 62:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 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:   PetscLogStage stage1, stage2;
 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:   PetscFunctionBeginUser;
 84:   PetscCall(PetscInitialize(&argc, &argv, "poweroptions", help));
 85:   PetscCallMPI(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:     PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
 93:     /* Register the components in the network */
 94:     PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &User.compkey_branch));
 95:     PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &User.compkey_bus));
 96:     PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &User.compkey_gen));
 97:     PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &User.compkey_load));

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

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

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

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

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

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

131:     /* Add network components only process 0 has any data to add */
132:     if (!crank) {
133:       genj  = 0;
134:       loadj = 0;
135:       PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
136:       for (i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_branch, &pfdata->branch[i - eStart], 0));
137:       PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
138:       for (i = vStart; i < vEnd; i++) {
139:         PetscCall(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++) PetscCall(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++) PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_load, &pfdata->load[loadj++], 0));
145:         }
146:       }
147:     }

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

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

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

163:     PetscCall(PetscLogStagePop());
164:     PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
165:     PetscCall(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:       PetscCall(DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge));
176:       PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents));
177:       PetscCall(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:       PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents));
182:       for (kk=0; kk < numComponents; kk++) {
183:         PetscCall(DMNetworkGetComponent(networkdm,i,kk,&key,&component));
184:         if (key == 1) {
185:           bus = (VERTEX_Power)(component);
186:           PetscCall(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:           PetscCall(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:           PetscCall(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:     PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));

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

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

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

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

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

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

222:     PetscCall(SNESDestroy(&snes));
223:     PetscCall(DMDestroy(&networkdm));
224:   }
225:   PetscCall(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*/