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: 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;
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;
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));
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;
110: numEdges = pfdata->nbranch;
111: PetscCall(PetscMalloc1(2 * numEdges, &edges));
112: PetscCall(GetListofEdges_Power(pfdata, edges));
113: }
115: /* If external option activated. Introduce error in jacobian */
116: PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &User.jac_error));
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));
126: /* Set up the network layout */
127: PetscCall(DMNetworkLayoutSetUp(networkdm));
129: if (!crank) PetscCall(PetscFree(edges));
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: }
149: /* Set up DM for use */
150: PetscCall(DMSetUp(networkdm));
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: }
160: /* Distribute networkdm to multiple processes */
161: PetscCall(DMNetworkDistribute(&networkdm, 0));
163: PetscCall(PetscLogStagePop());
164: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
165: PetscCall(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: 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: }
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));
200: PetscCall(DMCreateGlobalVector(networkdm, &X));
201: PetscCall(VecDuplicate(X, &F));
203: PetscCall(DMCreateMatrix(networkdm, &J));
204: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
206: PetscCall(SetInitialValues(networkdm, X, &User));
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));
215: PetscCall(SNESSolve(snes, NULL, X));
216: /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
218: PetscCall(VecDestroy(&X));
219: PetscCall(VecDestroy(&F));
220: PetscCall(MatDestroy(&J));
222: PetscCall(SNESDestroy(&snes));
223: PetscCall(DMDestroy(&networkdm));
224: }
225: PetscCall(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*/