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
FormFunction(SNES snes,Vec X,Vec F,void * appctx)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
SetInitialValues(DM networkdm,Vec X,void * appctx)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
main(int argc,char ** argv)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