Actual source code: ex9busopt.c

      1: static char help[] = "Application of adjoint sensitivity analysis for power grid stability analysis of WECC 9 bus system.\n\
      2: This example is based on the 9-bus (node) example given in the book Power\n\
      3: Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
      4: The power grid in this example consists of 9 buses (nodes), 3 generators,\n\
      5: 3 loads, and 9 transmission lines. The network equations are written\n\
      6: in current balance form using rectangular coordinates.\n\n";

      8: /*
      9:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
     10:   The objectivie is to find optimal parameter PG for each generator to minizie the frequency violations due to faults.
     11:   The problem features discontinuities and a cost function in integral form.
     12:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex9busadj.c for details.
     13: */

     15: #include <petsctao.h>
     16: #include <petscts.h>
     17: #include <petscdm.h>
     18: #include <petscdmda.h>
     19: #include <petscdmcomposite.h>
     20: #include <petsctime.h>

     22: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);

     24: #define freq 60
     25: #define w_s  (2 * PETSC_PI * freq)

     27: /* Sizes and indices */
     28: const PetscInt nbus    = 9;         /* Number of network buses */
     29: const PetscInt ngen    = 3;         /* Number of generators */
     30: const PetscInt nload   = 3;         /* Number of loads */
     31: const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
     32: const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */

     34: /* Generator real and reactive powers (found via loadflow) */
     35: PetscScalar PG[3] = {0.69, 1.59, 0.69};
     36: /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/

     38: const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
     39: /* Generator constants */
     40: const PetscScalar H[3]    = {23.64, 6.4, 3.01};       /* Inertia constant */
     41: const PetscScalar Rs[3]   = {0.0, 0.0, 0.0};          /* Stator Resistance */
     42: const PetscScalar Xd[3]   = {0.146, 0.8958, 1.3125};  /* d-axis reactance */
     43: const PetscScalar Xdp[3]  = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
     44: 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 */
     45: const PetscScalar Xqp[3]  = {0.0969, 0.1969, 0.25};   /* q-axis transient reactance */
     46: const PetscScalar Td0p[3] = {8.96, 6.0, 5.89};        /* d-axis open circuit time constant */
     47: const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6};       /* q-axis open circuit time constant */
     48: PetscScalar       M[3];                               /* M = 2*H/w_s */
     49: PetscScalar       D[3];                               /* D = 0.1*M */

     51: PetscScalar TM[3]; /* Mechanical Torque */
     52: /* Exciter system constants */
     53: const PetscScalar KA[3] = {20.0, 20.0, 20.0};    /* Voltage regulartor gain constant */
     54: const PetscScalar TA[3] = {0.2, 0.2, 0.2};       /* Voltage regulator time constant */
     55: const PetscScalar KE[3] = {1.0, 1.0, 1.0};       /* Exciter gain constant */
     56: const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
     57: const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
     58: const PetscScalar TF[3] = {0.35, 0.35, 0.35};    /* Feedback stabilizer time constant */
     59: const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
     60: const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

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

     76:     Note: All loads have the same characteristic currently.
     77: */
     78: const PetscScalar PD0[3]       = {1.25, 0.9, 1.0};
     79: const PetscScalar QD0[3]       = {0.5, 0.3, 0.35};
     80: const PetscInt    ld_nsegsp[3] = {3, 3, 3};
     81: const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
     82: const PetscScalar ld_betap[3]  = {2.0, 1.0, 0.0};
     83: const PetscInt    ld_nsegsq[3] = {3, 3, 3};
     84: const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
     85: const PetscScalar ld_betaq[3]  = {2.0, 1.0, 0.0};

     87: typedef struct {
     88:   DM          dmgen, dmnet;        /* DMs to manage generator and network subsystem */
     89:   DM          dmpgrid;             /* Composite DM to manage the entire power grid */
     90:   Mat         Ybus;                /* Network admittance matrix */
     91:   Vec         V0;                  /* Initial voltage vector (Power flow solution) */
     92:   PetscReal   tfaulton, tfaultoff; /* Fault on and off times */
     93:   PetscInt    faultbus;            /* Fault bus */
     94:   PetscScalar Rfault;
     95:   PetscReal   t0, tmax;
     96:   PetscInt    neqs_gen, neqs_net, neqs_pgrid;
     97:   Mat         Sol; /* Matrix to save solution at each time step */
     98:   PetscInt    stepnum;
     99:   PetscBool   alg_flg;
    100:   PetscReal   t;
    101:   IS          is_diff;        /* indices for differential equations */
    102:   IS          is_alg;         /* indices for algebraic equations */
    103:   PetscReal   freq_u, freq_l; /* upper and lower frequency limit */
    104:   PetscInt    pow;            /* power coefficient used in the cost function */
    105:   PetscBool   jacp_flg;
    106:   Mat         J, Jacp;
    107:   Mat         DRDU, DRDP;
    108: } Userctx;

    110: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
    111: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
    112: {
    113:   PetscFunctionBegin;
    114:   *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
    115:   *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
    120: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
    121: {
    122:   PetscFunctionBegin;
    123:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
    124:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
    125:   PetscFunctionReturn(PETSC_SUCCESS);
    126: }

    128: /* Saves the solution at each time to a matrix */
    129: PetscErrorCode SaveSolution(TS ts)
    130: {
    131:   Userctx           *user;
    132:   Vec                X;
    133:   PetscScalar       *mat;
    134:   const PetscScalar *x;
    135:   PetscInt           idx;
    136:   PetscReal          t;

    138:   PetscFunctionBegin;
    139:   PetscCall(TSGetApplicationContext(ts, &user));
    140:   PetscCall(TSGetTime(ts, &t));
    141:   PetscCall(TSGetSolution(ts, &X));
    142:   idx = user->stepnum * (user->neqs_pgrid + 1);
    143:   PetscCall(MatDenseGetArray(user->Sol, &mat));
    144:   PetscCall(VecGetArrayRead(X, &x));
    145:   mat[idx] = t;
    146:   PetscCall(PetscArraycpy(mat + idx + 1, x, user->neqs_pgrid));
    147:   PetscCall(MatDenseRestoreArray(user->Sol, &mat));
    148:   PetscCall(VecRestoreArrayRead(X, &x));
    149:   user->stepnum++;
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    151: }

    153: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
    154: {
    155:   Vec          Xgen, Xnet;
    156:   PetscScalar *xgen, *xnet;
    157:   PetscInt     i, idx = 0;
    158:   PetscScalar  Vr, Vi, IGr, IGi, Vm, Vm2;
    159:   PetscScalar  Eqp, Edp, delta;
    160:   PetscScalar  Efd, RF, VR; /* Exciter variables */
    161:   PetscScalar  Id, Iq;      /* Generator dq axis currents */
    162:   PetscScalar  theta, Vd, Vq, SE;

    164:   PetscFunctionBegin;
    165:   M[0] = 2 * H[0] / w_s;
    166:   M[1] = 2 * H[1] / w_s;
    167:   M[2] = 2 * H[2] / w_s;
    168:   D[0] = 0.1 * M[0];
    169:   D[1] = 0.1 * M[1];
    170:   D[2] = 0.1 * M[2];

    172:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));

    174:   /* Network subsystem initialization */
    175:   PetscCall(VecCopy(user->V0, Xnet));

    177:   /* Generator subsystem initialization */
    178:   PetscCall(VecGetArray(Xgen, &xgen));
    179:   PetscCall(VecGetArray(Xnet, &xnet));

    181:   for (i = 0; i < ngen; i++) {
    182:     Vr  = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
    183:     Vi  = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
    184:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    185:     Vm2 = Vm * Vm;
    186:     IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
    187:     IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;

    189:     delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */

    191:     theta = PETSC_PI / 2.0 - delta;

    193:     Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
    194:     Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */

    196:     Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
    197:     Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

    199:     Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
    200:     Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */

    202:     TM[i] = PG[i];

    204:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
    205:     xgen[idx]     = Eqp;
    206:     xgen[idx + 1] = Edp;
    207:     xgen[idx + 2] = delta;
    208:     xgen[idx + 3] = w_s;

    210:     idx = idx + 4;

    212:     xgen[idx]     = Id;
    213:     xgen[idx + 1] = Iq;

    215:     idx = idx + 2;

    217:     /* Exciter */
    218:     Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
    219:     SE  = k1[i] * PetscExpScalar(k2[i] * Efd);
    220:     VR  = KE[i] * Efd + SE;
    221:     RF  = KF[i] * Efd / TF[i];

    223:     xgen[idx]     = Efd;
    224:     xgen[idx + 1] = RF;
    225:     xgen[idx + 2] = VR;

    227:     Vref[i] = Vm + (VR / KA[i]);

    229:     idx = idx + 3;
    230:   }

    232:   PetscCall(VecRestoreArray(Xgen, &xgen));
    233:   PetscCall(VecRestoreArray(Xnet, &xnet));

    235:   /* PetscCall(VecView(Xgen,0)); */
    236:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
    237:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    238:   PetscFunctionReturn(PETSC_SUCCESS);
    239: }

    241: PetscErrorCode InitialGuess(Vec X, Userctx *user, const PetscScalar PGv[])
    242: {
    243:   Vec          Xgen, Xnet;
    244:   PetscScalar *xgen, *xnet;
    245:   PetscInt     i, idx = 0;
    246:   PetscScalar  Vr, Vi, IGr, IGi, Vm, Vm2;
    247:   PetscScalar  Eqp, Edp, delta;
    248:   PetscScalar  Efd, RF, VR; /* Exciter variables */
    249:   PetscScalar  Id, Iq;      /* Generator dq axis currents */
    250:   PetscScalar  theta, Vd, Vq, SE;

    252:   PetscFunctionBegin;
    253:   M[0] = 2 * H[0] / w_s;
    254:   M[1] = 2 * H[1] / w_s;
    255:   M[2] = 2 * H[2] / w_s;
    256:   D[0] = 0.1 * M[0];
    257:   D[1] = 0.1 * M[1];
    258:   D[2] = 0.1 * M[2];

    260:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));

    262:   /* Network subsystem initialization */
    263:   PetscCall(VecCopy(user->V0, Xnet));

    265:   /* Generator subsystem initialization */
    266:   PetscCall(VecGetArray(Xgen, &xgen));
    267:   PetscCall(VecGetArray(Xnet, &xnet));

    269:   for (i = 0; i < ngen; i++) {
    270:     Vr  = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
    271:     Vi  = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
    272:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    273:     Vm2 = Vm * Vm;
    274:     IGr = (Vr * PGv[i] + Vi * QG[i]) / Vm2;
    275:     IGi = (Vi * PGv[i] - Vr * QG[i]) / Vm2;

    277:     delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */

    279:     theta = PETSC_PI / 2.0 - delta;

    281:     Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
    282:     Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */

    284:     Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
    285:     Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

    287:     Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
    288:     Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */

    290:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
    291:     xgen[idx]     = Eqp;
    292:     xgen[idx + 1] = Edp;
    293:     xgen[idx + 2] = delta;
    294:     xgen[idx + 3] = w_s;

    296:     idx = idx + 4;

    298:     xgen[idx]     = Id;
    299:     xgen[idx + 1] = Iq;

    301:     idx = idx + 2;

    303:     /* Exciter */
    304:     Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
    305:     SE  = k1[i] * PetscExpScalar(k2[i] * Efd);
    306:     VR  = KE[i] * Efd + SE;
    307:     RF  = KF[i] * Efd / TF[i];

    309:     xgen[idx]     = Efd;
    310:     xgen[idx + 1] = RF;
    311:     xgen[idx + 2] = VR;

    313:     idx = idx + 3;
    314:   }

    316:   PetscCall(VecRestoreArray(Xgen, &xgen));
    317:   PetscCall(VecRestoreArray(Xnet, &xnet));

    319:   /* PetscCall(VecView(Xgen,0)); */
    320:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
    321:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: PetscErrorCode DICDPFiniteDifference(Vec X, Vec *DICDP, Userctx *user)
    326: {
    327:   Vec         Y;
    328:   PetscScalar PGv[3], eps;
    329:   PetscInt    i, j;

    331:   PetscFunctionBegin;
    332:   eps = 1.e-7;
    333:   PetscCall(VecDuplicate(X, &Y));

    335:   for (i = 0; i < ngen; i++) {
    336:     for (j = 0; j < 3; j++) PGv[j] = PG[j];
    337:     PGv[i] = PG[i] + eps;
    338:     PetscCall(InitialGuess(Y, user, PGv));
    339:     PetscCall(InitialGuess(X, user, PG));

    341:     PetscCall(VecAXPY(Y, -1.0, X));
    342:     PetscCall(VecScale(Y, 1. / eps));
    343:     PetscCall(VecCopy(Y, DICDP[i]));
    344:   }
    345:   PetscCall(VecDestroy(&Y));
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: /* Computes F = [-f(x,y);g(x,y)] */
    350: PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
    351: {
    352:   Vec          Xgen, Xnet, Fgen, Fnet;
    353:   PetscScalar *xgen, *xnet, *fgen, *fnet;
    354:   PetscInt     i, idx = 0;
    355:   PetscScalar  Vr, Vi, Vm, Vm2;
    356:   PetscScalar  Eqp, Edp, delta, w; /* Generator variables */
    357:   PetscScalar  Efd, RF, VR;        /* Exciter variables */
    358:   PetscScalar  Id, Iq;             /* Generator dq axis currents */
    359:   PetscScalar  Vd, Vq, SE;
    360:   PetscScalar  IGr, IGi, IDr, IDi;
    361:   PetscScalar  Zdq_inv[4], det;
    362:   PetscScalar  PD, QD, Vm0, *v0;
    363:   PetscInt     k;

    365:   PetscFunctionBegin;
    366:   PetscCall(VecZeroEntries(F));
    367:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    368:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
    369:   PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
    370:   PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));

    372:   /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
    373:      The generator current injection, IG, and load current injection, ID are added later
    374:   */
    375:   /* Note that the values in Ybus are stored assuming the imaginary current balance
    376:      equation is ordered first followed by real current balance equation for each bus.
    377:      Thus imaginary current contribution goes in location 2*i, and
    378:      real current contribution in 2*i+1
    379:   */
    380:   PetscCall(MatMult(user->Ybus, Xnet, Fnet));

    382:   PetscCall(VecGetArray(Xgen, &xgen));
    383:   PetscCall(VecGetArray(Xnet, &xnet));
    384:   PetscCall(VecGetArray(Fgen, &fgen));
    385:   PetscCall(VecGetArray(Fnet, &fnet));

    387:   /* Generator subsystem */
    388:   for (i = 0; i < ngen; i++) {
    389:     Eqp   = xgen[idx];
    390:     Edp   = xgen[idx + 1];
    391:     delta = xgen[idx + 2];
    392:     w     = xgen[idx + 3];
    393:     Id    = xgen[idx + 4];
    394:     Iq    = xgen[idx + 5];
    395:     Efd   = xgen[idx + 6];
    396:     RF    = xgen[idx + 7];
    397:     VR    = xgen[idx + 8];

    399:     /* Generator differential equations */
    400:     fgen[idx]     = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
    401:     fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
    402:     fgen[idx + 2] = -w + w_s;
    403:     fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];

    405:     Vr = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
    406:     Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */

    408:     PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
    409:     /* Algebraic equations for stator currents */
    410:     det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];

    412:     Zdq_inv[0] = Rs[i] / det;
    413:     Zdq_inv[1] = Xqp[i] / det;
    414:     Zdq_inv[2] = -Xdp[i] / det;
    415:     Zdq_inv[3] = Rs[i] / det;

    417:     fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
    418:     fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

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

    423:     fnet[2 * gbus[i]] -= IGi;
    424:     fnet[2 * gbus[i] + 1] -= IGr;

    426:     Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);

    428:     SE = k1[i] * PetscExpScalar(k2[i] * Efd);

    430:     /* Exciter differential equations */
    431:     fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
    432:     fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
    433:     fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];

    435:     idx = idx + 9;
    436:   }

    438:   PetscCall(VecGetArray(user->V0, &v0));
    439:   for (i = 0; i < nload; i++) {
    440:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
    441:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
    442:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    443:     Vm2 = Vm * Vm;
    444:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
    445:     PD = QD = 0.0;
    446:     for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
    447:     for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);

    449:     /* Load currents */
    450:     IDr = (PD * Vr + QD * Vi) / Vm2;
    451:     IDi = (-QD * Vr + PD * Vi) / Vm2;

    453:     fnet[2 * lbus[i]] += IDi;
    454:     fnet[2 * lbus[i] + 1] += IDr;
    455:   }
    456:   PetscCall(VecRestoreArray(user->V0, &v0));

    458:   PetscCall(VecRestoreArray(Xgen, &xgen));
    459:   PetscCall(VecRestoreArray(Xnet, &xnet));
    460:   PetscCall(VecRestoreArray(Fgen, &fgen));
    461:   PetscCall(VecRestoreArray(Fnet, &fnet));

    463:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
    464:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    465:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: /* \dot{x} - f(x,y)
    470:      g(x,y) = 0
    471:  */
    472: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
    473: {
    474:   SNES               snes;
    475:   PetscScalar       *f;
    476:   const PetscScalar *xdot;
    477:   PetscInt           i;

    479:   PetscFunctionBegin;
    480:   user->t = t;

    482:   PetscCall(TSGetSNES(ts, &snes));
    483:   PetscCall(ResidualFunction(snes, X, F, user));
    484:   PetscCall(VecGetArray(F, &f));
    485:   PetscCall(VecGetArrayRead(Xdot, &xdot));
    486:   for (i = 0; i < ngen; i++) {
    487:     f[9 * i] += xdot[9 * i];
    488:     f[9 * i + 1] += xdot[9 * i + 1];
    489:     f[9 * i + 2] += xdot[9 * i + 2];
    490:     f[9 * i + 3] += xdot[9 * i + 3];
    491:     f[9 * i + 6] += xdot[9 * i + 6];
    492:     f[9 * i + 7] += xdot[9 * i + 7];
    493:     f[9 * i + 8] += xdot[9 * i + 8];
    494:   }
    495:   PetscCall(VecRestoreArray(F, &f));
    496:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
    497:   PetscFunctionReturn(PETSC_SUCCESS);
    498: }

    500: /* This function is used for solving the algebraic system only during fault on and
    501:    off times. It computes the entire F and then zeros out the part corresponding to
    502:    differential equations
    503:  F = [0;g(y)];
    504: */
    505: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
    506: {
    507:   Userctx     *user = (Userctx *)ctx;
    508:   PetscScalar *f;
    509:   PetscInt     i;

    511:   PetscFunctionBegin;
    512:   PetscCall(ResidualFunction(snes, X, F, user));
    513:   PetscCall(VecGetArray(F, &f));
    514:   for (i = 0; i < ngen; i++) {
    515:     f[9 * i]     = 0;
    516:     f[9 * i + 1] = 0;
    517:     f[9 * i + 2] = 0;
    518:     f[9 * i + 3] = 0;
    519:     f[9 * i + 6] = 0;
    520:     f[9 * i + 7] = 0;
    521:     f[9 * i + 8] = 0;
    522:   }
    523:   PetscCall(VecRestoreArray(F, &f));
    524:   PetscFunctionReturn(PETSC_SUCCESS);
    525: }

    527: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
    528: {
    529:   PetscInt *d_nnz;
    530:   PetscInt  i, idx = 0, start = 0;
    531:   PetscInt  ncols;

    533:   PetscFunctionBegin;
    534:   PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
    535:   for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
    536:   /* Generator subsystem */
    537:   for (i = 0; i < ngen; i++) {
    538:     d_nnz[idx] += 3;
    539:     d_nnz[idx + 1] += 2;
    540:     d_nnz[idx + 2] += 2;
    541:     d_nnz[idx + 3] += 5;
    542:     d_nnz[idx + 4] += 6;
    543:     d_nnz[idx + 5] += 6;

    545:     d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
    546:     d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;

    548:     d_nnz[idx + 6] += 2;
    549:     d_nnz[idx + 7] += 2;
    550:     d_nnz[idx + 8] += 5;

    552:     idx = idx + 9;
    553:   }

    555:   start = user->neqs_gen;
    556:   for (i = 0; i < nbus; i++) {
    557:     PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
    558:     d_nnz[start + 2 * i] += ncols;
    559:     d_nnz[start + 2 * i + 1] += ncols;
    560:     PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
    561:   }

    563:   PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
    564:   PetscCall(PetscFree(d_nnz));
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: }

    568: /*
    569:    J = [-df_dx, -df_dy
    570:         dg_dx, dg_dy]
    571: */
    572: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
    573: {
    574:   Userctx           *user = (Userctx *)ctx;
    575:   Vec                Xgen, Xnet;
    576:   PetscScalar       *xgen, *xnet;
    577:   PetscInt           i, idx = 0;
    578:   PetscScalar        Vr, Vi, Vm, Vm2;
    579:   PetscScalar        Eqp, Edp, delta; /* Generator variables */
    580:   PetscScalar        Efd;             /* Exciter variables */
    581:   PetscScalar        Id, Iq;          /* Generator dq axis currents */
    582:   PetscScalar        Vd, Vq;
    583:   PetscScalar        val[10];
    584:   PetscInt           row[2], col[10];
    585:   PetscInt           net_start = user->neqs_gen;
    586:   PetscInt           ncols;
    587:   const PetscInt    *cols;
    588:   const PetscScalar *yvals;
    589:   PetscInt           k;
    590:   PetscScalar        Zdq_inv[4], det;
    591:   PetscScalar        dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
    592:   PetscScalar        dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
    593:   PetscScalar        dSE_dEfd;
    594:   PetscScalar        dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
    595:   PetscScalar        PD, QD, Vm0, *v0, Vm4;
    596:   PetscScalar        dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
    597:   PetscScalar        dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;

    599:   PetscFunctionBegin;
    600:   PetscCall(MatZeroEntries(B));
    601:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    602:   PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));

    604:   PetscCall(VecGetArray(Xgen, &xgen));
    605:   PetscCall(VecGetArray(Xnet, &xnet));

    607:   /* Generator subsystem */
    608:   for (i = 0; i < ngen; i++) {
    609:     Eqp   = xgen[idx];
    610:     Edp   = xgen[idx + 1];
    611:     delta = xgen[idx + 2];
    612:     Id    = xgen[idx + 4];
    613:     Iq    = xgen[idx + 5];
    614:     Efd   = xgen[idx + 6];

    616:     /*    fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
    617:     row[0] = idx;
    618:     col[0] = idx;
    619:     col[1] = idx + 4;
    620:     col[2] = idx + 6;
    621:     val[0] = 1 / Td0p[i];
    622:     val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
    623:     val[2] = -1 / Td0p[i];

    625:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

    627:     /*    fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
    628:     row[0] = idx + 1;
    629:     col[0] = idx + 1;
    630:     col[1] = idx + 5;
    631:     val[0] = 1 / Tq0p[i];
    632:     val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
    633:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

    635:     /*    fgen[idx+2] = - w + w_s; */
    636:     row[0] = idx + 2;
    637:     col[0] = idx + 2;
    638:     col[1] = idx + 3;
    639:     val[0] = 0;
    640:     val[1] = -1;
    641:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

    643:     /*    fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
    644:     row[0] = idx + 3;
    645:     col[0] = idx;
    646:     col[1] = idx + 1;
    647:     col[2] = idx + 3;
    648:     col[3] = idx + 4;
    649:     col[4] = idx + 5;
    650:     val[0] = Iq / M[i];
    651:     val[1] = Id / M[i];
    652:     val[2] = D[i] / M[i];
    653:     val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
    654:     val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
    655:     PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));

    657:     Vr = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
    658:     Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
    659:     PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

    661:     det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];

    663:     Zdq_inv[0] = Rs[i] / det;
    664:     Zdq_inv[1] = Xqp[i] / det;
    665:     Zdq_inv[2] = -Xdp[i] / det;
    666:     Zdq_inv[3] = Rs[i] / det;

    668:     dVd_dVr    = PetscSinScalar(delta);
    669:     dVd_dVi    = -PetscCosScalar(delta);
    670:     dVq_dVr    = PetscCosScalar(delta);
    671:     dVq_dVi    = PetscSinScalar(delta);
    672:     dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
    673:     dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);

    675:     /*    fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
    676:     row[0] = idx + 4;
    677:     col[0] = idx;
    678:     col[1] = idx + 1;
    679:     col[2] = idx + 2;
    680:     val[0] = -Zdq_inv[1];
    681:     val[1] = -Zdq_inv[0];
    682:     val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
    683:     col[3] = idx + 4;
    684:     col[4] = net_start + 2 * gbus[i];
    685:     col[5] = net_start + 2 * gbus[i] + 1;
    686:     val[3] = 1;
    687:     val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
    688:     val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
    689:     PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));

    691:     /*  fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
    692:     row[0] = idx + 5;
    693:     col[0] = idx;
    694:     col[1] = idx + 1;
    695:     col[2] = idx + 2;
    696:     val[0] = -Zdq_inv[3];
    697:     val[1] = -Zdq_inv[2];
    698:     val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
    699:     col[3] = idx + 5;
    700:     col[4] = net_start + 2 * gbus[i];
    701:     col[5] = net_start + 2 * gbus[i] + 1;
    702:     val[3] = 1;
    703:     val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
    704:     val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
    705:     PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));

    707:     dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
    708:     dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
    709:     dIGr_dId    = PetscSinScalar(delta);
    710:     dIGr_dIq    = PetscCosScalar(delta);
    711:     dIGi_dId    = -PetscCosScalar(delta);
    712:     dIGi_dIq    = PetscSinScalar(delta);

    714:     /* fnet[2*gbus[i]]   -= IGi; */
    715:     row[0] = net_start + 2 * gbus[i];
    716:     col[0] = idx + 2;
    717:     col[1] = idx + 4;
    718:     col[2] = idx + 5;
    719:     val[0] = -dIGi_ddelta;
    720:     val[1] = -dIGi_dId;
    721:     val[2] = -dIGi_dIq;
    722:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

    724:     /* fnet[2*gbus[i]+1]   -= IGr; */
    725:     row[0] = net_start + 2 * gbus[i] + 1;
    726:     col[0] = idx + 2;
    727:     col[1] = idx + 4;
    728:     col[2] = idx + 5;
    729:     val[0] = -dIGr_ddelta;
    730:     val[1] = -dIGr_dId;
    731:     val[2] = -dIGr_dIq;
    732:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

    734:     Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);

    736:     /*    fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
    737:     /*    SE  = k1[i]*PetscExpScalar(k2[i]*Efd); */
    738:     dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);

    740:     row[0] = idx + 6;
    741:     col[0] = idx + 6;
    742:     col[1] = idx + 8;
    743:     val[0] = (KE[i] + dSE_dEfd) / TE[i];
    744:     val[1] = -1 / TE[i];
    745:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

    747:     /* Exciter differential equations */

    749:     /*    fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
    750:     row[0] = idx + 7;
    751:     col[0] = idx + 6;
    752:     col[1] = idx + 7;
    753:     val[0] = (-KF[i] / TF[i]) / TF[i];
    754:     val[1] = 1 / TF[i];
    755:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

    757:     /*    fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
    758:     /* Vm = (Vd^2 + Vq^2)^0.5; */
    759:     dVm_dVd = Vd / Vm;
    760:     dVm_dVq = Vq / Vm;
    761:     dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
    762:     dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
    763:     row[0]  = idx + 8;
    764:     col[0]  = idx + 6;
    765:     col[1]  = idx + 7;
    766:     col[2]  = idx + 8;
    767:     val[0]  = (KA[i] * KF[i] / TF[i]) / TA[i];
    768:     val[1]  = -KA[i] / TA[i];
    769:     val[2]  = 1 / TA[i];
    770:     col[3]  = net_start + 2 * gbus[i];
    771:     col[4]  = net_start + 2 * gbus[i] + 1;
    772:     val[3]  = KA[i] * dVm_dVr / TA[i];
    773:     val[4]  = KA[i] * dVm_dVi / TA[i];
    774:     PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
    775:     idx = idx + 9;
    776:   }

    778:   for (i = 0; i < nbus; i++) {
    779:     PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
    780:     row[0] = net_start + 2 * i;
    781:     for (k = 0; k < ncols; k++) {
    782:       col[k] = net_start + cols[k];
    783:       val[k] = yvals[k];
    784:     }
    785:     PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
    786:     PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));

    788:     PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
    789:     row[0] = net_start + 2 * i + 1;
    790:     for (k = 0; k < ncols; k++) {
    791:       col[k] = net_start + cols[k];
    792:       val[k] = yvals[k];
    793:     }
    794:     PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
    795:     PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
    796:   }

    798:   PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
    799:   PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));

    801:   PetscCall(VecGetArray(user->V0, &v0));
    802:   for (i = 0; i < nload; i++) {
    803:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
    804:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
    805:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
    806:     Vm2 = Vm * Vm;
    807:     Vm4 = Vm2 * Vm2;
    808:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
    809:     PD = QD = 0.0;
    810:     dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
    811:     for (k = 0; k < ld_nsegsp[i]; k++) {
    812:       PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
    813:       dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
    814:       dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
    815:     }
    816:     for (k = 0; k < ld_nsegsq[i]; k++) {
    817:       QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
    818:       dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
    819:       dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
    820:     }

    822:     /*    IDr = (PD*Vr + QD*Vi)/Vm2; */
    823:     /*    IDi = (-QD*Vr + PD*Vi)/Vm2; */

    825:     dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
    826:     dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;

    828:     dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
    829:     dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;

    831:     /*    fnet[2*lbus[i]]   += IDi; */
    832:     row[0] = net_start + 2 * lbus[i];
    833:     col[0] = net_start + 2 * lbus[i];
    834:     col[1] = net_start + 2 * lbus[i] + 1;
    835:     val[0] = dIDi_dVr;
    836:     val[1] = dIDi_dVi;
    837:     PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
    838:     /*    fnet[2*lbus[i]+1] += IDr; */
    839:     row[0] = net_start + 2 * lbus[i] + 1;
    840:     col[0] = net_start + 2 * lbus[i];
    841:     col[1] = net_start + 2 * lbus[i] + 1;
    842:     val[0] = dIDr_dVr;
    843:     val[1] = dIDr_dVi;
    844:     PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
    845:   }
    846:   PetscCall(VecRestoreArray(user->V0, &v0));

    848:   PetscCall(VecRestoreArray(Xgen, &xgen));
    849:   PetscCall(VecRestoreArray(Xnet, &xnet));

    851:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));

    853:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    854:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    855:   PetscFunctionReturn(PETSC_SUCCESS);
    856: }

    858: /*
    859:    J = [I, 0
    860:         dg_dx, dg_dy]
    861: */
    862: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
    863: {
    864:   Userctx *user = (Userctx *)ctx;

    866:   PetscFunctionBegin;
    867:   PetscCall(ResidualJacobian(snes, X, A, B, ctx));
    868:   PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
    869:   PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
    870:   PetscFunctionReturn(PETSC_SUCCESS);
    871: }

    873: /*
    874:    J = [a*I-df_dx, -df_dy
    875:         dg_dx, dg_dy]
    876: */

    878: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
    879: {
    880:   SNES        snes;
    881:   PetscScalar atmp = (PetscScalar)a;
    882:   PetscInt    i, row;

    884:   PetscFunctionBegin;
    885:   user->t = t;

    887:   PetscCall(TSGetSNES(ts, &snes));
    888:   PetscCall(ResidualJacobian(snes, X, A, B, user));
    889:   for (i = 0; i < ngen; i++) {
    890:     row = 9 * i;
    891:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    892:     row = 9 * i + 1;
    893:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    894:     row = 9 * i + 2;
    895:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    896:     row = 9 * i + 3;
    897:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    898:     row = 9 * i + 6;
    899:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    900:     row = 9 * i + 7;
    901:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    902:     row = 9 * i + 8;
    903:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
    904:   }
    905:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    906:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    907:   PetscFunctionReturn(PETSC_SUCCESS);
    908: }

    910: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
    911: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
    912: {
    913:   PetscScalar a;
    914:   PetscInt    row, col;
    915:   Userctx    *ctx = (Userctx *)ctx0;

    917:   PetscFunctionBeginUser;
    918:   if (ctx->jacp_flg) {
    919:     PetscCall(MatZeroEntries(A));

    921:     for (col = 0; col < 3; col++) {
    922:       a   = 1.0 / M[col];
    923:       row = 9 * col + 3;
    924:       PetscCall(MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES));
    925:     }

    927:     ctx->jacp_flg = PETSC_FALSE;

    929:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    930:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    931:   }
    932:   PetscFunctionReturn(PETSC_SUCCESS);
    933: }

    935: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
    936: {
    937:   const PetscScalar *u;
    938:   PetscInt           idx;
    939:   Vec                Xgen, Xnet;
    940:   PetscScalar       *r, *xgen;
    941:   PetscInt           i;

    943:   PetscFunctionBegin;
    944:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    945:   PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));

    947:   PetscCall(VecGetArray(Xgen, &xgen));

    949:   PetscCall(VecGetArrayRead(U, &u));
    950:   PetscCall(VecGetArray(R, &r));
    951:   r[0] = 0.;
    952:   idx  = 0;
    953:   for (i = 0; i < ngen; i++) {
    954:     r[0] += PetscPowScalarInt(PetscMax(0., PetscMax(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->freq_l - xgen[idx + 3] / (2. * PETSC_PI))), user->pow);
    955:     idx += 9;
    956:   }
    957:   PetscCall(VecRestoreArrayRead(U, &u));
    958:   PetscCall(VecRestoreArray(R, &r));
    959:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    960:   PetscFunctionReturn(PETSC_SUCCESS);
    961: }

    963: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
    964: {
    965:   Vec          Xgen, Xnet, Dgen, Dnet;
    966:   PetscScalar *xgen, *dgen;
    967:   PetscInt     i;
    968:   PetscInt     idx;
    969:   Vec          drdu_col;
    970:   PetscScalar *xarr;

    972:   PetscFunctionBegin;
    973:   PetscCall(VecDuplicate(U, &drdu_col));
    974:   PetscCall(MatDenseGetColumn(DRDU, 0, &xarr));
    975:   PetscCall(VecPlaceArray(drdu_col, xarr));
    976:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    977:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet));
    978:   PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
    979:   PetscCall(DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet));

    981:   PetscCall(VecGetArray(Xgen, &xgen));
    982:   PetscCall(VecGetArray(Dgen, &dgen));

    984:   idx = 0;
    985:   for (i = 0; i < ngen; i++) {
    986:     dgen[idx + 3] = 0.;
    987:     if (xgen[idx + 3] / (2. * PETSC_PI) > user->freq_u) dgen[idx + 3] = user->pow * PetscPowScalarInt(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->pow - 1) / (2. * PETSC_PI);
    988:     if (xgen[idx + 3] / (2. * PETSC_PI) < user->freq_l) dgen[idx + 3] = user->pow * PetscPowScalarInt(user->freq_l - xgen[idx + 3] / (2. * PETSC_PI), user->pow - 1) / (-2. * PETSC_PI);
    989:     idx += 9;
    990:   }

    992:   PetscCall(VecRestoreArray(Dgen, &dgen));
    993:   PetscCall(VecRestoreArray(Xgen, &xgen));
    994:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet));
    995:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet));
    996:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
    997:   PetscCall(VecResetArray(drdu_col));
    998:   PetscCall(MatDenseRestoreColumn(DRDU, &xarr));
    999:   PetscCall(VecDestroy(&drdu_col));
    1000:   PetscFunctionReturn(PETSC_SUCCESS);
    1001: }

    1003: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
    1004: {
    1005:   PetscFunctionBegin;
    1006:   PetscFunctionReturn(PETSC_SUCCESS);
    1007: }

    1009: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
    1010: {
    1011:   PetscScalar *x, *y, sensip;
    1012:   PetscInt     i;

    1014:   PetscFunctionBegin;
    1015:   PetscCall(VecGetArray(lambda, &x));
    1016:   PetscCall(VecGetArray(mu, &y));

    1018:   for (i = 0; i < 3; i++) {
    1019:     PetscCall(VecDot(lambda, DICDP[i], &sensip));
    1020:     sensip = sensip + y[i];
    1021:     /* PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip)); */
    1022:     y[i] = sensip;
    1023:   }
    1024:   PetscCall(VecRestoreArray(mu, &y));
    1025:   PetscFunctionReturn(PETSC_SUCCESS);
    1026: }

    1028: int main(int argc, char **argv)
    1029: {
    1030:   Userctx      user;
    1031:   Vec          p;
    1032:   PetscScalar *x_ptr;
    1033:   PetscMPIInt  size;
    1034:   PetscInt     i;
    1035:   PetscViewer  Xview, Ybusview;
    1036:   PetscInt    *idx2;
    1037:   Tao          tao;
    1038:   KSP          ksp;
    1039:   PC           pc;
    1040:   Vec          lowerb, upperb;

    1042:   PetscFunctionBeginUser;
    1043:   PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
    1044:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    1045:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    1047:   user.jacp_flg   = PETSC_TRUE;
    1048:   user.neqs_gen   = 9 * ngen; /* # eqs. for generator subsystem */
    1049:   user.neqs_net   = 2 * nbus; /* # eqs. for network subsystem   */
    1050:   user.neqs_pgrid = user.neqs_gen + user.neqs_net;

    1052:   /* Create indices for differential and algebraic equations */
    1053:   PetscCall(PetscMalloc1(7 * ngen, &idx2));
    1054:   for (i = 0; i < ngen; i++) {
    1055:     idx2[7 * i]     = 9 * i;
    1056:     idx2[7 * i + 1] = 9 * i + 1;
    1057:     idx2[7 * i + 2] = 9 * i + 2;
    1058:     idx2[7 * i + 3] = 9 * i + 3;
    1059:     idx2[7 * i + 4] = 9 * i + 6;
    1060:     idx2[7 * i + 5] = 9 * i + 7;
    1061:     idx2[7 * i + 6] = 9 * i + 8;
    1062:   }
    1063:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
    1064:   PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
    1065:   PetscCall(PetscFree(idx2));

    1067:   /* Set run time options */
    1068:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
    1069:   {
    1070:     user.tfaulton  = 1.0;
    1071:     user.tfaultoff = 1.2;
    1072:     user.Rfault    = 0.0001;
    1073:     user.faultbus  = 8;
    1074:     PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
    1075:     PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
    1076:     PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
    1077:     user.t0   = 0.0;
    1078:     user.tmax = 1.3;
    1079:     PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
    1080:     PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
    1081:     user.freq_u = 61.0;
    1082:     user.freq_l = 59.0;
    1083:     user.pow    = 2;
    1084:     PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
    1085:     PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
    1086:     PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
    1087:   }
    1088:   PetscOptionsEnd();

    1090:   /* Create DMs for generator and network subsystems */
    1091:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
    1092:   PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
    1093:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
    1094:   PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
    1095:   PetscCall(DMSetFromOptions(user.dmnet));
    1096:   PetscCall(DMSetUp(user.dmnet));
    1097:   /* Create a composite DM packer and add the two DMs */
    1098:   PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
    1099:   PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
    1100:   PetscCall(DMSetFromOptions(user.dmgen));
    1101:   PetscCall(DMSetUp(user.dmgen));
    1102:   PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
    1103:   PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));

    1105:   /* Read initial voltage vector and Ybus */
    1106:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
    1107:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));

    1109:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
    1110:   PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
    1111:   PetscCall(VecLoad(user.V0, Xview));

    1113:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
    1114:   PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
    1115:   PetscCall(MatSetType(user.Ybus, MATBAIJ));
    1116:   /*  PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
    1117:   PetscCall(MatLoad(user.Ybus, Ybusview));

    1119:   PetscCall(PetscViewerDestroy(&Xview));
    1120:   PetscCall(PetscViewerDestroy(&Ybusview));

    1122:   /* Allocate space for Jacobians */
    1123:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.J));
    1124:   PetscCall(MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
    1125:   PetscCall(MatSetFromOptions(user.J));
    1126:   PetscCall(PreallocateJacobian(user.J, &user));

    1128:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
    1129:   PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3));
    1130:   PetscCall(MatSetFromOptions(user.Jacp));
    1131:   PetscCall(MatSetUp(user.Jacp));
    1132:   PetscCall(MatZeroEntries(user.Jacp)); /* initialize to zeros */

    1134:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP));
    1135:   PetscCall(MatSetUp(user.DRDP));
    1136:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU));
    1137:   PetscCall(MatSetUp(user.DRDU));

    1139:   /* Create TAO solver and set desired solution method */
    1140:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
    1141:   PetscCall(TaoSetType(tao, TAOBLMVM));
    1142:   /*
    1143:      Optimization starts
    1144:   */
    1145:   /* Set initial solution guess */
    1146:   PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
    1147:   PetscCall(VecGetArray(p, &x_ptr));
    1148:   x_ptr[0] = PG[0];
    1149:   x_ptr[1] = PG[1];
    1150:   x_ptr[2] = PG[2];
    1151:   PetscCall(VecRestoreArray(p, &x_ptr));

    1153:   PetscCall(TaoSetSolution(tao, p));
    1154:   /* Set routine for function and gradient evaluation */
    1155:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));

    1157:   /* Set bounds for the optimization */
    1158:   PetscCall(VecDuplicate(p, &lowerb));
    1159:   PetscCall(VecDuplicate(p, &upperb));
    1160:   PetscCall(VecGetArray(lowerb, &x_ptr));
    1161:   x_ptr[0] = 0.5;
    1162:   x_ptr[1] = 0.5;
    1163:   x_ptr[2] = 0.5;
    1164:   PetscCall(VecRestoreArray(lowerb, &x_ptr));
    1165:   PetscCall(VecGetArray(upperb, &x_ptr));
    1166:   x_ptr[0] = 2.0;
    1167:   x_ptr[1] = 2.0;
    1168:   x_ptr[2] = 2.0;
    1169:   PetscCall(VecRestoreArray(upperb, &x_ptr));
    1170:   PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));

    1172:   /* Check for any TAO command line options */
    1173:   PetscCall(TaoSetFromOptions(tao));
    1174:   PetscCall(TaoGetKSP(tao, &ksp));
    1175:   if (ksp) {
    1176:     PetscCall(KSPGetPC(ksp, &pc));
    1177:     PetscCall(PCSetType(pc, PCNONE));
    1178:   }

    1180:   /* SOLVE THE APPLICATION */
    1181:   PetscCall(TaoSolve(tao));

    1183:   PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
    1184:   /* Free TAO data structures */
    1185:   PetscCall(TaoDestroy(&tao));

    1187:   PetscCall(DMDestroy(&user.dmgen));
    1188:   PetscCall(DMDestroy(&user.dmnet));
    1189:   PetscCall(DMDestroy(&user.dmpgrid));
    1190:   PetscCall(ISDestroy(&user.is_diff));
    1191:   PetscCall(ISDestroy(&user.is_alg));

    1193:   PetscCall(MatDestroy(&user.J));
    1194:   PetscCall(MatDestroy(&user.Jacp));
    1195:   PetscCall(MatDestroy(&user.Ybus));
    1196:   /* PetscCall(MatDestroy(&user.Sol)); */
    1197:   PetscCall(VecDestroy(&user.V0));
    1198:   PetscCall(VecDestroy(&p));
    1199:   PetscCall(VecDestroy(&lowerb));
    1200:   PetscCall(VecDestroy(&upperb));
    1201:   PetscCall(MatDestroy(&user.DRDU));
    1202:   PetscCall(MatDestroy(&user.DRDP));
    1203:   PetscCall(PetscFinalize());
    1204:   return 0;
    1205: }

    1207: /* ------------------------------------------------------------------ */
    1208: /*
    1209:    FormFunction - Evaluates the function and corresponding gradient.

    1211:    Input Parameters:
    1212:    tao - the Tao context
    1213:    X   - the input vector
    1214:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

    1216:    Output Parameters:
    1217:    f   - the newly evaluated function
    1218:    G   - the newly evaluated gradient
    1219: */
    1220: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx0)
    1221: {
    1222:   TS       ts, quadts;
    1223:   SNES     snes_alg;
    1224:   Userctx *ctx = (Userctx *)ctx0;
    1225:   Vec      X;
    1226:   PetscInt i;
    1227:   /* sensitivity context */
    1228:   PetscScalar *x_ptr;
    1229:   Vec          lambda[1], q;
    1230:   Vec          mu[1];
    1231:   PetscInt     steps1, steps2, steps3;
    1232:   Vec          DICDP[3];
    1233:   Vec          F_alg;
    1234:   PetscInt     row_loc, col_loc;
    1235:   PetscScalar  val;
    1236:   Vec          Xdot;

    1238:   PetscFunctionBegin;
    1239:   PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
    1240:   PG[0] = x_ptr[0];
    1241:   PG[1] = x_ptr[1];
    1242:   PG[2] = x_ptr[2];
    1243:   PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));

    1245:   ctx->stepnum = 0;

    1247:   PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));

    1249:   /* Create matrix to save solutions at each time step */
    1250:   /* PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */
    1251:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1252:      Create timestepping solver context
    1253:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    1254:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    1255:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
    1256:   PetscCall(TSSetType(ts, TSCN));
    1257:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
    1258:   PetscCall(TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobianFn *)IJacobian, ctx));
    1259:   PetscCall(TSSetApplicationContext(ts, ctx));
    1260:   /*   Set RHS JacobianP */
    1261:   PetscCall(TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx));

    1263:   PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
    1264:   PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, ctx));
    1265:   PetscCall(TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, ctx));
    1266:   PetscCall(TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, ctx));

    1268:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1269:      Set initial conditions
    1270:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    1271:   PetscCall(SetInitialGuess(X, ctx));

    1273:   /* Approximate DICDP with finite difference, we want to zero out network variables */
    1274:   for (i = 0; i < 3; i++) PetscCall(VecDuplicate(X, &DICDP[i]));
    1275:   PetscCall(DICDPFiniteDifference(X, DICDP, ctx));

    1277:   PetscCall(VecDuplicate(X, &F_alg));
    1278:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
    1279:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
    1280:   PetscCall(MatZeroEntries(ctx->J));
    1281:   PetscCall(SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx));
    1282:   PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
    1283:   PetscCall(SNESSetFromOptions(snes_alg));
    1284:   ctx->alg_flg = PETSC_TRUE;
    1285:   /* Solve the algebraic equations */
    1286:   PetscCall(SNESSolve(snes_alg, NULL, X));

    1288:   /* Just to set up the Jacobian structure */
    1289:   PetscCall(VecDuplicate(X, &Xdot));
    1290:   PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx));
    1291:   PetscCall(VecDestroy(&Xdot));

    1293:   ctx->stepnum++;

    1295:   /*
    1296:     Save trajectory of solution so that TSAdjointSolve() may be used
    1297:   */
    1298:   PetscCall(TSSetSaveTrajectory(ts));

    1300:   PetscCall(TSSetTimeStep(ts, 0.01));
    1301:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    1302:   PetscCall(TSSetFromOptions(ts));
    1303:   /* PetscCall(TSSetPostStep(ts,SaveSolution)); */

    1305:   /* Prefault period */
    1306:   ctx->alg_flg = PETSC_FALSE;
    1307:   PetscCall(TSSetTime(ts, 0.0));
    1308:   PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
    1309:   PetscCall(TSSolve(ts, X));
    1310:   PetscCall(TSGetStepNumber(ts, &steps1));

    1312:   /* Create the nonlinear solver for solving the algebraic system */
    1313:   /* Note that although the algebraic system needs to be solved only for
    1314:      Idq and V, we reuse the entire system including xgen. The xgen
    1315:      variables are held constant by setting their residuals to 0 and
    1316:      putting a 1 on the Jacobian diagonal for xgen rows
    1317:   */
    1318:   PetscCall(MatZeroEntries(ctx->J));

    1320:   /* Apply disturbance - resistive fault at ctx->faultbus */
    1321:   /* This is done by adding shunt conductance to the diagonal location
    1322:      in the Ybus matrix */
    1323:   row_loc = 2 * ctx->faultbus;
    1324:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
    1325:   val     = 1 / ctx->Rfault;
    1326:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
    1327:   row_loc = 2 * ctx->faultbus + 1;
    1328:   col_loc = 2 * ctx->faultbus; /* Location for G */
    1329:   val     = 1 / ctx->Rfault;
    1330:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

    1332:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
    1333:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

    1335:   ctx->alg_flg = PETSC_TRUE;
    1336:   /* Solve the algebraic equations */
    1337:   PetscCall(SNESSolve(snes_alg, NULL, X));

    1339:   ctx->stepnum++;

    1341:   /* Disturbance period */
    1342:   ctx->alg_flg = PETSC_FALSE;
    1343:   PetscCall(TSSetTime(ts, ctx->tfaulton));
    1344:   PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
    1345:   PetscCall(TSSolve(ts, X));
    1346:   PetscCall(TSGetStepNumber(ts, &steps2));
    1347:   steps2 -= steps1;

    1349:   /* Remove the fault */
    1350:   row_loc = 2 * ctx->faultbus;
    1351:   col_loc = 2 * ctx->faultbus + 1;
    1352:   val     = -1 / ctx->Rfault;
    1353:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
    1354:   row_loc = 2 * ctx->faultbus + 1;
    1355:   col_loc = 2 * ctx->faultbus;
    1356:   val     = -1 / ctx->Rfault;
    1357:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

    1359:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
    1360:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

    1362:   PetscCall(MatZeroEntries(ctx->J));

    1364:   ctx->alg_flg = PETSC_TRUE;

    1366:   /* Solve the algebraic equations */
    1367:   PetscCall(SNESSolve(snes_alg, NULL, X));

    1369:   ctx->stepnum++;

    1371:   /* Post-disturbance period */
    1372:   ctx->alg_flg = PETSC_TRUE;
    1373:   PetscCall(TSSetTime(ts, ctx->tfaultoff));
    1374:   PetscCall(TSSetMaxTime(ts, ctx->tmax));
    1375:   PetscCall(TSSolve(ts, X));
    1376:   PetscCall(TSGetStepNumber(ts, &steps3));
    1377:   steps3 -= steps2;
    1378:   steps3 -= steps1;

    1380:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1381:      Adjoint model starts here
    1382:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    1383:   PetscCall(TSSetPostStep(ts, NULL));
    1384:   PetscCall(MatCreateVecs(ctx->J, &lambda[0], NULL));
    1385:   /*   Set initial conditions for the adjoint integration */
    1386:   PetscCall(VecZeroEntries(lambda[0]));

    1388:   PetscCall(MatCreateVecs(ctx->Jacp, &mu[0], NULL));
    1389:   PetscCall(VecZeroEntries(mu[0]));
    1390:   PetscCall(TSSetCostGradients(ts, 1, lambda, mu));

    1392:   PetscCall(TSAdjointSetSteps(ts, steps3));
    1393:   PetscCall(TSAdjointSolve(ts));

    1395:   PetscCall(MatZeroEntries(ctx->J));
    1396:   /* Applying disturbance - resistive fault at ctx->faultbus */
    1397:   /* This is done by deducting shunt conductance to the diagonal location
    1398:      in the Ybus matrix */
    1399:   row_loc = 2 * ctx->faultbus;
    1400:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
    1401:   val     = 1. / ctx->Rfault;
    1402:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
    1403:   row_loc = 2 * ctx->faultbus + 1;
    1404:   col_loc = 2 * ctx->faultbus; /* Location for G */
    1405:   val     = 1. / ctx->Rfault;
    1406:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

    1408:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
    1409:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

    1411:   /*   Set number of steps for the adjoint integration */
    1412:   PetscCall(TSAdjointSetSteps(ts, steps2));
    1413:   PetscCall(TSAdjointSolve(ts));

    1415:   PetscCall(MatZeroEntries(ctx->J));
    1416:   /* remove the fault */
    1417:   row_loc = 2 * ctx->faultbus;
    1418:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
    1419:   val     = -1. / ctx->Rfault;
    1420:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
    1421:   row_loc = 2 * ctx->faultbus + 1;
    1422:   col_loc = 2 * ctx->faultbus; /* Location for G */
    1423:   val     = -1. / ctx->Rfault;
    1424:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

    1426:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
    1427:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

    1429:   /*   Set number of steps for the adjoint integration */
    1430:   PetscCall(TSAdjointSetSteps(ts, steps1));
    1431:   PetscCall(TSAdjointSolve(ts));

    1433:   PetscCall(ComputeSensiP(lambda[0], mu[0], DICDP, ctx));
    1434:   PetscCall(VecCopy(mu[0], G));

    1436:   PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
    1437:   PetscCall(TSGetSolution(quadts, &q));
    1438:   PetscCall(VecGetArray(q, &x_ptr));
    1439:   *f       = x_ptr[0];
    1440:   x_ptr[0] = 0;
    1441:   PetscCall(VecRestoreArray(q, &x_ptr));

    1443:   PetscCall(VecDestroy(&lambda[0]));
    1444:   PetscCall(VecDestroy(&mu[0]));

    1446:   PetscCall(SNESDestroy(&snes_alg));
    1447:   PetscCall(VecDestroy(&F_alg));
    1448:   PetscCall(VecDestroy(&X));
    1449:   PetscCall(TSDestroy(&ts));
    1450:   for (i = 0; i < 3; i++) PetscCall(VecDestroy(&DICDP[i]));
    1451:   PetscFunctionReturn(PETSC_SUCCESS);
    1452: }

    1454: /*TEST

    1456:    build:
    1457:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

    1459:    test:
    1460:       args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
    1461:       localrunfiles: petscoptions X.bin Ybus.bin

    1463: TEST*/