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