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: SNESGetDM(snes, &networkdm);
21: DMGetLocalVector(networkdm, &localX);
22: DMGetLocalVector(networkdm, &localF);
23: VecSet(F, 0.0);
24: VecSet(localF, 0.0);
26: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
27: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
29: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
30: FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, User);
32: DMRestoreLocalVector(networkdm, &localX);
34: DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
35: DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
36: DMRestoreLocalVector(networkdm, &localF);
37: return 0;
38: }
40: PetscErrorCode SetInitialValues(DM networkdm, Vec X, void *appctx)
41: {
42: PetscInt vStart, vEnd, nv, ne;
43: const PetscInt *vtx, *edges;
44: Vec localX;
45: UserCtx_Power *user_power = (UserCtx_Power *)appctx;
47: DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);
49: DMGetLocalVector(networkdm, &localX);
51: VecSet(X, 0.0);
52: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
53: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
55: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
56: SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, user_power);
58: DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X);
59: DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X);
60: DMRestoreLocalVector(networkdm, &localX);
61: return 0;
62: }
64: int main(int argc, char **argv)
65: {
66: char pfdata_file[PETSC_MAX_PATH_LEN] = "case9.m";
67: PFDATA *pfdata;
68: PetscInt numEdges = 0;
69: PetscInt *edges = NULL;
70: PetscInt i;
71: DM networkdm;
72: UserCtx_Power User;
73: #if defined(PETSC_USE_LOG)
74: PetscLogStage stage1, stage2;
75: #endif
76: PetscMPIInt rank;
77: PetscInt eStart, eEnd, vStart, vEnd, j;
78: PetscInt genj, loadj;
79: Vec X, F;
80: Mat J;
81: SNES snes;
84: PetscInitialize(&argc, &argv, "poweroptions", help);
85: 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: DMNetworkCreate(PETSC_COMM_WORLD, &networkdm);
93: /* Register the components in the network */
94: DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &User.compkey_branch);
95: DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &User.compkey_bus);
96: DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &User.compkey_gen);
97: DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &User.compkey_load);
99: PetscLogStageRegister("Read Data", &stage1);
100: PetscLogStagePush(stage1);
101: /* READ THE DATA */
102: if (!crank) {
103: /* READ DATA */
104: /* Only rank 0 reads the data */
105: PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL);
106: PetscNew(&pfdata);
107: PFReadMatPowerData(pfdata, pfdata_file);
108: User.Sbase = pfdata->sbase;
110: numEdges = pfdata->nbranch;
111: PetscMalloc1(2 * numEdges, &edges);
112: GetListofEdges_Power(pfdata, edges);
113: }
115: /* If external option activated. Introduce error in jacobian */
116: PetscOptionsHasName(NULL, NULL, "-jac_error", &User.jac_error);
118: PetscLogStagePop();
119: MPI_Barrier(PETSC_COMM_WORLD);
120: PetscLogStageRegister("Create network", &stage2);
121: PetscLogStagePush(stage2);
122: /* Set number of nodes/edges */
123: DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1);
124: DMNetworkAddSubnetwork(networkdm, "", numEdges, edges, NULL);
126: /* Set up the network layout */
127: DMNetworkLayoutSetUp(networkdm);
129: if (!crank) PetscFree(edges);
131: /* Add network components only process 0 has any data to add */
132: if (!crank) {
133: genj = 0;
134: loadj = 0;
135: DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
136: for (i = eStart; i < eEnd; i++) DMNetworkAddComponent(networkdm, i, User.compkey_branch, &pfdata->branch[i - eStart], 0);
137: DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);
138: for (i = vStart; i < vEnd; i++) {
139: 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++) 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++) DMNetworkAddComponent(networkdm, i, User.compkey_load, &pfdata->load[loadj++], 0);
145: }
146: }
147: }
149: /* Set up DM for use */
150: DMSetUp(networkdm);
152: if (!crank) {
153: PetscFree(pfdata->bus);
154: PetscFree(pfdata->gen);
155: PetscFree(pfdata->branch);
156: PetscFree(pfdata->load);
157: PetscFree(pfdata);
158: }
160: /* Distribute networkdm to multiple processes */
161: DMNetworkDistribute(&networkdm, 0);
163: PetscLogStagePop();
164: DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
165: 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: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
176: DMNetworkGetNumComponents(networkdm,i,&numComponents);
177: 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: DMNetworkGetNumComponents(networkdm,i,&numComponents);
182: for (kk=0; kk < numComponents; kk++) {
183: DMNetworkGetComponent(networkdm,i,kk,&key,&component);
184: if (key == 1) {
185: bus = (VERTEX_Power)(component);
186: 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: 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: 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: MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD);
200: DMCreateGlobalVector(networkdm, &X);
201: VecDuplicate(X, &F);
203: DMCreateMatrix(networkdm, &J);
204: MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
206: SetInitialValues(networkdm, X, &User);
208: /* HOOK UP SOLVER */
209: SNESCreate(PETSC_COMM_WORLD, &snes);
210: SNESSetDM(snes, networkdm);
211: SNESSetFunction(snes, F, FormFunction, &User);
212: SNESSetJacobian(snes, J, J, FormJacobian_Power, &User);
213: SNESSetFromOptions(snes);
215: SNESSolve(snes, NULL, X);
216: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
218: VecDestroy(&X);
219: VecDestroy(&F);
220: MatDestroy(&J);
222: SNESDestroy(&snes);
223: DMDestroy(&networkdm);
224: }
225: 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*/