SWE
/import/home/rettenbs/src/SWE/submodules/swe_solvers/src/solver/Hybrid.hpp
00001 
00036 #ifndef HYBRID_HPP_
00037 #define HYBRID_HPP_
00038 
00039 #include "FWave.hpp"
00040 #include "AugRie.hpp"
00041 
00042 #ifndef NDEBUG
00043 #include <iostream>
00044 #endif
00045 
00046 namespace solver {
00047   template <typename T> class Hybrid;
00048 }
00049 
00056 template <typename T> class solver::Hybrid {
00057   private:
00059     const T g;
00061     const T dryTol;
00063     solver::FWave<T> fWaveSolver;
00065     solver::AugRie<T> augRieSolver;
00066 
00067 #ifndef NDEBUG
00068 
00069     long counterFWave;
00070     
00072     long counterAugRie;
00073 #endif
00074 
00075   public:
00085     Hybrid( T i_dryTolerance =                  (T) 0.01,
00086             T i_gravity =                       (T) 9.81,
00087             T i_newtonTolerance =               (T) 0.000001,
00088             int i_maxNumberOfNewtonIterations =     10,
00089             T zeroTolerance =                   (T) 0.000000001 ):
00090               g(i_gravity),
00091               dryTol(i_dryTolerance),
00092               fWaveSolver( i_dryTolerance,
00093                            i_gravity,
00094                            zeroTolerance ),
00095               augRieSolver( i_dryTolerance,
00096                             i_gravity,
00097                             i_newtonTolerance,
00098                             i_maxNumberOfNewtonIterations,
00099                             zeroTolerance ) {
00100 #ifndef NDEBUG
00101       //set the counters to zero.
00102       counterFWave = counterAugRie = 0;
00103 
00104 #ifndef SUPPRESS_SOLVER_DEBUG_OUTPUT
00105       //print some information about the used solver.
00106       std::cout << "  *** solver::Hybrid created" << std::endl
00107                 << "    zeroTolerance=" << zeroTolerance << std::endl
00108                 << "    gravity=" << i_gravity << std::endl
00109                 << "    dryTolerance=" << i_dryTolerance << std::endl
00110                 << "    newtonTolerance=" << i_newtonTolerance << std::endl
00111                 << "    maxNumberNumberOfNewtonIterations=" << i_maxNumberOfNewtonIterations << std::endl
00112                 << "  ***\n\n";
00113 #endif
00114 #endif
00115     };
00116 
00134     void computeNetUpdates ( const T &i_hLeft,  const T &i_hRight,
00135                              const T &i_huLeft, const T &i_huRight,
00136                              const T &i_bLeft,  const T &i_bRight,
00137 
00138                              T &o_hUpdateLeft,
00139                              T &o_hUpdateRight,
00140                              T &o_huUpdateLeft,
00141                              T &o_huUpdateRight,
00142                              T &o_maxWaveSpeed ) {
00143       if( useAugmentedRiemannSolver( i_hLeft, i_hRight,
00144                                      i_huLeft, i_huRight,
00145                                      i_bLeft, i_bRight,
00146                                      o_hUpdateLeft, o_hUpdateRight,
00147                                      o_huUpdateLeft, o_huUpdateRight,
00148                                      o_maxWaveSpeed ) == true ) {
00149 
00150         augRieSolver.computeNetUpdates( i_hLeft, i_hRight,
00151                                         i_huLeft, i_huRight,
00152                                         i_bLeft, i_bRight,
00153 
00154                                         o_hUpdateLeft, o_hUpdateRight,
00155                                         o_huUpdateLeft, o_huUpdateRight,
00156                                         o_maxWaveSpeed );
00157 
00158 #ifndef NDEBUG
00159         counterAugRie++;
00160 #endif
00161 
00162       }
00163 #ifndef NDEBUG
00164       else {
00165         counterFWave++;
00166       }
00167 #endif
00168     }
00169 
00190     bool useAugmentedRiemannSolver( const T &i_hLeft,  const T &i_hRight,
00191                                     const T &i_huLeft, const T &i_huRight,
00192                                     const T &i_bLeft,  const T &i_bRight,
00193 
00194                                     T &o_hUpdateLeft,
00195                                     T &o_hUpdateRight,
00196                                     T &o_huUpdateLeft,
00197                                     T &o_huUpdateRight,
00198                                     T &o_maxWaveSpeed ) {
00199       //use the augmented solver by default
00200       bool useAugRieSolver = true;
00201 
00202       if(i_hLeft > dryTol && i_hRight > dryTol) { //both cells wet?
00203 
00204         //compute the velocities
00205         T uLeft = i_huLeft / i_hLeft;
00206         T uRight = i_huRight / i_hRight;
00207 
00208         //definition of a "simple" Riemann-problem (GeoClaw)
00209         //TODO: Literature?
00210         if( std::fabs(i_hLeft - i_hRight) < (T)0.001*std::min(i_hLeft, i_hRight) && //"small" jump in height
00211             std::max(uLeft*uLeft, uRight*uRight)/std::max(g*i_hLeft, g*i_hRight) < (T)0.01 ) { //"small" jump in velocity
00212           useAugRieSolver = false;
00213 
00214           //compute simplified wave speeds (Geoclaw)
00215           //TODO: Literature?
00216           T waveSpeeds[2];
00217           waveSpeeds[0] = (uLeft + uRight)*(T).5 - std::sqrt(g*(i_hLeft + i_hRight)*(T).5);
00218           waveSpeeds[1] = (uLeft + uRight)*(T).5 + std::sqrt(g*(i_hLeft + i_hRight)*(T).5);
00219 
00220           //compute the f-waves
00221           fWaveSolver.computeNetUpdatesHybrid( i_hLeft, i_hRight,
00222                                                i_huLeft, i_huRight,
00223                                                i_bLeft, i_bRight,
00224                                                uLeft, uRight,
00225                                                waveSpeeds,
00226 
00227                                                o_hUpdateLeft, o_hUpdateRight,
00228                                                o_huUpdateLeft, o_huUpdateRight,
00229                                                o_maxWaveSpeed );
00230         }
00231       }
00232       //DryDry case: Don't use any solver and set everything to zero
00233       else if(i_hLeft < dryTol && i_hRight < dryTol) {
00234         o_hUpdateLeft = o_hUpdateRight = o_huUpdateLeft = o_huUpdateRight = o_maxWaveSpeed = 0;
00235         useAugRieSolver = false;
00236       }
00237 
00238       return useAugRieSolver;
00239     }
00240 
00241 #ifndef NDEBUG
00242 
00245     void printStats() {
00246       std::cout << "  *** solver::Hybrid printing statistics" << std::endl
00247                 << "    times the f-Wave solver was used: " << counterFWave << std::endl
00248                 << "    times the AugRie solver was used: " << counterAugRie << std::endl
00249                 << "  ***\n\n";
00250     }
00251 
00258     void getStats( long &o_counterFWave, long &o_counterAugRie ) {
00259       o_counterFWave = counterFWave;
00260       o_counterAugRie = counterAugRie;
00261     }
00262 #endif
00263 };
00264 
00265 #endif /* HYBRID_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends