SWE
/import/home/rettenbs/src/SWE/submodules/swe_solvers/src/solver/FWave.hpp
00001 
00053 #ifndef FWAVE_HPP_
00054 #define FWAVE_HPP_
00055 
00056 #include <cassert>
00057 #include <cmath>
00058 #include <algorithm>
00059 #include "WavePropagation.hpp"
00060 
00061 namespace solver {
00062   template <typename T> class FWave;
00063 }
00064 
00070 template <typename T> class solver::FWave: public WavePropagation<T> {
00071   private: //explicit for unit tests
00072     //use nondependent names (template base class)
00073     using solver::WavePropagation<T>::zeroTol;
00074     using solver::WavePropagation<T>::g;
00075     using solver::WavePropagation<T>::dryTol;
00076 
00077     using solver::WavePropagation<T>::hLeft;
00078     using solver::WavePropagation<T>::hRight;
00079     using solver::WavePropagation<T>::huLeft;
00080     using solver::WavePropagation<T>::huRight;
00081     using solver::WavePropagation<T>::bLeft;
00082     using solver::WavePropagation<T>::bRight;
00083     using solver::WavePropagation<T>::uLeft;
00084     using solver::WavePropagation<T>::uRight;
00085 
00086     using solver::WavePropagation<T>::wetDryState;
00087     using solver::WavePropagation<T>::DryDry;
00088     using solver::WavePropagation<T>::WetWet;
00089     using solver::WavePropagation<T>::WetDryWall;
00090     using solver::WavePropagation<T>::DryWetWall;
00091 
00092     using solver::WavePropagation<T>::storeParameters;
00093 
00094   public:
00102     FWave( T i_dryTolerance =  (T) 0.01,
00103            T i_gravity =       (T) 9.81,
00104            T i_zeroTolerance = (T) 0.000000001 ):
00105              WavePropagation<T>( i_dryTolerance, i_gravity, i_zeroTolerance ) {
00106 
00107 #ifndef NDEBUG
00108 #ifndef SUPPRESS_SOLVER_DEBUG_OUTPUT
00109       //print some information about the used solver
00110       std::cout << "  *** solver::FWave created" << std::endl
00111                 << "    zeroTolerance=" << zeroTol << std::endl
00112                 << "    gravity=" << g << std::endl
00113                 << "    dryTolerance=" << dryTol << std::endl
00114                 << "  ***\n\n";
00115 #endif
00116 #endif
00117     }
00118 
00141     void computeNetUpdates ( const T &i_hLeft,  const T &i_hRight,
00142                              const T &i_huLeft, const T &i_huRight,
00143                              const T &i_bLeft,  const T &i_bRight,
00144 
00145                              T &o_hUpdateLeft,
00146                              T &o_hUpdateRight,
00147                              T &o_huUpdateLeft,
00148                              T &o_huUpdateRight,
00149                              T &o_maxWaveSpeed ) {
00150 
00151       //set speeds to zero (will be determined later)
00152       uLeft = uRight = 0;
00153 
00154       //reset the maximum wave speed
00155       o_maxWaveSpeed = 0;
00156 
00158       T waveSpeeds[2];
00159 
00160       //store parameters to member variables
00161       storeParameters( i_hLeft, i_hRight,
00162                        i_huLeft, i_huRight,
00163                        i_bLeft, i_bRight );
00164 
00165       //determine the wet/dry state and compute local variables correspondingly
00166       determineWetDryState();
00167 
00168       //zero updates and return in the case of dry cells
00169       if(wetDryState == DryDry) {
00170         o_hUpdateLeft = o_hUpdateRight = o_huUpdateLeft = o_huUpdateRight  = (T)0.;
00171         return;
00172       }
00173 
00174       //compute the wave speeds
00175       computeWaveSpeeds(waveSpeeds);
00176 
00177       //use the wave speeds to compute the net-updates
00178       computeNetUpdates_WithWaveSpeeds( waveSpeeds,
00179                                         o_hUpdateLeft, o_hUpdateRight,
00180                                         o_huUpdateLeft, o_huUpdateRight,
00181                                         o_maxWaveSpeed );
00182 
00183       //zero ghost updates (wall boundary)
00184       if(wetDryState == WetDryWall) {
00185         o_hUpdateRight = 0;
00186         o_huUpdateRight = 0;
00187       }
00188       else if(wetDryState == DryWetWall) {
00189         o_hUpdateLeft = 0;
00190         o_huUpdateLeft = 0;
00191       }
00192     }
00193 
00219     void computeNetUpdatesHybrid ( const T &i_hLeft,  const T &i_hRight,
00220                                    const T &i_huLeft, const T &i_huRight,
00221                                    const T &i_bLeft,  const T &i_bRight,
00222                                    const T &i_uLeft,  const T &i_uRight,
00223 
00224                                    const T i_waveSpeeds[2],
00225 
00226                                    T &o_hUpdateLeft,
00227                                    T &o_hUpdateRight,
00228                                    T &o_huUpdateLeft,
00229                                    T &o_huUpdateRight,
00230                                    T &o_maxWaveSpeed ) {
00231       //store parameters to member variables
00232       storeParameters( i_hLeft,  i_hRight,
00233                        i_huLeft, i_huRight,
00234                        i_bLeft,  i_bRight,
00235                        i_uLeft,  i_uRight );
00236 
00237       computeNetUpdates_WithWaveSpeeds( i_waveSpeeds,
00238                                                                                   o_hUpdateLeft, o_hUpdateRight,
00239                                                                             o_huUpdateLeft, o_huUpdateRight,
00240                                                                                   o_maxWaveSpeed );
00241     }
00242 
00243   private:
00256     void computeNetUpdates_WithWaveSpeeds ( const T i_waveSpeeds[2],
00257 
00258                                                                                   T &o_hUpdateLeft,
00259                                                                                         T &o_hUpdateRight,
00260                                                                       T &o_huUpdateLeft,
00261                                                                                         T &o_huUpdateRight,
00262                                                                                   T &o_maxWaveSpeed ) {
00263       //reset net updates
00264       o_hUpdateLeft = o_hUpdateRight = o_huUpdateLeft = o_huUpdateRight  = (T)0.;
00265 
00267       T fWaves[2][2];
00268 
00269       //compute the decomposition into f-waves
00270       computeWaveDecomposition( i_waveSpeeds, fWaves );
00271 
00272       //compute the net-updates
00273       //1st wave family
00274       if(i_waveSpeeds[0] < -zeroTol) { //left going
00275         o_hUpdateLeft +=  fWaves[0][0];
00276         o_huUpdateLeft += fWaves[0][1];
00277       }
00278       else if(i_waveSpeeds[0] > zeroTol) { //right going
00279         o_hUpdateRight +=  fWaves[0][0];
00280         o_huUpdateRight += fWaves[0][1];
00281       }
00282       else { //split waves
00283         o_hUpdateLeft +=   (T).5*fWaves[0][0];
00284         o_huUpdateLeft +=  (T).5*fWaves[0][1];
00285         o_hUpdateRight +=  (T).5*fWaves[0][0];
00286         o_huUpdateRight += (T).5*fWaves[0][1];
00287       }
00288 
00289       //2nd wave family
00290       if(i_waveSpeeds[1] < -zeroTol) { //left going
00291           o_hUpdateLeft +=  fWaves[1][0];
00292           o_huUpdateLeft += fWaves[1][1];
00293       }
00294       else if(i_waveSpeeds[1] > zeroTol) { //right going
00295           o_hUpdateRight += fWaves[1][0];
00296           o_huUpdateRight += fWaves[1][1];
00297       }
00298       else { //split waves
00299         o_hUpdateLeft +=   (T).5*fWaves[1][0];
00300         o_huUpdateLeft +=  (T).5*fWaves[1][1];
00301         o_hUpdateRight +=  (T).5*fWaves[1][0];
00302         o_huUpdateRight += (T).5*fWaves[1][1];
00303       }
00304 
00305       //compute maximum wave speed (-> CFL-condition)
00306       o_maxWaveSpeed = std::max( std::fabs(i_waveSpeeds[0]) , std::fabs(i_waveSpeeds[1]) );
00307     }
00308 
00312     void determineWetDryState() {
00313       //determine the wet/dry state
00314       if(hLeft < dryTol && hRight < dryTol) { //both cells are dry
00315         wetDryState = DryDry;
00316       }
00317       else if(hLeft < dryTol) { // left cell dry, right cell wet
00318         uRight = huRight/hRight;
00319 
00320         //Set wall boundary conditions.
00321         //This is not correct in the case of inundation problems.
00322         hLeft = hRight;
00323         bLeft = bRight;
00324         huLeft = -huRight;
00325         uLeft = -uRight;
00326         wetDryState = DryWetWall;
00327       }
00328       else if(hRight < dryTol) { //left cell wet, right cell dry
00329         uLeft = huLeft/hLeft;
00330 
00331         //Set wall boundary conditions.
00332         //This is not correct in the case of inundation problems.
00333           hRight = hLeft;
00334           bRight = bLeft;
00335           huRight = -huLeft;
00336           uLeft = -uRight;
00337           wetDryState = WetDryWall;
00338       }
00339       else { //both cells wet
00340         uLeft = huLeft/hLeft;
00341         uRight = huRight/hRight;
00342 
00343         wetDryState = WetWet;
00344       }
00345     }
00346 
00353     void computeWaveDecomposition( const T i_waveSpeeds[2],
00354                                          T o_fWaves[2][2]
00355                                  ) const {
00356       //TODO: Update documentation.
00357       //Eigenvalues***********************************************************************************************
00358       //Computed somewhere before.
00359       //An option would be to use the char. speeds:
00360       //
00361       //lambda^1 = u_{i-1} - sqrt(g*h_{i-1})
00362       //lambda^2 = u_i     + sqrt(g*h_i)
00363       //Matrix of right eigenvectors******************************************************************************
00364       //     1                              1
00365       // R =
00366       //     u_{i-1} - sqrt(g * h_{i-1})    u_i + sqrt(g * h_i)
00367       //**********************************************************************************************************
00368       //                                                                      u_i + sqrt(g * h_i)              -1
00369       // R^{-1} = 1 / (u_i - sqrt(g * h_i) - u_{i-1} + sqrt(g * h_{i-1}) *
00370       //                                                                   -( u_{i-1} - sqrt(g * h_{i-1}) )     1
00371       //**********************************************************************************************************
00372       //                hu
00373       // f(q) =
00374       //         hu^2 + 1/2 g * h^2
00375       //
00376       //**********************************************************************************************************
00377       //                                    0
00378       //  \delta x \Psi =
00379       //                  -g * 1/2 * (h_i + h_{i-1}) * (b_i - b_{i+1})
00380       //**********************************************************************************************************
00381       // beta = R^{-1} * (f(Q_i) - f(Q_{i-1}) - \delta x \Psi)
00382       //**********************************************************************************************************
00383 
00384       //assert: wave speed of the 1st wave family should be less than the speed of the 2nd wave family.
00385       assert( i_waveSpeeds[0] < i_waveSpeeds[1] );
00386 
00387       T lambdaDif = i_waveSpeeds[1] - i_waveSpeeds[0];
00388 
00389       //assert: no division by zero
00390       assert(std::abs(lambdaDif) > zeroTol);
00391 
00392       //compute the inverse matrix R^{-1}
00393       T Rinv[2][2];
00394 
00395       T oneDivLambdaDif = (T)1. /  lambdaDif;
00396       Rinv[0][0] = oneDivLambdaDif *  i_waveSpeeds[1];
00397       Rinv[0][1] = -oneDivLambdaDif;
00398 
00399       Rinv[1][0] = oneDivLambdaDif * -i_waveSpeeds[0];
00400       Rinv[1][1] = oneDivLambdaDif;
00401 
00402       //right hand side
00403       T fDif[2];
00404 
00405       //calculate modified (bathymetry!) flux difference
00406       // f(Q_i) - f(Q_{i-1})
00407       fDif[0] = huRight - huLeft;
00408       fDif[1] = huRight * uRight + (T).5 * g * hRight * hRight
00409               -(huLeft  * uLeft  + (T).5 * g * hLeft  * hLeft);
00410 
00411       // \delta x \Psi[2]
00412       T psi = -g * (T).5 * (hRight + hLeft)*(bRight - bLeft);
00413       fDif[1] -= psi;
00414 
00415       //solve linear equations
00416       T beta[2];
00417       beta[0] = Rinv[0][0] * fDif[0] + Rinv[0][1] * fDif[1];
00418       beta[1] = Rinv[1][0] * fDif[0] + Rinv[1][1] * fDif[1];
00419 
00420       //return f-waves
00421       o_fWaves[0][0] = beta[0];
00422       o_fWaves[0][1] = beta[0] * i_waveSpeeds[0];
00423 
00424       o_fWaves[1][0] = beta[1];
00425       o_fWaves[1][1] = beta[1] * i_waveSpeeds[1];
00426     }
00427 
00433     void computeWaveSpeeds( T o_waveSpeeds[2]) const {
00434       //compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i}
00435       T characteristicSpeeds[2];
00436       characteristicSpeeds[0] = uLeft - std::sqrt(g*hLeft);
00437       characteristicSpeeds[1] = uRight + std::sqrt(g*hRight);
00438 
00439       //compute "Roe speeds"
00440       T hRoe = (T).5 * (hRight + hLeft);
00441       T uRoe = uLeft * std::sqrt(hLeft) + uRight * std::sqrt(hRight);
00442       uRoe /= std::sqrt(hLeft) + std::sqrt(hRight);
00443 
00444       T roeSpeeds[2];
00445       roeSpeeds[0] = uRoe - std::sqrt(g*hRoe);
00446       roeSpeeds[1] = uRoe + std::sqrt(g*hRoe);
00447 
00448       //computer eindfeldt speeds
00449       T einfeldtSpeeds[2];
00450       einfeldtSpeeds[0] = std::min(characteristicSpeeds[0], roeSpeeds[0]);
00451       einfeldtSpeeds[1] = std::max(characteristicSpeeds[1], roeSpeeds[1]);
00452 
00453       //set wave speeds
00454       o_waveSpeeds[0] = einfeldtSpeeds[0];
00455       o_waveSpeeds[1] = einfeldtSpeeds[1];
00456     }
00457 };
00458 
00459 
00460 
00461 
00462 #endif /* FWAVE_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends