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: }