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

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

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

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

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

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

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

128:     /* Set up the network layout */
129:     PetscCall(DMNetworkLayoutSetUp(networkdm));

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

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

151:     /* Set up DM for use */
152:     PetscCall(DMSetUp(networkdm));

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

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

165:     PetscCall(PetscLogStagePop());
166:     PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
167:     PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

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

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

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

202:     PetscCall(DMCreateGlobalVector(networkdm, &X));
203:     PetscCall(VecDuplicate(X, &F));

205:     PetscCall(DMCreateMatrix(networkdm, &J));
206:     PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

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

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

217:     PetscCall(SNESSolve(snes, NULL, X));
218:     /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */

220:     PetscCall(VecDestroy(&X));
221:     PetscCall(VecDestroy(&F));
222:     PetscCall(MatDestroy(&J));

224:     PetscCall(SNESDestroy(&snes));
225:     PetscCall(DMDestroy(&networkdm));
226:   }
227:   PetscCall(PetscFinalize());
228:   return 0;
229: }

231: /*TEST

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

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

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

249: TEST*/