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