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