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