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: }