Actual source code: ex9busdmnetwork.c

      1: static char help[] = "This example uses the same problem set up of ex9busdmnetwork.c. \n\
      2: It demonstrates setting and accessing of variables for individual components, instead of \n\
      3: the network vertices (as used in ex9busdmnetwork.c). This is especially useful where vertices \n\
      4: /edges have multiple-components associated with them and one or more components has physics \n\
      5: associated with it. \n\
      6: Input parameters include:\n\
      7:   -nc : number of copies of the base case\n\n";

      9: /*
     10:    This example was modified from ex9busdmnetwork.c.
     11: */

     13: #include <petscts.h>
     14: #include <petscdmnetwork.h>

     16: #define FREQ    60
     17: #define W_S     (2 * PETSC_PI * FREQ)
     18: #define NGEN    3 /* No. of generators in the 9 bus system */
     19: #define NLOAD   3 /* No. of loads in the 9 bus system */
     20: #define NBUS    9 /* No. of buses in the 9 bus system */
     21: #define NBRANCH 9 /* No. of branches in the 9 bus system */

     23: typedef struct {
     24:   PetscInt    id;    /* Bus Number or extended bus name*/
     25:   PetscScalar mbase; /* MVA base of the machine */
     26:   PetscScalar PG;    /* Generator active power output */
     27:   PetscScalar QG;    /* Generator reactive power output */

     29:   /* Generator constants */
     30:   PetscScalar H;    /* Inertia constant */
     31:   PetscScalar Rs;   /* Stator Resistance */
     32:   PetscScalar Xd;   /* d-axis reactance */
     33:   PetscScalar Xdp;  /* d-axis transient reactance */
     34:   PetscScalar Xq;   /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
     35:   PetscScalar Xqp;  /* q-axis transient reactance */
     36:   PetscScalar Td0p; /* d-axis open circuit time constant */
     37:   PetscScalar Tq0p; /* q-axis open circuit time constant */
     38:   PetscScalar M;    /* M = 2*H/W_S */
     39:   PetscScalar D;    /* D = 0.1*M */
     40:   PetscScalar TM;   /* Mechanical Torque */
     41: } Gen;

     43: typedef struct {
     44:   /* Exciter system constants */
     45:   PetscScalar KA;     /* Voltage regulator gain constant */
     46:   PetscScalar TA;     /* Voltage regulator time constant */
     47:   PetscScalar KE;     /* Exciter gain constant */
     48:   PetscScalar TE;     /* Exciter time constant */
     49:   PetscScalar KF;     /* Feedback stabilizer gain constant */
     50:   PetscScalar TF;     /* Feedback stabilizer time constant */
     51:   PetscScalar k1, k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
     52:   PetscScalar Vref;   /* Voltage regulator voltage setpoint */
     53: } Exc;

     55: typedef struct {
     56:   PetscInt    id;      /* node id */
     57:   PetscInt    nofgen;  /* Number of generators at the bus*/
     58:   PetscInt    nofload; /*  Number of load at the bus*/
     59:   PetscScalar yff[2];  /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/
     60:   PetscScalar vr;      /* Real component of bus voltage */
     61:   PetscScalar vi;      /* Imaginary component of bus voltage */
     62: } Bus;

     64: /* Load constants
     65:   We use a composite load model that describes the load and reactive powers at each time instant as follows
     66:   P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
     67:   Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
     68:   where
     69:     id                  - index of the load
     70:     ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
     71:     ld_alphap,ld_alphap - Percentage contribution (weights) or loads
     72:     P_D0                - Real power load
     73:     Q_D0                - Reactive power load
     74:     Vm(t)              - Voltage magnitude at time t
     75:     Vm0                - Voltage magnitude at t = 0
     76:     ld_betap, ld_betaq  - exponents describing the load model for real and reactive part

     78:     Note: All loads have the same characteristic currently.
     79:   */
     80: typedef struct {
     81:   PetscInt    id; /* bus id */
     82:   PetscInt    ld_nsegsp, ld_nsegsq;
     83:   PetscScalar PD0, QD0;
     84:   PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */
     85:   PetscScalar ld_betap[3], ld_alphaq[3], ld_betaq[3];
     86: } Load;

     88: typedef struct {
     89:   PetscInt    id;     /* node id */
     90:   PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/
     91: } Branch;

     93: typedef struct {
     94:   PetscReal    tfaulton, tfaultoff; /* Fault on and off times */
     95:   PetscReal    t;
     96:   PetscReal    t0, tmax; /* initial time and final time */
     97:   PetscInt     faultbus; /* Fault bus */
     98:   PetscScalar  Rfault;   /* Fault resistance (pu) */
     99:   PetscScalar *ybusfault;
    100:   PetscBool    alg_flg;
    101: } Userctx;

    103: /* Used to read data into the DMNetwork components */
    104: PetscErrorCode read_data(PetscInt nc, Gen **pgen, Exc **pexc, Load **pload, Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
    105: {
    106:   PetscInt           i, j, row[1], col[2];
    107:   PetscInt          *edgelist;
    108:   PetscInt           nofgen[9]  = {1, 1, 1, 0, 0, 0, 0, 0, 0}; /* Buses at which generators are incident */
    109:   PetscInt           nofload[9] = {0, 0, 0, 0, 1, 1, 0, 1, 0}; /* Buses at which loads are incident */
    110:   const PetscScalar *varr;
    111:   PetscScalar        M[3], D[3];
    112:   Bus               *bus;
    113:   Branch            *branch;
    114:   Gen               *gen;
    115:   Exc               *exc;
    116:   Load              *load;
    117:   Mat                Ybus;
    118:   Vec                V0;

    120:   /*10 parameters*/
    121:   /* Generator real and reactive powers (found via loadflow) */
    122:   static const PetscScalar PG[3] = {0.716786142395021, 1.630000000000000, 0.850000000000000};
    123:   static const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};

    125:   /* Generator constants */
    126:   static const PetscScalar H[3]    = {23.64, 6.4, 3.01};       /* Inertia constant */
    127:   static const PetscScalar Rs[3]   = {0.0, 0.0, 0.0};          /* Stator Resistance */
    128:   static const PetscScalar Xd[3]   = {0.146, 0.8958, 1.3125};  /* d-axis reactance */
    129:   static const PetscScalar Xdp[3]  = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
    130:   static const PetscScalar Xq[3]   = {0.4360, 0.8645, 1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
    131:   static const PetscScalar Xqp[3]  = {0.0969, 0.1969, 0.25};   /* q-axis transient reactance */
    132:   static const PetscScalar Td0p[3] = {8.96, 6.0, 5.89};        /* d-axis open circuit time constant */
    133:   static const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6};       /* q-axis open circuit time constant */

    135:   /* Exciter system constants (8 parameters)*/
    136:   static const PetscScalar KA[3] = {20.0, 20.0, 20.0};    /* Voltage regulartor gain constant */
    137:   static const PetscScalar TA[3] = {0.2, 0.2, 0.2};       /* Voltage regulator time constant */
    138:   static const PetscScalar KE[3] = {1.0, 1.0, 1.0};       /* Exciter gain constant */
    139:   static const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
    140:   static const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
    141:   static const PetscScalar TF[3] = {0.35, 0.35, 0.35};    /* Feedback stabilizer time constant */
    142:   static const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
    143:   static const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

    145:   /* Load constants */
    146:   static const PetscScalar PD0[3]       = {1.25, 0.9, 1.0};
    147:   static const PetscScalar QD0[3]       = {0.5, 0.3, 0.35};
    148:   static const PetscScalar ld_alphaq[3] = {1, 0, 0};
    149:   static const PetscScalar ld_betaq[3]  = {2, 1, 0};
    150:   static const PetscScalar ld_betap[3]  = {2, 1, 0};
    151:   static const PetscScalar ld_alphap[3] = {1, 0, 0};
    152:   PetscInt                 ld_nsegsp[3] = {3, 3, 3};
    153:   PetscInt                 ld_nsegsq[3] = {3, 3, 3};
    154:   PetscViewer              Xview, Ybusview;
    155:   PetscInt                 neqs_net, m, n;

    157:   PetscFunctionBeginUser;
    158:   /* Read V0 and Ybus from files */
    159:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "X.bin", FILE_MODE_READ, &Xview));
    160:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "Ybus.bin", FILE_MODE_READ, &Ybusview));
    161:   PetscCall(VecCreate(PETSC_COMM_SELF, &V0));
    162:   PetscCall(VecLoad(V0, Xview));

    164:   PetscCall(MatCreate(PETSC_COMM_SELF, &Ybus));
    165:   PetscCall(MatSetType(Ybus, MATBAIJ));
    166:   PetscCall(MatLoad(Ybus, Ybusview));

    168:   /* Destroy unnecessary stuff */
    169:   PetscCall(PetscViewerDestroy(&Xview));
    170:   PetscCall(PetscViewerDestroy(&Ybusview));

    172:   PetscCall(MatGetLocalSize(Ybus, &m, &n));
    173:   neqs_net = 2 * NBUS; /* # eqs. for network subsystem   */
    174:   PetscCheck(m == neqs_net && n == neqs_net, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix Ybus is in wrong sizes");

    176:   M[0] = 2 * H[0] / W_S;
    177:   M[1] = 2 * H[1] / W_S;
    178:   M[2] = 2 * H[2] / W_S;
    179:   D[0] = 0.1 * M[0];
    180:   D[1] = 0.1 * M[1];
    181:   D[2] = 0.1 * M[2];

    183:   /* Allocate memory for bus, generators, exciter, loads and branches */
    184:   PetscCall(PetscCalloc5(NBUS * nc, &bus, NGEN * nc, &gen, NLOAD * nc, &load, NBRANCH * nc + (nc - 1), &branch, NGEN * nc, &exc));

    186:   PetscCall(VecGetArrayRead(V0, &varr));

    188:   /* read bus data */
    189:   for (i = 0; i < nc; i++) {
    190:     for (j = 0; j < NBUS; j++) {
    191:       bus[i * 9 + j].id      = i * 9 + j;
    192:       bus[i * 9 + j].nofgen  = nofgen[j];
    193:       bus[i * 9 + j].nofload = nofload[j];
    194:       bus[i * 9 + j].vr      = varr[2 * j];
    195:       bus[i * 9 + j].vi      = varr[2 * j + 1];
    196:       row[0]                 = 2 * j;
    197:       col[0]                 = 2 * j;
    198:       col[1]                 = 2 * j + 1;
    199:       /* real and imaginary part of admittance from Ybus into yff */
    200:       PetscCall(MatGetValues(Ybus, 1, row, 2, col, bus[i * 9 + j].yff));
    201:     }
    202:   }

    204:   /* read generator data */
    205:   for (i = 0; i < nc; i++) {
    206:     for (j = 0; j < NGEN; j++) {
    207:       gen[i * 3 + j].id   = i * 3 + j;
    208:       gen[i * 3 + j].PG   = PG[j];
    209:       gen[i * 3 + j].QG   = QG[j];
    210:       gen[i * 3 + j].H    = H[j];
    211:       gen[i * 3 + j].Rs   = Rs[j];
    212:       gen[i * 3 + j].Xd   = Xd[j];
    213:       gen[i * 3 + j].Xdp  = Xdp[j];
    214:       gen[i * 3 + j].Xq   = Xq[j];
    215:       gen[i * 3 + j].Xqp  = Xqp[j];
    216:       gen[i * 3 + j].Td0p = Td0p[j];
    217:       gen[i * 3 + j].Tq0p = Tq0p[j];
    218:       gen[i * 3 + j].M    = M[j];
    219:       gen[i * 3 + j].D    = D[j];
    220:     }
    221:   }

    223:   for (i = 0; i < nc; i++) {
    224:     for (j = 0; j < NGEN; j++) {
    225:       /* exciter system */
    226:       exc[i * 3 + j].KA = KA[j];
    227:       exc[i * 3 + j].TA = TA[j];
    228:       exc[i * 3 + j].KE = KE[j];
    229:       exc[i * 3 + j].TE = TE[j];
    230:       exc[i * 3 + j].KF = KF[j];
    231:       exc[i * 3 + j].TF = TF[j];
    232:       exc[i * 3 + j].k1 = k1[j];
    233:       exc[i * 3 + j].k2 = k2[j];
    234:     }
    235:   }

    237:   /* read load data */
    238:   for (i = 0; i < nc; i++) {
    239:     for (j = 0; j < NLOAD; j++) {
    240:       load[i * 3 + j].id        = i * 3 + j;
    241:       load[i * 3 + j].PD0       = PD0[j];
    242:       load[i * 3 + j].QD0       = QD0[j];
    243:       load[i * 3 + j].ld_nsegsp = ld_nsegsp[j];

    245:       load[i * 3 + j].ld_alphap[0] = ld_alphap[0];
    246:       load[i * 3 + j].ld_alphap[1] = ld_alphap[1];
    247:       load[i * 3 + j].ld_alphap[2] = ld_alphap[2];

    249:       load[i * 3 + j].ld_alphaq[0] = ld_alphaq[0];
    250:       load[i * 3 + j].ld_alphaq[1] = ld_alphaq[1];
    251:       load[i * 3 + j].ld_alphaq[2] = ld_alphaq[2];

    253:       load[i * 3 + j].ld_betap[0] = ld_betap[0];
    254:       load[i * 3 + j].ld_betap[1] = ld_betap[1];
    255:       load[i * 3 + j].ld_betap[2] = ld_betap[2];
    256:       load[i * 3 + j].ld_nsegsq   = ld_nsegsq[j];

    258:       load[i * 3 + j].ld_betaq[0] = ld_betaq[0];
    259:       load[i * 3 + j].ld_betaq[1] = ld_betaq[1];
    260:       load[i * 3 + j].ld_betaq[2] = ld_betaq[2];
    261:     }
    262:   }
    263:   PetscCall(PetscCalloc1(2 * NBRANCH * nc + 2 * (nc - 1), &edgelist));

    265:   /* read edgelist */
    266:   for (i = 0; i < nc; i++) {
    267:     for (j = 0; j < NBRANCH; j++) {
    268:       switch (j) {
    269:       case 0:
    270:         edgelist[i * 18 + 2 * j]     = 0 + 9 * i;
    271:         edgelist[i * 18 + 2 * j + 1] = 3 + 9 * i;
    272:         break;
    273:       case 1:
    274:         edgelist[i * 18 + 2 * j]     = 1 + 9 * i;
    275:         edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
    276:         break;
    277:       case 2:
    278:         edgelist[i * 18 + 2 * j]     = 2 + 9 * i;
    279:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
    280:         break;
    281:       case 3:
    282:         edgelist[i * 18 + 2 * j]     = 3 + 9 * i;
    283:         edgelist[i * 18 + 2 * j + 1] = 4 + 9 * i;
    284:         break;
    285:       case 4:
    286:         edgelist[i * 18 + 2 * j]     = 3 + 9 * i;
    287:         edgelist[i * 18 + 2 * j + 1] = 5 + 9 * i;
    288:         break;
    289:       case 5:
    290:         edgelist[i * 18 + 2 * j]     = 4 + 9 * i;
    291:         edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
    292:         break;
    293:       case 6:
    294:         edgelist[i * 18 + 2 * j]     = 5 + 9 * i;
    295:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
    296:         break;
    297:       case 7:
    298:         edgelist[i * 18 + 2 * j]     = 6 + 9 * i;
    299:         edgelist[i * 18 + 2 * j + 1] = 7 + 9 * i;
    300:         break;
    301:       case 8:
    302:         edgelist[i * 18 + 2 * j]     = 7 + 9 * i;
    303:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
    304:         break;
    305:       default:
    306:         break;
    307:       }
    308:     }
    309:   }

    311:   /* for connecting last bus of previous network(9*i-1) to first bus of next network(9*i), the branch admittance=-0.0301407+j17.3611 */
    312:   for (i = 1; i < nc; i++) {
    313:     edgelist[18 * nc + 2 * (i - 1)]     = 8 + (i - 1) * 9;
    314:     edgelist[18 * nc + 2 * (i - 1) + 1] = 9 * i;

    316:     /* adding admittances to the off-diagonal elements */
    317:     branch[9 * nc + (i - 1)].id     = 9 * nc + (i - 1);
    318:     branch[9 * nc + (i - 1)].yft[0] = 17.3611;
    319:     branch[9 * nc + (i - 1)].yft[1] = -0.0301407;

    321:     /* subtracting admittances from the diagonal elements */
    322:     bus[9 * i - 1].yff[0] -= 17.3611;
    323:     bus[9 * i - 1].yff[1] -= -0.0301407;

    325:     bus[9 * i].yff[0] -= 17.3611;
    326:     bus[9 * i].yff[1] -= -0.0301407;
    327:   }

    329:   /* read branch data */
    330:   for (i = 0; i < nc; i++) {
    331:     for (j = 0; j < NBRANCH; j++) {
    332:       branch[i * 9 + j].id = i * 9 + j;

    334:       row[0] = edgelist[2 * j] * 2;
    335:       col[0] = edgelist[2 * j + 1] * 2;
    336:       col[1] = edgelist[2 * j + 1] * 2 + 1;
    337:       PetscCall(MatGetValues(Ybus, 1, row, 2, col, branch[i * 9 + j].yft)); /*imaginary part of admittance*/
    338:     }
    339:   }

    341:   *pgen      = gen;
    342:   *pexc      = exc;
    343:   *pload     = load;
    344:   *pbus      = bus;
    345:   *pbranch   = branch;
    346:   *pedgelist = edgelist;

    348:   PetscCall(VecRestoreArrayRead(V0, &varr));

    350:   /* Destroy unnecessary stuff */
    351:   PetscCall(MatDestroy(&Ybus));
    352:   PetscCall(VecDestroy(&V0));
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

    356: PetscErrorCode SetInitialGuess(DM networkdm, Vec X)
    357: {
    358:   Bus         *bus;
    359:   Gen         *gen;
    360:   Exc         *exc;
    361:   PetscInt     v, vStart, vEnd, offset;
    362:   PetscInt     key, numComps, j;
    363:   PetscBool    ghostvtex;
    364:   Vec          localX;
    365:   PetscScalar *xarr;
    366:   PetscScalar  Vr = 0, Vi = 0, Vm = 0, Vm2; /* Terminal voltage variables */
    367:   PetscScalar  IGr, IGi;                    /* Generator real and imaginary current */
    368:   PetscScalar  Eqp, Edp, delta;             /* Generator variables */
    369:   PetscScalar  Efd = 0, RF, VR;             /* Exciter variables */
    370:   PetscScalar  Vd, Vq;                      /* Generator dq axis voltages */
    371:   PetscScalar  Id, Iq;                      /* Generator dq axis currents */
    372:   PetscScalar  theta;                       /* Generator phase angle */
    373:   PetscScalar  SE;
    374:   void        *component;

    376:   PetscFunctionBegin;
    377:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
    378:   PetscCall(DMGetLocalVector(networkdm, &localX));

    380:   PetscCall(VecSet(X, 0.0));
    381:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
    382:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

    384:   PetscCall(VecGetArray(localX, &xarr));

    386:   for (v = vStart; v < vEnd; v++) {
    387:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
    388:     if (ghostvtex) continue;

    390:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
    391:     for (j = 0; j < numComps; j++) {
    392:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
    393:       if (key == 1) {
    394:         bus = (Bus *)(component);

    396:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
    397:         xarr[offset]     = bus->vr;
    398:         xarr[offset + 1] = bus->vi;

    400:         Vr = bus->vr;
    401:         Vi = bus->vi;
    402:       } else if (key == 2) {
    403:         gen = (Gen *)(component);
    404:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
    405:         Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    406:         Vm2 = Vm * Vm;
    407:         /* Real part of gen current */
    408:         IGr = (Vr * gen->PG + Vi * gen->QG) / Vm2;
    409:         /* Imaginary part of gen current */
    410:         IGi = (Vi * gen->PG - Vr * gen->QG) / Vm2;

    412:         /* Machine angle */
    413:         delta = atan2(Vi + gen->Xq * IGr, Vr - gen->Xq * IGi);
    414:         theta = PETSC_PI / 2.0 - delta;

    416:         /* d-axis stator current */
    417:         Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta);

    419:         /* q-axis stator current */
    420:         Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta);

    422:         Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
    423:         Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

    425:         /* d-axis transient EMF */
    426:         Edp = Vd + gen->Rs * Id - gen->Xqp * Iq;

    428:         /* q-axis transient EMF */
    429:         Eqp = Vq + gen->Rs * Iq + gen->Xdp * Id;

    431:         gen->TM = gen->PG;

    433:         xarr[offset]     = Eqp;
    434:         xarr[offset + 1] = Edp;
    435:         xarr[offset + 2] = delta;
    436:         xarr[offset + 3] = W_S;
    437:         xarr[offset + 4] = Id;
    438:         xarr[offset + 5] = Iq;

    440:         Efd = Eqp + (gen->Xd - gen->Xdp) * Id;

    442:       } else if (key == 3) {
    443:         exc = (Exc *)(component);
    444:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));

    446:         SE = exc->k1 * PetscExpScalar(exc->k2 * Efd);
    447:         VR = exc->KE * Efd + SE;
    448:         RF = exc->KF * Efd / exc->TF;

    450:         xarr[offset]     = Efd;
    451:         xarr[offset + 1] = RF;
    452:         xarr[offset + 2] = VR;

    454:         exc->Vref = Vm + (VR / exc->KA);
    455:       }
    456:     }
    457:   }
    458:   PetscCall(VecRestoreArray(localX, &xarr));
    459:   PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
    460:   PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
    461:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
    462:   PetscFunctionReturn(PETSC_SUCCESS);
    463: }

    465: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
    466: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
    467: {
    468:   PetscFunctionBegin;
    469:   *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
    470:   *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
    471:   PetscFunctionReturn(PETSC_SUCCESS);
    472: }

    474: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
    475: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
    476: {
    477:   PetscFunctionBegin;
    478:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
    479:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
    480:   PetscFunctionReturn(PETSC_SUCCESS);
    481: }

    483: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
    484: PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
    485: {
    486:   DM                 networkdm;
    487:   Vec                localX, localXdot, localF;
    488:   PetscInt           vfrom, vto, offsetfrom, offsetto;
    489:   PetscInt           v, vStart, vEnd, e;
    490:   PetscScalar       *farr;
    491:   PetscScalar        Vd = 0, Vq = 0, SE;
    492:   const PetscScalar *xarr, *xdotarr;
    493:   void              *component;
    494:   PetscScalar        Vr = 0, Vi = 0;

    496:   PetscFunctionBegin;
    497:   PetscCall(VecSet(F, 0.0));

    499:   PetscCall(TSGetDM(ts, &networkdm));
    500:   PetscCall(DMGetLocalVector(networkdm, &localF));
    501:   PetscCall(DMGetLocalVector(networkdm, &localX));
    502:   PetscCall(DMGetLocalVector(networkdm, &localXdot));
    503:   PetscCall(VecSet(localF, 0.0));

    505:   /* update ghost values of localX and localXdot */
    506:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
    507:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

    509:   PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
    510:   PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));

    512:   PetscCall(VecGetArrayRead(localX, &xarr));
    513:   PetscCall(VecGetArrayRead(localXdot, &xdotarr));
    514:   PetscCall(VecGetArray(localF, &farr));

    516:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

    518:   for (v = vStart; v < vEnd; v++) {
    519:     PetscInt    i, j, offsetbus, offsetgen, offsetexc, key;
    520:     Bus        *bus;
    521:     Gen        *gen;
    522:     Exc        *exc;
    523:     Load       *load;
    524:     PetscBool   ghostvtex;
    525:     PetscInt    numComps;
    526:     PetscScalar Yffr, Yffi; /* Real and imaginary fault admittances */
    527:     PetscScalar Vm, Vm2, Vm0;
    528:     PetscScalar Vr0 = 0, Vi0 = 0;
    529:     PetscScalar PD, QD;

    531:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
    532:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));

    534:     for (j = 0; j < numComps; j++) {
    535:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
    536:       if (key == 1) {
    537:         PetscInt        nconnedges;
    538:         const PetscInt *connedges;

    540:         bus = (Bus *)(component);
    541:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
    542:         if (!ghostvtex) {
    543:           Vr = xarr[offsetbus];
    544:           Vi = xarr[offsetbus + 1];

    546:           Yffr = bus->yff[1];
    547:           Yffi = bus->yff[0];

    549:           if (user->alg_flg) {
    550:             Yffr += user->ybusfault[bus->id * 2 + 1];
    551:             Yffi += user->ybusfault[bus->id * 2];
    552:           }
    553:           Vr0 = bus->vr;
    554:           Vi0 = bus->vi;

    556:           /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
    557:            The generator current injection, IG, and load current injection, ID are added later
    558:            */
    559:           farr[offsetbus] += Yffi * Vr + Yffr * Vi;     /* imaginary current due to diagonal elements */
    560:           farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi; /* real current due to diagonal elements */
    561:         }

    563:         PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));

    565:         for (i = 0; i < nconnedges; i++) {
    566:           Branch         *branch;
    567:           PetscInt        keye;
    568:           PetscScalar     Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
    569:           const PetscInt *cone;

    571:           e = connedges[i];
    572:           PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));

    574:           Yfti = branch->yft[0];
    575:           Yftr = branch->yft[1];

    577:           PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));

    579:           vfrom = cone[0];
    580:           vto   = cone[1];

    582:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
    583:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));

    585:           /* From bus and to bus real and imaginary voltages */
    586:           Vfr = xarr[offsetfrom];
    587:           Vfi = xarr[offsetfrom + 1];
    588:           Vtr = xarr[offsetto];
    589:           Vti = xarr[offsetto + 1];

    591:           if (vfrom == v) {
    592:             farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
    593:             farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
    594:           } else {
    595:             farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
    596:             farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
    597:           }
    598:         }
    599:       } else if (key == 2) {
    600:         if (!ghostvtex) {
    601:           PetscScalar Eqp, Edp, delta, w; /* Generator variables */
    602:           PetscScalar Efd;                /* Exciter field voltage */
    603:           PetscScalar Id, Iq;             /* Generator dq axis currents */
    604:           PetscScalar IGr, IGi, Zdq_inv[4], det;
    605:           PetscScalar Xd, Xdp, Td0p, Xq, Xqp, Tq0p, TM, D, M, Rs; /* Generator parameters */

    607:           gen = (Gen *)(component);
    608:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));

    610:           /* Generator state variables */
    611:           Eqp   = xarr[offsetgen];
    612:           Edp   = xarr[offsetgen + 1];
    613:           delta = xarr[offsetgen + 2];
    614:           w     = xarr[offsetgen + 3];
    615:           Id    = xarr[offsetgen + 4];
    616:           Iq    = xarr[offsetgen + 5];

    618:           /* Generator parameters */
    619:           Xd   = gen->Xd;
    620:           Xdp  = gen->Xdp;
    621:           Td0p = gen->Td0p;
    622:           Xq   = gen->Xq;
    623:           Xqp  = gen->Xqp;
    624:           Tq0p = gen->Tq0p;
    625:           TM   = gen->TM;
    626:           D    = gen->D;
    627:           M    = gen->M;
    628:           Rs   = gen->Rs;

    630:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, 2, &offsetexc));
    631:           Efd = xarr[offsetexc];

    633:           /* Generator differential equations */
    634:           farr[offsetgen]     = (Eqp + (Xd - Xdp) * Id - Efd) / Td0p + xdotarr[offsetgen];
    635:           farr[offsetgen + 1] = (Edp - (Xq - Xqp) * Iq) / Tq0p + xdotarr[offsetgen + 1];
    636:           farr[offsetgen + 2] = -w + W_S + xdotarr[offsetgen + 2];
    637:           farr[offsetgen + 3] = (-TM + Edp * Id + Eqp * Iq + (Xqp - Xdp) * Id * Iq + D * (w - W_S)) / M + xdotarr[offsetgen + 3];

    639:           PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

    641:           /* Algebraic equations for stator currents */
    642:           det = Rs * Rs + Xdp * Xqp;

    644:           Zdq_inv[0] = Rs / det;
    645:           Zdq_inv[1] = Xqp / det;
    646:           Zdq_inv[2] = -Xdp / det;
    647:           Zdq_inv[3] = Rs / det;

    649:           farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
    650:           farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

    652:           PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));

    654:           /* Add generator current injection to network */
    655:           farr[offsetbus] -= IGi;
    656:           farr[offsetbus + 1] -= IGr;
    657:         }
    658:       } else if (key == 3) {
    659:         if (!ghostvtex) {
    660:           PetscScalar k1, k2, KE, TE, TF, KA, KF, Vref, TA; /* Generator parameters */
    661:           PetscScalar Efd, RF, VR;                          /* Exciter variables */

    663:           exc = (Exc *)(component);
    664:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));

    666:           Efd = xarr[offsetexc];
    667:           RF  = xarr[offsetexc + 1];
    668:           VR  = xarr[offsetexc + 2];

    670:           k1   = exc->k1;
    671:           k2   = exc->k2;
    672:           KE   = exc->KE;
    673:           TE   = exc->TE;
    674:           TF   = exc->TF;
    675:           KA   = exc->KA;
    676:           KF   = exc->KF;
    677:           Vref = exc->Vref;
    678:           TA   = exc->TA;

    680:           Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
    681:           SE = k1 * PetscExpScalar(k2 * Efd);

    683:           /* Exciter differential equations */
    684:           farr[offsetexc]     = (KE * Efd + SE - VR) / TE + xdotarr[offsetexc];
    685:           farr[offsetexc + 1] = (RF - KF * Efd / TF) / TF + xdotarr[offsetexc + 1];
    686:           farr[offsetexc + 2] = (VR - KA * RF + KA * KF * Efd / TF - KA * (Vref - Vm)) / TA + xdotarr[offsetexc + 2];
    687:         }
    688:       } else if (key == 4) {
    689:         if (!ghostvtex) {
    690:           PetscInt     k;
    691:           PetscInt     ld_nsegsp;
    692:           PetscInt     ld_nsegsq;
    693:           PetscScalar *ld_alphap;
    694:           PetscScalar *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;

    696:           load = (Load *)(component);

    698:           /* Load Parameters */
    699:           ld_nsegsp = load->ld_nsegsp;
    700:           ld_alphap = load->ld_alphap;
    701:           ld_betap  = load->ld_betap;
    702:           ld_nsegsq = load->ld_nsegsq;
    703:           ld_alphaq = load->ld_alphaq;
    704:           ld_betaq  = load->ld_betaq;
    705:           PD0       = load->PD0;
    706:           QD0       = load->QD0;

    708:           Vr  = xarr[offsetbus];     /* Real part of generator terminal voltage */
    709:           Vi  = xarr[offsetbus + 1]; /* Imaginary part of the generator terminal voltage */
    710:           Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    711:           Vm2 = Vm * Vm;
    712:           Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
    713:           PD = QD = 0.0;
    714:           for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
    715:           for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);

    717:           /* Load currents */
    718:           IDr = (PD * Vr + QD * Vi) / Vm2;
    719:           IDi = (-QD * Vr + PD * Vi) / Vm2;

    721:           /* Load current contribution to the network */
    722:           farr[offsetbus] += IDi;
    723:           farr[offsetbus + 1] += IDr;
    724:         }
    725:       }
    726:     }
    727:   }

    729:   PetscCall(VecRestoreArrayRead(localX, &xarr));
    730:   PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
    731:   PetscCall(VecRestoreArray(localF, &farr));
    732:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
    733:   PetscCall(DMRestoreLocalVector(networkdm, &localXdot));

    735:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
    736:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
    737:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
    738:   PetscFunctionReturn(PETSC_SUCCESS);
    739: }

    741: /* This function is used for solving the algebraic system only during fault on and
    742:    off times. It computes the entire F and then zeros out the part corresponding to
    743:    differential equations
    744:  F = [0;g(y)];
    745: */
    746: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
    747: {
    748:   DM                 networkdm;
    749:   Vec                localX, localF;
    750:   PetscInt           vfrom, vto, offsetfrom, offsetto;
    751:   PetscInt           v, vStart, vEnd, e;
    752:   PetscScalar       *farr;
    753:   Userctx           *user = (Userctx *)ctx;
    754:   const PetscScalar *xarr;
    755:   void              *component;
    756:   PetscScalar        Vr = 0, Vi = 0;

    758:   PetscFunctionBegin;
    759:   PetscCall(VecSet(F, 0.0));
    760:   PetscCall(SNESGetDM(snes, &networkdm));
    761:   PetscCall(DMGetLocalVector(networkdm, &localF));
    762:   PetscCall(DMGetLocalVector(networkdm, &localX));
    763:   PetscCall(VecSet(localF, 0.0));

    765:   /* update ghost values of locaX and locaXdot */
    766:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
    767:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

    769:   PetscCall(VecGetArrayRead(localX, &xarr));
    770:   PetscCall(VecGetArray(localF, &farr));

    772:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

    774:   for (v = vStart; v < vEnd; v++) {
    775:     PetscInt    i, j, offsetbus, offsetgen, key, numComps;
    776:     PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0 = 0, Vi0 = 0, PD, QD;
    777:     Bus        *bus;
    778:     Gen        *gen;
    779:     Load       *load;
    780:     PetscBool   ghostvtex;

    782:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
    783:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));

    785:     for (j = 0; j < numComps; j++) {
    786:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
    787:       if (key == 1) {
    788:         PetscInt        nconnedges;
    789:         const PetscInt *connedges;

    791:         bus = (Bus *)(component);
    792:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
    793:         if (!ghostvtex) {
    794:           Vr = xarr[offsetbus];
    795:           Vi = xarr[offsetbus + 1];

    797:           Yffr = bus->yff[1];
    798:           Yffi = bus->yff[0];
    799:           if (user->alg_flg) {
    800:             Yffr += user->ybusfault[bus->id * 2 + 1];
    801:             Yffi += user->ybusfault[bus->id * 2];
    802:           }
    803:           Vr0 = bus->vr;
    804:           Vi0 = bus->vi;

    806:           farr[offsetbus] += Yffi * Vr + Yffr * Vi;
    807:           farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi;
    808:         }
    809:         PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));

    811:         for (i = 0; i < nconnedges; i++) {
    812:           Branch         *branch;
    813:           PetscInt        keye;
    814:           PetscScalar     Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
    815:           const PetscInt *cone;

    817:           e = connedges[i];
    818:           PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));

    820:           Yfti = branch->yft[0];
    821:           Yftr = branch->yft[1];

    823:           PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
    824:           vfrom = cone[0];
    825:           vto   = cone[1];

    827:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
    828:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));

    830:           /*From bus and to bus real and imaginary voltages */
    831:           Vfr = xarr[offsetfrom];
    832:           Vfi = xarr[offsetfrom + 1];
    833:           Vtr = xarr[offsetto];
    834:           Vti = xarr[offsetto + 1];

    836:           if (vfrom == v) {
    837:             farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
    838:             farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
    839:           } else {
    840:             farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
    841:             farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
    842:           }
    843:         }
    844:       } else if (key == 2) {
    845:         if (!ghostvtex) {
    846:           PetscScalar Eqp, Edp, delta; /* Generator variables */
    847:           PetscScalar Id, Iq;          /* Generator dq axis currents */
    848:           PetscScalar Vd, Vq, IGr, IGi, Zdq_inv[4], det;
    849:           PetscScalar Xdp, Xqp, Rs; /* Generator parameters */

    851:           gen = (Gen *)(component);
    852:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));

    854:           /* Generator state variables */
    855:           Eqp   = xarr[offsetgen];
    856:           Edp   = xarr[offsetgen + 1];
    857:           delta = xarr[offsetgen + 2];
    858:           /* w     = xarr[idx+3]; not being used */
    859:           Id = xarr[offsetgen + 4];
    860:           Iq = xarr[offsetgen + 5];

    862:           /* Generator parameters */
    863:           Xdp = gen->Xdp;
    864:           Xqp = gen->Xqp;
    865:           Rs  = gen->Rs;

    867:           /* Set generator differential equation residual functions to zero */
    868:           farr[offsetgen]     = 0;
    869:           farr[offsetgen + 1] = 0;
    870:           farr[offsetgen + 2] = 0;
    871:           farr[offsetgen + 3] = 0;

    873:           PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

    875:           /* Algebraic equations for stator currents */
    876:           det = Rs * Rs + Xdp * Xqp;

    878:           Zdq_inv[0] = Rs / det;
    879:           Zdq_inv[1] = Xqp / det;
    880:           Zdq_inv[2] = -Xdp / det;
    881:           Zdq_inv[3] = Rs / det;

    883:           farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
    884:           farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

    886:           /* Add generator current injection to network */
    887:           PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));

    889:           farr[offsetbus] -= IGi;
    890:           farr[offsetbus + 1] -= IGr;

    892:           /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
    893:         }
    894:       } else if (key == 3) {
    895:         if (!ghostvtex) {
    896:           PetscInt offsetexc;
    897:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
    898:           /* Set exciter differential equation residual functions equal to zero*/
    899:           farr[offsetexc]     = 0;
    900:           farr[offsetexc + 1] = 0;
    901:           farr[offsetexc + 2] = 0;
    902:         }
    903:       } else if (key == 4) {
    904:         if (!ghostvtex) {
    905:           PetscInt     k, ld_nsegsp, ld_nsegsq;
    906:           PetscScalar *ld_alphap, *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;

    908:           load = (Load *)(component);

    910:           /* Load Parameters */
    911:           ld_nsegsp = load->ld_nsegsp;
    912:           ld_alphap = load->ld_alphap;
    913:           ld_betap  = load->ld_betap;
    914:           ld_nsegsq = load->ld_nsegsq;
    915:           ld_alphaq = load->ld_alphaq;
    916:           ld_betaq  = load->ld_betaq;

    918:           PD0 = load->PD0;
    919:           QD0 = load->QD0;

    921:           Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    922:           Vm2 = Vm * Vm;
    923:           Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
    924:           PD = QD = 0.0;
    925:           for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
    926:           for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);

    928:           /* Load currents */
    929:           IDr = (PD * Vr + QD * Vi) / Vm2;
    930:           IDi = (-QD * Vr + PD * Vi) / Vm2;

    932:           farr[offsetbus] += IDi;
    933:           farr[offsetbus + 1] += IDr;
    934:         }
    935:       }
    936:     }
    937:   }

    939:   PetscCall(VecRestoreArrayRead(localX, &xarr));
    940:   PetscCall(VecRestoreArray(localF, &farr));
    941:   PetscCall(DMRestoreLocalVector(networkdm, &localX));

    943:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
    944:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
    945:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
    946:   PetscFunctionReturn(PETSC_SUCCESS);
    947: }

    949: int main(int argc, char **argv)
    950: {
    951:   PetscInt      i, j, *edgelist = NULL, eStart, eEnd, vStart, vEnd;
    952:   PetscInt      genj, excj, loadj, componentkey[5];
    953:   PetscInt      nc = 1; /* No. of copies (default = 1) */
    954:   PetscMPIInt   size, rank;
    955:   Vec           X, F_alg;
    956:   TS            ts;
    957:   SNES          snes_alg, snes;
    958:   Bus          *bus;
    959:   Branch       *branch;
    960:   Gen          *gen;
    961:   Exc          *exc;
    962:   Load         *load;
    963:   DM            networkdm;
    964:   PetscLogStage stage1;
    965:   Userctx       user;
    966:   KSP           ksp;
    967:   PC            pc;
    968:   PetscInt      numEdges = 0;

    970:   PetscFunctionBeginUser;
    971:   PetscCall(PetscInitialize(&argc, &argv, "ex9busnetworkops", help));
    972:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
    973:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    974:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

    976:   /* Read initial voltage vector and Ybus */
    977:   if (rank == 0) PetscCall(read_data(nc, &gen, &exc, &load, &bus, &branch, &edgelist));

    979:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
    980:   PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(Branch), &componentkey[0]));
    981:   PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(Bus), &componentkey[1]));
    982:   PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(Gen), &componentkey[2]));
    983:   PetscCall(DMNetworkRegisterComponent(networkdm, "excstruct", sizeof(Exc), &componentkey[3]));
    984:   PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(Load), &componentkey[4]));

    986:   PetscCall(PetscLogStageRegister("Create network", &stage1));
    987:   PetscCall(PetscLogStagePush(stage1));

    989:   /* Set local number of edges and edge connectivity */
    990:   if (rank == 0) numEdges = NBRANCH * nc + (nc - 1);
    991:   PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
    992:   PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, numEdges, edgelist, NULL));

    994:   /* Set up the network layout */
    995:   PetscCall(DMNetworkLayoutSetUp(networkdm));

    997:   if (rank == 0) PetscCall(PetscFree(edgelist));

    999:   /* Add network components (physical parameters of nodes and branches) and number of variables */
    1000:   if (rank == 0) {
    1001:     PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    1002:     genj  = 0;
    1003:     loadj = 0;
    1004:     excj  = 0;
    1005:     for (i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[0], &branch[i - eStart], 0));

    1007:     PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

    1009:     for (i = vStart; i < vEnd; i++) {
    1010:       PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[1], &bus[i - vStart], 2));
    1011:       if (bus[i - vStart].nofgen) {
    1012:         for (j = 0; j < bus[i - vStart].nofgen; j++) {
    1013:           /* Add generator */
    1014:           PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[2], &gen[genj++], 6));
    1015:           /* Add exciter */
    1016:           PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[3], &exc[excj++], 3));
    1017:         }
    1018:       }
    1019:       if (bus[i - vStart].nofload) {
    1020:         for (j = 0; j < bus[i - vStart].nofload; j++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[4], &load[loadj++], 0));
    1021:       }
    1022:     }
    1023:   }

    1025:   PetscCall(DMSetUp(networkdm));

    1027:   if (rank == 0) PetscCall(PetscFree5(bus, gen, load, branch, exc));

    1029:   /* for parallel options: Network partitioning and distribution of data */
    1030:   if (size > 1) PetscCall(DMNetworkDistribute(&networkdm, 0));
    1031:   PetscCall(PetscLogStagePop());

    1033:   PetscCall(DMCreateGlobalVector(networkdm, &X));

    1035:   PetscCall(SetInitialGuess(networkdm, X));

    1037:   /* Options for fault simulation */
    1038:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
    1039:   user.tfaulton  = 0.02;
    1040:   user.tfaultoff = 0.05;
    1041:   user.Rfault    = 0.0001;
    1042:   user.faultbus  = 8;
    1043:   PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
    1044:   PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
    1045:   PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
    1046:   user.t0   = 0.0;
    1047:   user.tmax = 0.1;
    1048:   PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
    1049:   PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));

    1051:   PetscCall(PetscMalloc1(18 * nc, &user.ybusfault));
    1052:   for (i = 0; i < 18 * nc; i++) user.ybusfault[i] = 0;
    1053:   user.ybusfault[user.faultbus * 2 + 1] = 1 / user.Rfault;
    1054:   PetscOptionsEnd();

    1056:   /* Setup TS solver                                           */
    1057:   /*--------------------------------------------------------*/
    1058:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    1059:   PetscCall(TSSetDM(ts, (DM)networkdm));
    1060:   PetscCall(TSSetType(ts, TSCN));

    1062:   PetscCall(TSGetSNES(ts, &snes));
    1063:   PetscCall(SNESGetKSP(snes, &ksp));
    1064:   PetscCall(KSPGetPC(ksp, &pc));
    1065:   PetscCall(PCSetType(pc, PCBJACOBI));

    1067:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
    1068:   PetscCall(TSSetMaxTime(ts, user.tfaulton));
    1069:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1070:   PetscCall(TSSetTimeStep(ts, 0.01));
    1071:   PetscCall(TSSetFromOptions(ts));

    1073:   /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
    1074:     eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
    1075:   user.alg_flg = PETSC_FALSE;

    1077:   /* Prefault period */
    1078:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "... (1) Prefault period ... \n"));

    1080:   PetscCall(TSSetSolution(ts, X));
    1081:   PetscCall(TSSetUp(ts));
    1082:   PetscCall(TSSolve(ts, X));

    1084:   /* Create the nonlinear solver for solving the algebraic system */
    1085:   PetscCall(VecDuplicate(X, &F_alg));

    1087:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
    1088:   PetscCall(SNESSetDM(snes_alg, (DM)networkdm));
    1089:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
    1090:   PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
    1091:   PetscCall(SNESSetFromOptions(snes_alg));

    1093:   /* Apply disturbance - resistive fault at user.faultbus */
    1094:   /* This is done by adding shunt conductance to the diagonal location
    1095:      in the Ybus matrix */
    1096:   user.alg_flg = PETSC_TRUE;

    1098:   /* Solve the algebraic equations */
    1099:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (2) Apply disturbance, solve algebraic equations ... \n"));
    1100:   PetscCall(SNESSolve(snes_alg, NULL, X));

    1102:   /* Disturbance period */
    1103:   PetscCall(TSSetTime(ts, user.tfaulton));
    1104:   PetscCall(TSSetMaxTime(ts, user.tfaultoff));
    1105:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1106:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));

    1108:   user.alg_flg = PETSC_TRUE;
    1109:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (3) Disturbance period ... \n"));
    1110:   PetscCall(TSSolve(ts, X));

    1112:   /* Remove the fault */
    1113:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));

    1115:   user.alg_flg = PETSC_FALSE;
    1116:   /* Solve the algebraic equations */
    1117:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (4) Remove fault, solve algebraic equations ... \n"));
    1118:   PetscCall(SNESSolve(snes_alg, NULL, X));
    1119:   PetscCall(SNESDestroy(&snes_alg));

    1121:   /* Post-disturbance period */
    1122:   PetscCall(TSSetTime(ts, user.tfaultoff));
    1123:   PetscCall(TSSetMaxTime(ts, user.tmax));
    1124:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    1125:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));

    1127:   user.alg_flg = PETSC_FALSE;
    1128:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (5) Post-disturbance period ... \n"));
    1129:   PetscCall(TSSolve(ts, X));

    1131:   PetscCall(PetscFree(user.ybusfault));
    1132:   PetscCall(VecDestroy(&F_alg));
    1133:   PetscCall(VecDestroy(&X));
    1134:   PetscCall(DMDestroy(&networkdm));
    1135:   PetscCall(TSDestroy(&ts));
    1136:   PetscCall(PetscFinalize());
    1137:   return 0;
    1138: }

    1140: /*TEST

    1142:    build:
    1143:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

    1145:    test:
    1146:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
    1147:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

    1149:    test:
    1150:       suffix: 2
    1151:       nsize: 2
    1152:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
    1153:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

    1155: TEST*/