Actual source code: pffunctions.c
1: /* function subroutines used by power.c */
3: #include "power.h"
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;
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: }
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;
46: PetscFunctionBegin;
47: PetscCall(VecGetArrayRead(localX, &xarr));
49: for (v = 0; v < nv; v++) {
50: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex));
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));
58: if (key == user_power->compkey_bus) {
59: PetscInt nconnedges;
60: const PetscInt *connedges;
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: }
80: Vm = xarr[offset + 1];
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: }
95: PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
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;
105: e = connedges[i];
106: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL));
107: if (!branch->status) continue;
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];
118: PetscCall(DMNetworkGetConnectedVertices(networkdm, edges[e], &cone));
119: vfrom = cone[0];
120: vto = cone[1];
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));
127: if (goffsetto < 0) goffsetto = -goffsetto - 1;
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;
136: PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL));
137: PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL));
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 */
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));
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 */
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));
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: }
210: PetscCall(VecRestoreArrayRead(localX, &xarr));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
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;
221: PetscFunctionBegin;
222: PetscCall(MatZeroEntries(J));
224: PetscCall(SNESGetDM(snes, &networkdm));
225: PetscCall(DMGetLocalVector(networkdm, &localX));
227: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
228: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
230: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
231: PetscCall(FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx));
233: PetscCall(DMRestoreLocalVector(networkdm, &localX));
235: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
236: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
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;
255: PetscFunctionBegin;
256: PetscCall(VecGetArrayRead(localX, &xarr));
257: PetscCall(VecGetArray(localF, &farr));
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));
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;
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: }
278: if (!ghostvtex) {
279: Vm = xarr[offset + 1];
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: }
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;
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];
306: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
307: vfrom = cone[0];
308: vto = cone[1];
310: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
311: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
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;
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: }
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;
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;
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: }