Chombo + EB + MF  3.2
GodunovUtilities.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _GODUNOVUTILITIES_H_
12 #define _GODUNOVUTILITIES_H_
13 
14 #include "Box.H"
15 #include "IntVectSet.H"
16 #include "Vector.H"
17 #include "REAL.H"
18 #include "FArrayBox.H"
19 #include "FluxBox.H"
20 #include "NamespaceHeader.H"
21 
22 class GodunovPhysics;
23 
24 ///
25 /**
26  Utility class for higher-order Godunov methods: slopes, parabolic
27  interpolants, limiters. Contains no physics-dependent methods, but
28  one of the member functions (PPMFaceValues()) may require a pointer
29  to a GodunovPhysics class in order to set boundary slopes.
30 */
32 {
33 public:
34  /// Constructor
35  /**
36  */
38 
39  /// Destructor
40  /**
41  */
43 
44  /// Define the object
45  /**
46  */
47  void define(const ProblemDomain& a_domain,
48  const Real& a_dx);
49 
50  /**
51  \name Slopes
52  */
53 
54  /**@{*/
55 
56  /// Compute componentwise van Leer slopes.
57  /**
58  Given cell averages W, compute van Leer slopes dW on a
59  component-by-component basis.
60  */
61  void vanLeerSlopes(FArrayBox& a_dW,
62  const FArrayBox& a_W,
63  const int& a_numSlopes,
64  const bool& a_useLimiting,
65  const int& a_dir,
66  const Box& a_box);
67 
68  /// Compute fourth-order slopes.
69  /**
70  Given cell averages W and van Leer slopes dWvL, compute fourth-order
71  slopes dW4. Limiting is performed in a separate pass.
72  */
73  void fourthOrderSlopes(FArrayBox& a_dW4,
74  const FArrayBox& a_W,
75  const FArrayBox& a_dWvL,
76  const int& a_numSlopes,
77  const int& a_dir,
78  const Box& a_box);
79 
80  /// Compute slopes (dW- and dW+) using one sided differences
81  /**
82  */
83  void oneSidedDifferences(FArrayBox& a_dWMinus,
84  FArrayBox& a_dWPlus,
85  const FArrayBox& a_W,
86  const int& a_dir,
87  const Box& a_box);
88 
89  /// Compute slopes (dW (center), dW-, and dW+)
90  /**
91  a_dwCent, a_dwMinus, a_dWPlus, and a_dW all live on
92  cell-centered a_entireBox.
93 
94  For i in a_centerBox:
95  a_dwCent[i] = (a_W[i+1] - a_W[i-1])/2;
96  a_dWMinus[i] = a_W[i] - a_W[i-1];
97  a_dWPlus[i] = a_W[i+1] - a_W[i].
98 
99  For i in a_loBox, set only a_dWPlus[i] = a_W[i+1] - a_W[i]
100  and copy it to a_dWCent[i] and a_dWMinus[i].
101 
102  For i in a_hiBox, set only a_dWMinus[i] = a_W[i] - a_W[i-1]
103  and copy it to a_dWCent[i] and a_dWPlus[i].
104  */
105  void slopes(FArrayBox& a_dWCent,
106  FArrayBox& a_dWMinus,
107  FArrayBox& a_dWPlus,
108  const FArrayBox& a_W,
109  const int& a_numSlopes,
110  const int& a_dir,
111  const Box& a_loBox,
112  const Box& a_hiBox,
113  const Box& a_centerBox,
114  const Box& a_entireBox,
115  const int& a_hasLo,
116  const int& a_hasHi);
117 
118  /**@}*/
119 
120  /**
121  \name Parabolic interpolants
122  */
123 
124  /**@{*/
125 
126  /// PLM normal predictor.
127  /**
128  Compute the increments in the characteristic amplitudes.
129  */
130  void PLMNormalPred(FArrayBox& a_dWCharMinus,
131  FArrayBox& a_dWCharPlus,
132  const FArrayBox& a_dWChar,
133  const FArrayBox& a_Lambda,
134  const Real& a_dtbydx,
135  const Box& a_box);
136 
137  /// PPM face-centered interpolant.
138  /**
139  Given the cell average a_W, compute fourth-order accurate
140  face-centered values WFace on a_box by differentiating the indefinite
141  integral. Limiting is performed in a separate pass.
142  */
143  void PPMFaceValues(FArrayBox& a_WFace,
144  const FArrayBox& a_W,
145  const int& a_numSlopes,
146  const bool& a_useLimiting,
147  const int& a_dir,
148  const Box& a_box,
149  const Real& a_time,
150  const GodunovPhysics* a_physPtr = NULL);
151 
152 
153  /// PPM normal predictor.
154  /**
155  On input, dW(Minus,Plus), contain the characteristic
156  expansions of the differences between the (minus, plus) face values
157  and the cell average. On output, dW(Minus,Plus) contain the
158  characteristic amplitudes of the corrections required to compute
159  the normal predictor.
160  */
161  void PPMNormalPred(FArrayBox& a_dWMinus,
162  FArrayBox& a_dWPlus,
163  const FArrayBox& a_Lambda,
164  const Real& a_dtbydx,
165  const int& a_numSlopes,
166  const Box& a_box);
167 
168  /**@}*/
169 
170  /**
171  \name Limiters
172  */
173 
174  /**@{*/
175 
176  /// van Leer slope limiter.
177  /**
178  On input, dW contains the centered, unlimited slopes, and
179  dW(Minus,Plus) contain the one-sided slopes from the minus, plus sides.
180  On output, dW contains the limited slopes.
181  slopes dW4. Limiting is performed in a separate pass.
182  */
183  void slopeLimiter(FArrayBox& a_dW,
184  const FArrayBox& a_dWMinus,
185  const FArrayBox& a_dWPlus,
186  const int& a_numSlopes,
187  const Box& a_box);
188 
189  /// extremum-preserving van Leer slope limiter.
190  /**
191  On input, dW contains the centered, unlimited slopes, and
192  dW(Minus,Plus) contain the one-sided slopes from the minus, plus sides.
193  On output, dW contains the limited slopes.
194  slopes dW4. Limiting is performed in a separate pass.
195  */
197  const FArrayBox& a_dWMinus,
198  const FArrayBox& a_dWPlus,
199  const int& a_numSlopes,
200  const Box& a_box,
201  const int& a_dir);
202 
203  /// PPM Limiter.
204  /**
205  a_dWMinus and a_dWPlus are the differences between the face values
206  on the minus and plus sides of cells and the average in the cell.
207  That is,
208  a_dWMinus[i] = WFace[i - e/2] - a_W[i]
209  a_dWPlus[i] = WFace[i + e/2] - a_W[i]
210  where e is the unit vector in dimension a_dir.
211  The PPM limiter is applied to these values to obtain a monotone
212  interpolant in the cell.
213  The function returns the limited a_dWMinus and a_dWPlus on a_box.
214  petermc, 4 Sep 2008: included a_W in argument list
215 
216  If m_highOrderLimiter,
217  then need a_dWMinus and a_dWPlus on a_box,
218  and need a_W on on a_box grown by 3 in dimension a_dir.
219  Returns limited a_dWMinus and a_dWPlus on a_box.
220  */
221  void PPMLimiter(FArrayBox& a_dWMinus,
222  FArrayBox& a_dWPlus,
223  const FArrayBox& a_W,
224  const int& a_numSlopes,
225  const int& a_dir,
226  const Box& a_box);
227 
228  /// Set whether to use high-order limiter.
229  void highOrderLimiter(bool a_highOrderLimiter);
230 
231  /// query whether we're using the high-order limiter
233 
234  /**@}*/
235 
236  /**
237  \name Dissipation mechanisms
238  */
239 
240  /**@{*/
241 
242  /// Compute the slope flattening coefficients
243  /**
244  Compute the slope flattening coefficients, a_flattening, using the
245  primitive variables, a_W, within a_box.
246  */
247  void computeFlattening(/// flattening coeffs, 1 component on a_box
248  FArrayBox& a_flattening,
249  /// primitive variables, on a_box grown by 3 within m_domain
250  const FArrayBox& a_W,
251  /// interval of a_W with velocity; length SpaceDim
252  const Interval& a_velInt,
253  /// index of a_W with pressure
254  const int& a_presInd,
255  /// minimum pressure to use in ratio
256  const Real& a_smallPres,
257  /// index of a_W with bulk modulus
258  const int& a_bulkModulusInd,
259  /// cell-centered box on which a_flattening lives
260  const Box& a_box);
261 
262  /// Apply the flattening to slopes.
263  /**
264  Multiply every component of a_dW by a_flat, on a_box.
265  */
266  void applyFlattening(/// slopes to be flattened, on at least a_box
267  FArrayBox& a_dW,
268  /// cell-centered flattening coefficients, on at least a_box
269  const FArrayBox& a_flat,
270  /// cell-centered box on which to flatten
271  const Box& a_box);
272 
273  /// Compute face-centered artificial viscosity flux.
274  /**
275  Increments face-centered flux in the a_dir direction with quadratic
276  artificial viscosity.
277  */
278  void artificialViscosity(FArrayBox& a_F,
279  const FArrayBox& a_U,
280  const FArrayBox& a_divVel,
281  const Real& a_scale,
282  const int& a_dir,
283  const Box& a_box);
284 
285  /**@}*/
286 
287 
288  /**
289  \name Divergence
290  */
291 
292  /**@{*/
293 
294  /// Compute face-centered velocity divergence.
295  /**
296  Returns face-centered velocity divergence on faces in the
297  direction a_dir. The velocities are the components a_velInterval
298  of a_W.
299  */
300  void divVel(FArrayBox& a_divVel,
301  const FArrayBox& a_W,
302  const Interval& a_velInt,
303  const int& a_dir,
304  const Box& a_box);
305 
306  /// Compute high-order face-centered velocity divergence for artificial viscosity
307  /**
308  Computes a face-centered nonlinear function of the divergence
309  suitable for use as an artificial viscosity coefficient for
310  a fourth-order method.
311 
312  F_visc = -dx*K*dU/dx
313  K = max(-dx*div(u)*min(1,h^2*|div(u)/(c^2)/a_M0sq|),0)
314  */
315  void divVelHO(FArrayBox& a_divVel,
316  const FArrayBox& a_W,
317  const int& a_dir,
318  const Box& a_box,
319  GodunovPhysics* a_physPtr);
320 
321  /**@}*/
322 
323  /**
324  \name Deconvolution and convolution
325  */
326 
327  /**@{*/
328 
329  /// Deconvolve from cell-centered to cell-averaged data (or vice versa)
330  /**
331  Deconvolves (or convolves, if a_sign = -1) by the formula
332 
333  a_avgFab += a_sign * laplacian(a_cellFab) * m_dx^2 / 24, on a_cellFab.box()
334  */
335  void deconvolve(FArrayBox& a_avgFab,
336  const FArrayBox& a_cellFab,
337  int a_sign = 1);
338 
339  /// Deconvolve from cell-centered to cell-averaged data (or vice versa) on a specified box
340  /**
341  Deconvolves (or convolves, if a_sign = -1) by the formula
342 
343  a_avgFab += a_sign * laplacian(a_cellFab) * m_dx^2 / 24, on a_box
344  */
345  void deconvolve(FArrayBox& a_avgFab,
346  const FArrayBox& a_cellFab,
347  const Box& a_box,
348  int a_sign = 1);
349 
350  /// Deconvolve from face-centered to face-averaged data (or vice versa)
351  /**
352  Deconvolves (or convolves, if a_sign = -1) by the formula
353 
354  a_avgFlux += a_sign * laplacian(a_cellFlux) * m_dx^2 / 24, on faces of a_cellFlux.box()
355  */
356  void deconvolveFace(FluxBox& a_avgFlux,
357  const FluxBox& a_cellFlux,
358  int a_sign = 1);
359 
360  /// Deconvolve from face-centered to face-averaged data (or vice versa) on a specified box
361  /**
362  Deconvolves (or convolves, if a_sign = -1) by the formula
363 
364  a_avgFlux += a_sign * laplacian(a_cellFlux) * m_dx^2 / 24, on faces of a_box
365  */
366  void deconvolveFace(FluxBox& a_avgFlux,
367  const FluxBox& a_cellFlux,
368  const Box& a_box,
369  int a_sign = 1);
370 
371  /**@}*/
372 
373 protected:
374  // Problem domain and grid spacing
377 
378  // Has this object been defined
380 
381  // Use a high-order limiter? (default false)
383 
384 private:
385  // We may allow copying and assignment later.
386  // Disallowed for all the usual reasons
387  void operator=(const GodunovUtilities&);
389 };
390 
391 #include "NamespaceFooter.H"
392 #endif
void deconvolve(FArrayBox &a_avgFab, const FArrayBox &a_cellFab, int a_sign=1)
Deconvolve from cell-centered to cell-averaged data (or vice versa)
void slopeLimiterExtPreserving(FArrayBox &a_dW, const FArrayBox &a_dWMinus, const FArrayBox &a_dWPlus, const int &a_numSlopes, const Box &a_box, const int &a_dir)
extremum-preserving van Leer slope limiter.
bool useHighOrderLimiter() const
query whether we're using the high-order limiter
Definition: GodunovUtilities.H:232
void PPMLimiter(FArrayBox &a_dWMinus, FArrayBox &a_dWPlus, const FArrayBox &a_W, const int &a_numSlopes, const int &a_dir, const Box &a_box)
PPM Limiter.
void deconvolveFace(FluxBox &a_avgFlux, const FluxBox &a_cellFlux, int a_sign=1)
Deconvolve from face-centered to face-averaged data (or vice versa)
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
Definition: GodunovUtilities.H:31
void PPMNormalPred(FArrayBox &a_dWMinus, FArrayBox &a_dWPlus, const FArrayBox &a_Lambda, const Real &a_dtbydx, const int &a_numSlopes, const Box &a_box)
PPM normal predictor.
void applyFlattening(FArrayBox &a_dW, const FArrayBox &a_flat, const Box &a_box)
Apply the flattening to slopes.
GodunovUtilities()
Constructor.
void oneSidedDifferences(FArrayBox &a_dWMinus, FArrayBox &a_dWPlus, const FArrayBox &a_W, const int &a_dir, const Box &a_box)
Compute slopes (dW- and dW+) using one sided differences.
bool m_highOrderLimiter
Definition: GodunovUtilities.H:382
void highOrderLimiter(bool a_highOrderLimiter)
Set whether to use high-order limiter.
A FArrayBox-like container for face-centered fluxes.
Definition: FluxBox.H:22
Structure for passing component ranges in code.
Definition: Interval.H:23
void divVelHO(FArrayBox &a_divVel, const FArrayBox &a_W, const int &a_dir, const Box &a_box, GodunovPhysics *a_physPtr)
Compute high-order face-centered velocity divergence for artificial viscosity.
ProblemDomain m_domain
Definition: GodunovUtilities.H:375
void operator=(const GodunovUtilities &)
void PPMFaceValues(FArrayBox &a_WFace, const FArrayBox &a_W, const int &a_numSlopes, const bool &a_useLimiting, const int &a_dir, const Box &a_box, const Real &a_time, const GodunovPhysics *a_physPtr=NULL)
PPM face-centered interpolant.
double Real
Definition: REAL.H:33
bool m_isDefined
Definition: GodunovUtilities.H:379
Real m_dx
Definition: GodunovUtilities.H:376
void PLMNormalPred(FArrayBox &a_dWCharMinus, FArrayBox &a_dWCharPlus, const FArrayBox &a_dWChar, const FArrayBox &a_Lambda, const Real &a_dtbydx, const Box &a_box)
PLM normal predictor.
Definition: GodunovPhysics.H:40
void divVel(FArrayBox &a_divVel, const FArrayBox &a_W, const Interval &a_velInt, const int &a_dir, const Box &a_box)
Compute face-centered velocity divergence.
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
void vanLeerSlopes(FArrayBox &a_dW, const FArrayBox &a_W, const int &a_numSlopes, const bool &a_useLimiting, const int &a_dir, const Box &a_box)
Compute componentwise van Leer slopes.
void slopeLimiter(FArrayBox &a_dW, const FArrayBox &a_dWMinus, const FArrayBox &a_dWPlus, const int &a_numSlopes, const Box &a_box)
van Leer slope limiter.
void fourthOrderSlopes(FArrayBox &a_dW4, const FArrayBox &a_W, const FArrayBox &a_dWvL, const int &a_numSlopes, const int &a_dir, const Box &a_box)
Compute fourth-order slopes.
~GodunovUtilities()
Destructor.
void computeFlattening(FArrayBox &a_flattening, const FArrayBox &a_W, const Interval &a_velInt, const int &a_presInd, const Real &a_smallPres, const int &a_bulkModulusInd, const Box &a_box)
Compute the slope flattening coefficients.
Definition: FArrayBox.H:45
void artificialViscosity(FArrayBox &a_F, const FArrayBox &a_U, const FArrayBox &a_divVel, const Real &a_scale, const int &a_dir, const Box &a_box)
Compute face-centered artificial viscosity flux.
void slopes(FArrayBox &a_dWCent, FArrayBox &a_dWMinus, FArrayBox &a_dWPlus, const FArrayBox &a_W, const int &a_numSlopes, const int &a_dir, const Box &a_loBox, const Box &a_hiBox, const Box &a_centerBox, const Box &a_entireBox, const int &a_hasLo, const int &a_hasHi)
Compute slopes (dW (center), dW-, and dW+)
void define(const ProblemDomain &a_domain, const Real &a_dx)
Define the object.