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