Actual source code: ex1.c
petsc-3.13.4 2020-08-01
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: PetscErrorCode ierr;
36: PetscScalar *f,r;
37: const PetscScalar *u,*udot;
38: static PetscScalar R = .4;
41: PetscRandomGetValue(ctx->rand,&r);
42: if (r > .9) R = .5;
43: if (r < .1) R = .4;
44: R = .4;
45: /* The next three lines allow us to access the entries of the vectors directly */
46: VecGetArrayRead(U,&u);
47: VecGetArrayRead(Udot,&udot);
48: VecGetArray(F,&f);
49: f[0] = 2.0*ctx->H*udot[0]/ctx->omega_s + ctx->E*ctx->V*PetscSinScalar(u[1])/ctx->X - R;
50: f[1] = udot[1] - u[0] + ctx->omega_s;
52: VecRestoreArrayRead(U,&u);
53: VecRestoreArrayRead(Udot,&udot);
54: VecRestoreArray(F,&f);
55: return(0);
56: }
58: /*
59: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
60: */
61: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
62: {
63: PetscErrorCode ierr;
64: PetscInt rowcol[] = {0,1};
65: PetscScalar J[2][2];
66: const PetscScalar *u,*udot;
69: VecGetArrayRead(U,&u);
70: VecGetArrayRead(Udot,&udot);
71: J[0][0] = 2.0*ctx->H*a/ctx->omega_s; J[0][1] = -ctx->E*ctx->V*PetscCosScalar(u[1])/ctx->X;
72: J[1][0] = -1.0; J[1][1] = a;
73: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
74: VecRestoreArrayRead(U,&u);
75: VecRestoreArrayRead(Udot,&udot);
77: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
79: if (A != B) {
80: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
82: }
83: return(0);
84: }
86: int main(int argc,char **argv)
87: {
88: TS ts; /* ODE integrator */
89: Vec U; /* solution will be stored here */
90: Mat A; /* Jacobian matrix */
92: PetscMPIInt size;
93: PetscInt n = 2;
94: AppCtx ctx;
95: PetscScalar *u;
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Initialize program
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
101: MPI_Comm_size(PETSC_COMM_WORLD,&size);
102: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Create necessary matrix and vectors
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: MatCreate(PETSC_COMM_WORLD,&A);
108: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
109: MatSetFromOptions(A);
110: MatSetUp(A);
112: MatCreateVecs(A,&U,NULL);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Set runtime options
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Reaction options","");
118: {
119: ctx.omega_s = 1.0;
120: PetscOptionsScalar("-omega_s","","",ctx.omega_s,&ctx.omega_s,NULL);
121: ctx.H = 1.0;
122: PetscOptionsScalar("-H","","",ctx.H,&ctx.H,NULL);
123: ctx.E = 1.0;
124: PetscOptionsScalar("-E","","",ctx.E,&ctx.E,NULL);
125: ctx.V = 1.0;
126: PetscOptionsScalar("-V","","",ctx.V,&ctx.V,NULL);
127: ctx.X = 1.0;
128: PetscOptionsScalar("-X","","",ctx.X,&ctx.X,NULL);
130: VecGetArray(U,&u);
131: u[0] = 1;
132: u[1] = .7;
133: VecRestoreArray(U,&u);
134: PetscOptionsGetVec(NULL,NULL,"-initial",U,NULL);
135: }
136: PetscOptionsEnd();
138: PetscRandomCreate(PETSC_COMM_WORLD,&ctx.rand);
139: PetscRandomSetFromOptions(ctx.rand);
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Create timestepping solver context
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: TSCreate(PETSC_COMM_WORLD,&ts);
145: TSSetProblemType(ts,TS_NONLINEAR);
146: TSSetType(ts,TSROSW);
147: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);
148: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Set initial conditions
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: TSSetSolution(ts,U);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Set solver options
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: TSSetMaxTime(ts,2000.0);
159: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
160: TSSetTimeStep(ts,.001);
161: TSSetFromOptions(ts);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Solve nonlinear system
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: TSSolve(ts,U);
168: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: Free work space. All PETSc objects should be destroyed when they are no longer needed.
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: MatDestroy(&A);
172: VecDestroy(&U);
173: TSDestroy(&ts);
174: PetscRandomDestroy(&ctx.rand);
175: PetscFinalize();
176: return ierr;
177: }
180: /*TEST
182: build:
183: requires: !complex !single
185: test:
186: args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_max_steps 10
188: test:
189: suffix: 2
190: args: -ts_max_steps 10
192: TEST*/