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 12 13 Ensemble of initial conditions 14 ./ex2 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly 15 16 Fault at .1 seconds 17 ./ex2 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly 18 19 Initial conditions same as when fault is ended 20 ./ex2 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly 21 22 23 F*/ 24 25 /* 26 Include "petscts.h" so that we can use TS solvers. Note that this 27 file automatically includes: 28 petscsys.h - base PETSc routines petscvec.h - vectors 29 petscmat.h - matrices 30 petscis.h - index sets petscksp.h - Krylov subspace methods 31 petscviewer.h - viewers petscpc.h - preconditioners 32 petscksp.h - linear solvers 33 */ 34 35 #include <petscts.h> 36 37 typedef struct { 38 PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X; 39 PetscReal tf,tcl; 40 } AppCtx; 41 42 /* 43 Defines the ODE passed to the ODE solver 44 */ 45 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 46 { 47 PetscErrorCode ierr; 48 const PetscScalar *u; 49 PetscScalar *f,Pmax; 50 51 PetscFunctionBegin; 52 /* The next three lines allow us to access the entries of the vectors directly */ 53 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 54 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 55 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 */ 56 else Pmax = ctx->Pmax; 57 58 f[0] = ctx->omega_b*(u[1] - ctx->omega_s); 59 f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H); 60 61 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 62 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 /* 67 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 68 */ 69 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx) 70 { 71 PetscErrorCode ierr; 72 PetscInt rowcol[] = {0,1}; 73 PetscScalar J[2][2],Pmax; 74 const PetscScalar *u; 75 76 PetscFunctionBegin; 77 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 78 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 */ 79 else Pmax = ctx->Pmax; 80 81 J[0][0] = 0; J[0][1] = ctx->omega_b; 82 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); 83 84 ierr = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 85 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 86 87 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 88 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 89 if (A != B) { 90 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 91 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 92 } 93 PetscFunctionReturn(0); 94 } 95 96 int main(int argc,char **argv) 97 { 98 TS ts; /* ODE integrator */ 99 Vec U; /* solution will be stored here */ 100 Mat A; /* Jacobian matrix */ 101 PetscErrorCode ierr; 102 PetscMPIInt size; 103 PetscInt n = 2; 104 AppCtx ctx; 105 PetscScalar *u; 106 PetscReal du[2] = {0.0,0.0}; 107 PetscBool ensemble = PETSC_FALSE,flg1,flg2; 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Initialize program 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 113 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 114 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 115 116 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 117 Create necessary matrix and vectors 118 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 119 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 120 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 121 ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr); 122 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 123 ierr = MatSetUp(A);CHKERRQ(ierr); 124 125 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Set runtime options 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");CHKERRQ(ierr); 131 { 132 ctx.omega_b = 1.0; 133 ctx.omega_s = 2.0*PETSC_PI*60.0; 134 ctx.H = 5.0; 135 ierr = PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);CHKERRQ(ierr); 136 ctx.D = 5.0; 137 ierr = PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);CHKERRQ(ierr); 138 ctx.E = 1.1378; 139 ctx.V = 1.0; 140 ctx.X = 0.545; 141 ctx.Pmax = ctx.E*ctx.V/ctx.X; 142 ierr = PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);CHKERRQ(ierr); 143 ctx.Pm = 0.9; 144 ierr = PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);CHKERRQ(ierr); 145 ctx.tf = 1.0; 146 ctx.tcl = 1.05; 147 ierr = PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);CHKERRQ(ierr); 148 ierr = PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);CHKERRQ(ierr); 149 ierr = PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL);CHKERRQ(ierr); 150 if (ensemble) { 151 ctx.tf = -1; 152 ctx.tcl = -1; 153 } 154 155 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 156 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 157 u[1] = 1.0; 158 ierr = PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1);CHKERRQ(ierr); 159 n = 2; 160 ierr = PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2);CHKERRQ(ierr); 161 u[0] += du[0]; 162 u[1] += du[1]; 163 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 164 if (flg1 || flg2) { 165 ctx.tf = -1; 166 ctx.tcl = -1; 167 } 168 } 169 ierr = PetscOptionsEnd();CHKERRQ(ierr); 170 171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172 Create timestepping solver context 173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 174 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 175 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 176 ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr); 177 ierr = TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);CHKERRQ(ierr); 178 ierr = TSSetRHSJacobian(ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);CHKERRQ(ierr); 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Set initial conditions 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 184 185 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186 Set solver options 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 ierr = TSSetMaxTime(ts,35.0);CHKERRQ(ierr); 189 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 190 ierr = TSSetTimeStep(ts,.01);CHKERRQ(ierr); 191 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 192 193 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 Solve nonlinear system 195 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196 if (ensemble) { 197 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) { 198 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 199 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax); 200 u[1] = ctx.omega_s; 201 u[0] += du[0]; 202 u[1] += du[1]; 203 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 204 ierr = TSSetTimeStep(ts,.01);CHKERRQ(ierr); 205 ierr = TSSolve(ts,U);CHKERRQ(ierr); 206 } 207 } else { 208 ierr = TSSolve(ts,U);CHKERRQ(ierr); 209 } 210 ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 211 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 212 Free work space. All PETSc objects should be destroyed when they are no longer needed. 213 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 214 ierr = MatDestroy(&A);CHKERRQ(ierr); 215 ierr = VecDestroy(&U);CHKERRQ(ierr); 216 ierr = TSDestroy(&ts);CHKERRQ(ierr); 217 ierr = PetscFinalize(); 218 return ierr; 219 } 220 221 222 /*TEST 223 224 build: 225 requires: !complex 226 227 test: 228 229 TEST*/ 230