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 _GODUNOVPHYSICS_H_ 00012 #define _GODUNOVPHYSICS_H_ 00013 00014 #include <string> 00015 using std::string; 00016 00017 #include "Box.H" 00018 #include "IntVectSet.H" 00019 #include "Vector.H" 00020 #include "PhysIBC.H" 00021 #include "FluxBox.H" 00022 00023 #include "GodunovUtilities.H" 00024 #include "NamespaceHeader.H" 00025 00026 class HDF5HeaderData; 00027 00028 /// 00029 /** 00030 The base class GodunovPhysics provides the physics-dependent components 00031 for a higher-order Godunov method for a single patch: characteristic 00032 analysis, Riemann solver, quasilinear update, conservative update, 00033 and transformations between conserved, primitive, and flux variables. 00034 This class is essentially pure virtual; i.e., all of its member functions 00035 are virtual; and the ones that have default implementations are ones 00036 that are optionally defined; i.e., the default definition is to send 00037 an error message. Physics-dependent versions of this class that are 00038 required in real applications are derived from this class by inheritance. 00039 */ 00040 class GodunovPhysics 00041 { 00042 public: 00043 /// Constructor 00044 /** 00045 */ 00046 GodunovPhysics(); 00047 00048 /// Get the initial and boundary condition object 00049 /** 00050 */ 00051 PhysIBC* getPhysIBC() const; 00052 00053 /// Set the initial and boundary condition object 00054 /** 00055 */ 00056 void setPhysIBC(PhysIBC* a_bc); 00057 00058 /// Destructor 00059 /** 00060 */ 00061 virtual ~GodunovPhysics(); 00062 00063 /// Define the object 00064 /** 00065 */ 00066 virtual void define(const ProblemDomain& a_domain, 00067 const Real& a_dx); 00068 00069 /// Set the current box (default implementation - do nothing) 00070 /** 00071 */ 00072 virtual void setCurrentBox(const Box& a_currentBox); 00073 00074 /// Compute the maximum wave speed. 00075 /** 00076 */ 00077 virtual Real getMaxWaveSpeed(const FArrayBox& a_U, 00078 const Box& a_box) = 0; 00079 00080 /// Compute the maximum wave speed 00081 /** 00082 */ 00083 virtual void soundSpeed(FArrayBox& a_speed, 00084 const FArrayBox& a_U, 00085 const Box& a_box); 00086 00087 /// Object factory for this class 00088 /** 00089 */ 00090 virtual GodunovPhysics* new_godunovPhysics() const = 0; 00091 00092 /// Compute the increment in the conserved variables from face variables. 00093 /** 00094 Compute dU = dt*dUdt, the change in the conserved variables over 00095 the time step. The fluxes are returned are suitable for use in refluxing. 00096 This has a default implementation but can be redefined as needed. 00097 */ 00098 virtual void computeUpdate(FArrayBox& a_dU, 00099 FluxBox& a_F, 00100 const FArrayBox& a_U, 00101 const FluxBox& a_WHalf, 00102 const bool& a_useArtificialViscosity, 00103 const Real& a_artificialViscosity, 00104 const Real& a_currentTime, 00105 const Real& a_dx, 00106 const Real& a_dt, 00107 const Box& a_box); 00108 00109 /// Compute the fluxes from primitive variable values on a face 00110 /** 00111 This has a default implementation which throws an error. The method is 00112 here so that the default implementation of "computeUpdate" can use it 00113 and the user can supply it. It has an implementation so if the user 00114 redefines "computeUpdate" they aren't force to implement "getFlux" - 00115 which is only used by the default implementation of "computeUpdate". 00116 */ 00117 virtual void getFlux(FArrayBox& a_flux, 00118 const FArrayBox& a_WHalf, 00119 const int& a_dir, 00120 const Box& a_box); 00121 00122 /// Transform a_dW from primitive to characteristic variables 00123 /** 00124 On input, a_dW contains the increments of the primitive variables. On 00125 output, it contains the increments in the characteristic variables. 00126 00127 IMPORTANT NOTE: It is assumed that the characteristic analysis puts the 00128 smallest eigenvalue first, the largest eigenvalue last, and orders the 00129 characteristic variables accordingly. 00130 */ 00131 virtual void charAnalysis(FArrayBox& a_dW, 00132 const FArrayBox& a_W, 00133 const int& a_dir, 00134 const Box& a_box) = 0; 00135 00136 /// Transform a_dW from characteristic to primitive variables 00137 /** 00138 On input, a_dW contains the increments of the characteristic variables. 00139 On output, it contains the increments in the primitive variables. 00140 00141 IMPORTANT NOTE: It is assumed that the characteristic analysis puts the 00142 smallest eigenvalue first, the largest eigenvalue last, and orders the 00143 characteristic variables accordingly. 00144 */ 00145 virtual void charSynthesis(FArrayBox& a_dW, 00146 const FArrayBox& a_W, 00147 const int& a_dir, 00148 const Box& a_box) = 0; 00149 00150 /// Compute the characteristic values (eigenvalues) 00151 /** 00152 Compute the characteristic values (eigenvalues). 00153 00154 IMPORTANT NOTE: It is assumed that the characteristic analysis puts the 00155 smallest eigenvalue first, the largest eigenvalue last, and orders the 00156 characteristic variables accordingly. 00157 */ 00158 virtual void charValues(FArrayBox& a_lambda, 00159 const FArrayBox& a_W, 00160 const int& a_dir, 00161 const Box& a_box) = 0; 00162 00163 /// Add to (increment) the source terms given the current state 00164 /** 00165 On input, a_S contains the current source terms. On output, a_S has 00166 had any additional source terms (based on the current state, a_W) 00167 added to it. This should all be done on the region defined by a_box. 00168 */ 00169 virtual void incrementSource(FArrayBox& a_S, 00170 const FArrayBox& a_W, 00171 const Box& a_box) = 0; 00172 00173 /// Compute the solution to the Riemann problem. 00174 /** 00175 Given input left and right states in a direction, a_dir, compute a 00176 Riemann problem and generate fluxes at the faces within a_box. 00177 */ 00178 virtual void riemann(/// face-centered solution to Riemann problem 00179 FArrayBox& a_WStar, 00180 /// left state, on cells to left of each face 00181 const FArrayBox& a_WLeft, 00182 /// right state, on cells to right of each face 00183 const FArrayBox& a_WRight, 00184 /// state on cells, used to set boundary conditions 00185 const FArrayBox& a_W, 00186 /// current time 00187 const Real& a_time, 00188 /// direction of faces 00189 const int& a_dir, 00190 /// face-centered box on which to set a_WStar 00191 const Box& a_box) = 0; 00192 00193 /// Post-normal predictor calculation. 00194 /** 00195 Add increment to normal predictor, e.g. to account for source terms due to 00196 spatially-varying coefficients, to bound primitive variable ranges. 00197 */ 00198 virtual void postNormalPred(FArrayBox& a_dWMinus, 00199 FArrayBox& a_dWPlus, 00200 const FArrayBox& a_W, 00201 const Real& a_dt, 00202 const Real& a_dx, 00203 const int& a_dir, 00204 const Box& a_box) = 0; 00205 00206 /// Compute the quasilinear update A*dW/dx. 00207 /** 00208 */ 00209 virtual void quasilinearUpdate(FArrayBox& a_AdWdx, 00210 const FArrayBox& a_wHalf, 00211 const FArrayBox& a_W, 00212 const Real& a_scale, 00213 const int& a_dir, 00214 const Box& a_box) = 0; 00215 00216 /// Compute primitive variables from conserved variables. 00217 /** 00218 */ 00219 virtual void consToPrim(FArrayBox& a_W, 00220 const FArrayBox& a_U, 00221 const Box& a_box) = 0; 00222 00223 /// Compute the artificial viscosity contribution to the flux 00224 /** 00225 Compute the artificial viscosity contribution to the flux. This has 00226 a default implementation but this can be overridded as needed. 00227 */ 00228 virtual void artVisc(FArrayBox& a_F, 00229 const FArrayBox& a_U, 00230 const Real& a_artificialViscosity, 00231 const Real& a_currentTime, 00232 const int& a_dir, 00233 const Box& a_box); 00234 00235 #ifdef CH_USE_HDF5 00236 virtual void expressions(HDF5HeaderData& a_holder) const 00237 { 00238 } 00239 #endif 00240 00241 /** 00242 \name Access functions 00243 */ 00244 /**@{*/ 00245 00246 /// Number of conserved variables 00247 /** 00248 Return the number of conserved variables. 00249 */ 00250 virtual int numConserved() = 0; 00251 00252 /// Names of the conserved variables 00253 /** 00254 Return the names of the conserved variables. 00255 */ 00256 virtual Vector<string> stateNames() = 0; 00257 00258 /// Returns true if 4th-order artificial viscosity is defined. 00259 /** 00260 */ 00261 virtual bool fourthOrderArtificialViscosityIsDefined() const; 00262 00263 /// Defines fourth-order artifical viscosity strong shock threshold. 00264 /** 00265 */ 00266 virtual void setFourthOrderArtificialViscosityParameter(const Real& M0sq); 00267 /// Returns fourth-order artifical viscosity strong shock threshold. 00268 /** 00269 */ 00270 virtual Real getFourthOrderArtificialViscosityParameter() const; 00271 00272 /// Number of flux variables 00273 /** 00274 Return the number of flux variables. This can be greater than the number 00275 of conserved variables if addition fluxes/face-centered quantities are 00276 computed. 00277 */ 00278 virtual int numFluxes() = 0; 00279 00280 /// Is the object completely defined 00281 /** 00282 Return true if the object is completely defined. 00283 */ 00284 virtual bool isDefined() const; 00285 00286 /// Number of primitive variables 00287 /** 00288 Return the number of primitive variables. This may be greater than the 00289 number of conserved variables if derived/redundant quantities are also 00290 stored for convenience. 00291 */ 00292 virtual int numPrimitives() = 0; 00293 00294 /// Interval within the primitive variables corresponding to the velocities 00295 /** 00296 Return the interval of component indices within the primitive variable 00297 of the velocities. Used for slope flattening (slope computation) and 00298 computing the divergence of the velocity (artificial viscosity). 00299 */ 00300 virtual Interval velocityInterval() = 0; 00301 00302 /// Component index within the primitive variables of the pressure 00303 /** 00304 Return the component index withn the primitive variables for the 00305 pressure. Used for slope flattening (slope computation). 00306 */ 00307 virtual int pressureIndex() = 0; 00308 00309 /// Used to limit the absolute value of a "pressure" difference 00310 /** 00311 Return a value that is used by slope flattening to limit (away from 00312 zero) the absolute value of a slope in the pressureIndex() component 00313 (slope computation). 00314 */ 00315 virtual Real smallPressure() = 0; 00316 00317 /// Component index within the primitive variables of the bulk modulus 00318 /** 00319 Return the component index withn the primitive variables for the 00320 bulk modulus. Used for slope flattening (slope computation) used 00321 as a normalization to measure shock strength. 00322 */ 00323 virtual int bulkModulusIndex() = 0; 00324 00325 /// Component index within the primitive variables of the density. 00326 /** 00327 Return the component index within the primitive variables for the 00328 density. Used for fourth-order accurate artificial viscosity. 00329 */ 00330 virtual int densityIndex(); 00331 00332 /**@}*/ 00333 00334 protected: 00335 // Has this object been defined 00336 bool m_isDefined; 00337 00338 // Problem domain and grid spacing 00339 ProblemDomain m_domain; 00340 Real m_dx; 00341 00342 // Object containing various methods for the Godunov computation 00343 GodunovUtilities m_util; 00344 00345 // Flag to use fourth-order accurate artificial viscosity, strong shock 00346 // threshhold M0sq (looks like the shock Mach number squared). 00347 bool m_useFourthOrderArtificialViscosity; 00348 Real m_M0sq; 00349 00350 // Initial and boundary condition object and has it been set 00351 PhysIBC* m_bc; 00352 bool m_isBCSet; 00353 00354 00355 private: 00356 00357 // Disallowed for all the usual reasons 00358 void operator=(const GodunovPhysics&); 00359 GodunovPhysics(const GodunovPhysics&); 00360 }; 00361 00362 #include "NamespaceFooter.H" 00363 #endif