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 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *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 */ 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 */ 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 */ 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 */ 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 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; 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 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 */ 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 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 } 211 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