Actual source code: ex3.h

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

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

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

 25:   /* Event for fault-on time */
 26:   fvalue[0] = t - user->tf;
 27:   /* Event for fault-off time */
 28:   fvalue[1] = t - user->tcl;

 30:   return 0;
 31: }

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

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

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

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

 63:   VecRestoreArrayRead(U, &u);
 64:   VecRestoreArray(F, &f);
 65:   return 0;
 66: }

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

 77:   VecGetArrayRead(U, &u);
 78:   Pmax    = ctx->Pmax;
 79:   J[0][0] = 0;
 80:   J[0][1] = ctx->omega_b;
 81:   J[1][0] = -ctx->omega_s / (2.0 * ctx->H) * Pmax * PetscCosScalar(u[0]);
 82:   J[1][1] = -ctx->omega_s / (2.0 * ctx->H) * ctx->D;
 83:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
 84:   VecRestoreArrayRead(U, &u);

 86:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 87:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 88:   if (A != B) {
 89:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 90:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 91:   }
 92:   return 0;
 93: }

 95: /*
 96:      Defines the ODE passed to the ODE solver
 97: */
 98: PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
 99: {
100:   PetscScalar       *f, Pmax;
101:   const PetscScalar *u, *udot;

103:   /*  The next three lines allow us to access the entries of the vectors directly */
104:   VecGetArrayRead(U, &u);
105:   VecGetArrayRead(Udot, &udot);
106:   VecGetArray(F, &f);
107:   Pmax = ctx->Pmax;
108:   f[0] = udot[0] - ctx->omega_b * (u[1] - ctx->omega_s);
109:   f[1] = 2.0 * ctx->H / ctx->omega_s * udot[1] + Pmax * PetscSinScalar(u[0]) + ctx->D * (u[1] - ctx->omega_s) - ctx->Pm;

111:   VecRestoreArrayRead(U, &u);
112:   VecRestoreArrayRead(Udot, &udot);
113:   VecRestoreArray(F, &f);
114:   return 0;
115: }

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

126:   VecGetArrayRead(U, &u);
127:   VecGetArrayRead(Udot, &udot);
128:   Pmax    = ctx->Pmax;
129:   J[0][0] = a;
130:   J[0][1] = -ctx->omega_b;
131:   J[1][1] = 2.0 * ctx->H / ctx->omega_s * a + ctx->D;
132:   J[1][0] = Pmax * PetscCosScalar(u[0]);

134:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
135:   VecRestoreArrayRead(U, &u);
136:   VecRestoreArrayRead(Udot, &udot);

138:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
139:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
140:   if (A != B) {
141:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
142:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
143:   }
144:   return 0;
145: }

147: PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
148: {
149:   PetscInt     row[] = {0, 1}, col[] = {0};
150:   PetscScalar *x, J[2][1];
151:   AppCtx      *ctx = (AppCtx *)ctx0;

154:   VecGetArray(X, &x);
155:   J[0][0] = 0;
156:   J[1][0] = ctx->omega_s / (2.0 * ctx->H);
157:   MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES);

159:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
160:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
161:   return 0;
162: }

164: PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
165: {
166:   PetscScalar       *r;
167:   const PetscScalar *u;

169:   VecGetArrayRead(U, &u);
170:   VecGetArray(R, &r);
171:   r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
172:   VecRestoreArray(R, &r);
173:   VecRestoreArrayRead(U, &u);
174:   return 0;
175: }

177: /* Transpose of DRDU */
178: PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
179: {
180:   PetscScalar        ru[2];
181:   PetscInt           row[] = {0, 1}, col[] = {0};
182:   const PetscScalar *u;

184:   VecGetArrayRead(U, &u);
185:   ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
186:   ru[1] = 0;
187:   MatSetValues(DRDU, 2, row, 1, col, ru, INSERT_VALUES);
188:   VecRestoreArrayRead(U, &u);
189:   MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY);
190:   MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY);
191:   return 0;
192: }

194: PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, void *ctx)
195: {
196:   MatZeroEntries(DRDP);
197:   MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY);
198:   MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY);
199:   return 0;
200: }

202: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
203: {
204:   PetscScalar       *y, sensip;
205:   const PetscScalar *x;

207:   VecGetArrayRead(lambda, &x);
208:   VecGetArray(mu, &y);
209:   sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
210:   y[0]   = sensip;
211:   VecRestoreArray(mu, &y);
212:   VecRestoreArrayRead(lambda, &x);
213:   return 0;
214: }