Actual source code: ex6.c
1: static char help[] = "Model Equations for Advection \n";
3: /*
4: Modified from ex3.c
5: Page 9, Section 1.2 Model Equations for Advection-Diffusion
7: u_t + a u_x = 0, 0<= x <= 1.0
9: The initial conditions used here different from the book.
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: */
16: #include <petscts.h>
17: #include <petscdm.h>
18: #include <petscdmda.h>
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;
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 *);
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;
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));
50: appctx.a = -1.0;
51: PetscCall(PetscOptionsGetReal(NULL, NULL, "-a", &appctx.a, NULL));
53: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 60, 1, 1, NULL, &da));
54: PetscCall(DMSetFromOptions(da));
55: PetscCall(DMSetUp(da));
57: /* Create vector data structures for approximate and exact solutions */
58: PetscCall(DMCreateGlobalVector(da, &U));
60: /* Create timestepping solver context */
61: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
62: PetscCall(TSSetDM(ts, da));
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: }
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));
84: /* Evaluate initial conditions */
85: PetscCall(InitialConditions(ts, U, &appctx));
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));
90: /* Run the timestepping solver */
91: PetscCall(TSSolve(ts, U));
93: /* Free work space */
94: PetscCall(TSDestroy(&ts));
95: PetscCall(VecDestroy(&U));
96: PetscCall(DMDestroy(&da));
98: PetscCall(PetscFinalize());
99: return 0;
100: }
101: /* --------------------------------------------------------------------- */
102: /*
103: InitialConditions - Computes the solution at the initial time.
105: Input Parameter:
106: u - uninitialized solution vector (global)
107: appctx - user-defined application context
109: Output Parameter:
110: u - vector with solution at initial time (global)
111: */
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;
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));
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);
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
151: Input Parameters:
152: t - current time
153: solution - vector in which exact solution will be computed
154: appctx - user-defined application context
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: */
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;
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;
174: /* Get a pointer to vector data. */
175: PetscCall(DMDAVecGetArray(da, U, &u));
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));
182: /* Restore vector */
183: PetscCall(DMDAVecRestoreArray(da, U, &u));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /* --------------------------------------------------------------------- */
188: /*
189: Use Lax-Friedrichs method to evaluate F(u,t) = du/dt + a * du/dx
191: See https://en.wikipedia.org/wiki/Lax%E2%80%93Friedrichs_method
192: */
193: PetscErrorCode IFunction_LaxFriedrichs(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *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;
202: PetscFunctionBegin;
203: PetscCall(TSGetTimeStep(ts, &dt));
204: PetscCall(TSGetSolution(ts, &Uold));
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); */
213: PetscCall(DMGetLocalVector(da, &localUold));
214: PetscCall(DMGlobalToLocalBegin(da, Uold, INSERT_VALUES, localUold));
215: PetscCall(DMGlobalToLocalEnd(da, Uold, INSERT_VALUES, localUold));
217: /* Get pointers to vector data */
218: PetscCall(DMDAVecGetArrayRead(da, U, &uarray));
219: PetscCall(DMDAVecGetArrayRead(da, localUold, &uoldarray));
220: PetscCall(DMDAVecGetArray(da, F, &f));
222: /* advection */
223: c = appctx->a * dt / h; /* Courant-Friedrichs-Lewy number (CFL number) */
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: }
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: }
238: /*
239: Use Lax-Wendroff method to evaluate F(u,t) = du/dt + a * du/dx
240: */
241: PetscErrorCode IFunction_LaxWendroff(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *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;
250: PetscFunctionBegin;
251: PetscCall(TSGetTimeStep(ts, &dt));
252: PetscCall(TSGetSolution(ts, &Uold));
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); */
261: PetscCall(DMGetLocalVector(da, &localUold));
262: PetscCall(DMGlobalToLocalBegin(da, Uold, INSERT_VALUES, localUold));
263: PetscCall(DMGlobalToLocalEnd(da, Uold, INSERT_VALUES, localUold));
265: /* Get pointers to vector data */
266: PetscCall(DMDAVecGetArrayRead(da, U, &uarray));
267: PetscCall(DMDAVecGetArrayRead(da, localUold, &uoldarray));
268: PetscCall(DMDAVecGetArray(da, F, &f));
270: /* advection -- finite volume (appctx->a < 0 -- can be relaxed?) */
271: lambda = dt / h;
272: a = appctx->a;
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: }
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: }
288: /*TEST
290: test:
291: args: -ts_max_steps 10 -ts_monitor
293: test:
294: suffix: 2
295: nsize: 3
296: args: -ts_max_steps 10 -ts_monitor
297: output_file: output/ex6_1.out
299: test:
300: suffix: 3
301: args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false
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
309: TEST*/