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