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