Actual source code: ex9adj.c
2: static char help[] = "Basic equation for generator stability analysis.\n";
\begin{eqnarray}
\frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
\frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
\end{eqnarray}
Ensemble of initial conditions
./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
Fault at .1 seconds
./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
Initial conditions same as when fault is ended
./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
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: */
32: #include <petscts.h>
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;
40: PetscErrorCode PostStepFunction(TS ts)
41: {
42: Vec U;
43: PetscReal t;
44: const PetscScalar *u;
46: TSGetTime(ts,&t);
47: TSGetSolution(ts,&U);
48: VecGetArrayRead(U,&u);
49: PetscPrintf(PETSC_COMM_SELF,"delta(%3.2f) = %8.7f\n",(double)t,(double)u[0]);
50: VecRestoreArrayRead(U,&u);
51: return 0;
52: }
54: /*
55: Defines the ODE passed to the ODE solver
56: */
57: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
58: {
59: PetscScalar *f,Pmax;
60: const PetscScalar *u;
62: /* The next three lines allow us to access the entries of the vectors directly */
63: VecGetArrayRead(U,&u);
64: VecGetArray(F,&f);
65: 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 */
66: else Pmax = ctx->Pmax;
68: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
69: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
71: VecRestoreArrayRead(U,&u);
72: VecRestoreArray(F,&f);
73: return 0;
74: }
76: /*
77: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
78: */
79: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
80: {
81: PetscInt rowcol[] = {0,1};
82: PetscScalar J[2][2],Pmax;
83: const PetscScalar *u;
85: VecGetArrayRead(U,&u);
86: 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 */
87: else Pmax = ctx->Pmax;
89: J[0][0] = 0; J[0][1] = ctx->omega_b;
90: 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);
92: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
93: VecRestoreArrayRead(U,&u);
95: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
97: if (A != B) {
98: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
100: }
101: return 0;
102: }
104: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
105: {
106: PetscInt row[] = {0,1},col[]={0};
107: PetscScalar J[2][1];
108: AppCtx *ctx=(AppCtx*)ctx0;
111: J[0][0] = 0;
112: J[1][0] = ctx->omega_s/(2.0*ctx->H);
113: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
114: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
116: return 0;
117: }
119: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
120: {
121: PetscScalar *r;
122: const PetscScalar *u;
124: VecGetArrayRead(U,&u);
125: VecGetArray(R,&r);
126: r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta);
127: VecRestoreArray(R,&r);
128: VecRestoreArrayRead(U,&u);
129: return 0;
130: }
132: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx)
133: {
134: PetscScalar ru[1];
135: const PetscScalar *u;
136: PetscInt row[] = {0},col[] = {0};
138: VecGetArrayRead(U,&u);
139: ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1);
140: VecRestoreArrayRead(U,&u);
141: MatSetValues(DRDU,1,row,1,col,ru,INSERT_VALUES);
142: MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
143: MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
144: return 0;
145: }
147: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx)
148: {
149: MatZeroEntries(DRDP);
150: MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
152: return 0;
153: }
155: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
156: {
157: PetscScalar sensip;
158: const PetscScalar *x,*y;
160: VecGetArrayRead(lambda,&x);
161: VecGetArrayRead(mu,&y);
162: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
163: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameter pm: %.7f \n",(double)sensip);
164: VecRestoreArrayRead(lambda,&x);
165: VecRestoreArrayRead(mu,&y);
166: return 0;
167: }
169: int main(int argc,char **argv)
170: {
171: TS ts,quadts; /* ODE integrator */
172: Vec U; /* solution will be stored here */
173: Mat A; /* Jacobian matrix */
174: Mat Jacp; /* Jacobian matrix */
175: Mat DRDU,DRDP;
177: PetscMPIInt size;
178: PetscInt n = 2;
179: AppCtx ctx;
180: PetscScalar *u;
181: PetscReal du[2] = {0.0,0.0};
182: PetscBool ensemble = PETSC_FALSE,flg1,flg2;
183: PetscReal ftime;
184: PetscInt steps;
185: PetscScalar *x_ptr,*y_ptr;
186: Vec lambda[1],q,mu[1];
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Initialize program
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: PetscInitialize(&argc,&argv,(char*)0,help);
192: MPI_Comm_size(PETSC_COMM_WORLD,&size);
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Create necessary matrix and vectors
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: MatCreate(PETSC_COMM_WORLD,&A);
199: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
200: MatSetType(A,MATDENSE);
201: MatSetFromOptions(A);
202: MatSetUp(A);
204: MatCreateVecs(A,&U,NULL);
206: MatCreate(PETSC_COMM_WORLD,&Jacp);
207: MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
208: MatSetFromOptions(Jacp);
209: MatSetUp(Jacp);
211: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&DRDP);
212: MatSetUp(DRDP);
213: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,2,NULL,&DRDU);
214: MatSetUp(DRDU);
216: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: Set runtime options
218: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
219: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
220: {
221: ctx.beta = 2;
222: ctx.c = 10000.0;
223: ctx.u_s = 1.0;
224: ctx.omega_s = 1.0;
225: ctx.omega_b = 120.0*PETSC_PI;
226: ctx.H = 5.0;
227: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
228: ctx.D = 5.0;
229: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
230: ctx.E = 1.1378;
231: ctx.V = 1.0;
232: ctx.X = 0.545;
233: ctx.Pmax = ctx.E*ctx.V/ctx.X;
234: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
235: ctx.Pm = 1.1;
236: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
237: ctx.tf = 0.1;
238: ctx.tcl = 0.2;
239: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
240: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
241: PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL);
242: if (ensemble) {
243: ctx.tf = -1;
244: ctx.tcl = -1;
245: }
247: VecGetArray(U,&u);
248: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
249: u[1] = 1.0;
250: PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1);
251: n = 2;
252: PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2);
253: u[0] += du[0];
254: u[1] += du[1];
255: VecRestoreArray(U,&u);
256: if (flg1 || flg2) {
257: ctx.tf = -1;
258: ctx.tcl = -1;
259: }
260: }
261: PetscOptionsEnd();
263: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: Create timestepping solver context
265: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
266: TSCreate(PETSC_COMM_WORLD,&ts);
267: TSSetProblemType(ts,TS_NONLINEAR);
268: TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
269: TSSetType(ts,TSRK);
270: TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
271: TSSetRHSJacobian(ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
272: TSCreateQuadratureTS(ts,PETSC_TRUE,&quadts);
273: TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
274: TSSetRHSJacobian(quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
275: TSSetRHSJacobianP(quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx);
277: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: Set initial conditions
279: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
280: TSSetSolution(ts,U);
282: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
283: Save trajectory of solution so that TSAdjointSolve() may be used
284: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
285: TSSetSaveTrajectory(ts);
287: MatCreateVecs(A,&lambda[0],NULL);
288: /* Set initial conditions for the adjoint integration */
289: VecGetArray(lambda[0],&y_ptr);
290: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
291: VecRestoreArray(lambda[0],&y_ptr);
293: MatCreateVecs(Jacp,&mu[0],NULL);
294: VecGetArray(mu[0],&x_ptr);
295: x_ptr[0] = -1.0;
296: VecRestoreArray(mu[0],&x_ptr);
297: TSSetCostGradients(ts,1,lambda,mu);
299: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300: Set solver options
301: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
302: TSSetMaxTime(ts,10.0);
303: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
304: TSSetTimeStep(ts,.01);
305: TSSetFromOptions(ts);
307: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308: Solve nonlinear system
309: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
310: if (ensemble) {
311: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
312: VecGetArray(U,&u);
313: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
314: u[1] = ctx.omega_s;
315: u[0] += du[0];
316: u[1] += du[1];
317: VecRestoreArray(U,&u);
318: TSSetTimeStep(ts,.01);
319: TSSolve(ts,U);
320: }
321: } else {
322: TSSolve(ts,U);
323: }
324: VecView(U,PETSC_VIEWER_STDOUT_WORLD);
325: TSGetSolveTime(ts,&ftime);
326: TSGetStepNumber(ts,&steps);
328: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
329: Adjoint model starts here
330: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
331: /* Set initial conditions for the adjoint integration */
332: VecGetArray(lambda[0],&y_ptr);
333: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
334: VecRestoreArray(lambda[0],&y_ptr);
336: VecGetArray(mu[0],&x_ptr);
337: x_ptr[0] = -1.0;
338: VecRestoreArray(mu[0],&x_ptr);
340: /* Set RHS JacobianP */
341: TSSetRHSJacobianP(ts,Jacp,RHSJacobianP,&ctx);
343: TSAdjointSolve(ts);
345: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n");
346: VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);
347: VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD);
348: TSGetCostIntegral(ts,&q);
349: VecView(q,PETSC_VIEWER_STDOUT_WORLD);
350: VecGetArray(q,&x_ptr);
351: PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm));
352: VecRestoreArray(q,&x_ptr);
354: ComputeSensiP(lambda[0],mu[0],&ctx);
356: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: Free work space. All PETSc objects should be destroyed when they are no longer needed.
358: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
359: MatDestroy(&A);
360: MatDestroy(&Jacp);
361: MatDestroy(&DRDU);
362: MatDestroy(&DRDP);
363: VecDestroy(&U);
364: VecDestroy(&lambda[0]);
365: VecDestroy(&mu[0]);
366: TSDestroy(&ts);
367: PetscFinalize();
368: return 0;
369: }
371: /*TEST
373: build:
374: requires: !complex
376: test:
377: args: -viewer_binary_skip_info -ts_adapt_type none
379: TEST*/