Actual source code: ex19.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Solves the van der Pol DAE.\n\
3: Input parameters include:\n";
5: /*
6: Concepts: TS^time-dependent nonlinear problems
7: Concepts: TS^van der Pol DAE
8: Processors: 1
9: */
10: /* ------------------------------------------------------------------------
12: This program solves the van der Pol DAE
13: y' = -z = f(y,z) (1)
14: 0 = y-(z^3/3 - z) = g(y,z)
15: on the domain 0 <= x <= 1, with the boundary conditions
16: y(0) = -2, y'(0) = -2.355301397608119909925287735864250951918
17: This is a nonlinear equation.
19: Notes:
20: This code demonstrates the TS solver interface with the Van der Pol DAE,
21: namely it is the case when there is no RHS (meaning the RHS == 0), and the
22: equations are converted to two variants of linear problems, u_t = f(u,t),
23: namely turning (1) into a vector equation in terms of u,
25: [ y' + z ] = [ 0 ]
26: [ (z^3/3 - z) - y ] [ 0 ]
28: which then we can write as a vector equation
30: [ u_1' + u_2 ] = [ 0 ] (2)
31: [ (u_2^3/3 - u_2) - u_1 ] [ 0 ]
33: which is now in the desired form of u_t = f(u,t). As this is a DAE, and
34: there is no u_2', there is no need for a split,
36: so
38: [ G(u',u,t) ] = [ u_1' ] + [ u_2 ]
39: [ 0 ] [ (u_2^3/3 - u_2) - u_1 ]
41: Using the definition of the Jacobian of G (from the PETSc user manual),
42: in the equation G(u',u,t) = F(u,t),
44: dG dG
45: J(G) = a * -- - --
46: du' du
48: where d is the partial derivative. In this example,
50: dG [ 1 ; 0 ]
51: -- = [ ]
52: du' [ 0 ; 0 ]
54: dG [ 0 ; 1 ]
55: -- = [ ]
56: du [ -1 ; 1 - u_2^2 ]
58: Hence,
60: [ a ; -1 ]
61: J(G) = [ ]
62: [ 1 ; u_2^2 - 1 ]
64: ------------------------------------------------------------------------- */
66: #include <petscts.h>
68: typedef struct _n_User *User;
69: struct _n_User {
70: PetscReal next_output;
71: };
73: /*
74: * User-defined routines
75: */
77: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
78: {
79: PetscErrorCode ierr;
80: PetscScalar *f;
81: const PetscScalar *x,*xdot;
84: VecGetArrayRead(X,&x);
85: VecGetArrayRead(Xdot,&xdot);
86: VecGetArray(F,&f);
87: f[0] = xdot[0] + x[1];
88: f[1] = (x[1]*x[1]*x[1]/3.0 - x[1])-x[0];
89: VecRestoreArrayRead(X,&x);
90: VecRestoreArrayRead(Xdot,&xdot);
91: VecRestoreArray(F,&f);
92: return(0);
93: }
95: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
96: {
97: PetscErrorCode ierr;
98: PetscInt rowcol[] = {0,1};
99: PetscScalar J[2][2];
100: const PetscScalar *x;
103: VecGetArrayRead(X,&x);
104: J[0][0] = a; J[0][1] = -1.;
105: J[1][0] = 1.; J[1][1] = -1. + x[1]*x[1];
106: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
107: VecRestoreArrayRead(X,&x);
109: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
111: if (A != B) {
112: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
114: }
115: return(0);
116: }
118: static PetscErrorCode RegisterMyARK2(void)
119: {
123: {
124: const PetscReal
125: A[3][3] = {{0,0,0},
126: {0.41421356237309504880,0,0},
127: {0.75,0.25,0}},
128: At[3][3] = {{0,0,0},
129: {0.12132034355964257320,0.29289321881345247560,0},
130: {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
131: *bembedt = NULL,*bembed = NULL;
132: TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);
133: }
134: return(0);
135: }
137: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
138: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
139: {
140: PetscErrorCode ierr;
141: const PetscScalar *x;
142: PetscReal tfinal, dt;
143: User user = (User)ctx;
144: Vec interpolatedX;
147: TSGetTimeStep(ts,&dt);
148: TSGetMaxTime(ts,&tfinal);
150: while (user->next_output <= t && user->next_output <= tfinal) {
151: VecDuplicate(X,&interpolatedX);
152: TSInterpolate(ts,user->next_output,interpolatedX);
153: VecGetArrayRead(interpolatedX,&x);
154: 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]));
155: VecRestoreArrayRead(interpolatedX,&x);
156: VecDestroy(&interpolatedX);
157: user->next_output += PetscRealConstant(0.1);
158: }
159: return(0);
160: }
162: int main(int argc,char **argv)
163: {
164: TS ts; /* nonlinear solver */
165: Vec x; /* solution, residual vectors */
166: Mat A; /* Jacobian matrix */
167: PetscInt steps;
168: PetscReal ftime = 0.5;
169: PetscBool monitor = PETSC_FALSE;
170: PetscScalar *x_ptr;
171: PetscMPIInt size;
172: struct _n_User user;
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Initialize program
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178: PetscInitialize(&argc,&argv,NULL,help);
179: MPI_Comm_size(PETSC_COMM_WORLD,&size);
180: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
182: RegisterMyARK2();
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Set runtime options
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188: user.next_output = 0.0;
189: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
191: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: Create necessary matrix and vectors, solve same ODE on every process
193: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194: MatCreate(PETSC_COMM_WORLD,&A);
195: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
196: MatSetFromOptions(A);
197: MatSetUp(A);
198: MatCreateVecs(A,&x,NULL);
200: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: Create timestepping solver context
202: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203: TSCreate(PETSC_COMM_WORLD,&ts);
204: TSSetType(ts,TSBEULER);
205: TSSetIFunction(ts,NULL,IFunction,&user);
206: TSSetIJacobian(ts,A,A,IJacobian,&user);
207: TSSetMaxTime(ts,ftime);
208: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
209: if (monitor) {
210: TSMonitorSet(ts,Monitor,&user,NULL);
211: }
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Set initial conditions
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: VecGetArray(x,&x_ptr);
217: x_ptr[0] = -2; x_ptr[1] = -2.355301397608119909925287735864250951918;
218: VecRestoreArray(x,&x_ptr);
219: TSSetTimeStep(ts,.001);
221: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: Set runtime options
223: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: TSSetFromOptions(ts);
226: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: Solve nonlinear system
228: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229: TSSolve(ts,x);
230: TSGetSolveTime(ts,&ftime);
231: TSGetStepNumber(ts,&steps);
232: PetscPrintf(PETSC_COMM_WORLD,"steps %3D, ftime %g\n",steps,(double)ftime);
233: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: Free work space. All PETSc objects should be destroyed when they
237: are no longer needed.
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: MatDestroy(&A);
240: VecDestroy(&x);
241: TSDestroy(&ts);
243: PetscFinalize();
244: return ierr;
245: }
247: /*TEST
249: test:
250: requires: !single
251: suffix: a
252: args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp
253: output_file: output/ex19_pi42.out
255: test:
256: requires: !single
257: suffix: b
258: args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_filter PI42
259: output_file: output/ex19_pi42.out
261: test:
262: requires: !single
263: suffix: c
264: 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
265: output_file: output/ex19_pi42.out
267: TEST*/