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