xref: /petsc/src/snes/tutorials/network/power/power2.c (revision 6a98f8dc3f2c9149905a87dc2e9d0fedaf64e09a) 
1 static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks 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                       The data file format used is from the MatPower package (http://www.pserc.cornell.edu//matpower/).\n\
4                       This example shows the use of subnetwork feature in DMNetwork. It creates duplicates of the same network which are treated as subnetworks.\n\
5                       Run this program: mpiexec -n <n> ./pf2\n\
6                       mpiexec -n <n> ./pf2 \n";
7 
8 /* T
9    Concepts: DMNetwork
10    Concepts: PETSc SNES solver
11 */
12 
13 #include "power.h"
14 #include <petscdmnetwork.h>
15 
16 PetscErrorCode FormFunction_Subnet(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
17 {
18   PetscErrorCode ierr;
19   UserCtx_Power  *User=(UserCtx_Power*)appctx;
20   PetscInt       e,v,vfrom,vto;
21   const PetscScalar *xarr;
22   PetscScalar    *farr;
23   PetscInt       offsetfrom,offsetto,offset;
24 
25   PetscFunctionBegin;
26   ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
27   ierr = VecGetArray(localF,&farr);CHKERRQ(ierr);
28 
29   for (v=0; v<nv; v++) {
30     PetscInt    i,j,key;
31     PetscScalar Vm;
32     PetscScalar Sbase=User->Sbase;
33     VERTEX_Power  bus=NULL;
34     GEN         gen;
35     LOAD        load;
36     PetscBool   ghostvtex;
37     PetscInt    numComps;
38     void*       component;
39 
40     ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr);
41     ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr);
42     ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr);
43     for (j = 0; j < numComps; j++) {
44       ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr);
45       if (key == 1) {
46         PetscInt       nconnedges;
47 	const PetscInt *connedges;
48 
49 	bus = (VERTEX_Power)(component);
50 	/* Handle reference bus constrained dofs */
51 	if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
52 	  farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
53 	  farr[offset+1] = xarr[offset+1] - bus->vm;
54 	  break;
55 	}
56 
57 	if (!ghostvtex) {
58 	  Vm = xarr[offset+1];
59 
60 	  /* Shunt injections */
61 	  farr[offset] += Vm*Vm*bus->gl/Sbase;
62 	  if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
63 	}
64 
65 	ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr);
66 	for (i=0; i < nconnedges; i++) {
67 	  EDGE_Power       branch;
68 	  PetscInt       keye;
69           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
70           const PetscInt *cone;
71           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
72 
73 	  e = connedges[i];
74 	  ierr = DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);CHKERRQ(ierr);
75 	  if (!branch->status) continue;
76 	  Gff = branch->yff[0];
77 	  Bff = branch->yff[1];
78 	  Gft = branch->yft[0];
79 	  Bft = branch->yft[1];
80 	  Gtf = branch->ytf[0];
81 	  Btf = branch->ytf[1];
82 	  Gtt = branch->ytt[0];
83 	  Btt = branch->ytt[1];
84 
85 	  ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr);
86 	  vfrom = cone[0];
87 	  vto   = cone[1];
88 
89 	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);
90 	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);
91 
92 	  thetaf = xarr[offsetfrom];
93 	  Vmf     = xarr[offsetfrom+1];
94 	  thetat = xarr[offsetto];
95 	  Vmt     = xarr[offsetto+1];
96 	  thetaft = thetaf - thetat;
97 	  thetatf = thetat - thetaf;
98 
99 	  if (vfrom == vtx[v]) {
100 	    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
101 	    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
102 	  } else {
103 	    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
104 	    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
105 	  }
106 	}
107       } else if (key == 2) {
108 	if (!ghostvtex) {
109 	  gen = (GEN)(component);
110 	  if (!gen->status) continue;
111 	  farr[offset] += -gen->pg/Sbase;
112 	  farr[offset+1] += -gen->qg/Sbase;
113 	}
114       } else if (key == 3) {
115 	if (!ghostvtex) {
116 	  load = (LOAD)(component);
117 	  farr[offset] += load->pl/Sbase;
118 	  farr[offset+1] += load->ql/Sbase;
119 	}
120       }
121     }
122     if (bus && bus->ide == PV_BUS) {
123       farr[offset+1] = xarr[offset+1] - bus->vm;
124     }
125   }
126   ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
127   ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr);
128   PetscFunctionReturn(0);
129 }
130 
131 
132 PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
133 {
134   PetscErrorCode ierr;
135   DM             networkdm;
136   Vec           localX,localF;
137   PetscInt      nv,ne;
138   const PetscInt *vtx,*edges;
139 
140   PetscFunctionBegin;
141   ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);
142   ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);
143   ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr);
144   ierr = VecSet(F,0.0);CHKERRQ(ierr);
145 
146   ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
147   ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
148 
149   ierr = DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);
150   ierr = DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);
151 
152   /* Form Function for first subnetwork */
153   ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
154   ierr = FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
155 
156   /* Form Function for second subnetwork */
157   ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
158   ierr = FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
159 
160   ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);
161 
162   ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr);
163   ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr);
164   ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr);
165   PetscFunctionReturn(0);
166 }
167 
168 PetscErrorCode FormJacobian_Subnet(DM networkdm,Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
169 {
170   PetscErrorCode ierr;
171   UserCtx_Power  *User=(UserCtx_Power*)appctx;
172   PetscInt       e,v,vfrom,vto;
173   const PetscScalar *xarr;
174   PetscInt       offsetfrom,offsetto,goffsetfrom,goffsetto;
175   PetscInt       row[2],col[8];
176   PetscScalar    values[8];
177 
178   PetscFunctionBegin;
179   ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
180 
181   for (v=0; v<nv; v++) {
182     PetscInt    i,j,key;
183     PetscInt    offset,goffset;
184     PetscScalar Vm;
185     PetscScalar Sbase=User->Sbase;
186     VERTEX_Power  bus;
187     PetscBool   ghostvtex;
188     PetscInt    numComps;
189     void*       component;
190 
191     ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr);
192     ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr);
193     for (j = 0; j < numComps; j++) {
194       ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr);
195       ierr = DMNetworkGetVariableGlobalOffset(networkdm,vtx[v],&goffset);CHKERRQ(ierr);
196       ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr);
197       if (key == 1) {
198         PetscInt       nconnedges;
199 	const PetscInt *connedges;
200 
201 	bus = (VERTEX_Power)(component);
202 	if (!ghostvtex) {
203 	  /* Handle reference bus constrained dofs */
204 	  if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
205 	    row[0] = goffset; row[1] = goffset+1;
206 	    col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
207 	    values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
208 	    ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);
209 	    break;
210 	  }
211 
212 	  Vm = xarr[offset+1];
213 
214 	  /* Shunt injections */
215           row[0] = goffset; row[1] = goffset+1;
216           col[0] = goffset; col[1] = goffset+1;
217           values[0] = values[1] = values[2] = values[3] = 0.0;
218           if (bus->ide != PV_BUS) {
219             values[1] = 2.0*Vm*bus->gl/Sbase;
220             values[3] = -2.0*Vm*bus->bl/Sbase;
221           }
222           ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);
223 	}
224 
225 	ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr);
226 	for (i=0; i < nconnedges; i++) {
227 	  EDGE_Power       branch;
228 	  VERTEX_Power     busf,bust;
229 	  PetscInt       keyf,keyt;
230           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
231           const PetscInt *cone;
232           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
233 
234 	  e = connedges[i];
235 	  ierr = DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch);CHKERRQ(ierr);
236 	  if (!branch->status) continue;
237 
238 	  Gff = branch->yff[0];
239 	  Bff = branch->yff[1];
240 	  Gft = branch->yft[0];
241 	  Bft = branch->yft[1];
242 	  Gtf = branch->ytf[0];
243 	  Btf = branch->ytf[1];
244 	  Gtt = branch->ytt[0];
245 	  Btt = branch->ytt[1];
246 
247 	  ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr);
248 	  vfrom = cone[0];
249 	  vto   = cone[1];
250 
251 	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);
252 	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);
253 	  ierr = DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);CHKERRQ(ierr);
254 	  ierr = DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);CHKERRQ(ierr);
255 
256 	  if (goffsetto < 0) goffsetto = -goffsetto - 1;
257 
258 	  thetaf = xarr[offsetfrom];
259 	  Vmf     = xarr[offsetfrom+1];
260 	  thetat = xarr[offsetto];
261 	  Vmt     = xarr[offsetto+1];
262 	  thetaft = thetaf - thetat;
263 	  thetatf = thetat - thetaf;
264 
265 	  ierr = DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf);CHKERRQ(ierr);
266 	  ierr = DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust);CHKERRQ(ierr);
267 
268 	  if (vfrom == vtx[v]) {
269 	    if (busf->ide != REF_BUS) {
270 	      /*    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));  */
271 	      row[0]  = goffsetfrom;
272 	      col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
273 	      values[0] =  Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
274 	      values[1] =  2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
275 	      values[2] =  Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
276 	      values[3] =  Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
277 
278 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
279 	    }
280 	    if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
281 	      row[0] = goffsetfrom+1;
282 	      col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
283 	      /*    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
284 	      values[0] =  Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
285 	      values[1] =  -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
286 	      values[2] =  Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
287 	      values[3] =  Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
288 
289 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
290 	    }
291 	  } else {
292 	    if (bust->ide != REF_BUS) {
293 	      row[0] = goffsetto;
294 	      col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
295 	      /*    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
296 	      values[0] =  Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
297 	      values[1] =  2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
298 	      values[2] =  Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
299 	      values[3] =  Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
300 
301 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
302 	    }
303 	    if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
304 	      row[0] = goffsetto+1;
305 	      col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
306 	      /*    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
307 	      values[0] =  Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
308 	      values[1] =  -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
309 	      values[2] =  Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
310 	      values[3] =  Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
311 
312 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
313 	    }
314 	  }
315 	}
316 	if (!ghostvtex && bus->ide == PV_BUS) {
317 	  row[0] = goffset+1; col[0] = goffset+1;
318 	  values[0]  = 1.0;
319 	  ierr = MatSetValues(J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
320 	}
321       }
322     }
323   }
324   ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
325   PetscFunctionReturn(0);
326 }
327 
328 PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
329 {
330   PetscErrorCode ierr;
331   DM            networkdm;
332   Vec           localX;
333   PetscInt      ne,nv;
334   const PetscInt *vtx,*edges;
335 
336   PetscFunctionBegin;
337   ierr = MatZeroEntries(J);CHKERRQ(ierr);
338 
339   ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);
340   ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);
341 
342   ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
343   ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
344 
345   /* Form Jacobian for first subnetwork */
346   ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
347   ierr = FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
348 
349   /* Form Jacobian for second subnetwork */
350   ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
351   ierr = FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
352 
353   ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);
354 
355   ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
356   ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
357   PetscFunctionReturn(0);
358 }
359 
360 PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
361 {
362   PetscErrorCode ierr;
363   VERTEX_Power     bus;
364   PetscInt       i;
365   GEN            gen;
366   PetscBool      ghostvtex;
367   PetscScalar    *xarr;
368   PetscInt       key,numComps,j,offset;
369   void*          component;
370 
371   PetscFunctionBegin;
372   ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr);
373   for (i = 0; i < nv; i++) {
374     ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr);
375     if (ghostvtex) continue;
376 
377     ierr = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr);
378     ierr = DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);CHKERRQ(ierr);
379     for (j=0; j < numComps; j++) {
380       ierr = DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component);CHKERRQ(ierr);
381       if (key == 1) {
382 	bus = (VERTEX_Power)(component);
383 	xarr[offset] = bus->va*PETSC_PI/180.0;
384 	xarr[offset+1] = bus->vm;
385       } else if(key == 2) {
386 	gen = (GEN)(component);
387 	if (!gen->status) continue;
388 	xarr[offset+1] = gen->vs;
389 	break;
390       }
391     }
392   }
393   ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr);
394   PetscFunctionReturn(0);
395 }
396 
397 PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
398 {
399   PetscErrorCode ierr;
400   PetscInt       nv,ne;
401   const PetscInt *vtx,*edges;
402   Vec            localX;
403 
404   PetscFunctionBegin;
405   ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);
406 
407   ierr = VecSet(X,0.0);CHKERRQ(ierr);
408   ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
409   ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
410 
411   /* Set initial guess for first subnetwork */
412   ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
413   ierr = SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
414 
415   /* Set initial guess for second subnetwork */
416   ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
417   ierr = SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
418 
419   ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr);
420   ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr);
421   ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);
422   PetscFunctionReturn(0);
423 }
424 
425 int main(int argc,char ** argv)
426 {
427   PetscErrorCode   ierr;
428   char             pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
429   PFDATA           *pfdata1,*pfdata2;
430   PetscInt         numEdges1=0,numVertices1=0,numEdges2=0,numVertices2=0;
431   PetscInt         *edgelist1 = NULL,*edgelist2 = NULL;
432   DM               networkdm;
433   PetscInt         componentkey[4];
434   UserCtx_Power    User;
435   PetscLogStage    stage1,stage2;
436   PetscMPIInt      rank;
437   PetscInt         nsubnet = 2;
438   PetscInt         numVertices[2],numEdges[2];
439   PetscInt         *edgelist[2];
440   PetscInt         nv,ne;
441   const PetscInt   *vtx;
442   const PetscInt   *edges;
443   PetscInt         i,j;
444   PetscInt         genj,loadj;
445   Vec              X,F;
446   Mat              J;
447   SNES             snes;
448 
449 
450   ierr = PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr;
451   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
452   {
453     /* 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 */
454     /* this is an experiment to see how the analyzer reacts */
455     const PetscMPIInt crank = rank;
456 
457     /* Create an empty network object */
458     ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr);
459 
460     /* Register the components in the network */
461     ierr = DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);CHKERRQ(ierr);
462     ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);CHKERRQ(ierr);
463     ierr = DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);CHKERRQ(ierr);
464     ierr = DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);CHKERRQ(ierr);
465 
466     ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr);
467     PetscLogStagePush(stage1);
468     /* READ THE DATA */
469     if (!crank) {
470       /* Only rank 0 reads the data */
471       ierr = PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);CHKERRQ(ierr);
472       /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
473 
474       /*    READ DATA FOR THE FIRST SUBNETWORK */
475       ierr = PetscNew(&pfdata1);CHKERRQ(ierr);
476       ierr = PFReadMatPowerData(pfdata1,pfdata_file);CHKERRQ(ierr);
477       User.Sbase = pfdata1->sbase;
478 
479       numEdges1 = pfdata1->nbranch;
480       numVertices1 = pfdata1->nbus;
481 
482       ierr = PetscMalloc1(2*numEdges1,&edgelist1);CHKERRQ(ierr);
483       ierr = GetListofEdges_Power(pfdata1,edgelist1);CHKERRQ(ierr);
484 
485       /*    READ DATA FOR THE SECOND SUBNETWORK */
486       ierr = PetscNew(&pfdata2);CHKERRQ(ierr);
487       ierr = PFReadMatPowerData(pfdata2,pfdata_file);CHKERRQ(ierr);
488       User.Sbase = pfdata2->sbase;
489 
490       numEdges2 = pfdata2->nbranch;
491       numVertices2 = pfdata2->nbus;
492 
493       ierr = PetscMalloc1(2*numEdges2,&edgelist2);CHKERRQ(ierr);
494       ierr = GetListofEdges_Power(pfdata2,edgelist2);CHKERRQ(ierr);
495     }
496 
497     PetscLogStagePop();
498     ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
499     ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr);
500     PetscLogStagePush(stage2);
501 
502     /* Set number of nodes/edges */
503     numVertices[0] = numVertices1; numVertices[1] = numVertices2;
504     numEdges[0] = numEdges1; numEdges[1] = numEdges2;
505     ierr = DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,0,NULL);CHKERRQ(ierr);
506 
507     edgelist[0] = edgelist1; edgelist[1] = edgelist2;
508 
509     /* Add edge connectivity */
510     ierr = DMNetworkSetEdgeList(networkdm,edgelist,NULL);CHKERRQ(ierr);
511 
512     /* Set up the network layout */
513     ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr);
514 
515     /* Add network components only process 0 has any data to add*/
516     if (!crank) {
517       genj=0; loadj=0;
518 
519       /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
520       ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
521 
522       for (i = 0; i < ne; i++) {
523         ierr = DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i]);CHKERRQ(ierr);
524       }
525 
526       for (i = 0; i < nv; i++) {
527         ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i]);CHKERRQ(ierr);
528         if (pfdata1->bus[i].ngen) {
529           for (j = 0; j < pfdata1->bus[i].ngen; j++) {
530             ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++]);CHKERRQ(ierr);
531           }
532         }
533         if (pfdata1->bus[i].nload) {
534           for (j=0; j < pfdata1->bus[i].nload; j++) {
535             ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++]);CHKERRQ(ierr);
536           }
537         }
538         /* Add number of variables */
539         ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr);
540       }
541 
542       genj=0; loadj=0;
543 
544       /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
545       ierr = DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
546 
547       for (i = 0; i < ne; i++) {
548         ierr = DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i]);CHKERRQ(ierr);
549       }
550 
551       for (i = 0; i < nv; i++) {
552         ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i]);CHKERRQ(ierr);
553         if (pfdata2->bus[i].ngen) {
554           for (j = 0; j < pfdata2->bus[i].ngen; j++) {
555             ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++]);CHKERRQ(ierr);
556           }
557         }
558         if (pfdata2->bus[i].nload) {
559           for (j=0; j < pfdata2->bus[i].nload; j++) {
560             ierr = DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++]);CHKERRQ(ierr);
561           }
562         }
563         /* Add number of variables */
564         ierr = DMNetworkAddNumVariables(networkdm,vtx[i],2);CHKERRQ(ierr);
565       }
566     }
567 
568     /* Set up DM for use */
569     ierr = DMSetUp(networkdm);CHKERRQ(ierr);
570 
571     if (!crank) {
572       ierr = PetscFree(edgelist1);CHKERRQ(ierr);
573       ierr = PetscFree(edgelist2);CHKERRQ(ierr);
574     }
575 
576     if (!crank) {
577       ierr = PetscFree(pfdata1->bus);CHKERRQ(ierr);
578       ierr = PetscFree(pfdata1->gen);CHKERRQ(ierr);
579       ierr = PetscFree(pfdata1->branch);CHKERRQ(ierr);
580       ierr = PetscFree(pfdata1->load);CHKERRQ(ierr);
581       ierr = PetscFree(pfdata1);CHKERRQ(ierr);
582 
583       ierr = PetscFree(pfdata2->bus);CHKERRQ(ierr);
584       ierr = PetscFree(pfdata2->gen);CHKERRQ(ierr);
585       ierr = PetscFree(pfdata2->branch);CHKERRQ(ierr);
586       ierr = PetscFree(pfdata2->load);CHKERRQ(ierr);
587       ierr = PetscFree(pfdata2);CHKERRQ(ierr);
588     }
589 
590     /* Distribute networkdm to multiple processes */
591     ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr);
592 
593     PetscLogStagePop();
594 
595     /* Broadcast Sbase to all processors */
596     ierr = MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
597 
598     ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr);
599     ierr = VecDuplicate(X,&F);CHKERRQ(ierr);
600 
601     ierr = DMCreateMatrix(networkdm,&J);CHKERRQ(ierr);
602     ierr = MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
603 
604     ierr = SetInitialValues(networkdm,X,&User);CHKERRQ(ierr);
605 
606     /* HOOK UP SOLVER */
607     ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
608     ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr);
609     ierr = SNESSetFunction(snes,F,FormFunction,&User);CHKERRQ(ierr);
610     ierr = SNESSetJacobian(snes,J,J,FormJacobian,&User);CHKERRQ(ierr);
611     ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
612 
613     ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr);
614     /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
615 
616     ierr = VecDestroy(&X);CHKERRQ(ierr);
617     ierr = VecDestroy(&F);CHKERRQ(ierr);
618     ierr = MatDestroy(&J);CHKERRQ(ierr);
619 
620     ierr = SNESDestroy(&snes);CHKERRQ(ierr);
621     ierr = DMDestroy(&networkdm);CHKERRQ(ierr);
622   }
623   ierr = PetscFinalize();
624   return ierr;
625 }
626 
627 /*TEST
628 
629    build:
630      depends: PFReadData.c pffunctions.c
631      requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
632 
633 
634    test:
635      args: -snes_rtol 1.e-3
636      localrunfiles: poweroptions case9.m
637      output_file: output/power2_1.out
638      requires: double !complex
639 
640    test:
641      suffix: 2
642      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
643      nsize: 4
644      localrunfiles: poweroptions case9.m
645      output_file: output/power2_1.out
646      requires: double !complex
647 
648 TEST*/
649