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