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*/