1 #pragma once
2
3 typedef enum {
4 SA_ADJ,
5 SA_TLM
6 } SAMethod;
7 static const char *const SAMethods[] = {"ADJ", "TLM", "SAMethod", "SA_", 0};
8
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;
21
22 /* Event check */
EventFunction(TS ts,PetscReal t,Vec X,PetscReal * fvalue,PetscCtx ctx)23 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec X, PetscReal *fvalue, PetscCtx ctx)
24 {
25 AppCtx *user = (AppCtx *)ctx;
26
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 }
34
PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec X,PetscBool forwardsolve,PetscCtx ctx)35 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, PetscCtx ctx)
36 {
37 AppCtx *user = (AppCtx *)ctx;
38
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 }
50
51 /*
52 Defines the ODE passed to the ODE solver
53 */
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)54 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
55 {
56 PetscScalar *f, Pmax;
57 const PetscScalar *u;
58
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));
66
67 PetscCall(VecRestoreArrayRead(U, &u));
68 PetscCall(VecRestoreArray(F, &f));
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71
72 /*
73 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of a and the Jacobian.
74 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx * ctx)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;
80
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));
90
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 }
99
100 /*
101 Defines the ODE passed to the ODE solver
102 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)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;
107
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;
116
117 PetscCall(VecRestoreArrayRead(U, &u));
118 PetscCall(VecRestoreArrayRead(Udot, &udot));
119 PetscCall(VecRestoreArray(F, &f));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
123 /*
124 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
125 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)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;
131
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]);
140
141 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
142 PetscCall(VecRestoreArrayRead(U, &u));
143 PetscCall(VecRestoreArrayRead(Udot, &udot));
144
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 }
153
RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,PetscCtx ctx0)154 PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx0)
155 {
156 PetscInt row[] = {0, 1}, col[] = {0};
157 PetscScalar *x, J[2][1];
158 AppCtx *ctx = (AppCtx *)ctx0;
159
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));
165
166 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
167 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
168 PetscFunctionReturn(PETSC_SUCCESS);
169 }
170
CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx * ctx)171 PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
172 {
173 PetscScalar *r;
174 const PetscScalar *u;
175
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 }
184
185 /* Transpose of DRDU */
DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx * ctx)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;
191
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 }
202
DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,PetscCtx ctx)203 PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, PetscCtx 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 }
211
ComputeSensiP(Vec lambda,Vec mu,AppCtx * ctx)212 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
213 {
214 PetscScalar *y, sensip;
215 const PetscScalar *x;
216
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 }
226