Actual source code: ex1.c
2: static char help[] = "Basic equation for generator stability analysis.\n";
4: /*F
6: \begin{eqnarray}
7: \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - \frac{EV}{X} \sin(\theta) \\
8: \frac{d \theta}{dt} = \omega - \omega_s
9: \end{eqnarray}
11: F*/
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;
68: J[0][1] = -ctx->E * ctx->V * PetscCosScalar(u[1]) / ctx->X;
69: J[1][0] = -1.0;
70: J[1][1] = a;
71: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
72: VecRestoreArrayRead(U, &u);
73: VecRestoreArrayRead(Udot, &udot);
75: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
77: if (A != B) {
78: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
80: }
81: return 0;
82: }
84: int main(int argc, char **argv)
85: {
86: TS ts; /* ODE integrator */
87: Vec U; /* solution will be stored here */
88: Mat A; /* Jacobian matrix */
89: PetscMPIInt size;
90: PetscInt n = 2;
91: AppCtx ctx;
92: PetscScalar *u;
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Initialize program
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscInitialize(&argc, &argv, (char *)0, help);
99: MPI_Comm_size(PETSC_COMM_WORLD, &size);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create necessary matrix and vectors
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: MatCreate(PETSC_COMM_WORLD, &A);
106: MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
107: MatSetFromOptions(A);
108: MatSetUp(A);
110: MatCreateVecs(A, &U, NULL);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Set runtime options
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Reaction options", "");
116: {
117: ctx.omega_s = 1.0;
118: PetscOptionsScalar("-omega_s", "", "", ctx.omega_s, &ctx.omega_s, NULL);
119: ctx.H = 1.0;
120: PetscOptionsScalar("-H", "", "", ctx.H, &ctx.H, NULL);
121: ctx.E = 1.0;
122: PetscOptionsScalar("-E", "", "", ctx.E, &ctx.E, NULL);
123: ctx.V = 1.0;
124: PetscOptionsScalar("-V", "", "", ctx.V, &ctx.V, NULL);
125: ctx.X = 1.0;
126: PetscOptionsScalar("-X", "", "", ctx.X, &ctx.X, NULL);
128: VecGetArray(U, &u);
129: u[0] = 1;
130: u[1] = .7;
131: VecRestoreArray(U, &u);
132: PetscOptionsGetVec(NULL, NULL, "-initial", U, NULL);
133: }
134: PetscOptionsEnd();
136: PetscRandomCreate(PETSC_COMM_WORLD, &ctx.rand);
137: PetscRandomSetFromOptions(ctx.rand);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Create timestepping solver context
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: TSCreate(PETSC_COMM_WORLD, &ts);
143: TSSetProblemType(ts, TS_NONLINEAR);
144: TSSetType(ts, TSROSW);
145: TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx);
146: TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Set initial conditions
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: TSSetSolution(ts, U);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Set solver options
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: TSSetMaxTime(ts, 2000.0);
157: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
158: TSSetTimeStep(ts, .001);
159: TSSetFromOptions(ts);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Solve nonlinear system
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: TSSolve(ts, U);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Free work space. All PETSc objects should be destroyed when they are no longer needed.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: MatDestroy(&A);
170: VecDestroy(&U);
171: TSDestroy(&ts);
172: PetscRandomDestroy(&ctx.rand);
173: PetscFinalize();
174: return 0;
175: }
177: /*TEST
179: build:
180: requires: !complex !single
182: test:
183: args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
185: test:
186: suffix: 2
187: args: -ts_max_steps 10
189: TEST*/