Actual source code: power.c
petsc-3.14.4 2021-02-03
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: #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: numVertices = pfdata->nbus;
122: PetscMalloc1(2*numEdges,&edges);
123: GetListofEdges_Power(pfdata,edges);
124: }
126: /* If external option activated. Introduce error in jacobian */
127: PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);
129: PetscLogStagePop();
130: MPI_Barrier(PETSC_COMM_WORLD);
131: PetscLogStageRegister("Create network",&stage2);
132: PetscLogStagePush(stage2);
133: /* Set number of nodes/edges */
134: DMNetworkSetSizes(networkdm,1,&numVertices,&numEdges,0,NULL);
135: /* Add edge connectivity */
136: DMNetworkSetEdgeList(networkdm,&edges,NULL);
137: /* Set up the network layout */
138: DMNetworkLayoutSetUp(networkdm);
140: if (!crank) {
141: PetscFree(edges);
142: }
144: /* Add network components only process 0 has any data to add*/
145: if (!crank) {
146: genj=0; loadj=0;
147: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
148: for (i = eStart; i < eEnd; i++) {
149: DMNetworkAddComponent(networkdm,i,User.compkey_branch,&pfdata->branch[i-eStart]);
150: }
151: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
152: for (i = vStart; i < vEnd; i++) {
153: DMNetworkAddComponent(networkdm,i,User.compkey_bus,&pfdata->bus[i-vStart]);
154: if (pfdata->bus[i-vStart].ngen) {
155: for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) {
156: DMNetworkAddComponent(networkdm,i,User.compkey_gen,&pfdata->gen[genj++]);
157: }
158: }
159: if (pfdata->bus[i-vStart].nload) {
160: for (j=0; j < pfdata->bus[i-vStart].nload; j++) {
161: DMNetworkAddComponent(networkdm,i,User.compkey_load,&pfdata->load[loadj++]);
162: }
163: }
164: /* Add number of variables */
165: DMNetworkAddNumVariables(networkdm,i,2);
166: }
167: }
169: /* Set up DM for use */
170: DMSetUp(networkdm);
172: if (!crank) {
173: PetscFree(pfdata->bus);
174: PetscFree(pfdata->gen);
175: PetscFree(pfdata->branch);
176: PetscFree(pfdata->load);
177: PetscFree(pfdata);
178: }
180: /* Distribute networkdm to multiple processes */
181: DMNetworkDistribute(&networkdm,0);
183: PetscLogStagePop();
184: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
185: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
187: #if 0
188: EDGE_Power edge;
189: PetscInt key,kk,numComponents;
190: VERTEX_Power bus;
191: GEN gen;
192: LOAD load;
194: for (i = eStart; i < eEnd; i++) {
195: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&edge);
196: DMNetworkGetNumComponents(networkdm,i,&numComponents);
197: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
198: }
200: for (i = vStart; i < vEnd; i++) {
201: DMNetworkGetNumComponents(networkdm,i,&numComponents);
202: for (kk=0; kk < numComponents; kk++) {
203: DMNetworkGetComponent(networkdm,i,kk,&key,&component);
204: if (key == 1) {
205: bus = (VERTEX_Power)(component);
206: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
207: } else if (key == 2) {
208: gen = (GEN)(component);
209: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
210: } else if (key == 3) {
211: load = (LOAD)(component);
212: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
213: }
214: }
215: }
216: #endif
217: /* Broadcast Sbase to all processors */
218: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
220: DMCreateGlobalVector(networkdm,&X);
221: VecDuplicate(X,&F);
223: DMCreateMatrix(networkdm,&J);
224: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
226: SetInitialValues(networkdm,X,&User);
228: /* HOOK UP SOLVER */
229: SNESCreate(PETSC_COMM_WORLD,&snes);
230: SNESSetDM(snes,networkdm);
231: SNESSetFunction(snes,F,FormFunction,&User);
232: SNESSetJacobian(snes,J,J,FormJacobian_Power,&User);
233: SNESSetFromOptions(snes);
235: SNESSolve(snes,NULL,X);
236: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
238: VecDestroy(&X);
239: VecDestroy(&F);
240: MatDestroy(&J);
242: SNESDestroy(&snes);
243: DMDestroy(&networkdm);
244: }
245: PetscFinalize();
246: return ierr;
247: }
249: /*TEST
251: build:
252: depends: PFReadData.c pffunctions.c
253: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
256: test:
257: args: -snes_rtol 1.e-3
258: localrunfiles: poweroptions case9.m
259: output_file: output/power_1.out
260: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
262: test:
263: suffix: 2
264: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
265: nsize: 4
266: localrunfiles: poweroptions case9.m
267: output_file: output/power_1.out
268: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
270: TEST*/