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: DM networkdm;
73: UserCtx_Power User;
74: PetscLogStage stage1, stage2;
75: PetscMPIInt rank;
76: PetscInt eStart, eEnd, vStart, vEnd, j;
77: PetscInt genj, loadj;
78: Vec X, F;
79: Mat J;
80: SNES snes;
82: PetscFunctionBeginUser;
83: PetscCall(PetscInitialize(&argc, &argv, "poweroptions", help));
84: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
85: {
86: /* 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 */
87: /* this is an experiment to see how the analyzer reacts */
88: const PetscMPIInt crank = rank;
90: /* Create an empty network object */
91: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
92: /* Register the components in the network */
93: PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &User.compkey_branch));
94: PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &User.compkey_bus));
95: PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &User.compkey_gen));
96: PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &User.compkey_load));
98: PetscCall(PetscLogStageRegister("Read Data", &stage1));
99: PetscCall(PetscLogStagePush(stage1));
100: /* READ THE DATA */
101: if (!crank) {
102: /* READ DATA */
103: /* Only rank 0 reads the data */
104: PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL));
105: PetscCall(PetscNew(&pfdata));
106: PetscCall(PFReadMatPowerData(pfdata, pfdata_file));
107: User.Sbase = pfdata->sbase;
109: numEdges = pfdata->nbranch;
110: PetscCall(PetscMalloc1(2 * numEdges, &edges));
111: PetscCall(GetListofEdges_Power(pfdata, edges));
112: }
114: /* If external option activated. Introduce error in jacobian */
115: PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &User.jac_error));
117: PetscCall(PetscLogStagePop());
118: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
119: PetscCall(PetscLogStageRegister("Create network", &stage2));
120: PetscCall(PetscLogStagePush(stage2));
121: /* Set number of nodes/edges */
122: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
123: PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges, edges, NULL));
125: /* Set up the network layout */
126: PetscCall(DMNetworkLayoutSetUp(networkdm));
128: if (!crank) PetscCall(PetscFree(edges));
130: /* Add network components only process 0 has any data to add */
131: if (!crank) {
132: genj = 0;
133: loadj = 0;
134: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
135: for (PetscInt i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_branch, &pfdata->branch[i - eStart], 0));
136: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
137: for (PetscInt i = vStart; i < vEnd; i++) {
138: PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_bus, &pfdata->bus[i - vStart], 2));
139: if (pfdata->bus[i - vStart].ngen) {
140: for (j = 0; j < pfdata->bus[i - vStart].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_gen, &pfdata->gen[genj++], 0));
141: }
142: if (pfdata->bus[i - vStart].nload) {
143: for (j = 0; j < pfdata->bus[i - vStart].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, i, User.compkey_load, &pfdata->load[loadj++], 0));
144: }
145: }
146: }
148: /* Set up DM for use */
149: PetscCall(DMSetUp(networkdm));
151: if (!crank) {
152: PetscCall(PetscFree(pfdata->bus));
153: PetscCall(PetscFree(pfdata->gen));
154: PetscCall(PetscFree(pfdata->branch));
155: PetscCall(PetscFree(pfdata->load));
156: PetscCall(PetscFree(pfdata));
157: }
159: /* Distribute networkdm to multiple processes */
160: PetscCall(DMNetworkDistribute(&networkdm, 0));
162: PetscCall(PetscLogStagePop());
163: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
164: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
166: #if 0
167: EDGE_Power edge;
168: PetscInt key,kk,numComponents;
169: VERTEX_Power bus;
170: GEN gen;
171: LOAD load;
173: for (PetscInt i = eStart; i < eEnd; i++) {
174: PetscCall(DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge));
175: PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents));
176: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j));
177: }
179: for (PetscInt i = vStart; i < vEnd; i++) {
180: PetscCall(DMNetworkGetNumComponents(networkdm,i,&numComponents));
181: for (kk=0; kk < numComponents; kk++) {
182: PetscCall(DMNetworkGetComponent(networkdm,i,kk,&key,&component));
183: if (key == 1) {
184: bus = (VERTEX_Power)component;
185: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i));
186: } else if (key == 2) {
187: gen = (GEN)component;
188: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,(double)gen->pg,(double)gen->qg));
189: } else if (key == 3) {
190: load = (LOAD)component;
191: PetscCall(PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,(double)load->pl,(double)load->ql));
192: }
193: }
194: }
195: #endif
196: /* Broadcast Sbase to all processors */
197: PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
199: PetscCall(DMCreateGlobalVector(networkdm, &X));
200: PetscCall(VecDuplicate(X, &F));
202: PetscCall(DMCreateMatrix(networkdm, &J));
203: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
205: PetscCall(SetInitialValues(networkdm, X, &User));
207: /* HOOK UP SOLVER */
208: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
209: PetscCall(SNESSetDM(snes, networkdm));
210: PetscCall(SNESSetFunction(snes, F, FormFunction, &User));
211: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian_Power, &User));
212: PetscCall(SNESSetFromOptions(snes));
214: PetscCall(SNESSolve(snes, NULL, X));
215: /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
217: PetscCall(VecDestroy(&X));
218: PetscCall(VecDestroy(&F));
219: PetscCall(MatDestroy(&J));
221: PetscCall(SNESDestroy(&snes));
222: PetscCall(DMDestroy(&networkdm));
223: }
224: PetscCall(PetscFinalize());
225: return 0;
226: }
228: /*TEST
230: build:
231: depends: PFReadData.c pffunctions.c
232: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
234: test:
235: args: -snes_rtol 1.e-3
236: localrunfiles: poweroptions case9.m
237: output_file: output/power_1.out
239: test:
240: suffix: 2
241: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
242: nsize: 4
243: localrunfiles: poweroptions case9.m
244: output_file: output/power_1.out
246: TEST*/