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: /* 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:   DMNetworkGetSubnetwork(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:   DMNetworkGetSubnetwork(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:     DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
    135:     DMNetworkAddSubnetwork(networkdm,"",numVertices,numEdges,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],0);
    150:       }
    151:       DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    152:       for (i = vStart; i < vEnd; i++) {
    153:         DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart],2);
    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++],0);
    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++],0);
    162:           }
    163:         }
    164:       }
    165:     }

    167:     /* Set up DM for use */
    168:     DMSetUp(networkdm);

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

    178:     /* Distribute networkdm to multiple processes */
    179:     DMNetworkDistribute(&networkdm,0);

    181:     PetscLogStagePop();
    182:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    183:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

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

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

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

    218:     DMCreateGlobalVector(networkdm,&X);
    219:     VecDuplicate(X,&F);

    221:     DMCreateMatrix(networkdm,&J);
    222:     MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

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

    233:     SNESSolve(snes,NULL,X);
    234:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

    236:     VecDestroy(&X);
    237:     VecDestroy(&F);
    238:     MatDestroy(&J);

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

    247: /*TEST

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

    253:    test:
    254:      args: -snes_rtol 1.e-3
    255:      localrunfiles: poweroptions case9.m
    256:      output_file: output/power_1.out

    258:    test:
    259:      suffix: 2
    260:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
    261:      nsize: 4
    262:      localrunfiles: poweroptions case9.m
    263:      output_file: output/power_1.out

    265: TEST*/