Actual source code: ct_vdp_imex.c
1: /*
2: * ex_vdp.c
3: *
4: * Created on: Jun 1, 2012
5: * Author: Hong Zhang
6: */
7: static char help[] = "Solves the van der Pol equation. \n Input parameters include:\n";
9: /*
10: * This program solves the van der Pol equation
11: * y' = z (1)
12: * z' = (((1-y^2)*z-y)/eps (2)
13: * on the domain 0<=x<=0.5, with the initial conditions
14: * y(0) = 2,
15: * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3
16: * IMEX schemes are applied to the splitted equation
17: * [y'] = [z] + [0 ]
18: * [z'] [0] [(((1-y^2)*z-y)/eps]
19: *
20: * F(x)= [z]
21: * [0]
22: *
23: * G(x)= [y'] - [0 ]
24: * [z'] [(((1-y^2)*z-y)/eps]
25: *
26: * JG(x) = G_x + a G_xdot
27: */
29: #include <petscdmda.h>
30: #include <petscts.h>
32: typedef struct _User *User;
33: struct _User {
34: PetscReal mu; /*stiffness control coefficient: epsilon*/
35: };
37: static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
38: static PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
39: static PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
41: int main(int argc, char **argv)
42: {
43: TS ts;
44: Vec x; /*solution vector*/
45: Mat A; /*Jacobian*/
46: PetscInt steps,mx,eimex_rowcol[2],two;
47: PetscErrorCode ierr;
48: PetscScalar *x_ptr;
49: PetscReal ftime,dt,norm;
50: Vec ref;
51: struct _User user; /* user-defined work context */
52: PetscViewer viewer;
54: PetscInitialize(&argc,&argv,NULL,help);
55: /* Initialize user application context */
56: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options","");
57: user.mu = 1e0;
58: PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL);
59: PetscOptionsEnd();
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Set runtime options
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
66: */
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create necessary matrix and vectors, solve same ODE on every process
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: MatCreate(PETSC_COMM_WORLD,&A);
72: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
73: MatSetFromOptions(A);
74: MatSetUp(A);
75: MatCreateVecs(A,&x,NULL);
77: MatCreateVecs(A,&ref,NULL);
78: VecGetArray(ref,&x_ptr);
79: /*
80: * [0,1], mu=10^-3
81: */
82: /*
83: x_ptr[0] = -1.8881254106283;
84: x_ptr[1] = 0.7359074233370;*/
86: /*
87: * [0,0.5],mu=10^-3
88: */
89: /*
90: x_ptr[0] = 1.596980778659137;
91: x_ptr[1] = -1.029103015879544;
92: */
93: /*
94: * [0,0.5],mu=1
95: */
96: x_ptr[0] = 1.619084329683235;
97: x_ptr[1] = -0.803530465176385;
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create timestepping solver context
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: TSCreate(PETSC_COMM_WORLD,&ts);
103: TSSetType(ts,TSEIMEX);
104: TSSetRHSFunction(ts,NULL,RHSFunction,&user);
105: TSSetIFunction(ts,NULL,IFunction,&user);
106: TSSetIJacobian(ts,A,A,IJacobian,&user);
108: dt = 0.00001;
109: ftime = 1.1;
110: TSSetTimeStep(ts,dt);
111: TSSetMaxTime(ts,ftime);
112: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Set initial conditions
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: VecGetArray(x,&x_ptr);
117: x_ptr[0] = 2.;
118: x_ptr[1] = -2./3. + 10./81.*(user.mu) - 292./2187.* (user.mu) * (user.mu)
119: -1814./19683.*(user.mu)*(user.mu)*(user.mu);
120: TSSetSolution(ts,x);
121: VecGetSize(x,&mx);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Set runtime options
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: TSSetFromOptions(ts);
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Solve nonlinear system
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: TSSolve(ts,x);
132: TSGetTime(ts,&ftime);
133: TSGetStepNumber(ts,&steps);
135: VecAXPY(x,-1.0,ref);
136: VecNorm(x,NORM_2,&norm);
137: TSGetTimeStep(ts,&dt);
139: eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2;
140: PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL);
141: PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1");
142: VecGetArray(x,&x_ptr);
143: PetscPrintf(PETSC_COMM_WORLD,"(%D,%D) %10.8f %18.15f %18.15f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm,(double)PetscRealPart(x_ptr[0]),(double)PetscRealPart(x_ptr[1]));
144: VecRestoreArray(x,&x_ptr);
146: /* Write line in convergence log */
147: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
148: PetscViewerSetType(viewer,PETSCVIEWERASCII);
149: PetscViewerFileSetMode(viewer,FILE_MODE_APPEND);
150: PetscViewerFileSetName(viewer,"eimex_nonstiff_vdp.txt");
151: PetscViewerASCIIPrintf(viewer,"%D %D %10.8f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm);
152: PetscViewerDestroy(&viewer);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Free work space.
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: MatDestroy(&A);
158: VecDestroy(&x);
159: VecDestroy(&ref);
160: TSDestroy(&ts);
161: PetscFinalize();
162: return 0;
163: }
165: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
166: {
167: PetscScalar *f;
168: const PetscScalar *x;
170: VecGetArrayRead(X,&x);
171: VecGetArray(F,&f);
172: f[0] = x[1];
173: f[1] = 0.0;
174: VecRestoreArrayRead(X,&x);
175: VecRestoreArray(F,&f);
176: return 0;
177: }
179: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
180: {
181: User user = (User)ptr;
182: PetscScalar *f;
183: const PetscScalar *x,*xdot;
185: VecGetArrayRead(X,&x);
186: VecGetArrayRead(Xdot,&xdot);
187: VecGetArray(F,&f);
188: f[0] = xdot[0];
189: f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu;
190: VecRestoreArrayRead(X,&x);
191: VecRestoreArrayRead(Xdot,&xdot);
192: VecRestoreArray(F,&f);
193: return 0;
194: }
196: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr)
197: {
198: User user = (User)ptr;
199: PetscReal mu = user->mu;
200: PetscInt rowcol[] = {0,1};
201: PetscScalar J[2][2];
202: const PetscScalar *x;
204: VecGetArrayRead(X,&x);
205: J[0][0] = a;
206: J[0][1] = 0;
207: J[1][0] = (2.*x[0]*x[1]+1.)/mu;
208: J[1][1] = a - (1. - x[0]*x[0])/mu;
209: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
210: VecRestoreArrayRead(X,&x);
212: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
213: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
214: if (A != B) {
215: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
216: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
217: }
218: return 0;
219: }
221: /*TEST
223: test:
224: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution
225: requires: x
227: test:
228: suffix: adapt
229: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution
230: requires: x
232: test:
233: suffix: loop
234: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt {{0.005 0.001 0.0005}separate output} -ts_max_steps {{100 500 1000}separate output} -ts_eimex_row_col {{1,1 2,1 3,1 2,2 3,2 3,3}separate output}
235: requires: x
237: TEST*/