Chombo + EB  3.0
GodunovPhysics.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 _GODUNOVPHYSICS_H_
12 #define _GODUNOVPHYSICS_H_
13 
14 #include <string>
15 using std::string;
16 
17 #include "Box.H"
18 #include "IntVectSet.H"
19 #include "Vector.H"
20 #include "PhysIBC.H"
21 #include "FluxBox.H"
22 
23 #include "GodunovUtilities.H"
24 #include "NamespaceHeader.H"
25 
26 class HDF5HeaderData;
27 
28 ///
29 /**
30  The base class GodunovPhysics provides the physics-dependent components
31  for a higher-order Godunov method for a single patch: characteristic
32  analysis, Riemann solver, quasilinear update, conservative update,
33  and transformations between conserved, primitive, and flux variables.
34  This class is essentially pure virtual; i.e., all of its member functions
35  are virtual; and the ones that have default implementations are ones
36  that are optionally defined; i.e., the default definition is to send
37  an error message. Physics-dependent versions of this class that are
38  required in real applications are derived from this class by inheritance.
39 */
41 {
42 public:
43  /// Constructor
44  /**
45  */
47 
48  /// Get the initial and boundary condition object
49  /**
50  */
51  PhysIBC* getPhysIBC() const;
52 
53  /// Set the initial and boundary condition object
54  /**
55  */
56  void setPhysIBC(PhysIBC* a_bc);
57 
58  /// Destructor
59  /**
60  */
61  virtual ~GodunovPhysics();
62 
63  /// Define the object
64  /**
65  */
66  virtual void define(const ProblemDomain& a_domain,
67  const Real& a_dx);
68 
69  /// Set the current box (default implementation - do nothing)
70  /**
71  */
72  virtual void setCurrentBox(const Box& a_currentBox);
73 
74  /// Compute the maximum wave speed.
75  /**
76  */
77  virtual Real getMaxWaveSpeed(const FArrayBox& a_U,
78  const Box& a_box) = 0;
79 
80  /// Compute the maximum wave speed
81  /**
82  */
83  virtual void soundSpeed(FArrayBox& a_speed,
84  const FArrayBox& a_U,
85  const Box& a_box);
86 
87  /// Object factory for this class
88  /**
89  */
90  virtual GodunovPhysics* new_godunovPhysics() const = 0;
91 
92  /// Compute the increment in the conserved variables from face variables.
93  /**
94  Compute dU = dt*dUdt, the change in the conserved variables over
95  the time step. The fluxes are returned are suitable for use in refluxing.
96  This has a default implementation but can be redefined as needed.
97  */
98  virtual void computeUpdate(FArrayBox& a_dU,
99  FluxBox& a_F,
100  const FArrayBox& a_U,
101  const FluxBox& a_WHalf,
102  const bool& a_useArtificialViscosity,
103  const Real& a_artificialViscosity,
104  const Real& a_currentTime,
105  const Real& a_dx,
106  const Real& a_dt,
107  const Box& a_box);
108 
109  /// Compute the fluxes from primitive variable values on a face
110  /**
111  This has a default implementation which throws an error. The method is
112  here so that the default implementation of "computeUpdate" can use it
113  and the user can supply it. It has an implementation so if the user
114  redefines "computeUpdate" they aren't force to implement "getFlux" -
115  which is only used by the default implementation of "computeUpdate".
116  */
117  virtual void getFlux(FArrayBox& a_flux,
118  const FArrayBox& a_WHalf,
119  const int& a_dir,
120  const Box& a_box);
121 
122  /// Transform a_dW from primitive to characteristic variables
123  /**
124  On input, a_dW contains the increments of the primitive variables. On
125  output, it contains the increments in the characteristic variables.
126 
127  IMPORTANT NOTE: It is assumed that the characteristic analysis puts the
128  smallest eigenvalue first, the largest eigenvalue last, and orders the
129  characteristic variables accordingly.
130  */
131  virtual void charAnalysis(FArrayBox& a_dW,
132  const FArrayBox& a_W,
133  const int& a_dir,
134  const Box& a_box) = 0;
135 
136  /// Transform a_dW from characteristic to primitive variables
137  /**
138  On input, a_dW contains the increments of the characteristic variables.
139  On output, it contains the increments in the primitive variables.
140 
141  IMPORTANT NOTE: It is assumed that the characteristic analysis puts the
142  smallest eigenvalue first, the largest eigenvalue last, and orders the
143  characteristic variables accordingly.
144  */
145  virtual void charSynthesis(FArrayBox& a_dW,
146  const FArrayBox& a_W,
147  const int& a_dir,
148  const Box& a_box) = 0;
149 
150  /// Compute the characteristic values (eigenvalues)
151  /**
152  Compute the characteristic values (eigenvalues).
153 
154  IMPORTANT NOTE: It is assumed that the characteristic analysis puts the
155  smallest eigenvalue first, the largest eigenvalue last, and orders the
156  characteristic variables accordingly.
157  */
158  virtual void charValues(FArrayBox& a_lambda,
159  const FArrayBox& a_W,
160  const int& a_dir,
161  const Box& a_box) = 0;
162 
163  /// Add to (increment) the source terms given the current state
164  /**
165  On input, a_S contains the current source terms. On output, a_S has
166  had any additional source terms (based on the current state, a_W)
167  added to it. This should all be done on the region defined by a_box.
168  */
169  virtual void incrementSource(FArrayBox& a_S,
170  const FArrayBox& a_W,
171  const Box& a_box) = 0;
172 
173  /// Compute the solution to the Riemann problem.
174  /**
175  Given input left and right states in a direction, a_dir, compute a
176  Riemann problem and generate fluxes at the faces within a_box.
177  */
178  virtual void riemann(/// face-centered solution to Riemann problem
179  FArrayBox& a_WStar,
180  /// left state, on cells to left of each face
181  const FArrayBox& a_WLeft,
182  /// right state, on cells to right of each face
183  const FArrayBox& a_WRight,
184  /// state on cells, used to set boundary conditions
185  const FArrayBox& a_W,
186  /// current time
187  const Real& a_time,
188  /// direction of faces
189  const int& a_dir,
190  /// face-centered box on which to set a_WStar
191  const Box& a_box) = 0;
192 
193  /// Post-normal predictor calculation.
194  /**
195  Add increment to normal predictor, e.g. to account for source terms due to
196  spatially-varying coefficients, to bound primitive variable ranges.
197  */
198  virtual void postNormalPred(FArrayBox& a_dWMinus,
199  FArrayBox& a_dWPlus,
200  const FArrayBox& a_W,
201  const Real& a_dt,
202  const Real& a_dx,
203  const int& a_dir,
204  const Box& a_box) = 0;
205 
206  /// Compute the quasilinear update A*dW/dx.
207  /**
208  */
209  virtual void quasilinearUpdate(FArrayBox& a_AdWdx,
210  const FArrayBox& a_wHalf,
211  const FArrayBox& a_W,
212  const Real& a_scale,
213  const int& a_dir,
214  const Box& a_box) = 0;
215 
216  /// Compute primitive variables from conserved variables.
217  /**
218  */
219  virtual void consToPrim(FArrayBox& a_W,
220  const FArrayBox& a_U,
221  const Box& a_box) = 0;
222 
223  /// Compute the artificial viscosity contribution to the flux
224  /**
225  Compute the artificial viscosity contribution to the flux. This has
226  a default implementation but this can be overridded as needed.
227  */
228  virtual void artVisc(FArrayBox& a_F,
229  const FArrayBox& a_U,
230  const Real& a_artificialViscosity,
231  const Real& a_currentTime,
232  const int& a_dir,
233  const Box& a_box);
234 
235 #ifdef CH_USE_HDF5
236  virtual void expressions(HDF5HeaderData& a_holder) const
237  {
238  }
239 #endif
240 
241  /**
242  \name Access functions
243  */
244  /**@{*/
245 
246  /// Number of conserved variables
247  /**
248  Return the number of conserved variables.
249  */
250  virtual int numConserved() = 0;
251 
252  /// Names of the conserved variables
253  /**
254  Return the names of the conserved variables.
255  */
256  virtual Vector<string> stateNames() = 0;
257 
258  /// Returns true if 4th-order artificial viscosity is defined.
259  /**
260  */
261  virtual bool fourthOrderArtificialViscosityIsDefined() const;
262 
263  /// Defines fourth-order artifical viscosity strong shock threshold.
264  /**
265  */
266  virtual void setFourthOrderArtificialViscosityParameter(const Real& M0sq);
267  /// Returns fourth-order artifical viscosity strong shock threshold.
268  /**
269  */
271 
272  /// Number of flux variables
273  /**
274  Return the number of flux variables. This can be greater than the number
275  of conserved variables if addition fluxes/face-centered quantities are
276  computed.
277  */
278  virtual int numFluxes() = 0;
279 
280  /// Is the object completely defined
281  /**
282  Return true if the object is completely defined.
283  */
284  virtual bool isDefined() const;
285 
286  /// Number of primitive variables
287  /**
288  Return the number of primitive variables. This may be greater than the
289  number of conserved variables if derived/redundant quantities are also
290  stored for convenience.
291  */
292  virtual int numPrimitives() = 0;
293 
294  /// Interval within the primitive variables corresponding to the velocities
295  /**
296  Return the interval of component indices within the primitive variable
297  of the velocities. Used for slope flattening (slope computation) and
298  computing the divergence of the velocity (artificial viscosity).
299  */
300  virtual Interval velocityInterval() = 0;
301 
302  /// Component index within the primitive variables of the pressure
303  /**
304  Return the component index withn the primitive variables for the
305  pressure. Used for slope flattening (slope computation).
306  */
307  virtual int pressureIndex() = 0;
308 
309  /// Used to limit the absolute value of a "pressure" difference
310  /**
311  Return a value that is used by slope flattening to limit (away from
312  zero) the absolute value of a slope in the pressureIndex() component
313  (slope computation).
314  */
315  virtual Real smallPressure() = 0;
316 
317  /// Component index within the primitive variables of the bulk modulus
318  /**
319  Return the component index withn the primitive variables for the
320  bulk modulus. Used for slope flattening (slope computation) used
321  as a normalization to measure shock strength.
322  */
323  virtual int bulkModulusIndex() = 0;
324 
325  /// Component index within the primitive variables of the density.
326  /**
327  Return the component index within the primitive variables for the
328  density. Used for fourth-order accurate artificial viscosity.
329  */
330  virtual int densityIndex();
331 
332  /**@}*/
333 
334 protected:
335  // Has this object been defined
337 
338  // Problem domain and grid spacing
341 
342  // Object containing various methods for the Godunov computation
344 
345  // Flag to use fourth-order accurate artificial viscosity, strong shock
346  // threshhold M0sq (looks like the shock Mach number squared).
349 
350  // Initial and boundary condition object and has it been set
352  bool m_isBCSet;
353 
354 
355 private:
356 
357  // Disallowed for all the usual reasons
358  void operator=(const GodunovPhysics&);
360 };
361 
362 #include "NamespaceFooter.H"
363 #endif
virtual void riemann(FArrayBox &a_WStar, const FArrayBox &a_WLeft, const FArrayBox &a_WRight, const FArrayBox &a_W, const Real &a_time, const int &a_dir, const Box &a_box)=0
Compute the solution to the Riemann problem.
virtual void charValues(FArrayBox &a_lambda, const FArrayBox &a_W, const int &a_dir, const Box &a_box)=0
Compute the characteristic values (eigenvalues)
Real m_M0sq
Definition: GodunovPhysics.H:348
virtual void quasilinearUpdate(FArrayBox &a_AdWdx, const FArrayBox &a_wHalf, const FArrayBox &a_W, const Real &a_scale, const int &a_dir, const Box &a_box)=0
Compute the quasilinear update A*dW/dx.
virtual void setCurrentBox(const Box &a_currentBox)
Set the current box (default implementation - do nothing)
virtual void define(const ProblemDomain &a_domain, const Real &a_dx)
Define the object.
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:130
data to be added to HDF5 files.
Definition: CH_HDF5.H:490
Definition: GodunovUtilities.H:31
bool m_isBCSet
Definition: GodunovPhysics.H:352
PhysIBC * getPhysIBC() const
Get the initial and boundary condition object.
virtual void incrementSource(FArrayBox &a_S, const FArrayBox &a_W, const Box &a_box)=0
Add to (increment) the source terms given the current state.
GodunovPhysics()
Constructor.
virtual Vector< string > stateNames()=0
Names of the conserved variables.
virtual void postNormalPred(FArrayBox &a_dWMinus, FArrayBox &a_dWPlus, const FArrayBox &a_W, const Real &a_dt, const Real &a_dx, const int &a_dir, const Box &a_box)=0
Post-normal predictor calculation.
virtual Interval velocityInterval()=0
Interval within the primitive variables corresponding to the velocities.
PhysIBC * m_bc
Definition: GodunovPhysics.H:351
A FArrayBox-like container for face-centered fluxes.
Definition: FluxBox.H:22
Structure for passing component ranges in code.
Definition: Interval.H:23
virtual void charAnalysis(FArrayBox &a_dW, const FArrayBox &a_W, const int &a_dir, const Box &a_box)=0
Transform a_dW from primitive to characteristic variables.
double Real
Definition: REAL.H:33
virtual int numConserved()=0
Number of conserved variables.
virtual void consToPrim(FArrayBox &a_W, const FArrayBox &a_U, const Box &a_box)=0
Compute primitive variables from conserved variables.
void setPhysIBC(PhysIBC *a_bc)
Set the initial and boundary condition object.
virtual void expressions(HDF5HeaderData &a_holder) const
Definition: GodunovPhysics.H:236
Definition: GodunovPhysics.H:40
virtual int numPrimitives()=0
Number of primitive variables.
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
virtual int pressureIndex()=0
Component index within the primitive variables of the pressure.
Real m_dx
Definition: GodunovPhysics.H:340
virtual void soundSpeed(FArrayBox &a_speed, const FArrayBox &a_U, const Box &a_box)
Compute the maximum wave speed.
virtual void getFlux(FArrayBox &a_flux, const FArrayBox &a_WHalf, const int &a_dir, const Box &a_box)
Compute the fluxes from primitive variable values on a face.
GodunovUtilities m_util
Definition: GodunovPhysics.H:343
bool m_isDefined
Definition: GodunovPhysics.H:336
virtual bool fourthOrderArtificialViscosityIsDefined() const
Returns true if 4th-order artificial viscosity is defined.
Definition: FArrayBox.H:44
virtual Real getFourthOrderArtificialViscosityParameter() const
Returns fourth-order artifical viscosity strong shock threshold.
virtual int numFluxes()=0
Number of flux variables.
virtual GodunovPhysics * new_godunovPhysics() const =0
Object factory for this class.
virtual void computeUpdate(FArrayBox &a_dU, FluxBox &a_F, const FArrayBox &a_U, const FluxBox &a_WHalf, const bool &a_useArtificialViscosity, const Real &a_artificialViscosity, const Real &a_currentTime, const Real &a_dx, const Real &a_dt, const Box &a_box)
Compute the increment in the conserved variables from face variables.
Physical/domain initial and boundary conditions.
Definition: PhysIBC.H:33
virtual void setFourthOrderArtificialViscosityParameter(const Real &M0sq)
Defines fourth-order artifical viscosity strong shock threshold.
virtual void artVisc(FArrayBox &a_F, const FArrayBox &a_U, const Real &a_artificialViscosity, const Real &a_currentTime, const int &a_dir, const Box &a_box)
Compute the artificial viscosity contribution to the flux.
ProblemDomain m_domain
Definition: GodunovPhysics.H:339
void operator=(const GodunovPhysics &)
virtual Real smallPressure()=0
Used to limit the absolute value of a "pressure" difference.
virtual Real getMaxWaveSpeed(const FArrayBox &a_U, const Box &a_box)=0
Compute the maximum wave speed.
virtual int bulkModulusIndex()=0
Component index within the primitive variables of the bulk modulus.
bool m_useFourthOrderArtificialViscosity
Definition: GodunovPhysics.H:347
virtual bool isDefined() const
Is the object completely defined.
virtual int densityIndex()
Component index within the primitive variables of the density.
virtual void charSynthesis(FArrayBox &a_dW, const FArrayBox &a_W, const int &a_dir, const Box &a_box)=0
Transform a_dW from characteristic to primitive variables.
virtual ~GodunovPhysics()
Destructor.