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