00001 #ifdef CH_LANG_CC 00002 /* 00003 * _______ __ 00004 * / ___/ / ___ __ _ / / ___ 00005 * / /__/ _ \/ _ \/ V \/ _ \/ _ \ 00006 * \___/_//_/\___/_/_/_/_.__/\___/ 00007 * Please refer to Copyright.txt, in Chombo's root directory. 00008 */ 00009 #endif 00010 00011 #ifndef _GODUNOVUTILITIES_H_ 00012 #define _GODUNOVUTILITIES_H_ 00013 00014 #include "Box.H" 00015 #include "IntVectSet.H" 00016 #include "Vector.H" 00017 #include "REAL.H" 00018 #include "FArrayBox.H" 00019 #include "FluxBox.H" 00020 #include "NamespaceHeader.H" 00021 00022 class GodunovPhysics; 00023 00024 /// 00025 /** 00026 Utility class for higher-order Godunov methods: slopes, parabolic 00027 interpolants, limiters. Contains no physics-dependent methods, but 00028 one of the member functions (PPMFaceValues()) may require a pointer 00029 to a GodunovPhysics class in order to set boundary slopes. 00030 */ 00031 class GodunovUtilities 00032 { 00033 public: 00034 /// Constructor 00035 /** 00036 */ 00037 GodunovUtilities(); 00038 00039 /// Destructor 00040 /** 00041 */ 00042 ~GodunovUtilities(); 00043 00044 /// Define the object 00045 /** 00046 */ 00047 void define(const ProblemDomain& a_domain, 00048 const Real& a_dx); 00049 00050 /** 00051 \name Slopes 00052 */ 00053 00054 /**@{*/ 00055 00056 /// Compute componentwise van Leer slopes. 00057 /** 00058 Given cell averages W, compute van Leer slopes dW on a 00059 component-by-component basis. 00060 */ 00061 void vanLeerSlopes(FArrayBox& a_dW, 00062 const FArrayBox& a_W, 00063 const int& a_numSlopes, 00064 const bool& a_useLimiting, 00065 const int& a_dir, 00066 const Box& a_box); 00067 00068 /// Compute fourth-order slopes. 00069 /** 00070 Given cell averages W and van Leer slopes dWvL, compute fourth-order 00071 slopes dW4. Limiting is performed in a separate pass. 00072 */ 00073 void fourthOrderSlopes(FArrayBox& a_dW4, 00074 const FArrayBox& a_W, 00075 const FArrayBox& a_dWvL, 00076 const int& a_numSlopes, 00077 const int& a_dir, 00078 const Box& a_box); 00079 00080 /// Compute slopes (dW- and dW+) using one sided differences 00081 /** 00082 */ 00083 void oneSidedDifferences(FArrayBox& a_dWMinus, 00084 FArrayBox& a_dWPlus, 00085 const FArrayBox& a_W, 00086 const int& a_dir, 00087 const Box& a_box); 00088 00089 /// Compute slopes (dW (center), dW-, and dW+) 00090 /** 00091 a_dwCent, a_dwMinus, a_dWPlus, and a_dW all live on 00092 cell-centered a_entireBox. 00093 00094 For i in a_centerBox: 00095 a_dwCent[i] = (a_W[i+1] - a_W[i-1])/2; 00096 a_dWMinus[i] = a_W[i] - a_W[i-1]; 00097 a_dWPlus[i] = a_W[i+1] - a_W[i]. 00098 00099 For i in a_loBox, set only a_dWPlus[i] = a_W[i+1] - a_W[i] 00100 and copy it to a_dWCent[i] and a_dWMinus[i]. 00101 00102 For i in a_hiBox, set only a_dWMinus[i] = a_W[i] - a_W[i-1] 00103 and copy it to a_dWCent[i] and a_dWPlus[i]. 00104 */ 00105 void slopes(FArrayBox& a_dWCent, 00106 FArrayBox& a_dWMinus, 00107 FArrayBox& a_dWPlus, 00108 const FArrayBox& a_W, 00109 const int& a_numSlopes, 00110 const int& a_dir, 00111 const Box& a_loBox, 00112 const Box& a_hiBox, 00113 const Box& a_centerBox, 00114 const Box& a_entireBox, 00115 const int& a_hasLo, 00116 const int& a_hasHi); 00117 00118 /**@}*/ 00119 00120 /** 00121 \name Parabolic interpolants 00122 */ 00123 00124 /**@{*/ 00125 00126 /// PLM normal predictor. 00127 /** 00128 Compute the increments in the characteristic amplitudes. 00129 */ 00130 void PLMNormalPred(FArrayBox& a_dWCharMinus, 00131 FArrayBox& a_dWCharPlus, 00132 const FArrayBox& a_dWChar, 00133 const FArrayBox& a_Lambda, 00134 const Real& a_dtbydx, 00135 const Box& a_box); 00136 00137 /// PPM face-centered interpolant. 00138 /** 00139 Given the cell average a_W, compute fourth-order accurate 00140 face-centered values WFace on a_box by differentiating the indefinite 00141 integral. Limiting is performed in a separate pass. 00142 */ 00143 void PPMFaceValues(FArrayBox& a_WFace, 00144 const FArrayBox& a_W, 00145 const int& a_numSlopes, 00146 const bool& a_useLimiting, 00147 const int& a_dir, 00148 const Box& a_box, 00149 const Real& a_time, 00150 const GodunovPhysics* a_physPtr = NULL); 00151 00152 00153 /// PPM normal predictor. 00154 /** 00155 On input, dW(Minus,Plus), contain the characteristic 00156 expansions of the differences between the (minus, plus) face values 00157 and the cell average. On output, dW(Minus,Plus) contain the 00158 characteristic amplitudes of the corrections required to compute 00159 the normal predictor. 00160 */ 00161 void PPMNormalPred(FArrayBox& a_dWMinus, 00162 FArrayBox& a_dWPlus, 00163 const FArrayBox& a_Lambda, 00164 const Real& a_dtbydx, 00165 const int& a_numSlopes, 00166 const Box& a_box); 00167 00168 /**@}*/ 00169 00170 /** 00171 \name Limiters 00172 */ 00173 00174 /**@{*/ 00175 00176 /// van Leer slope limiter. 00177 /** 00178 On input, dW contains the centered, unlimited slopes, and 00179 dW(Minus,Plus) contain the one-sided slopes from the minus, plus sides. 00180 On output, dW contains the limited slopes. 00181 slopes dW4. Limiting is performed in a separate pass. 00182 */ 00183 void slopeLimiter(FArrayBox& a_dW, 00184 const FArrayBox& a_dWMinus, 00185 const FArrayBox& a_dWPlus, 00186 const int& a_numSlopes, 00187 const Box& a_box); 00188 00189 /// extremum-preserving van Leer slope limiter. 00190 /** 00191 On input, dW contains the centered, unlimited slopes, and 00192 dW(Minus,Plus) contain the one-sided slopes from the minus, plus sides. 00193 On output, dW contains the limited slopes. 00194 slopes dW4. Limiting is performed in a separate pass. 00195 */ 00196 void slopeLimiterExtPreserving(FArrayBox& a_dW, 00197 const FArrayBox& a_dWMinus, 00198 const FArrayBox& a_dWPlus, 00199 const int& a_numSlopes, 00200 const Box& a_box, 00201 const int& a_dir); 00202 00203 /// PPM Limiter. 00204 /** 00205 a_dWMinus and a_dWPlus are the differences between the face values 00206 on the minus and plus sides of cells and the average in the cell. 00207 That is, 00208 a_dWMinus[i] = WFace[i - e/2] - a_W[i] 00209 a_dWPlus[i] = WFace[i + e/2] - a_W[i] 00210 where e is the unit vector in dimension a_dir. 00211 The PPM limiter is applied to these values to obtain a monotone 00212 interpolant in the cell. 00213 The function returns the limited a_dWMinus and a_dWPlus on a_box. 00214 petermc, 4 Sep 2008: included a_W in argument list 00215 00216 If m_highOrderLimiter, 00217 then need a_dWMinus and a_dWPlus on a_box, 00218 and need a_W on on a_box grown by 3 in dimension a_dir. 00219 Returns limited a_dWMinus and a_dWPlus on a_box. 00220 */ 00221 void PPMLimiter(FArrayBox& a_dWMinus, 00222 FArrayBox& a_dWPlus, 00223 const FArrayBox& a_W, 00224 const int& a_numSlopes, 00225 const int& a_dir, 00226 const Box& a_box); 00227 00228 /// Set whether to use high-order limiter. 00229 void highOrderLimiter(bool a_highOrderLimiter); 00230 00231 /// query whether we're using the high-order limiter 00232 bool useHighOrderLimiter() const {return m_highOrderLimiter;} 00233 00234 /**@}*/ 00235 00236 /** 00237 \name Dissipation mechanisms 00238 */ 00239 00240 /**@{*/ 00241 00242 /// Compute the slope flattening coefficients 00243 /** 00244 Compute the slope flattening coefficients, a_flattening, using the 00245 primitive variables, a_W, within a_box. 00246 */ 00247 void computeFlattening(/// flattening coeffs, 1 component on a_box 00248 FArrayBox& a_flattening, 00249 /// primitive variables, on a_box grown by 3 within m_domain 00250 const FArrayBox& a_W, 00251 /// interval of a_W with velocity; length SpaceDim 00252 const Interval& a_velInt, 00253 /// index of a_W with pressure 00254 const int& a_presInd, 00255 /// minimum pressure to use in ratio 00256 const Real& a_smallPres, 00257 /// index of a_W with bulk modulus 00258 const int& a_bulkModulusInd, 00259 /// cell-centered box on which a_flattening lives 00260 const Box& a_box); 00261 00262 /// Apply the flattening to slopes. 00263 /** 00264 Multiply every component of a_dW by a_flat, on a_box. 00265 */ 00266 void applyFlattening(/// slopes to be flattened, on at least a_box 00267 FArrayBox& a_dW, 00268 /// cell-centered flattening coefficients, on at least a_box 00269 const FArrayBox& a_flat, 00270 /// cell-centered box on which to flatten 00271 const Box& a_box); 00272 00273 /// Compute face-centered artificial viscosity flux. 00274 /** 00275 Increments face-centered flux in the a_dir direction with quadratic 00276 artificial viscosity. 00277 */ 00278 void artificialViscosity(FArrayBox& a_F, 00279 const FArrayBox& a_U, 00280 const FArrayBox& a_divVel, 00281 const Real& a_scale, 00282 const int& a_dir, 00283 const Box& a_box); 00284 00285 /**@}*/ 00286 00287 00288 /** 00289 \name Divergence 00290 */ 00291 00292 /**@{*/ 00293 00294 /// Compute face-centered velocity divergence. 00295 /** 00296 Returns face-centered velocity divergence on faces in the 00297 direction a_dir. The velocities are the components a_velInterval 00298 of a_W. 00299 */ 00300 void divVel(FArrayBox& a_divVel, 00301 const FArrayBox& a_W, 00302 const Interval& a_velInt, 00303 const int& a_dir, 00304 const Box& a_box); 00305 00306 /// Compute high-order face-centered velocity divergence for artificial viscosity 00307 /** 00308 Computes a face-centered nonlinear function of the divergence 00309 suitable for use as an artificial viscosity coefficient for 00310 a fourth-order method. 00311 00312 F_visc = -dx*K*dU/dx 00313 K = max(-dx*div(u)*min(1,h^2*|div(u)/(c^2)/a_M0sq|),0) 00314 */ 00315 void divVelHO(FArrayBox& a_divVel, 00316 const FArrayBox& a_W, 00317 const int& a_dir, 00318 const Box& a_box, 00319 GodunovPhysics* a_physPtr); 00320 00321 /**@}*/ 00322 00323 /** 00324 \name Deconvolution and convolution 00325 */ 00326 00327 /**@{*/ 00328 00329 /// Deconvolve from cell-centered to cell-averaged data (or vice versa) 00330 /** 00331 Deconvolves (or convolves, if a_sign = -1) by the formula 00332 00333 a_avgFab += a_sign * laplacian(a_cellFab) * m_dx^2 / 24, on a_cellFab.box() 00334 */ 00335 void deconvolve(FArrayBox& a_avgFab, 00336 const FArrayBox& a_cellFab, 00337 int a_sign = 1); 00338 00339 /// Deconvolve from cell-centered to cell-averaged data (or vice versa) on a specified box 00340 /** 00341 Deconvolves (or convolves, if a_sign = -1) by the formula 00342 00343 a_avgFab += a_sign * laplacian(a_cellFab) * m_dx^2 / 24, on a_box 00344 */ 00345 void deconvolve(FArrayBox& a_avgFab, 00346 const FArrayBox& a_cellFab, 00347 const Box& a_box, 00348 int a_sign = 1); 00349 00350 /// Deconvolve from face-centered to face-averaged data (or vice versa) 00351 /** 00352 Deconvolves (or convolves, if a_sign = -1) by the formula 00353 00354 a_avgFlux += a_sign * laplacian(a_cellFlux) * m_dx^2 / 24, on faces of a_cellFlux.box() 00355 */ 00356 void deconvolveFace(FluxBox& a_avgFlux, 00357 const FluxBox& a_cellFlux, 00358 int a_sign = 1); 00359 00360 /// Deconvolve from face-centered to face-averaged data (or vice versa) on a specified box 00361 /** 00362 Deconvolves (or convolves, if a_sign = -1) by the formula 00363 00364 a_avgFlux += a_sign * laplacian(a_cellFlux) * m_dx^2 / 24, on faces of a_box 00365 */ 00366 void deconvolveFace(FluxBox& a_avgFlux, 00367 const FluxBox& a_cellFlux, 00368 const Box& a_box, 00369 int a_sign = 1); 00370 00371 /**@}*/ 00372 00373 protected: 00374 // Problem domain and grid spacing 00375 ProblemDomain m_domain; 00376 Real m_dx; 00377 00378 // Has this object been defined 00379 bool m_isDefined; 00380 00381 // Use a high-order limiter? (default false) 00382 bool m_highOrderLimiter; 00383 00384 private: 00385 // We may allow copying and assignment later. 00386 // Disallowed for all the usual reasons 00387 void operator=(const GodunovUtilities&); 00388 GodunovUtilities(const GodunovUtilities&); 00389 }; 00390 00391 #include "NamespaceFooter.H" 00392 #endif