Actual source code: ex3.h

      1: #pragma once

      3: typedef enum {
      4:   SA_ADJ,
      5:   SA_TLM
      6: } SAMethod;
      7: static const char *const SAMethods[] = {"ADJ", "TLM", "SAMethod", "SA_", 0};

      9: typedef struct {
     10:   PetscScalar H, D, omega_b, omega_s, Pmax, Pmax_ini, Pm, E, V, X, u_s, c;
     11:   PetscInt    beta;
     12:   PetscReal   tf, tcl;
     13:   /* Solver context */
     14:   TS       ts, quadts;
     15:   Vec      U;    /* solution will be stored here */
     16:   Mat      Jac;  /* Jacobian matrix */
     17:   Mat      Jacp; /* Jacobianp matrix */
     18:   Mat      DRDU, DRDP;
     19:   SAMethod sa;
     20: } AppCtx;

     22: /* Event check */
     23: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec X, PetscReal *fvalue, void *ctx)
     24: {
     25:   AppCtx *user = (AppCtx *)ctx;

     27:   PetscFunctionBegin;
     28:   /* Event for fault-on time */
     29:   fvalue[0] = t - user->tf;
     30:   /* Event for fault-off time */
     31:   fvalue[1] = t - user->tcl;
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *ctx)
     36: {
     37:   AppCtx *user = (AppCtx *)ctx;

     39:   PetscFunctionBegin;
     40:   if (event_list[0] == 0) {
     41:     if (forwardsolve) user->Pmax = 0.0; /* Apply disturbance - this is done by setting Pmax = 0 */
     42:     else user->Pmax = user->Pmax_ini;   /* Going backward, reversal of event */
     43:   } else if (event_list[0] == 1) {
     44:     if (forwardsolve) user->Pmax = user->Pmax_ini; /* Remove the fault  - this is done by setting Pmax = Pmax_ini */
     45:     else user->Pmax = 0.0;                         /* Going backward, reversal of event */
     46:   }
     47:   PetscCall(TSRestartStep(ts)); /* Must set restart flag to true, otherwise methods with FSAL will fail */
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: /*
     52:      Defines the ODE passed to the ODE solver
     53: */
     54: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
     55: {
     56:   PetscScalar       *f, Pmax;
     57:   const PetscScalar *u;

     59:   PetscFunctionBegin;
     60:   /*  The next three lines allow us to access the entries of the vectors directly */
     61:   PetscCall(VecGetArrayRead(U, &u));
     62:   PetscCall(VecGetArray(F, &f));
     63:   Pmax = ctx->Pmax;
     64:   f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
     65:   f[1] = ctx->omega_s / (2.0 * ctx->H) * (ctx->Pm - Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s));

     67:   PetscCall(VecRestoreArrayRead(U, &u));
     68:   PetscCall(VecRestoreArray(F, &f));
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: /*
     73:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of a and the Jacobian.
     74: */
     75: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
     76: {
     77:   PetscInt           rowcol[] = {0, 1};
     78:   PetscScalar        J[2][2], Pmax;
     79:   const PetscScalar *u;

     81:   PetscFunctionBegin;
     82:   PetscCall(VecGetArrayRead(U, &u));
     83:   Pmax    = ctx->Pmax;
     84:   J[0][0] = 0;
     85:   J[0][1] = ctx->omega_b;
     86:   J[1][0] = -ctx->omega_s / (2.0 * ctx->H) * Pmax * PetscCosScalar(u[0]);
     87:   J[1][1] = -ctx->omega_s / (2.0 * ctx->H) * ctx->D;
     88:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
     89:   PetscCall(VecRestoreArrayRead(U, &u));

     91:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     92:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     93:   if (A != B) {
     94:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     95:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     96:   }
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: /*
    101:      Defines the ODE passed to the ODE solver
    102: */
    103: PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
    104: {
    105:   PetscScalar       *f, Pmax;
    106:   const PetscScalar *u, *udot;

    108:   PetscFunctionBegin;
    109:   /*  The next three lines allow us to access the entries of the vectors directly */
    110:   PetscCall(VecGetArrayRead(U, &u));
    111:   PetscCall(VecGetArrayRead(Udot, &udot));
    112:   PetscCall(VecGetArray(F, &f));
    113:   Pmax = ctx->Pmax;
    114:   f[0] = udot[0] - ctx->omega_b * (u[1] - ctx->omega_s);
    115:   f[1] = 2.0 * ctx->H / ctx->omega_s * udot[1] + Pmax * PetscSinScalar(u[0]) + ctx->D * (u[1] - ctx->omega_s) - ctx->Pm;

    117:   PetscCall(VecRestoreArrayRead(U, &u));
    118:   PetscCall(VecRestoreArrayRead(Udot, &udot));
    119:   PetscCall(VecRestoreArray(F, &f));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /*
    124:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
    125: */
    126: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
    127: {
    128:   PetscInt           rowcol[] = {0, 1};
    129:   PetscScalar        J[2][2], Pmax;
    130:   const PetscScalar *u, *udot;

    132:   PetscFunctionBegin;
    133:   PetscCall(VecGetArrayRead(U, &u));
    134:   PetscCall(VecGetArrayRead(Udot, &udot));
    135:   Pmax    = ctx->Pmax;
    136:   J[0][0] = a;
    137:   J[0][1] = -ctx->omega_b;
    138:   J[1][1] = 2.0 * ctx->H / ctx->omega_s * a + ctx->D;
    139:   J[1][0] = Pmax * PetscCosScalar(u[0]);

    141:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
    142:   PetscCall(VecRestoreArrayRead(U, &u));
    143:   PetscCall(VecRestoreArrayRead(Udot, &udot));

    145:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    146:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    147:   if (A != B) {
    148:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    149:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    150:   }
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
    155: {
    156:   PetscInt     row[] = {0, 1}, col[] = {0};
    157:   PetscScalar *x, J[2][1];
    158:   AppCtx      *ctx = (AppCtx *)ctx0;

    160:   PetscFunctionBeginUser;
    161:   PetscCall(VecGetArray(X, &x));
    162:   J[0][0] = 0;
    163:   J[1][0] = ctx->omega_s / (2.0 * ctx->H);
    164:   PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));

    166:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    167:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
    172: {
    173:   PetscScalar       *r;
    174:   const PetscScalar *u;

    176:   PetscFunctionBegin;
    177:   PetscCall(VecGetArrayRead(U, &u));
    178:   PetscCall(VecGetArray(R, &r));
    179:   r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
    180:   PetscCall(VecRestoreArray(R, &r));
    181:   PetscCall(VecRestoreArrayRead(U, &u));
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: /* Transpose of DRDU */
    186: PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
    187: {
    188:   PetscScalar        ru[2];
    189:   PetscInt           row[] = {0, 1}, col[] = {0};
    190:   const PetscScalar *u;

    192:   PetscFunctionBegin;
    193:   PetscCall(VecGetArrayRead(U, &u));
    194:   ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
    195:   ru[1] = 0;
    196:   PetscCall(MatSetValues(DRDU, 2, row, 1, col, ru, INSERT_VALUES));
    197:   PetscCall(VecRestoreArrayRead(U, &u));
    198:   PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
    199:   PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
    200:   PetscFunctionReturn(PETSC_SUCCESS);
    201: }

    203: PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, void *ctx)
    204: {
    205:   PetscFunctionBegin;
    206:   PetscCall(MatZeroEntries(DRDP));
    207:   PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
    208:   PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
    213: {
    214:   PetscScalar       *y, sensip;
    215:   const PetscScalar *x;

    217:   PetscFunctionBegin;
    218:   PetscCall(VecGetArrayRead(lambda, &x));
    219:   PetscCall(VecGetArray(mu, &y));
    220:   sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
    221:   y[0]   = sensip;
    222:   PetscCall(VecRestoreArray(mu, &y));
    223:   PetscCall(VecRestoreArrayRead(lambda, &x));
    224:   PetscFunctionReturn(PETSC_SUCCESS);
    225: }