1 typedef enum { 2 SA_ADJ, 3 SA_TLM 4 } SAMethod; 5 static const char *const SAMethods[] = {"ADJ", "TLM", "SAMethod", "SA_", 0}; 6 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; 19 20 /* Event check */ 21 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec X, PetscScalar *fvalue, void *ctx) 22 { 23 AppCtx *user = (AppCtx *)ctx; 24 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; 30 31 PetscFunctionReturn(0); 32 } 33 34 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *ctx) 35 { 36 AppCtx *user = (AppCtx *)ctx; 37 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 ture, otherwise methods with FSAL will fail */ 47 PetscFunctionReturn(0); 48 } 49 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; 57 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)); 65 66 PetscCall(VecRestoreArrayRead(U, &u)); 67 PetscCall(VecRestoreArray(F, &f)); 68 PetscFunctionReturn(0); 69 } 70 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; 79 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)); 89 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(0); 97 } 98 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; 106 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; 115 116 PetscCall(VecRestoreArrayRead(U, &u)); 117 PetscCall(VecRestoreArrayRead(Udot, &udot)); 118 PetscCall(VecRestoreArray(F, &f)); 119 PetscFunctionReturn(0); 120 } 121 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; 130 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]); 139 140 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 141 PetscCall(VecRestoreArrayRead(U, &u)); 142 PetscCall(VecRestoreArrayRead(Udot, &udot)); 143 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(0); 151 } 152 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; 158 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)); 164 165 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 166 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 167 PetscFunctionReturn(0); 168 } 169 170 PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx) 171 { 172 PetscScalar *r; 173 const PetscScalar *u; 174 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(0); 182 } 183 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; 190 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(0); 200 } 201 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(0); 209 } 210 211 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx) 212 { 213 PetscScalar *y, sensip; 214 const PetscScalar *x; 215 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(0); 224 } 225