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