Actual source code: ex1.c
1: static char help[] = "Basic equation for generator stability analysis.\n";
3: /*F
5: \begin{eqnarray}
6: \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - \frac{EV}{X} \sin(\theta) \\
7: \frac{d \theta}{dt} = \omega - \omega_s
8: \end{eqnarray}
10: F*/
12: /*
13: Include "petscts.h" so that we can use TS solvers. Note that this
14: file automatically includes:
15: petscsys.h - base PETSc routines petscvec.h - vectors
16: petscmat.h - matrices
17: petscis.h - index sets petscksp.h - Krylov subspace methods
18: petscviewer.h - viewers petscpc.h - preconditioners
19: petscksp.h - linear solvers
20: */
22: #include <petscts.h>
24: typedef struct {
25: PetscScalar H, omega_s, E, V, X;
26: PetscRandom rand;
27: } AppCtx;
29: /*
30: Defines the ODE passed to the ODE solver
31: */
32: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
33: {
34: PetscScalar *f, r;
35: const PetscScalar *u, *udot;
36: static PetscScalar R = .4;
38: PetscFunctionBegin;
39: PetscCall(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: PetscCall(VecGetArrayRead(U, &u));
45: PetscCall(VecGetArrayRead(Udot, &udot));
46: PetscCall(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: PetscCall(VecRestoreArrayRead(U, &u));
51: PetscCall(VecRestoreArrayRead(Udot, &udot));
52: PetscCall(VecRestoreArray(F, &f));
53: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
66: PetscCall(VecGetArrayRead(U, &u));
67: PetscCall(VecGetArrayRead(Udot, &udot));
68: J[0][0] = 2.0 * ctx->H * a / ctx->omega_s;
69: J[0][1] = -ctx->E * ctx->V * PetscCosScalar(u[1]) / ctx->X;
70: J[1][0] = -1.0;
71: J[1][1] = a;
72: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
73: PetscCall(VecRestoreArrayRead(U, &u));
74: PetscCall(VecRestoreArrayRead(Udot, &udot));
76: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
77: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
78: if (A != B) {
79: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
81: }
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: int main(int argc, char **argv)
86: {
87: TS ts; /* ODE integrator */
88: Vec U; /* solution will be stored here */
89: Mat A; /* Jacobian matrix */
90: PetscMPIInt size;
91: PetscInt n = 2;
92: AppCtx ctx;
93: PetscScalar *u;
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Initialize program
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscFunctionBeginUser;
99: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
100: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
101: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create necessary matrix and vectors
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
107: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
108: PetscCall(MatSetFromOptions(A));
109: PetscCall(MatSetUp(A));
111: PetscCall(MatCreateVecs(A, &U, NULL));
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Set runtime options
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Reaction options", "");
117: {
118: ctx.omega_s = 1.0;
119: PetscCall(PetscOptionsScalar("-omega_s", "", "", ctx.omega_s, &ctx.omega_s, NULL));
120: ctx.H = 1.0;
121: PetscCall(PetscOptionsScalar("-H", "", "", ctx.H, &ctx.H, NULL));
122: ctx.E = 1.0;
123: PetscCall(PetscOptionsScalar("-E", "", "", ctx.E, &ctx.E, NULL));
124: ctx.V = 1.0;
125: PetscCall(PetscOptionsScalar("-V", "", "", ctx.V, &ctx.V, NULL));
126: ctx.X = 1.0;
127: PetscCall(PetscOptionsScalar("-X", "", "", ctx.X, &ctx.X, NULL));
129: PetscCall(VecGetArray(U, &u));
130: u[0] = 1;
131: u[1] = .7;
132: PetscCall(VecRestoreArray(U, &u));
133: PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", U, NULL));
134: }
135: PetscOptionsEnd();
137: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ctx.rand));
138: PetscCall(PetscRandomSetFromOptions(ctx.rand));
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create timestepping solver context
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
144: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
145: PetscCall(TSSetType(ts, TSROSW));
146: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx));
147: PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Set initial conditions
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: PetscCall(TSSetSolution(ts, U));
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Set solver options
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: PetscCall(TSSetMaxTime(ts, 2000.0));
158: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
159: PetscCall(TSSetTimeStep(ts, .001));
160: PetscCall(TSSetFromOptions(ts));
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Solve nonlinear system
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: PetscCall(TSSolve(ts, U));
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Free work space. All PETSc objects should be destroyed when they are no longer needed.
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(MatDestroy(&A));
171: PetscCall(VecDestroy(&U));
172: PetscCall(TSDestroy(&ts));
173: PetscCall(PetscRandomDestroy(&ctx.rand));
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: /*TEST
180: build:
181: requires: !complex !single
183: test:
184: args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
186: test:
187: suffix: 2
188: args: -ts_max_steps 10
190: TEST*/