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:   PetscFunctionBegin;
 26:   /* Event for fault-on time */
 27:   fvalue[0] = t - user->tf;
 28:   /* Event for fault-off time */
 29:   fvalue[1] = t - user->tcl;

 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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