Actual source code: ex1.c
2: static char help[] = "Basic equation for generator stability analysis.\n";
\begin{eqnarray}
\frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - \frac{EV}{X} \sin(\theta) \\
\frac{d \theta}{dt} = \omega - \omega_s
\end{eqnarray}
13: /*
14: Include "petscts.h" so that we can use TS solvers. Note that this
15: file automatically includes:
16: petscsys.h - base PETSc routines petscvec.h - vectors
17: petscmat.h - matrices
18: petscis.h - index sets petscksp.h - Krylov subspace methods
19: petscviewer.h - viewers petscpc.h - preconditioners
20: petscksp.h - linear solvers
21: */
23: #include <petscts.h>
25: typedef struct {
26: PetscScalar H,omega_s,E,V,X;
27: PetscRandom rand;
28: } AppCtx;
30: /*
31: Defines the ODE passed to the ODE solver
32: */
33: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
34: {
35: PetscScalar *f,r;
36: const PetscScalar *u,*udot;
37: static PetscScalar R = .4;
39: PetscRandomGetValue(ctx->rand,&r);
40: if (r > .9) R = .5;
41: if (r < .1) R = .4;
42: R = .4;
43: /* The next three lines allow us to access the entries of the vectors directly */
44: VecGetArrayRead(U,&u);
45: VecGetArrayRead(Udot,&udot);
46: VecGetArray(F,&f);
47: f[0] = 2.0*ctx->H*udot[0]/ctx->omega_s + ctx->E*ctx->V*PetscSinScalar(u[1])/ctx->X - R;
48: f[1] = udot[1] - u[0] + ctx->omega_s;
50: VecRestoreArrayRead(U,&u);
51: VecRestoreArrayRead(Udot,&udot);
52: VecRestoreArray(F,&f);
53: return 0;
54: }
56: /*
57: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
58: */
59: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
60: {
61: PetscInt rowcol[] = {0,1};
62: PetscScalar J[2][2];
63: const PetscScalar *u,*udot;
65: VecGetArrayRead(U,&u);
66: VecGetArrayRead(Udot,&udot);
67: J[0][0] = 2.0*ctx->H*a/ctx->omega_s; J[0][1] = -ctx->E*ctx->V*PetscCosScalar(u[1])/ctx->X;
68: J[1][0] = -1.0; J[1][1] = a;
69: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
70: VecRestoreArrayRead(U,&u);
71: VecRestoreArrayRead(Udot,&udot);
73: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
75: if (A != B) {
76: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
78: }
79: return 0;
80: }
82: int main(int argc,char **argv)
83: {
84: TS ts; /* ODE integrator */
85: Vec U; /* solution will be stored here */
86: Mat A; /* Jacobian matrix */
88: PetscMPIInt size;
89: PetscInt n = 2;
90: AppCtx ctx;
91: PetscScalar *u;
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Initialize program
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: PetscInitialize(&argc,&argv,(char*)0,help);
97: MPI_Comm_size(PETSC_COMM_WORLD,&size);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create necessary matrix and vectors
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: MatCreate(PETSC_COMM_WORLD,&A);
104: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
105: MatSetFromOptions(A);
106: MatSetUp(A);
108: MatCreateVecs(A,&U,NULL);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Set runtime options
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Reaction options","");
114: {
115: ctx.omega_s = 1.0;
116: PetscOptionsScalar("-omega_s","","",ctx.omega_s,&ctx.omega_s,NULL);
117: ctx.H = 1.0;
118: PetscOptionsScalar("-H","","",ctx.H,&ctx.H,NULL);
119: ctx.E = 1.0;
120: PetscOptionsScalar("-E","","",ctx.E,&ctx.E,NULL);
121: ctx.V = 1.0;
122: PetscOptionsScalar("-V","","",ctx.V,&ctx.V,NULL);
123: ctx.X = 1.0;
124: PetscOptionsScalar("-X","","",ctx.X,&ctx.X,NULL);
126: VecGetArray(U,&u);
127: u[0] = 1;
128: u[1] = .7;
129: VecRestoreArray(U,&u);
130: PetscOptionsGetVec(NULL,NULL,"-initial",U,NULL);
131: }
132: PetscOptionsEnd();
134: PetscRandomCreate(PETSC_COMM_WORLD,&ctx.rand);
135: PetscRandomSetFromOptions(ctx.rand);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create timestepping solver context
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: TSCreate(PETSC_COMM_WORLD,&ts);
141: TSSetProblemType(ts,TS_NONLINEAR);
142: TSSetType(ts,TSROSW);
143: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);
144: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Set initial conditions
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: TSSetSolution(ts,U);
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Set solver options
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: TSSetMaxTime(ts,2000.0);
155: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
156: TSSetTimeStep(ts,.001);
157: TSSetFromOptions(ts);
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Solve nonlinear system
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: TSSolve(ts,U);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Free work space. All PETSc objects should be destroyed when they are no longer needed.
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: MatDestroy(&A);
168: VecDestroy(&U);
169: TSDestroy(&ts);
170: PetscRandomDestroy(&ctx.rand);
171: PetscFinalize();
172: return 0;
173: }
175: /*TEST
177: build:
178: requires: !complex !single
180: test:
181: args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
183: test:
184: suffix: 2
185: args: -ts_max_steps 10
187: TEST*/