1
2 static char help[] ="Model Equations for Advection \n";
3
4 /*
5 Modified from ex3.c
6 Page 9, Section 1.2 Model Equations for Advection-Diffusion
7
8 u_t + a u_x = 0, 0<= x <= 1.0
9
10 The initial conditions used here different from the book.
11
12 Example:
13 ./ex6 -ts_monitor -ts_view_solution -ts_max_steps 100 -ts_monitor_solution draw -draw_pause .1
14 ./ex6 -ts_monitor -ts_max_steps 100 -ts_monitor_lg_error -draw_pause .1
15 */
16
17 #include <petscts.h>
18 #include <petscdm.h>
19 #include <petscdmda.h>
20
21 /*
22 User-defined application context - contains data needed by the
23 application-provided call-back routines.
24 */
25 typedef struct {
26 PetscReal a; /* advection strength */
27 } AppCtx;
28
29 /* User-defined routines */
30 extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*);
31 extern PetscErrorCode Solution(TS,PetscReal,Vec,AppCtx*);
32 extern PetscErrorCode IFunction_LaxFriedrichs(TS,PetscReal,Vec,Vec,Vec,void*);
33 extern PetscErrorCode IFunction_LaxWendroff(TS,PetscReal,Vec,Vec,Vec,void*);
34
main(int argc,char ** argv)35 int main(int argc,char **argv)
36 {
37 AppCtx appctx; /* user-defined application context */
38 TS ts; /* timestepping context */
39 Vec U; /* approximate solution vector */
40 PetscErrorCode ierr;
41 PetscReal dt;
42 DM da;
43 PetscInt M;
44 PetscMPIInt rank;
45 PetscBool useLaxWendroff = PETSC_TRUE;
46
47 /* Initialize program and set problem parameters */
48 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
49 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
50
51 appctx.a = -1.0;
52 ierr = PetscOptionsGetReal(NULL,NULL,"-a",&appctx.a,NULL);CHKERRQ(ierr);
53
54 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);CHKERRQ(ierr);
55 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
56 ierr = DMSetUp(da);CHKERRQ(ierr);
57
58 /* Create vector data structures for approximate and exact solutions */
59 ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr);
60
61 /* Create timestepping solver context */
62 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
63 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
64
65 /* Function evaluation */
66 ierr = PetscOptionsGetBool(NULL,NULL,"-useLaxWendroff",&useLaxWendroff,NULL);CHKERRQ(ierr);
67 if (useLaxWendroff) {
68 if (!rank) {
69 ierr = PetscPrintf(PETSC_COMM_SELF,"... Use Lax-Wendroff finite volume\n");CHKERRQ(ierr);
70 }
71 ierr = TSSetIFunction(ts,NULL,IFunction_LaxWendroff,&appctx);CHKERRQ(ierr);
72 } else {
73 if (!rank) {
74 ierr = PetscPrintf(PETSC_COMM_SELF,"... Use Lax-LaxFriedrichs finite difference\n");CHKERRQ(ierr);
75 }
76 ierr = TSSetIFunction(ts,NULL,IFunction_LaxFriedrichs,&appctx);CHKERRQ(ierr);
77 }
78
79 /* Customize timestepping solver */
80 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
81 dt = 1.0/(PetscAbsReal(appctx.a)*M);
82 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
83 ierr = TSSetMaxSteps(ts,100);CHKERRQ(ierr);
84 ierr = TSSetMaxTime(ts,100.0);CHKERRQ(ierr);
85 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
86 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
87 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
88
89 /* Evaluate initial conditions */
90 ierr = InitialConditions(ts,U,&appctx);CHKERRQ(ierr);
91
92 /* For testing accuracy of TS with already known solution, e.g., '-ts_monitor_lg_error' */
93 ierr = TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);CHKERRQ(ierr);
94
95 /* Run the timestepping solver */
96 ierr = TSSolve(ts,U);CHKERRQ(ierr);
97
98 /* Free work space */
99 ierr = TSDestroy(&ts);CHKERRQ(ierr);
100 ierr = VecDestroy(&U);CHKERRQ(ierr);
101 ierr = DMDestroy(&da);CHKERRQ(ierr);
102
103 ierr = PetscFinalize();
104 return ierr;
105 }
106 /* --------------------------------------------------------------------- */
107 /*
108 InitialConditions - Computes the solution at the initial time.
109
110 Input Parameter:
111 u - uninitialized solution vector (global)
112 appctx - user-defined application context
113
114 Output Parameter:
115 u - vector with solution at initial time (global)
116 */
InitialConditions(TS ts,Vec U,AppCtx * appctx)117 PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
118 {
119 PetscScalar *u;
120 PetscErrorCode ierr;
121 PetscInt i,mstart,mend,um,M;
122 DM da;
123 PetscReal h;
124
125 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
126 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
127 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
128 h = 1.0/M;
129 mend = mstart + um;
130 /*
131 Get a pointer to vector data.
132 - For default PETSc vectors, VecGetArray() returns a pointer to
133 the data array. Otherwise, the routine is implementation dependent.
134 - You MUST call VecRestoreArray() when you no longer need access to
135 the array.
136 - Note that the Fortran interface to VecGetArray() differs from the
137 C version. See the users manual for details.
138 */
139 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
140
141 /*
142 We initialize the solution array by simply writing the solution
143 directly into the array locations. Alternatively, we could use
144 VecSetValues() or VecSetValuesLocal().
145 */
146 for (i=mstart; i<mend; i++) u[i] = PetscSinReal(PETSC_PI*i*6.*h) + 3.*PetscSinReal(PETSC_PI*i*2.*h);
147
148 /* Restore vector */
149 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
150 return 0;
151 }
152 /* --------------------------------------------------------------------- */
153 /*
154 Solution - Computes the exact solution at a given time
155
156 Input Parameters:
157 t - current time
158 solution - vector in which exact solution will be computed
159 appctx - user-defined application context
160
161 Output Parameter:
162 solution - vector with the newly computed exact solution
163 u(x,t) = sin(6*PI*(x - a*t)) + 3 * sin(2*PI*(x - a*t))
164 */
Solution(TS ts,PetscReal t,Vec U,AppCtx * appctx)165 PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
166 {
167 PetscScalar *u;
168 PetscReal a=appctx->a,h,PI6,PI2;
169 PetscErrorCode ierr;
170 PetscInt i,mstart,mend,um,M;
171 DM da;
172
173 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
174 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
175 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
176 h = 1.0/M;
177 mend = mstart + um;
178
179 /* Get a pointer to vector data. */
180 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
181
182 /* u[i] = sin(6*PI*(x[i] - a*t)) + 3 * sin(2*PI*(x[i] - a*t)) */
183 PI6 = PETSC_PI*6.;
184 PI2 = PETSC_PI*2.;
185 for (i=mstart; i<mend; i++) {
186 u[i] = PetscSinReal(PI6*(i*h - a*t)) + 3.*PetscSinReal(PI2*(i*h - a*t));
187 }
188
189 /* Restore vector */
190 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
191 return 0;
192 }
193
194 /* --------------------------------------------------------------------- */
195 /*
196 Use Lax-Friedrichs method to evaluate F(u,t) = du/dt + a * du/dx
197
198 See https://en.wikipedia.org/wiki/Lax%E2%80%93Friedrichs_method
199 */
IFunction_LaxFriedrichs(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void * ctx)200 PetscErrorCode IFunction_LaxFriedrichs(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void* ctx)
201 {
202 PetscErrorCode ierr;
203 AppCtx *appctx=(AppCtx*)ctx;
204 PetscInt mstart,mend,M,i,um;
205 DM da;
206 Vec Uold,localUold;
207 PetscScalar *uarray,*f,*uoldarray,h,uave,c;
208 PetscReal dt;
209
210 PetscFunctionBegin;
211 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
212 ierr = TSGetSolution(ts,&Uold);CHKERRQ(ierr);
213
214 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
215 ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
216 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
217 h = 1.0/M;
218 mend = mstart + um;
219 /* printf(" mstart %d, um %d\n",mstart,um); */
220
221 ierr = DMGetLocalVector(da,&localUold);CHKERRQ(ierr);
222 ierr = DMGlobalToLocalBegin(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);
223 ierr = DMGlobalToLocalEnd(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);
224
225 /* Get pointers to vector data */
226 ierr = DMDAVecGetArrayRead(da,U,&uarray);CHKERRQ(ierr);
227 ierr = DMDAVecGetArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);
228 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
229
230 /* advection */
231 c = appctx->a*dt/h; /* Courant-Friedrichs-Lewy number (CFL number) */
232
233 for (i=mstart; i<mend; i++) {
234 uave = 0.5*(uoldarray[i-1] + uoldarray[i+1]);
235 f[i] = uarray[i] - uave + c*0.5*(uoldarray[i+1] - uoldarray[i-1]);
236 }
237
238 /* Restore vectors */
239 ierr = DMDAVecRestoreArrayRead(da,U,&uarray);CHKERRQ(ierr);
240 ierr = DMDAVecRestoreArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);
241 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
242 ierr = DMRestoreLocalVector(da,&localUold);CHKERRQ(ierr);
243 PetscFunctionReturn(0);
244 }
245
246 /*
247 Use Lax-Wendroff method to evaluate F(u,t) = du/dt + a * du/dx
248 */
IFunction_LaxWendroff(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void * ctx)249 PetscErrorCode IFunction_LaxWendroff(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void* ctx)
250 {
251 PetscErrorCode ierr;
252 AppCtx *appctx=(AppCtx*)ctx;
253 PetscInt mstart,mend,M,i,um;
254 DM da;
255 Vec Uold,localUold;
256 PetscScalar *uarray,*f,*uoldarray,h,RFlux,LFlux,lambda;
257 PetscReal dt,a;
258
259 PetscFunctionBegin;
260 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
261 ierr = TSGetSolution(ts,&Uold);CHKERRQ(ierr);
262
263 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
264 ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
265 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
266 h = 1.0/M;
267 mend = mstart + um;
268 /* printf(" mstart %d, um %d\n",mstart,um); */
269
270 ierr = DMGetLocalVector(da,&localUold);CHKERRQ(ierr);
271 ierr = DMGlobalToLocalBegin(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);
272 ierr = DMGlobalToLocalEnd(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);
273
274 /* Get pointers to vector data */
275 ierr = DMDAVecGetArrayRead(da,U,&uarray);CHKERRQ(ierr);
276 ierr = DMDAVecGetArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);
277 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
278
279 /* advection -- finite volume (appctx->a < 0 -- can be relaxed?) */
280 lambda = dt/h;
281 a = appctx->a;
282
283 for (i=mstart; i<mend; i++) {
284 RFlux = 0.5 * a * (uoldarray[i+1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i+1] - uoldarray[i]);
285 LFlux = 0.5 * a * (uoldarray[i-1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i] - uoldarray[i-1]);
286 f[i] = uarray[i] - uoldarray[i] + lambda * (RFlux - LFlux);
287 }
288
289 /* Restore vectors */
290 ierr = DMDAVecRestoreArrayRead(da,U,&uarray);CHKERRQ(ierr);
291 ierr = DMDAVecRestoreArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);
292 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
293 ierr = DMRestoreLocalVector(da,&localUold);CHKERRQ(ierr);
294 PetscFunctionReturn(0);
295 }
296
297
298 /*TEST
299
300 test:
301 args: -ts_max_steps 10 -ts_monitor
302
303 test:
304 suffix: 2
305 nsize: 3
306 args: -ts_max_steps 10 -ts_monitor
307 output_file: output/ex6_1.out
308
309 test:
310 suffix: 3
311 args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false
312
313 test:
314 suffix: 4
315 nsize: 3
316 args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false
317 output_file: output/ex6_3.out
318
319 TEST*/
320