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 PetscErrorCode ierr; 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 ierr = TSRestartStep(ts);CHKERRQ(ierr); /* 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 { 53 PetscErrorCode ierr; 54 PetscScalar *f,Pmax; 55 const PetscScalar *u; 56 57 PetscFunctionBegin; 58 /* The next three lines allow us to access the entries of the vectors directly */ 59 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 60 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 61 Pmax = ctx->Pmax; 62 f[0] = ctx->omega_b*(u[1] - ctx->omega_s); 63 f[1] = ctx->omega_s/(2.0*ctx->H)*(ctx->Pm - Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s)); 64 65 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 66 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 /* 71 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of a and the Jacobian. 72 */ 73 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx) 74 { 75 PetscErrorCode ierr; 76 PetscInt rowcol[] = {0,1}; 77 PetscScalar J[2][2],Pmax; 78 const PetscScalar *u; 79 80 PetscFunctionBegin; 81 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 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 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 88 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 89 90 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 91 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 92 if (A != B) { 93 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 94 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 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 PetscErrorCode ierr; 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 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 111 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 112 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 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 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 118 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 119 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 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 PetscErrorCode ierr; 129 PetscInt rowcol[] = {0,1}; 130 PetscScalar J[2][2],Pmax; 131 const PetscScalar *u,*udot; 132 133 PetscFunctionBegin; 134 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 135 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 136 Pmax = ctx->Pmax; 137 J[0][0] = a; J[0][1] = -ctx->omega_b; 138 J[1][1] = 2.0*ctx->H/ctx->omega_s*a + ctx->D; J[1][0] = Pmax*PetscCosScalar(u[0]); 139 140 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 141 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 142 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 143 144 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 146 if (A != B) { 147 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 148 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0) 154 { 155 PetscErrorCode ierr; 156 PetscInt row[] = {0,1},col[] = {0}; 157 PetscScalar *x,J[2][1]; 158 AppCtx *ctx = (AppCtx*)ctx0; 159 160 PetscFunctionBeginUser; 161 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 162 J[0][0] = 0; 163 J[1][0] = ctx->omega_s/(2.0*ctx->H); 164 ierr = MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 165 166 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 167 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx) 172 { 173 PetscErrorCode ierr; 174 PetscScalar *r; 175 const PetscScalar *u; 176 177 PetscFunctionBegin; 178 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 179 ierr = VecGetArray(R,&r);CHKERRQ(ierr); 180 r[0] = ctx->c*PetscPowScalarInt(PetscMax(0.,u[0]-ctx->u_s),ctx->beta);CHKERRQ(ierr); 181 ierr = VecRestoreArray(R,&r);CHKERRQ(ierr); 182 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 /* Transpose of DRDU */ 187 PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx) 188 { 189 PetscScalar ru[2]; 190 PetscInt row[] = {0,1},col[] = {0}; 191 const PetscScalar *u; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 196 ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0.,u[0]-ctx->u_s),ctx->beta-1);CHKERRQ(ierr); 197 ru[1] = 0; 198 ierr = MatSetValues(DRDU,2,row,1,col,ru,INSERT_VALUES);CHKERRQ(ierr); 199 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 200 ierr = MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 201 ierr = MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 205 PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,void *ctx) 206 { 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = MatZeroEntries(DRDP);CHKERRQ(ierr); 211 ierr = MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 212 ierr = MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx) 217 { 218 PetscErrorCode ierr; 219 PetscScalar *y,sensip; 220 const PetscScalar *x; 221 222 PetscFunctionBegin; 223 ierr = VecGetArrayRead(lambda,&x);CHKERRQ(ierr); 224 ierr = VecGetArray(mu,&y);CHKERRQ(ierr); 225 sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0]; 226 y[0] = sensip; 227 ierr = VecRestoreArray(mu,&y);CHKERRQ(ierr); 228 ierr = VecRestoreArrayRead(lambda,&x);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231