Actual source code: ex19.c
2: static char help[] = "Solves the van der Pol DAE.\n\
3: Input parameters include:\n";
5: /* ------------------------------------------------------------------------
7: This program solves the van der Pol DAE
8: y' = -z = f(y,z) (1)
9: 0 = y-(z^3/3 - z) = g(y,z)
10: on the domain 0 <= x <= 1, with the boundary conditions
11: y(0) = -2, y'(0) = -2.355301397608119909925287735864250951918
12: This is a nonlinear equation.
14: Notes:
15: This code demonstrates the TS solver interface with the Van der Pol DAE,
16: namely it is the case when there is no RHS (meaning the RHS == 0), and the
17: equations are converted to two variants of linear problems, u_t = f(u,t),
18: namely turning (1) into a vector equation in terms of u,
20: [ y' + z ] = [ 0 ]
21: [ (z^3/3 - z) - y ] [ 0 ]
23: which then we can write as a vector equation
25: [ u_1' + u_2 ] = [ 0 ] (2)
26: [ (u_2^3/3 - u_2) - u_1 ] [ 0 ]
28: which is now in the desired form of u_t = f(u,t). As this is a DAE, and
29: there is no u_2', there is no need for a split,
31: so
33: [ F(u',u,t) ] = [ u_1' ] + [ u_2 ]
34: [ 0 ] [ (u_2^3/3 - u_2) - u_1 ]
36: Using the definition of the Jacobian of F (from the PETSc user manual),
37: in the equation F(u',u,t) = G(u,t),
39: dF dF
40: J(F) = a * -- - --
41: du' du
43: where d is the partial derivative. In this example,
45: dF [ 1 ; 0 ]
46: -- = [ ]
47: du' [ 0 ; 0 ]
49: dF [ 0 ; 1 ]
50: -- = [ ]
51: du [ -1 ; 1 - u_2^2 ]
53: Hence,
55: [ a ; -1 ]
56: J(F) = [ ]
57: [ 1 ; u_2^2 - 1 ]
59: ------------------------------------------------------------------------- */
61: #include <petscts.h>
63: typedef struct _n_User *User;
64: struct _n_User {
65: PetscReal next_output;
66: };
68: /*
69: User-defined routines
70: */
72: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
73: {
74: PetscScalar *f;
75: const PetscScalar *x,*xdot;
78: VecGetArrayRead(X,&x);
79: VecGetArrayRead(Xdot,&xdot);
80: VecGetArray(F,&f);
81: f[0] = xdot[0] + x[1];
82: f[1] = (x[1]*x[1]*x[1]/3.0 - x[1])-x[0];
83: VecRestoreArrayRead(X,&x);
84: VecRestoreArrayRead(Xdot,&xdot);
85: VecRestoreArray(F,&f);
86: return 0;
87: }
89: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
90: {
91: PetscInt rowcol[] = {0,1};
92: PetscScalar J[2][2];
93: const PetscScalar *x;
96: VecGetArrayRead(X,&x);
97: J[0][0] = a; J[0][1] = -1.;
98: J[1][0] = 1.; J[1][1] = -1. + x[1]*x[1];
99: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
100: VecRestoreArrayRead(X,&x);
102: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
103: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
104: if (A != B) {
105: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
107: }
108: return 0;
109: }
111: static PetscErrorCode RegisterMyARK2(void)
112: {
114: {
115: const PetscReal
116: A[3][3] = {{0,0,0},
117: {0.41421356237309504880,0,0},
118: {0.75,0.25,0}},
119: At[3][3] = {{0,0,0},
120: {0.12132034355964257320,0.29289321881345247560,0},
121: {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
122: *bembedt = NULL,*bembed = NULL;
123: TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);
124: }
125: return 0;
126: }
128: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
129: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
130: {
131: const PetscScalar *x;
132: PetscReal tfinal, dt;
133: User user = (User)ctx;
134: Vec interpolatedX;
137: TSGetTimeStep(ts,&dt);
138: TSGetMaxTime(ts,&tfinal);
140: while (user->next_output <= t && user->next_output <= tfinal) {
141: VecDuplicate(X,&interpolatedX);
142: TSInterpolate(ts,user->next_output,interpolatedX);
143: VecGetArrayRead(interpolatedX,&x);
144: PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %3D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",(double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),(double)PetscRealPart(x[1]));
145: VecRestoreArrayRead(interpolatedX,&x);
146: VecDestroy(&interpolatedX);
147: user->next_output += PetscRealConstant(0.1);
148: }
149: return 0;
150: }
152: int main(int argc,char **argv)
153: {
154: TS ts; /* nonlinear solver */
155: Vec x; /* solution, residual vectors */
156: Mat A; /* Jacobian matrix */
157: PetscInt steps;
158: PetscReal ftime = 0.5;
159: PetscBool monitor = PETSC_FALSE;
160: PetscScalar *x_ptr;
161: PetscMPIInt size;
162: struct _n_User user;
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Initialize program
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: PetscInitialize(&argc,&argv,NULL,help);
168: MPI_Comm_size(PETSC_COMM_WORLD,&size);
171: RegisterMyARK2();
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Set runtime options
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: user.next_output = 0.0;
178: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
180: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: Create necessary matrix and vectors, solve same ODE on every process
182: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: MatCreate(PETSC_COMM_WORLD,&A);
184: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
185: MatSetFromOptions(A);
186: MatSetUp(A);
187: MatCreateVecs(A,&x,NULL);
189: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: Create timestepping solver context
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192: TSCreate(PETSC_COMM_WORLD,&ts);
193: TSSetType(ts,TSBEULER);
194: TSSetIFunction(ts,NULL,IFunction,&user);
195: TSSetIJacobian(ts,A,A,IJacobian,&user);
196: TSSetMaxTime(ts,ftime);
197: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
198: if (monitor) {
199: TSMonitorSet(ts,Monitor,&user,NULL);
200: }
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Set initial conditions
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: VecGetArray(x,&x_ptr);
206: x_ptr[0] = -2; x_ptr[1] = -2.355301397608119909925287735864250951918;
207: VecRestoreArray(x,&x_ptr);
208: TSSetTimeStep(ts,.001);
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: Set runtime options
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213: TSSetFromOptions(ts);
215: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: Solve nonlinear system
217: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218: TSSolve(ts,x);
219: TSGetSolveTime(ts,&ftime);
220: TSGetStepNumber(ts,&steps);
221: PetscPrintf(PETSC_COMM_WORLD,"steps %3D, ftime %g\n",steps,(double)ftime);
222: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
224: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225: Free work space. All PETSc objects should be destroyed when they
226: are no longer needed.
227: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
228: MatDestroy(&A);
229: VecDestroy(&x);
230: TSDestroy(&ts);
232: PetscFinalize();
233: return 0;
234: }
236: /*TEST
238: test:
239: requires: !single
240: suffix: a
241: args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp
242: output_file: output/ex19_pi42.out
244: test:
245: requires: !single
246: suffix: b
247: args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_filter PI42
248: output_file: output/ex19_pi42.out
250: test:
251: requires: !single
252: suffix: c
253: args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_pid 0.4,0.2
254: output_file: output/ex19_pi42.out
256: test:
257: requires: !single
258: suffix: bdf_reject
259: args: -ts_type bdf -ts_dt 0.5 -ts_max_steps 1 -ts_max_reject {{0 1 2}separate_output} -ts_error_if_step_fails false -ts_adapt_monitor
261: TEST*/