xref: /petsc/src/snes/tutorials/network/power/pffunctions.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
1 /* function subroutines used by power.c */
2 
3 #include "power.h"
4 
5 PetscErrorCode GetListofEdges_Power(PFDATA *pfdata,PetscInt *edgelist)
6 {
7   PetscErrorCode ierr;
8   PetscInt       i,fbus,tbus,nbranches=pfdata->nbranch;
9   EDGE_Power     branch=pfdata->branch;
10   PetscBool      netview=PETSC_FALSE;
11 
12   PetscFunctionBegin;
13   ierr = PetscOptionsHasName(NULL,NULL, "-powernet_view",&netview);CHKERRQ(ierr);
14   for (i=0; i<nbranches; i++) {
15     fbus = branch[i].internal_i;
16     tbus = branch[i].internal_j;
17     edgelist[2*i]   = fbus;
18     edgelist[2*i+1] = tbus;
19     if (netview) {
20       ierr = PetscPrintf(PETSC_COMM_SELF,"branch %d, bus[%d] -> bus[%d]\n",i,fbus,tbus);CHKERRQ(ierr);
21     }
22   }
23   if (netview) {
24     for (i=0; i<pfdata->nbus; i++) {
25       if (pfdata->bus[i].ngen) {
26         ierr = PetscPrintf(PETSC_COMM_SELF," bus %D: gen\n",i);CHKERRQ(ierr);
27       } else if (pfdata->bus[i].nload) {
28         ierr = PetscPrintf(PETSC_COMM_SELF," bus %D: load\n",i);CHKERRQ(ierr);
29       }
30     }
31   }
32   PetscFunctionReturn(0);
33 }
34 
35 PetscErrorCode FormJacobian_Power_private(DM networkdm,Vec localX,Mat J,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
36 {
37   PetscErrorCode    ierr;
38   const PetscScalar *xarr;
39   PetscInt          i,v,row[2],col[8],e,vfrom,vto;
40   PetscInt          offsetfrom,offsetto,goffsetfrom,goffsetto,numComps;
41   PetscScalar       values[8];
42   PetscInt          j,key,offset,goffset;
43   PetscScalar       Vm;
44   UserCtx_Power     *user_power=(UserCtx_Power*)appctx;
45   PetscScalar       Sbase=user_power->Sbase;
46   VERTEX_Power      bus;
47   PetscBool         ghostvtex;
48   void*             component;
49 
50   PetscFunctionBegin;
51   ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
52 
53   for (v=0; v<nv; v++) {
54     ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr);
55 
56     ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr);
57     for (j = 0; j < numComps; j++) {
58       ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr);
59       ierr = DMNetworkGetVariableGlobalOffset(networkdm,vtx[v],&goffset);CHKERRQ(ierr);
60       ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr);
61 
62       if (key == user_power->compkey_bus) {
63         PetscInt       nconnedges;
64 	const PetscInt *connedges;
65 
66 	bus = (VERTEX_Power)(component);
67         if (!ghostvtex) {
68           /* Handle reference bus constrained dofs */
69           if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
70             row[0] = goffset; row[1] = goffset+1;
71             col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
72             values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
73             ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);
74             break;
75           }
76 
77           Vm = xarr[offset+1];
78 
79           /* Shunt injections */
80           row[0] = goffset; row[1] = goffset+1;
81           col[0] = goffset; col[1] = goffset+1;
82           values[0] = values[1] = values[2] = values[3] = 0.0;
83           if (bus->ide != PV_BUS) {
84             values[1] = 2.0*Vm*bus->gl/Sbase;
85             values[3] = -2.0*Vm*bus->bl/Sbase;
86           }
87           ierr = MatSetValues(J,2,row,2,col,values,ADD_VALUES);CHKERRQ(ierr);
88         }
89 
90 	ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr);
91 
92 	for (i=0; i < nconnedges; i++) {
93 	  EDGE_Power     branch;
94 	  VERTEX_Power   busf,bust;
95 	  PetscInt       keyf,keyt;
96           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
97           const PetscInt *cone;
98           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
99 
100 	  e = connedges[i];
101 	  ierr = DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch);CHKERRQ(ierr);
102 	  if (!branch->status) continue;
103 
104 	  Gff = branch->yff[0];
105 	  Bff = branch->yff[1];
106 	  Gft = branch->yft[0];
107 	  Bft = branch->yft[1];
108 	  Gtf = branch->ytf[0];
109 	  Btf = branch->ytf[1];
110 	  Gtt = branch->ytt[0];
111 	  Btt = branch->ytt[1];
112 
113 	  ierr = DMNetworkGetConnectedVertices(networkdm,edges[e],&cone);CHKERRQ(ierr);
114 	  vfrom = cone[0];
115 	  vto   = cone[1];
116 
117 	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);
118 	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);
119 	  ierr = DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);CHKERRQ(ierr);
120 	  ierr = DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);CHKERRQ(ierr);
121 
122 	  if (goffsetto < 0) goffsetto = -goffsetto - 1;
123 
124 	  thetaf  = xarr[offsetfrom];
125 	  Vmf     = xarr[offsetfrom+1];
126 	  thetat  = xarr[offsetto];
127 	  Vmt     = xarr[offsetto+1];
128 	  thetaft = thetaf - thetat;
129 	  thetatf = thetat - thetaf;
130 
131 	  ierr = DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf);CHKERRQ(ierr);
132 	  ierr = DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust);CHKERRQ(ierr);
133 
134 	  if (vfrom == vtx[v]) {
135 	    if (busf->ide != REF_BUS) {
136 	      /*    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));  */
137 	      row[0]  = goffsetfrom;
138 	      col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
139 	      values[0] =  Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
140 	      values[1] =  2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
141 	      values[2] =  Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
142 	      values[3] =  Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
143 
144 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
145 	    }
146 	    if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
147 	      row[0] = goffsetfrom+1;
148 	      col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
149 	      /*    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
150 	      values[0] =  Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
151 	      values[1] =  -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
152 	      values[2] =  Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
153 	      values[3] =  Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
154 
155 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
156 	    }
157 	  } else {
158 	    if (bust->ide != REF_BUS) {
159 	      row[0] = goffsetto;
160 	      col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
161 	      /*    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
162 	      values[0] =  Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
163 	      values[1] =  2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
164 	      values[2] =  Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
165 	      values[3] =  Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
166 
167 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
168 	    }
169 	    if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
170 	      row[0] = goffsetto+1;
171 	      col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
172 	      /*    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
173 	      values[0] =  Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
174 	      values[1] =  -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
175 	      values[2] =  Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
176 	      values[3] =  Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
177 
178 	      ierr = MatSetValues(J,1,row,4,col,values,ADD_VALUES);CHKERRQ(ierr);
179 	    }
180 	  }
181 	}
182 	if (!ghostvtex && bus->ide == PV_BUS) {
183 	  row[0] = goffset+1; col[0] = goffset+1; values[0]  = 1.0;
184           if (user_power->jac_error) values[0] = 50.0;
185           ierr = MatSetValues(J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
186 	}
187       }
188     }
189   }
190 
191   ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
192   PetscFunctionReturn(0);
193 }
194 
195 PetscErrorCode FormJacobian_Power(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
196 {
197   PetscErrorCode ierr;
198   DM             networkdm;
199   Vec            localX;
200   PetscInt       nv,ne;
201   const PetscInt *vtx,*edges;
202 
203   PetscFunctionBegin;
204   ierr = MatZeroEntries(J);CHKERRQ(ierr);
205 
206   ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);
207   ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);
208 
209   ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
210   ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
211 
212   ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
213   ierr = FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);CHKERRQ(ierr);
214 
215   ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);
216 
217   ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
218   ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
219   PetscFunctionReturn(0);
220 }
221 
222 PetscErrorCode FormFunction_Power(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
223 {
224   PetscErrorCode    ierr;
225   UserCtx_Power     *User=(UserCtx_Power*)appctx;
226   PetscInt          e,v,vfrom,vto;
227   const PetscScalar *xarr;
228   PetscScalar       *farr;
229   PetscInt          offsetfrom,offsetto,offset,i,j,key,numComps;
230   PetscScalar       Vm;
231   PetscScalar       Sbase=User->Sbase;
232   VERTEX_Power      bus=NULL;
233   GEN               gen;
234   LOAD              load;
235   PetscBool         ghostvtex;
236   void*             component;
237 
238   PetscFunctionBegin;
239   ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
240   ierr = VecGetArray(localF,&farr);CHKERRQ(ierr);
241 
242   for (v=0; v<nv; v++) {
243     ierr = DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);CHKERRQ(ierr);
244     ierr = DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);CHKERRQ(ierr);
245     ierr = DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);CHKERRQ(ierr);
246 
247     for (j = 0; j < numComps; j++) {
248       ierr = DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);CHKERRQ(ierr);
249       if (key == User->compkey_bus) {
250         PetscInt       nconnedges;
251 	const PetscInt *connedges;
252 
253 	bus = (VERTEX_Power)(component);
254 	/* Handle reference bus constrained dofs */
255 	if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
256 	  farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
257 	  farr[offset+1] = xarr[offset+1] - bus->vm;
258 	  break;
259 	}
260 
261 	if (!ghostvtex) {
262 	  Vm = xarr[offset+1];
263 
264 	  /* Shunt injections */
265 	  farr[offset] += Vm*Vm*bus->gl/Sbase;
266 	  if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
267 	}
268 
269 	ierr = DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);CHKERRQ(ierr);
270 	for (i=0; i < nconnedges; i++) {
271 	  EDGE_Power       branch;
272 	  PetscInt       keye;
273           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
274           const PetscInt *cone;
275           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
276 
277 	  e = connedges[i];
278 	  ierr = DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);CHKERRQ(ierr);
279 	  if (!branch->status) continue;
280 	  Gff = branch->yff[0];
281 	  Bff = branch->yff[1];
282 	  Gft = branch->yft[0];
283 	  Bft = branch->yft[1];
284 	  Gtf = branch->ytf[0];
285 	  Btf = branch->ytf[1];
286 	  Gtt = branch->ytt[0];
287 	  Btt = branch->ytt[1];
288 
289 	  ierr = DMNetworkGetConnectedVertices(networkdm,e,&cone);CHKERRQ(ierr);
290 	  vfrom = cone[0];
291 	  vto   = cone[1];
292 
293 	  ierr = DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);CHKERRQ(ierr);
294 	  ierr = DMNetworkGetVariableOffset(networkdm,vto,&offsetto);CHKERRQ(ierr);
295 
296 	  thetaf = xarr[offsetfrom];
297 	  Vmf     = xarr[offsetfrom+1];
298 	  thetat = xarr[offsetto];
299 	  Vmt     = xarr[offsetto+1];
300 	  thetaft = thetaf - thetat;
301 	  thetatf = thetat - thetaf;
302 
303 	  if (vfrom == vtx[v]) {
304 	    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
305 	    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
306 	  } else {
307 	    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
308 	    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
309 	  }
310 	}
311       } else if (key == User->compkey_gen) {
312 	if (!ghostvtex) {
313 	  gen = (GEN)(component);
314 	  if (!gen->status) continue;
315 	  farr[offset] += -gen->pg/Sbase;
316 	  farr[offset+1] += -gen->qg/Sbase;
317 	}
318       } else if (key == User->compkey_load) {
319 	if (!ghostvtex) {
320 	  load = (LOAD)(component);
321 	  farr[offset] += load->pl/Sbase;
322 	  farr[offset+1] += load->ql/Sbase;
323 	}
324       }
325     }
326     if (bus && bus->ide == PV_BUS) {
327       farr[offset+1] = xarr[offset+1] - bus->vm;
328     }
329   }
330   ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
331   ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr);
332   PetscFunctionReturn(0);
333 }
334 
335 PetscErrorCode SetInitialGuess_Power(DM networkdm,Vec localX,PetscInt nv,PetscInt ne,const PetscInt *vtx,const PetscInt *edges,void* appctx)
336 {
337   PetscErrorCode ierr;
338   VERTEX_Power   bus;
339   PetscInt       i;
340   GEN            gen;
341   PetscBool      ghostvtex;
342   PetscScalar    *xarr;
343   PetscInt       key,numComps,j,offset;
344   void*          component;
345   PetscMPIInt    rank;
346   MPI_Comm       comm;
347   UserCtx_Power  *User=(UserCtx_Power*)appctx;
348 
349   PetscFunctionBegin;
350   ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr);
351   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
352   ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr);
353   for (i = 0; i < nv; i++) {
354     ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr);
355     if (ghostvtex) continue;
356 
357     ierr = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr);
358     ierr = DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);CHKERRQ(ierr);
359     for (j=0; j < numComps; j++) {
360       ierr = DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component);CHKERRQ(ierr);
361       if (key == User->compkey_bus) {
362 	bus = (VERTEX_Power)(component);
363 	xarr[offset] = bus->va*PETSC_PI/180.0;
364 	xarr[offset+1] = bus->vm;
365       } else if(key == User->compkey_gen) {
366 	gen = (GEN)(component);
367 	if (!gen->status) continue;
368 	xarr[offset+1] = gen->vs;
369 	break;
370       }
371     }
372   }
373   ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr);
374   PetscFunctionReturn(0);
375 }
376