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"; 8 9 #include "power.h" 10 #include <petscdmnetwork.h> 11 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; 19 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)); 26 27 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 28 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 29 30 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 31 PetscCall(FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, User)); 32 33 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 34 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 } 40 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; 47 48 PetscFunctionBegin; 49 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd)); 50 51 PetscCall(DMGetLocalVector(networkdm, &localX)); 52 53 PetscCall(VecSet(X, 0.0)); 54 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 55 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 56 57 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges)); 58 PetscCall(SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, user_power)); 59 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 } 65 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; 82 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; 90 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)); 98 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; 109 110 numEdges = pfdata->nbranch; 111 PetscCall(PetscMalloc1(2 * numEdges, &edges)); 112 PetscCall(GetListofEdges_Power(pfdata, edges)); 113 } 114 115 /* If external option activated. Introduce error in jacobian */ 116 PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &User.jac_error)); 117 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)); 125 126 /* Set up the network layout */ 127 PetscCall(DMNetworkLayoutSetUp(networkdm)); 128 129 if (!crank) PetscCall(PetscFree(edges)); 130 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 } 148 149 /* Set up DM for use */ 150 PetscCall(DMSetUp(networkdm)); 151 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 } 159 160 /* Distribute networkdm to multiple processes */ 161 PetscCall(DMNetworkDistribute(&networkdm, 0)); 162 163 PetscCall(PetscLogStagePop()); 164 PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd)); 165 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd)); 166 167 #if 0 168 EDGE_Power edge; 169 PetscInt key,kk,numComponents; 170 VERTEX_Power bus; 171 GEN gen; 172 LOAD load; 173 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 } 179 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)); 199 200 PetscCall(DMCreateGlobalVector(networkdm, &X)); 201 PetscCall(VecDuplicate(X, &F)); 202 203 PetscCall(DMCreateMatrix(networkdm, &J)); 204 PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 205 206 PetscCall(SetInitialValues(networkdm, X, &User)); 207 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)); 214 215 PetscCall(SNESSolve(snes, NULL, X)); 216 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 217 218 PetscCall(VecDestroy(&X)); 219 PetscCall(VecDestroy(&F)); 220 PetscCall(MatDestroy(&J)); 221 222 PetscCall(SNESDestroy(&snes)); 223 PetscCall(DMDestroy(&networkdm)); 224 } 225 PetscCall(PetscFinalize()); 226 return 0; 227 } 228 229 /*TEST 230 231 build: 232 depends: PFReadData.c pffunctions.c 233 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED) 234 235 test: 236 args: -snes_rtol 1.e-3 237 localrunfiles: poweroptions case9.m 238 output_file: output/power_1.out 239 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 246 247 TEST*/ 248