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, PetscScalar *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;
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *ctx)
37: {
38: AppCtx *user = (AppCtx *)ctx;
40: PetscFunctionBegin;
41: if (event_list[0] == 0) {
42: if (forwardsolve) user->Pmax = 0.0; /* Apply disturbance - this is done by setting Pmax = 0 */
43: else user->Pmax = user->Pmax_ini; /* Going backward, reversal of event */
44: } else if (event_list[0] == 1) {
45: if (forwardsolve) user->Pmax = user->Pmax_ini; /* Remove the fault - this is done by setting Pmax = Pmax_ini */
46: else user->Pmax = 0.0; /* Going backward, reversal of event */
47: }
48: PetscCall(TSRestartStep(ts)); /* Must set restart flag to true, otherwise methods with FSAL will fail */
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*
53: Defines the ODE passed to the ODE solver
54: */
55: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
56: {
57: PetscScalar *f, Pmax;
58: const PetscScalar *u;
60: PetscFunctionBegin;
61: /* The next three lines allow us to access the entries of the vectors directly */
62: PetscCall(VecGetArrayRead(U, &u));
63: PetscCall(VecGetArray(F, &f));
64: Pmax = ctx->Pmax;
65: f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
66: f[1] = ctx->omega_s / (2.0 * ctx->H) * (ctx->Pm - Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s));
68: PetscCall(VecRestoreArrayRead(U, &u));
69: PetscCall(VecRestoreArray(F, &f));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*
74: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of a and the Jacobian.
75: */
76: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
77: {
78: PetscInt rowcol[] = {0, 1};
79: PetscScalar J[2][2], Pmax;
80: const PetscScalar *u;
82: PetscFunctionBegin;
83: PetscCall(VecGetArrayRead(U, &u));
84: Pmax = ctx->Pmax;
85: J[0][0] = 0;
86: J[0][1] = ctx->omega_b;
87: J[1][0] = -ctx->omega_s / (2.0 * ctx->H) * Pmax * PetscCosScalar(u[0]);
88: J[1][1] = -ctx->omega_s / (2.0 * ctx->H) * ctx->D;
89: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
90: PetscCall(VecRestoreArrayRead(U, &u));
92: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
93: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
94: if (A != B) {
95: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
96: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
97: }
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*
102: Defines the ODE passed to the ODE solver
103: */
104: PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
105: {
106: PetscScalar *f, Pmax;
107: const PetscScalar *u, *udot;
109: PetscFunctionBegin;
110: /* The next three lines allow us to access the entries of the vectors directly */
111: PetscCall(VecGetArrayRead(U, &u));
112: PetscCall(VecGetArrayRead(Udot, &udot));
113: PetscCall(VecGetArray(F, &f));
114: Pmax = ctx->Pmax;
115: f[0] = udot[0] - ctx->omega_b * (u[1] - ctx->omega_s);
116: f[1] = 2.0 * ctx->H / ctx->omega_s * udot[1] + Pmax * PetscSinScalar(u[0]) + ctx->D * (u[1] - ctx->omega_s) - ctx->Pm;
118: PetscCall(VecRestoreArrayRead(U, &u));
119: PetscCall(VecRestoreArrayRead(Udot, &udot));
120: PetscCall(VecRestoreArray(F, &f));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /*
125: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
126: */
127: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
128: {
129: PetscInt rowcol[] = {0, 1};
130: PetscScalar J[2][2], Pmax;
131: const PetscScalar *u, *udot;
133: PetscFunctionBegin;
134: PetscCall(VecGetArrayRead(U, &u));
135: PetscCall(VecGetArrayRead(Udot, &udot));
136: Pmax = ctx->Pmax;
137: J[0][0] = a;
138: J[0][1] = -ctx->omega_b;
139: J[1][1] = 2.0 * ctx->H / ctx->omega_s * a + ctx->D;
140: J[1][0] = Pmax * PetscCosScalar(u[0]);
142: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
143: PetscCall(VecRestoreArrayRead(U, &u));
144: PetscCall(VecRestoreArrayRead(Udot, &udot));
146: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
147: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
148: if (A != B) {
149: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
150: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
156: {
157: PetscInt row[] = {0, 1}, col[] = {0};
158: PetscScalar *x, J[2][1];
159: AppCtx *ctx = (AppCtx *)ctx0;
161: PetscFunctionBeginUser;
162: PetscCall(VecGetArray(X, &x));
163: J[0][0] = 0;
164: J[1][0] = ctx->omega_s / (2.0 * ctx->H);
165: PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
167: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
168: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
173: {
174: PetscScalar *r;
175: const PetscScalar *u;
177: PetscFunctionBegin;
178: PetscCall(VecGetArrayRead(U, &u));
179: PetscCall(VecGetArray(R, &r));
180: r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
181: PetscCall(VecRestoreArray(R, &r));
182: PetscCall(VecRestoreArrayRead(U, &u));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /* Transpose of DRDU */
187: PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
188: {
189: PetscScalar ru[2];
190: PetscInt row[] = {0, 1}, col[] = {0};
191: const PetscScalar *u;
193: PetscFunctionBegin;
194: PetscCall(VecGetArrayRead(U, &u));
195: ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
196: ru[1] = 0;
197: PetscCall(MatSetValues(DRDU, 2, row, 1, col, ru, INSERT_VALUES));
198: PetscCall(VecRestoreArrayRead(U, &u));
199: PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
200: PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, void *ctx)
205: {
206: PetscFunctionBegin;
207: PetscCall(MatZeroEntries(DRDP));
208: PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
209: PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
214: {
215: PetscScalar *y, sensip;
216: const PetscScalar *x;
218: PetscFunctionBegin;
219: PetscCall(VecGetArrayRead(lambda, &x));
220: PetscCall(VecGetArray(mu, &y));
221: sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
222: y[0] = sensip;
223: PetscCall(VecRestoreArray(mu, &y));
224: PetscCall(VecRestoreArrayRead(lambda, &x));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }