1 2 static char help[] = "Basic equation for generator stability analysis.\n"; 3 4 /*F 5 6 \begin{eqnarray} 7 \frac{d \theta}{dt} = \omega_b (\omega - \omega_s) 8 \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\ 9 \end{eqnarray} 10 11 Ensemble of initial conditions 12 ./ex9 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly 13 14 Fault at .1 seconds 15 ./ex9 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly 16 17 Initial conditions same as when fault is ended 18 ./ex9 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly 19 20 F*/ 21 22 /* 23 Include "petscts.h" so that we can use TS solvers. Note that this 24 file automatically includes: 25 petscsys.h - base PETSc routines petscvec.h - vectors 26 petscmat.h - matrices 27 petscis.h - index sets petscksp.h - Krylov subspace methods 28 petscviewer.h - viewers petscpc.h - preconditioners 29 petscksp.h - linear solvers 30 */ 31 32 #include <petscts.h> 33 34 typedef struct { 35 PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X,u_s,c; 36 PetscInt beta; 37 PetscReal tf,tcl; 38 } AppCtx; 39 40 PetscErrorCode PostStepFunction(TS ts) 41 { 42 Vec U; 43 PetscReal t; 44 const PetscScalar *u; 45 46 PetscFunctionBegin; 47 PetscCall(TSGetTime(ts,&t)); 48 PetscCall(TSGetSolution(ts,&U)); 49 PetscCall(VecGetArrayRead(U,&u)); 50 PetscCall(PetscPrintf(PETSC_COMM_SELF,"delta(%3.2f) = %8.7f\n",(double)t,(double)u[0])); 51 PetscCall(VecRestoreArrayRead(U,&u)); 52 PetscFunctionReturn(0); 53 } 54 55 /* 56 Defines the ODE passed to the ODE solver 57 */ 58 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 59 { 60 PetscScalar *f,Pmax; 61 const PetscScalar *u; 62 63 PetscFunctionBegin; 64 /* The next three lines allow us to access the entries of the vectors directly */ 65 PetscCall(VecGetArrayRead(U,&u)); 66 PetscCall(VecGetArray(F,&f)); 67 if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */ 68 else Pmax = ctx->Pmax; 69 70 f[0] = ctx->omega_b*(u[1] - ctx->omega_s); 71 f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H); 72 73 PetscCall(VecRestoreArrayRead(U,&u)); 74 PetscCall(VecRestoreArray(F,&f)); 75 PetscFunctionReturn(0); 76 } 77 78 /* 79 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 80 */ 81 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx) 82 { 83 PetscInt rowcol[] = {0,1}; 84 PetscScalar J[2][2],Pmax; 85 const PetscScalar *u; 86 87 PetscFunctionBegin; 88 PetscCall(VecGetArrayRead(U,&u)); 89 if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */ 90 else Pmax = ctx->Pmax; 91 92 J[0][0] = 0; J[0][1] = ctx->omega_b; 93 J[1][1] = -ctx->D*ctx->omega_s/(2.0*ctx->H); J[1][0] = -Pmax*PetscCosScalar(u[0])*ctx->omega_s/(2.0*ctx->H); 94 95 PetscCall(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 96 PetscCall(VecRestoreArrayRead(U,&u)); 97 98 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 99 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 100 if (A != B) { 101 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 102 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0) 108 { 109 PetscInt row[] = {0,1},col[]={0}; 110 PetscScalar J[2][1]; 111 AppCtx *ctx=(AppCtx*)ctx0; 112 113 PetscFunctionBeginUser; 114 J[0][0] = 0; 115 J[1][0] = ctx->omega_s/(2.0*ctx->H); 116 PetscCall(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 117 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 118 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 119 PetscFunctionReturn(0); 120 } 121 122 static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx) 123 { 124 PetscScalar *r; 125 const PetscScalar *u; 126 127 PetscFunctionBegin; 128 PetscCall(VecGetArrayRead(U,&u)); 129 PetscCall(VecGetArray(R,&r)); 130 r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta); 131 PetscCall(VecRestoreArray(R,&r)); 132 PetscCall(VecRestoreArrayRead(U,&u)); 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx) 137 { 138 PetscScalar ru[1]; 139 const PetscScalar *u; 140 PetscInt row[] = {0},col[] = {0}; 141 142 PetscFunctionBegin; 143 PetscCall(VecGetArrayRead(U,&u)); 144 ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1); 145 PetscCall(VecRestoreArrayRead(U,&u)); 146 PetscCall(MatSetValues(DRDU,1,row,1,col,ru,INSERT_VALUES)); 147 PetscCall(MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY)); 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx) 153 { 154 PetscFunctionBegin; 155 PetscCall(MatZeroEntries(DRDP)); 156 PetscCall(MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY)); 157 PetscCall(MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY)); 158 PetscFunctionReturn(0); 159 } 160 161 PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx) 162 { 163 PetscScalar sensip; 164 const PetscScalar *x,*y; 165 166 PetscFunctionBegin; 167 PetscCall(VecGetArrayRead(lambda,&x)); 168 PetscCall(VecGetArrayRead(mu,&y)); 169 sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0]; 170 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameter pm: %.7f \n",(double)sensip)); 171 PetscCall(VecRestoreArrayRead(lambda,&x)); 172 PetscCall(VecRestoreArrayRead(mu,&y)); 173 PetscFunctionReturn(0); 174 } 175 176 int main(int argc,char **argv) 177 { 178 TS ts,quadts; /* ODE integrator */ 179 Vec U; /* solution will be stored here */ 180 Mat A; /* Jacobian matrix */ 181 Mat Jacp; /* Jacobian matrix */ 182 Mat DRDU,DRDP; 183 PetscErrorCode ierr; 184 PetscMPIInt size; 185 PetscInt n = 2; 186 AppCtx ctx; 187 PetscScalar *u; 188 PetscReal du[2] = {0.0,0.0}; 189 PetscBool ensemble = PETSC_FALSE,flg1,flg2; 190 PetscReal ftime; 191 PetscInt steps; 192 PetscScalar *x_ptr,*y_ptr; 193 Vec lambda[1],q,mu[1]; 194 195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196 Initialize program 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 199 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 200 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 201 202 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 203 Create necessary matrix and vectors 204 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 205 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 206 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 207 PetscCall(MatSetType(A,MATDENSE)); 208 PetscCall(MatSetFromOptions(A)); 209 PetscCall(MatSetUp(A)); 210 211 PetscCall(MatCreateVecs(A,&U,NULL)); 212 213 PetscCall(MatCreate(PETSC_COMM_WORLD,&Jacp)); 214 PetscCall(MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1)); 215 PetscCall(MatSetFromOptions(Jacp)); 216 PetscCall(MatSetUp(Jacp)); 217 218 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&DRDP)); 219 PetscCall(MatSetUp(DRDP)); 220 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,2,NULL,&DRDU)); 221 PetscCall(MatSetUp(DRDU)); 222 223 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 224 Set runtime options 225 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 226 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");PetscCall(ierr); 227 { 228 ctx.beta = 2; 229 ctx.c = 10000.0; 230 ctx.u_s = 1.0; 231 ctx.omega_s = 1.0; 232 ctx.omega_b = 120.0*PETSC_PI; 233 ctx.H = 5.0; 234 PetscCall(PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL)); 235 ctx.D = 5.0; 236 PetscCall(PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL)); 237 ctx.E = 1.1378; 238 ctx.V = 1.0; 239 ctx.X = 0.545; 240 ctx.Pmax = ctx.E*ctx.V/ctx.X; 241 PetscCall(PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL)); 242 ctx.Pm = 1.1; 243 PetscCall(PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL)); 244 ctx.tf = 0.1; 245 ctx.tcl = 0.2; 246 PetscCall(PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL)); 247 PetscCall(PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL)); 248 PetscCall(PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL)); 249 if (ensemble) { 250 ctx.tf = -1; 251 ctx.tcl = -1; 252 } 253 254 PetscCall(VecGetArray(U,&u)); 255 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 256 u[1] = 1.0; 257 PetscCall(PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1)); 258 n = 2; 259 PetscCall(PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2)); 260 u[0] += du[0]; 261 u[1] += du[1]; 262 PetscCall(VecRestoreArray(U,&u)); 263 if (flg1 || flg2) { 264 ctx.tf = -1; 265 ctx.tcl = -1; 266 } 267 } 268 ierr = PetscOptionsEnd();PetscCall(ierr); 269 270 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 271 Create timestepping solver context 272 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 273 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 274 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 275 PetscCall(TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 276 PetscCall(TSSetType(ts,TSRK)); 277 PetscCall(TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx)); 278 PetscCall(TSSetRHSJacobian(ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx)); 279 PetscCall(TSCreateQuadratureTS(ts,PETSC_TRUE,&quadts)); 280 PetscCall(TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx)); 281 PetscCall(TSSetRHSJacobian(quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx)); 282 PetscCall(TSSetRHSJacobianP(quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx)); 283 284 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 285 Set initial conditions 286 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 287 PetscCall(TSSetSolution(ts,U)); 288 289 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 290 Save trajectory of solution so that TSAdjointSolve() may be used 291 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 292 PetscCall(TSSetSaveTrajectory(ts)); 293 294 PetscCall(MatCreateVecs(A,&lambda[0],NULL)); 295 /* Set initial conditions for the adjoint integration */ 296 PetscCall(VecGetArray(lambda[0],&y_ptr)); 297 y_ptr[0] = 0.0; y_ptr[1] = 0.0; 298 PetscCall(VecRestoreArray(lambda[0],&y_ptr)); 299 300 PetscCall(MatCreateVecs(Jacp,&mu[0],NULL)); 301 PetscCall(VecGetArray(mu[0],&x_ptr)); 302 x_ptr[0] = -1.0; 303 PetscCall(VecRestoreArray(mu[0],&x_ptr)); 304 PetscCall(TSSetCostGradients(ts,1,lambda,mu)); 305 306 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 307 Set solver options 308 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 309 PetscCall(TSSetMaxTime(ts,10.0)); 310 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 311 PetscCall(TSSetTimeStep(ts,.01)); 312 PetscCall(TSSetFromOptions(ts)); 313 314 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 315 Solve nonlinear system 316 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 317 if (ensemble) { 318 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) { 319 PetscCall(VecGetArray(U,&u)); 320 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 321 u[1] = ctx.omega_s; 322 u[0] += du[0]; 323 u[1] += du[1]; 324 PetscCall(VecRestoreArray(U,&u)); 325 PetscCall(TSSetTimeStep(ts,.01)); 326 PetscCall(TSSolve(ts,U)); 327 } 328 } else { 329 PetscCall(TSSolve(ts,U)); 330 } 331 PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 332 PetscCall(TSGetSolveTime(ts,&ftime)); 333 PetscCall(TSGetStepNumber(ts,&steps)); 334 335 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 336 Adjoint model starts here 337 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 338 /* Set initial conditions for the adjoint integration */ 339 PetscCall(VecGetArray(lambda[0],&y_ptr)); 340 y_ptr[0] = 0.0; y_ptr[1] = 0.0; 341 PetscCall(VecRestoreArray(lambda[0],&y_ptr)); 342 343 PetscCall(VecGetArray(mu[0],&x_ptr)); 344 x_ptr[0] = -1.0; 345 PetscCall(VecRestoreArray(mu[0],&x_ptr)); 346 347 /* Set RHS JacobianP */ 348 PetscCall(TSSetRHSJacobianP(ts,Jacp,RHSJacobianP,&ctx)); 349 350 PetscCall(TSAdjointSolve(ts)); 351 352 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n")); 353 PetscCall(VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD)); 354 PetscCall(VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD)); 355 PetscCall(TSGetCostIntegral(ts,&q)); 356 PetscCall(VecView(q,PETSC_VIEWER_STDOUT_WORLD)); 357 PetscCall(VecGetArray(q,&x_ptr)); 358 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm))); 359 PetscCall(VecRestoreArray(q,&x_ptr)); 360 361 PetscCall(ComputeSensiP(lambda[0],mu[0],&ctx)); 362 363 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 364 Free work space. All PETSc objects should be destroyed when they are no longer needed. 365 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 366 PetscCall(MatDestroy(&A)); 367 PetscCall(MatDestroy(&Jacp)); 368 PetscCall(MatDestroy(&DRDU)); 369 PetscCall(MatDestroy(&DRDP)); 370 PetscCall(VecDestroy(&U)); 371 PetscCall(VecDestroy(&lambda[0])); 372 PetscCall(VecDestroy(&mu[0])); 373 PetscCall(TSDestroy(&ts)); 374 PetscCall(PetscFinalize()); 375 return 0; 376 } 377 378 /*TEST 379 380 build: 381 requires: !complex 382 383 test: 384 args: -viewer_binary_skip_info -ts_adapt_type none 385 386 TEST*/ 387