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 */
 48:   PetscInitialize(&argc, &argv, (char *)0, help);
 49:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 51:   appctx.a = -1.0;
 52:   PetscOptionsGetReal(NULL, NULL, "-a", &appctx.a, NULL);

 54:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 60, 1, 1, NULL, &da);
 55:   DMSetFromOptions(da);
 56:   DMSetUp(da);

 58:   /* Create vector data structures for approximate and exact solutions */
 59:   DMCreateGlobalVector(da, &U);

 61:   /* Create timestepping solver context */
 62:   TSCreate(PETSC_COMM_WORLD, &ts);
 63:   TSSetDM(ts, da);

 65:   /* Function evaluation */
 66:   PetscOptionsGetBool(NULL, NULL, "-useLaxWendroff", &useLaxWendroff, NULL);
 67:   if (useLaxWendroff) {
 68:     if (rank == 0) PetscPrintf(PETSC_COMM_SELF, "... Use Lax-Wendroff finite volume\n");
 69:     TSSetIFunction(ts, NULL, IFunction_LaxWendroff, &appctx);
 70:   } else {
 71:     if (rank == 0) PetscPrintf(PETSC_COMM_SELF, "... Use Lax-LaxFriedrichs finite difference\n");
 72:     TSSetIFunction(ts, NULL, IFunction_LaxFriedrichs, &appctx);
 73:   }

 75:   /* Customize timestepping solver */
 76:   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:   TSSetTimeStep(ts, dt);
 79:   TSSetMaxSteps(ts, 100);
 80:   TSSetMaxTime(ts, 100.0);
 81:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 82:   TSSetType(ts, TSBEULER);
 83:   TSSetFromOptions(ts);

 85:   /* Evaluate initial conditions */
 86:   InitialConditions(ts, U, &appctx);

 88:   /* For testing accuracy of TS with already known solution, e.g., '-ts_monitor_lg_error' */
 89:   TSSetSolutionFunction(ts, (PetscErrorCode(*)(TS, PetscReal, Vec, void *))Solution, &appctx);

 91:   /* Run the timestepping solver */
 92:   TSSolve(ts, U);

 94:   /* Free work space */
 95:   TSDestroy(&ts);
 96:   VecDestroy(&U);
 97:   DMDestroy(&da);

 99:   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:   TSGetDM(ts, &da);
121:   DMDAGetCorners(da, &mstart, 0, 0, &um, 0, 0);
122:   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:   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:   DMDAVecRestoreArray(da, U, &u);
145:   return 0;
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:   TSGetDM(ts, &da);
168:   DMDAGetCorners(da, &mstart, 0, 0, &um, 0, 0);
169:   DMDAGetInfo(da, PETSC_IGNORE, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
170:   h    = 1.0 / M;
171:   mend = mstart + um;

173:   /* Get a pointer to vector data. */
174:   DMDAVecGetArray(da, U, &u);

176:   /* u[i] = sin(6*PI*(x[i] - a*t)) + 3 * sin(2*PI*(x[i] - a*t)) */
177:   PI6 = PETSC_PI * 6.;
178:   PI2 = PETSC_PI * 2.;
179:   for (i = mstart; i < mend; i++) u[i] = PetscSinReal(PI6 * (i * h - a * t)) + 3. * PetscSinReal(PI2 * (i * h - a * t));

181:   /* Restore vector */
182:   DMDAVecRestoreArray(da, U, &u);
183:   return 0;
184: }

186: /* --------------------------------------------------------------------- */
187: /*
188:  Use Lax-Friedrichs method to evaluate F(u,t) = du/dt + a *  du/dx

190:  See https://en.wikipedia.org/wiki/Lax%E2%80%93Friedrichs_method
191:  */
192: PetscErrorCode IFunction_LaxFriedrichs(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
193: {
194:   AppCtx      *appctx = (AppCtx *)ctx;
195:   PetscInt     mstart, mend, M, i, um;
196:   DM           da;
197:   Vec          Uold, localUold;
198:   PetscScalar *uarray, *f, *uoldarray, h, uave, c;
199:   PetscReal    dt;

201:   TSGetTimeStep(ts, &dt);
202:   TSGetSolution(ts, &Uold);

204:   TSGetDM(ts, &da);
205:   DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
206:   DMDAGetCorners(da, &mstart, 0, 0, &um, 0, 0);
207:   h    = 1.0 / M;
208:   mend = mstart + um;
209:   /* printf(" mstart %d, um %d\n",mstart,um); */

211:   DMGetLocalVector(da, &localUold);
212:   DMGlobalToLocalBegin(da, Uold, INSERT_VALUES, localUold);
213:   DMGlobalToLocalEnd(da, Uold, INSERT_VALUES, localUold);

215:   /* Get pointers to vector data */
216:   DMDAVecGetArrayRead(da, U, &uarray);
217:   DMDAVecGetArrayRead(da, localUold, &uoldarray);
218:   DMDAVecGetArray(da, F, &f);

220:   /* advection */
221:   c = appctx->a * dt / h; /* Courant-Friedrichs-Lewy number (CFL number) */

223:   for (i = mstart; i < mend; i++) {
224:     uave = 0.5 * (uoldarray[i - 1] + uoldarray[i + 1]);
225:     f[i] = uarray[i] - uave + c * 0.5 * (uoldarray[i + 1] - uoldarray[i - 1]);
226:   }

228:   /* Restore vectors */
229:   DMDAVecRestoreArrayRead(da, U, &uarray);
230:   DMDAVecRestoreArrayRead(da, localUold, &uoldarray);
231:   DMDAVecRestoreArray(da, F, &f);
232:   DMRestoreLocalVector(da, &localUold);
233:   return 0;
234: }

236: /*
237:  Use Lax-Wendroff method to evaluate F(u,t) = du/dt + a *  du/dx
238: */
239: PetscErrorCode IFunction_LaxWendroff(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
240: {
241:   AppCtx      *appctx = (AppCtx *)ctx;
242:   PetscInt     mstart, mend, M, i, um;
243:   DM           da;
244:   Vec          Uold, localUold;
245:   PetscScalar *uarray, *f, *uoldarray, h, RFlux, LFlux, lambda;
246:   PetscReal    dt, a;

248:   TSGetTimeStep(ts, &dt);
249:   TSGetSolution(ts, &Uold);

251:   TSGetDM(ts, &da);
252:   DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
253:   DMDAGetCorners(da, &mstart, 0, 0, &um, 0, 0);
254:   h    = 1.0 / M;
255:   mend = mstart + um;
256:   /* printf(" mstart %d, um %d\n",mstart,um); */

258:   DMGetLocalVector(da, &localUold);
259:   DMGlobalToLocalBegin(da, Uold, INSERT_VALUES, localUold);
260:   DMGlobalToLocalEnd(da, Uold, INSERT_VALUES, localUold);

262:   /* Get pointers to vector data */
263:   DMDAVecGetArrayRead(da, U, &uarray);
264:   DMDAVecGetArrayRead(da, localUold, &uoldarray);
265:   DMDAVecGetArray(da, F, &f);

267:   /* advection -- finite volume (appctx->a < 0 -- can be relaxed?) */
268:   lambda = dt / h;
269:   a      = appctx->a;

271:   for (i = mstart; i < mend; i++) {
272:     RFlux = 0.5 * a * (uoldarray[i + 1] + uoldarray[i]) - a * a * 0.5 * lambda * (uoldarray[i + 1] - uoldarray[i]);
273:     LFlux = 0.5 * a * (uoldarray[i - 1] + uoldarray[i]) - a * a * 0.5 * lambda * (uoldarray[i] - uoldarray[i - 1]);
274:     f[i]  = uarray[i] - uoldarray[i] + lambda * (RFlux - LFlux);
275:   }

277:   /* Restore vectors */
278:   DMDAVecRestoreArrayRead(da, U, &uarray);
279:   DMDAVecRestoreArrayRead(da, localUold, &uoldarray);
280:   DMDAVecRestoreArray(da, F, &f);
281:   DMRestoreLocalVector(da, &localUold);
282:   return 0;
283: }

285: /*TEST

287:    test:
288:       args: -ts_max_steps 10 -ts_monitor

290:    test:
291:       suffix: 2
292:       nsize: 3
293:       args: -ts_max_steps 10 -ts_monitor
294:       output_file: output/ex6_1.out

296:    test:
297:       suffix: 3
298:       args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false

300:    test:
301:       suffix: 4
302:       nsize: 3
303:       args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false
304:       output_file: output/ex6_3.out

306: TEST*/