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 */ 23 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec X, PetscReal *fvalue, void *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 33 PetscFunctionReturn(PETSC_SUCCESS); 34 } 35 36 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *ctx) 37 { 38 AppCtx *user = (AppCtx *)ctx; 39 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 } 51 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; 59 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)); 67 68 PetscCall(VecRestoreArrayRead(U, &u)); 69 PetscCall(VecRestoreArray(F, &f)); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 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; 81 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)); 91 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 } 100 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; 108 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; 117 118 PetscCall(VecRestoreArrayRead(U, &u)); 119 PetscCall(VecRestoreArrayRead(Udot, &udot)); 120 PetscCall(VecRestoreArray(F, &f)); 121 PetscFunctionReturn(PETSC_SUCCESS); 122 } 123 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; 132 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]); 141 142 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 143 PetscCall(VecRestoreArrayRead(U, &u)); 144 PetscCall(VecRestoreArrayRead(Udot, &udot)); 145 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 } 154 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; 160 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)); 166 167 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 168 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx) 173 { 174 PetscScalar *r; 175 const PetscScalar *u; 176 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 } 185 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; 192 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 } 203 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 } 212 213 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx) 214 { 215 PetscScalar *y, sensip; 216 const PetscScalar *x; 217 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 } 227