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*/