Actual source code: power.c
petsc-3.12.0 2019-09-29
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: DMNetworkGetSubnetworkInfo(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: DMNetworkGetSubnetworkInfo(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,numVertices=0;
79: PetscInt *edges = NULL;
80: PetscInt i;
81: DM networkdm;
82: UserCtx_Power User;
83: PetscLogStage stage1,stage2;
84: PetscMPIInt rank;
85: PetscInt eStart, eEnd, vStart, vEnd,j;
86: PetscInt genj,loadj;
87: Vec X,F;
88: Mat J;
89: SNES snes;
91: PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr;
92: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
93: {
94: /* 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 */
95: /* this is an experiment to see how the analyzer reacts */
96: const PetscMPIInt crank = rank;
98: /* Create an empty network object */
99: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
100: /* Register the components in the network */
101: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&User.compkey_branch);
102: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&User.compkey_bus);
103: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&User.compkey_gen);
104: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&User.compkey_load);
106: PetscLogStageRegister("Read Data",&stage1);
107: PetscLogStagePush(stage1);
108: /* READ THE DATA */
109: if (!crank) {
110: /* READ DATA */
111: /* Only rank 0 reads the data */
112: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
113: PetscNew(&pfdata);
114: PFReadMatPowerData(pfdata,pfdata_file);
115: User.Sbase = pfdata->sbase;
117: numEdges = pfdata->nbranch;
118: numVertices = pfdata->nbus;
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: DMNetworkSetSizes(networkdm,1,&numVertices,&numEdges,0,NULL);
133: /* Add edge connectivity */
134: DMNetworkSetEdgeList(networkdm,&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]);
148: }
149: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
150: for (i = vStart; i < vEnd; i++) {
151: DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart]);
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++]);
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++]);
160: }
161: }
162: /* Add number of variables */
163: DMNetworkAddNumVariables(networkdm,i,2);
164: }
165: }
167: /* Set up DM for use */
168: DMSetUp(networkdm);
170: if (!crank) {
171: PetscFree(pfdata->bus);
172: PetscFree(pfdata->gen);
173: PetscFree(pfdata->branch);
174: PetscFree(pfdata->load);
175: PetscFree(pfdata);
176: }
178: /* Distribute networkdm to multiple processes */
179: DMNetworkDistribute(&networkdm,0);
181: PetscLogStagePop();
182: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
183: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
185: #if 0
186: EDGE_Power edge;
187: PetscInt key,kk,numComponents;
188: VERTEX_Power bus;
189: GEN gen;
190: LOAD load;
192: for (i = eStart; i < eEnd; i++) {
193: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
194: DMNetworkGetNumComponents(networkdm,i,&numComponents);
195: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
196: }
198: for (i = vStart; i < vEnd; i++) {
199: DMNetworkGetNumComponents(networkdm,i,&numComponents);
200: for (kk=0; kk < numComponents; kk++) {
201: DMNetworkGetComponent(networkdm,i,kk,&key,&component);
202: if (key == 1) {
203: bus = (VERTEX_Power)(component);
204: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
205: } else if (key == 2) {
206: gen = (GEN)(component);
207: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
208: } else if (key == 3) {
209: load = (LOAD)(component);
210: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
211: }
212: }
213: }
214: #endif
215: /* Broadcast Sbase to all processors */
216: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
218: DMCreateGlobalVector(networkdm,&X);
219: VecDuplicate(X,&F);
221: DMCreateMatrix(networkdm,&J);
222: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
224: SetInitialValues(networkdm,X,&User);
226: /* HOOK UP SOLVER */
227: SNESCreate(PETSC_COMM_WORLD,&snes);
228: SNESSetDM(snes,networkdm);
229: SNESSetFunction(snes,F,FormFunction,&User);
230: SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);
231: SNESSetFromOptions(snes);
233: SNESSolve(snes,NULL,X);
234: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
236: VecDestroy(&X);
237: VecDestroy(&F);
238: MatDestroy(&J);
240: SNESDestroy(&snes);
241: DMDestroy(&networkdm);
242: }
243: PetscFinalize();
244: return ierr;
245: }
247: /*TEST
249: build:
250: depends: PFReadData.c pffunctions.c
251: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
254: test:
255: args: -snes_rtol 1.e-3
256: localrunfiles: poweroptions case9.m
257: output_file: output/power_1.out
258: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
260: test:
261: suffix: 2
262: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
263: nsize: 4
264: localrunfiles: poweroptions case9.m
265: output_file: output/power_1.out
266: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
268: TEST*/