Actual source code: ex22.c
1: static const char help[] = "Time-dependent advection-reaction PDE in 1d. Demonstrates IMEX methods.\n";
2: /*
3: u_t + a1*u_x = -k1*u + k2*v + s1
4: v_t + a2*v_x = k1*u - k2*v + s2
5: 0 < x < 1;
6: a1 = 1, k1 = 10^6, s1 = 0,
7: a2 = 0, k2 = 2*k1, s2 = 1
9: Initial conditions:
10: u(x,0) = 1 * s2*x
11: v(x,0) = k0/k1*u(x,0) + s1/k1
13: Upstream boundary conditions:
14: u(0,t) = 1-sin(12*t)^4
15: */
17: #include <petscdmda.h>
18: #include <petscts.h>
20: typedef PetscScalar Field[2];
22: typedef struct _User *User;
23: struct _User {
24: PetscReal a[2]; /* Advection speeds */
25: PetscReal k[2]; /* Reaction coefficients */
26: PetscReal s[2]; /* Source terms */
27: };
29: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
30: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
31: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
32: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
36: int main(int argc,char **argv)
37: {
38: TS ts; /* nonlinear solver */
39: Vec X; /* solution, residual vectors */
40: Mat J; /* Jacobian matrix */
41: PetscInt steps,maxsteps,mx;
43: DM da;
44: PetscReal ftime,dt;
45: struct _User user; /* user-defined work context */
47: PetscInitialize(&argc,&argv,(char *)0,help);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create distributed array (DMDA) to manage parallel grid and vectors
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-11,2,2,PETSC_NULL,&da);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Extract global vectors from DMDA;
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: DMCreateGlobalVector(da,&X);
59: /* Initialize user application context */
60: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Advection-reaction options",""); {
61: PetscOptionsReal("-a0","Advection rate 0","",user.a[0]=1,&user.a[0],PETSC_NULL);
62: PetscOptionsReal("-a1","Advection rate 1","",user.a[1]=0,&user.a[1],PETSC_NULL);
63: PetscOptionsReal("-k0","Reaction rate 0","",user.k[0]=1e6,&user.k[0],PETSC_NULL);
64: PetscOptionsReal("-k1","Reaction rate 1","",user.k[1]=2*user.k[0],&user.k[1],PETSC_NULL);
65: PetscOptionsReal("-s0","Source 0","",user.s[0]=0,&user.s[0],PETSC_NULL);
66: PetscOptionsReal("-s1","Source 1","",user.s[1]=1,&user.s[1],PETSC_NULL);
67: } PetscOptionsEnd();
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create timestepping solver context
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: TSCreate(PETSC_COMM_WORLD,&ts);
73: TSSetDM(ts,da);
74: TSSetType(ts,TSARKIMEX);
75: TSSetRHSFunction(ts,PETSC_NULL,FormRHSFunction,&user);
76: TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);
77: DMGetMatrix(da,MATAIJ,&J);
78: TSSetIJacobian(ts,J,J,FormIJacobian,&user);
80: ftime = 1.0;
81: maxsteps = 10000;
82: TSSetDuration(ts,maxsteps,ftime);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Set initial conditions
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: FormInitialSolution(ts,X,&user);
88: TSSetSolution(ts,X);
89: VecGetSize(X,&mx);
90: dt = .1 * PetscMax(user.a[0],user.a[1]) / mx; /* Advective CFL, I don't know why it needs so much safety factor. */
91: TSSetInitialTimeStep(ts,0.0,dt);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Set runtime options
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: TSSetFromOptions(ts);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Solve nonlinear system
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: TSSolve(ts,X,&ftime);
102: TSGetTimeStepNumber(ts,&steps);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Free work space.
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: MatDestroy(&J);
108: VecDestroy(&X);
109: TSDestroy(&ts);
110: DMDestroy(&da);
111: PetscFinalize();
112: return 0;
113: }
117: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
118: {
119: User user = (User)ptr;
120: DM da;
121: DMDALocalInfo info;
122: PetscInt i;
123: Field *x,*xdot,*f;
127: TSGetDM(ts,&da);
128: DMDAGetLocalInfo(da,&info);
130: /* Get pointers to vector data */
131: DMDAVecGetArray(da,X,&x);
132: DMDAVecGetArray(da,Xdot,&xdot);
133: DMDAVecGetArray(da,F,&f);
135: /* Compute function over the locally owned part of the grid */
136: for (i=info.xs; i<info.xs+info.xm; i++) {
137: f[i][0] = xdot[i][0] + user->k[0]*x[i][0] - user->k[1]*x[i][1] - user->s[0];
138: f[i][1] = xdot[i][1] - user->k[0]*x[i][0] + user->k[1]*x[i][1] - user->s[1];
139: }
141: /* Restore vectors */
142: DMDAVecRestoreArray(da,X,&x);
143: DMDAVecRestoreArray(da,Xdot,&xdot);
144: DMDAVecRestoreArray(da,F,&f);
145: return(0);
146: }
150: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
151: {
152: User user = (User)ptr;
153: DM da;
154: Vec Xloc;
155: DMDALocalInfo info;
156: PetscInt i,j;
157: PetscReal hx;
158: Field *x,*f;
162: TSGetDM(ts,&da);
163: DMDAGetLocalInfo(da,&info);
164: hx = 1.0/(PetscReal)info.mx;
166: /*
167: Scatter ghost points to local vector,using the 2-step process
168: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
169: By placing code between these two statements, computations can be
170: done while messages are in transition.
171: */
172: DMGetLocalVector(da,&Xloc);
173: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
174: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
176: /* Get pointers to vector data */
177: DMDAVecGetArray(da,Xloc,&x);
178: DMDAVecGetArray(da,F,&f);
180: /* Compute function over the locally owned part of the grid */
181: for (i=info.xs; i<info.xs+info.xm; i++) {
182: const PetscReal *a = user->a;
183: PetscReal u0t[2] = {1. - PetscPowScalar(sin(12*t),4.),0};
184: for (j=0; j<2; j++) {
185: if (i == 0) {
186: f[i][j] = a[j]/hx*(1./3*u0t[j] + 0.5*x[i][j] - x[i+1][j] + 1./6*x[i+2][j]);
187: } else if (i == 1) {
188: f[i][j] = a[j]/hx*(-1./12*u0t[j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);
189: } else if (i == info.mx-2) {
190: f[i][j] = a[j]/hx*(-1./6*x[i-2][j] + x[i-1][j] - 0.5*x[i][j] - 1./3*x[i+1][j]);
191: } else if (i == info.mx-1) {
192: f[i][j] = a[j]/hx*(-x[i][j] + x[i-1][j]);
193: } else {
194: f[i][j] = a[j]/hx*(-1./12*x[i-2][j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);
195: }
196: }
197: }
199: /* Restore vectors */
200: DMDAVecRestoreArray(da,Xloc,&x);
201: DMDAVecRestoreArray(da,F,&f);
202: DMRestoreLocalVector(da,&Xloc);
203: return(0);
204: }
206: /* --------------------------------------------------------------------- */
207: /*
208: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
209: */
212: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ptr)
213: {
214: User user = (User)ptr;
216: DMDALocalInfo info;
217: PetscInt i;
218: DM da;
219: Field *x,*xdot;
222: TSGetDM(ts,&da);
223: DMDAGetLocalInfo(da,&info);
225: /* Get pointers to vector data */
226: DMDAVecGetArray(da,X,&x);
227: DMDAVecGetArray(da,Xdot,&xdot);
229: /* Compute function over the locally owned part of the grid */
230: for (i=info.xs; i<info.xs+info.xm; i++) {
231: const PetscReal *k = user->k;
232: PetscScalar v[2][2] = {{a + k[0], -k[1]},
233: {-k[0], a+k[1]}};
234: MatSetValuesBlocked(*Jpre,1,&i,1,&i,&v[0][0],INSERT_VALUES);
235: }
237: /* Restore vectors */
238: DMDAVecRestoreArray(da,X,&x);
239: DMDAVecRestoreArray(da,Xdot,&xdot);
241: MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);
242: MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);
243: if (*J != *Jpre) {
244: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
245: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
246: }
247: return(0);
248: }
252: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
253: {
254: User user = (User)ctx;
255: DM da;
256: PetscInt i;
257: DMDALocalInfo info;
258: Field *x;
259: PetscReal hx;
263: TSGetDM(ts,&da);
264: DMDAGetLocalInfo(da,&info);
265: hx = 1.0/(PetscReal)info.mx;
267: /* Get pointers to vector data */
268: DMDAVecGetArray(da,X,&x);
270: /* Compute function over the locally owned part of the grid */
271: for (i=info.xs; i<info.xs+info.xm; i++) {
272: PetscReal r = (i+1)*hx,ik = user->k[1] != 0.0 ? 1.0/user->k[1] : 1.0;
273: x[i][0] = 1 + user->s[1]*r;
274: x[i][1] = user->k[0]*ik*x[i][0] + user->s[1]*ik;
275: }
276: DMDAVecRestoreArray(da,X,&x);
277: return(0);
278: }