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: /* T
10: Concepts: DMNetwork
11: Concepts: PETSc SNES solver
12: */
14: #include "power.h"
15: #include <petscdmnetwork.h>
17: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
18: {
19: DM networkdm;
20: UserCtx_Power *User=(UserCtx_Power*)appctx;
21: Vec localX,localF;
22: PetscInt nv,ne;
23: const PetscInt *vtx,*edges;
25: SNESGetDM(snes,&networkdm);
26: DMGetLocalVector(networkdm,&localX);
27: DMGetLocalVector(networkdm,&localF);
28: VecSet(F,0.0);
29: VecSet(localF,0.0);
31: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
32: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
34: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
35: FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,User);
37: DMRestoreLocalVector(networkdm,&localX);
39: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
40: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
41: DMRestoreLocalVector(networkdm,&localF);
42: return 0;
43: }
45: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
46: {
47: PetscInt vStart,vEnd,nv,ne;
48: const PetscInt *vtx,*edges;
49: Vec localX;
50: UserCtx_Power *user_power=(UserCtx_Power*)appctx;
52: DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);
54: DMGetLocalVector(networkdm,&localX);
56: VecSet(X,0.0);
57: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
58: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
60: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
61: SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,user_power);
63: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
64: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
65: DMRestoreLocalVector(networkdm,&localX);
66: return 0;
67: }
69: int main(int argc,char ** argv)
70: {
71: char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
72: PFDATA *pfdata;
73: PetscInt numEdges=0;
74: PetscInt *edges = NULL;
75: PetscInt i;
76: DM networkdm;
77: UserCtx_Power User;
78: #if defined(PETSC_USE_LOG)
79: PetscLogStage stage1,stage2;
80: #endif
81: PetscMPIInt rank;
82: PetscInt eStart, eEnd, vStart, vEnd,j;
83: PetscInt genj,loadj;
84: Vec X,F;
85: Mat J;
86: SNES snes;
88: PetscInitialize(&argc,&argv,"poweroptions",help);
89: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
90: {
91: /* 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 */
92: /* this is an experiment to see how the analyzer reacts */
93: const PetscMPIInt crank = rank;
95: /* Create an empty network object */
96: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
97: /* Register the components in the network */
98: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);
99: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);
100: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);
101: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);
103: PetscLogStageRegister("Read Data",&stage1);
104: PetscLogStagePush(stage1);
105: /* READ THE DATA */
106: if (!crank) {
107: /* READ DATA */
108: /* Only rank 0 reads the data */
109: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
110: PetscNew(&pfdata);
111: PFReadMatPowerData(pfdata,pfdata_file);
112: User.Sbase = pfdata->sbase;
114: numEdges = pfdata->nbranch;
115: PetscMalloc1(2*numEdges,&edges);
116: GetListofEdges_Power(pfdata,edges);
117: }
119: /* If external option activated. Introduce error in jacobian */
120: PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);
122: PetscLogStagePop();
123: MPI_Barrier(PETSC_COMM_WORLD);
124: PetscLogStageRegister("Create network",&stage2);
125: PetscLogStagePush(stage2);
126: /* Set number of nodes/edges */
127: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
128: DMNetworkAddSubnetwork(networkdm,"",numEdges,edges,NULL);
130: /* Set up the network layout */
131: DMNetworkLayoutSetUp(networkdm);
133: if (!crank) {
134: PetscFree(edges);
135: }
137: /* Add network components only process 0 has any data to add */
138: if (!crank) {
139: genj=0; loadj=0;
140: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
141: for (i = eStart; i < eEnd; i++) {
142: DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart],0);
143: }
144: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
145: for (i = vStart; i < vEnd; i++) {
146: DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart],2);
147: if (pfdata->bus[i-vStart].ngen) {
148: for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) {
149: DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++],0);
150: }
151: }
152: if (pfdata->bus[i-vStart].nload) {
153: for (j=0; j < pfdata->bus[i-vStart].nload; j++) {
154: DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++],0);
155: }
156: }
157: }
158: }
160: /* Set up DM for use */
161: DMSetUp(networkdm);
163: if (!crank) {
164: PetscFree(pfdata->bus);
165: PetscFree(pfdata->gen);
166: PetscFree(pfdata->branch);
167: PetscFree(pfdata->load);
168: PetscFree(pfdata);
169: }
171: /* Distribute networkdm to multiple processes */
172: DMNetworkDistribute(&networkdm,0);
174: PetscLogStagePop();
175: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
176: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
178: #if 0
179: EDGE_Power edge;
180: PetscInt key,kk,numComponents;
181: VERTEX_Power bus;
182: GEN gen;
183: LOAD load;
185: for (i = eStart; i < eEnd; i++) {
186: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
187: DMNetworkGetNumComponents(networkdm,i,&numComponents);
188: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
189: }
191: for (i = vStart; i < vEnd; i++) {
192: DMNetworkGetNumComponents(networkdm,i,&numComponents);
193: for (kk=0; kk < numComponents; kk++) {
194: DMNetworkGetComponent(networkdm,i,kk,&key,&component);
195: if (key == 1) {
196: bus = (VERTEX_Power)(component);
197: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
198: } else if (key == 2) {
199: gen = (GEN)(component);
200: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
201: } else if (key == 3) {
202: load = (LOAD)(component);
203: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
204: }
205: }
206: }
207: #endif
208: /* Broadcast Sbase to all processors */
209: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
211: DMCreateGlobalVector(networkdm,&X);
212: VecDuplicate(X,&F);
214: DMCreateMatrix(networkdm,&J);
215: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
217: SetInitialValues(networkdm,X,&User);
219: /* HOOK UP SOLVER */
220: SNESCreate(PETSC_COMM_WORLD,&snes);
221: SNESSetDM(snes,networkdm);
222: SNESSetFunction(snes,F,FormFunction,&User);
223: SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);
224: SNESSetFromOptions(snes);
226: SNESSolve(snes,NULL,X);
227: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
229: VecDestroy(&X);
230: VecDestroy(&F);
231: MatDestroy(&J);
233: SNESDestroy(&snes);
234: DMDestroy(&networkdm);
235: }
236: PetscFinalize();
237: return 0;
238: }
240: /*TEST
242: build:
243: depends: PFReadData.c pffunctions.c
244: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
246: test:
247: args: -snes_rtol 1.e-3
248: localrunfiles: poweroptions case9.m
249: output_file: output/power_1.out
251: test:
252: suffix: 2
253: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
254: nsize: 4
255: localrunfiles: poweroptions case9.m
256: output_file: output/power_1.out
258: TEST*/