1 static char help[] = "Transistor amplifier.\n";
2
3 /*F
4 ` This example illustrates the implementation of an implicit DAE index-1 of form M y'=f(t,y) with singular mass matrix, where
5
6 [ -C1 C1 ]
7 [ C1 -C1 ]
8 M =[ -C2 ]; Ck = k * 1e-06
9 [ -C3 C3]
10 [ C3 -C3]
11
12 [ -(U(t) - y[0])/1000 ]
13 [ -6/R + y[1]/4500 + 0.01 * h(y[1]-y[2]) ]
14 f(t,y)= [ y[2]/R - h(y[1]-y[2]) ]
15 [ (y[3]-6)/9000 + 0.99 * h([y1]-y[2]) ]
16 [ y[4]/9000 ]
17
18 U(t) = 0.4 * Sin(200 Pi t); h[V] = 1e-06 * Exp(V/0.026 - 1) `
19
20 Useful options: -ts_monitor_lg_solution -ts_monitor_lg_timestep -lg_indicate_data_points 0
21 F*/
22
23 /*
24 Include "petscts.h" so that we can use TS solvers. Note that this
25 file automatically includes:
26 petscsys.h - base PETSc routines petscvec.h - vectors
27 petscmat.h - matrices
28 petscis.h - index sets petscksp.h - Krylov subspace methods
29 petscviewer.h - viewers petscpc.h - preconditioners
30 petscksp.h - linear solvers
31 */
32 #include <petscts.h>
33
34 FILE *gfilepointer_data, *gfilepointer_info;
35
36 /* Defines the source */
Ue(PetscScalar t,PetscScalar * U)37 PetscErrorCode Ue(PetscScalar t, PetscScalar *U)
38 {
39 PetscFunctionBeginUser;
40 *U = 0.4 * PetscSinReal(200 * PETSC_PI * t);
41 PetscFunctionReturn(PETSC_SUCCESS);
42 }
43
44 /*
45 Defines the DAE passed to the time solver
46 */
IFunctionImplicit(TS ts,PetscReal t,Vec Y,Vec Ydot,Vec F,PetscCtx ctx)47 static PetscErrorCode IFunctionImplicit(TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, PetscCtx ctx)
48 {
49 const PetscScalar *y, *ydot;
50 PetscScalar *f;
51
52 PetscFunctionBeginUser;
53 /* The next three lines allow us to access the entries of the vectors directly */
54 PetscCall(VecGetArrayRead(Y, &y));
55 PetscCall(VecGetArrayRead(Ydot, &ydot));
56 PetscCall(VecGetArrayWrite(F, &f));
57
58 f[0] = ydot[0] / 1.e6 - ydot[1] / 1.e6 - PetscSinReal(200 * PETSC_PI * t) / 2500. + y[0] / 1000.;
59 f[1] = -ydot[0] / 1.e6 + ydot[1] / 1.e6 - 0.0006666766666666667 + PetscExpReal((500 * (y[1] - y[2])) / 13.) / 1.e8 + y[1] / 4500.;
60 f[2] = ydot[2] / 500000. + 1.e-6 - PetscExpReal((500 * (y[1] - y[2])) / 13.) / 1.e6 + y[2] / 9000.;
61 f[3] = (3 * ydot[3]) / 1.e6 - (3 * ydot[4]) / 1.e6 - 0.0006676566666666666 + (99 * PetscExpReal((500 * (y[1] - y[2])) / 13.)) / 1.e8 + y[3] / 9000.;
62 f[4] = (3 * ydot[4]) / 1.e6 - (3 * ydot[3]) / 1.e6 + y[4] / 9000.;
63
64 PetscCall(VecRestoreArrayRead(Y, &y));
65 PetscCall(VecRestoreArrayRead(Ydot, &ydot));
66 PetscCall(VecRestoreArrayWrite(F, &f));
67 PetscFunctionReturn(PETSC_SUCCESS);
68 }
69
70 /*
71 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
72 */
IJacobianImplicit(TS ts,PetscReal t,Vec Y,Vec Ydot,PetscReal a,Mat A,Mat B,PetscCtx ctx)73 static PetscErrorCode IJacobianImplicit(TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
74 {
75 PetscInt rowcol[] = {0, 1, 2, 3, 4};
76 const PetscScalar *y, *ydot;
77 PetscScalar J[5][5];
78
79 PetscFunctionBeginUser;
80 PetscCall(VecGetArrayRead(Y, &y));
81 PetscCall(VecGetArrayRead(Ydot, &ydot));
82
83 PetscCall(PetscMemzero(J, sizeof(J)));
84
85 J[0][0] = a / 1.e6 + 0.001;
86 J[0][1] = -a / 1.e6;
87 J[1][0] = -a / 1.e6;
88 J[1][1] = a / 1.e6 + 0.00022222222222222223 + PetscExpReal((500 * (y[1] - y[2])) / 13.) / 2.6e6;
89 J[1][2] = -PetscExpReal((500 * (y[1] - y[2])) / 13.) / 2.6e6;
90 J[2][1] = -PetscExpReal((500 * (y[1] - y[2])) / 13.) / 26000.;
91 J[2][2] = a / 500000 + 0.00011111111111111112 + PetscExpReal((500 * (y[1] - y[2])) / 13.) / 26000.;
92 J[3][1] = (99 * PetscExpReal((500 * (y[1] - y[2])) / 13.)) / 2.6e6;
93 J[3][2] = (-99 * PetscExpReal((500 * (y[1] - y[2])) / 13.)) / 2.6e6;
94 J[3][3] = (3 * a) / 1.e6 + 0.00011111111111111112;
95 J[3][4] = -(3 * a) / 1.e6;
96 J[4][3] = -(3 * a) / 1.e6;
97 J[4][4] = (3 * a) / 1.e6 + 0.00011111111111111112;
98
99 PetscCall(MatSetValues(B, 5, rowcol, 5, rowcol, &J[0][0], INSERT_VALUES));
100
101 PetscCall(VecRestoreArrayRead(Y, &y));
102 PetscCall(VecRestoreArrayRead(Ydot, &ydot));
103
104 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
105 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
106 if (A != B) {
107 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
108 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
109 }
110 PetscFunctionReturn(PETSC_SUCCESS);
111 }
112
main(int argc,char ** argv)113 int main(int argc, char **argv)
114 {
115 TS ts; /* ODE integrator */
116 Vec Y; /* solution will be stored here */
117 Mat A; /* Jacobian matrix */
118 PetscMPIInt size;
119 PetscInt n = 5;
120 PetscScalar *y;
121
122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123 Initialize program
124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125 PetscFunctionBeginUser;
126 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
127 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
128 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
129
130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 Create necessary matrix and vectors
132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
134 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
135 PetscCall(MatSetFromOptions(A));
136 PetscCall(MatSetUp(A));
137
138 PetscCall(MatCreateVecs(A, &Y, NULL));
139
140 PetscCall(VecGetArray(Y, &y));
141 y[0] = 0.0;
142 y[1] = 3.0;
143 y[2] = y[1];
144 y[3] = 6.0;
145 y[4] = 0.0;
146 PetscCall(VecRestoreArray(Y, &y));
147
148 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149 Create timestepping solver context
150 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
152 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
153 PetscCall(TSSetType(ts, TSARKIMEX));
154 /* Must use ARKIMEX with fully implicit stages since mass matrix is not the identity */
155 PetscCall(TSARKIMEXSetType(ts, TSARKIMEXPRSSP2));
156 PetscCall(TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1));
157 /*PetscCall(TSSetType(ts,TSROSW));*/
158 PetscCall(TSSetIFunction(ts, NULL, IFunctionImplicit, NULL));
159 PetscCall(TSSetIJacobian(ts, A, A, IJacobianImplicit, NULL));
160
161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162 Set initial conditions
163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164 PetscCall(TSSetSolution(ts, Y));
165
166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167 Set solver options
168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169 PetscCall(TSSetMaxTime(ts, 0.15));
170 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
171 PetscCall(TSSetTimeStep(ts, .001));
172 PetscCall(TSSetFromOptions(ts));
173
174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175 Do time stepping
176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177 PetscCall(TSSolve(ts, Y));
178
179 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180 Free work space. All PETSc objects should be destroyed when they are no longer needed.
181 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182 PetscCall(MatDestroy(&A));
183 PetscCall(VecDestroy(&Y));
184 PetscCall(TSDestroy(&ts));
185 PetscCall(PetscFinalize());
186 return 0;
187 }
188
189 /*TEST
190 build:
191 requires: !single !complex
192 test:
193 args: -ts_monitor
194
195 TEST*/
196