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