SWE
/import/home/rettenbs/src/SWE/submodules/swe_solvers/src/solver/FWaveVec.hpp
00001 
00049 #ifndef FWAVEVEC_HPP_
00050 #define FWAVEVEC_HPP_
00051 
00052 #include <cmath>
00053 
00054 namespace solver
00055 {
00056 
00060 template<typename T>
00061 class FWaveVec
00062 {
00063 private:
00064         const T dryTol;
00065         const T gravity;
00066         const T zeroTol;
00067 
00068 public:
00069         FWaveVec(T i_dryTol = (T) 100,
00070                          T i_gravity = (T) 9.81,
00071                          T i_zeroTol = (T) 0.0000001)
00072                 : dryTol(i_dryTol),
00073                   gravity(i_gravity),
00074                   zeroTol(i_zeroTol)
00075         {
00076         }
00077 
00078         void computeNetUpdates ( T i_hLeft,  T i_hRight,
00079                              T i_huLeft, T i_huRight,
00080                              T i_bLeft,  T i_bRight,
00081 
00082                              T &o_hUpdateLeft,
00083                              T &o_hUpdateRight,
00084                              T &o_huUpdateLeft,
00085                              T &o_huUpdateRight,
00086                              T &o_maxWaveSpeed ) const
00087         {
00088                   // determine the wet dry state and corr. values, if necessary.
00089                   if( i_hLeft < dryTol && i_hRight < dryTol ) {
00090                       // Dry/Dry case
00091                       // Set dummy values such that the result is zero
00092                       i_hLeft = dryTol;
00093                       i_huLeft = 0.; i_bLeft = 0.;
00094                       i_hRight = dryTol;
00095                       i_huRight = 0.; i_bRight = 0.;
00096                   } else if ( i_hLeft < dryTol ) {
00097                       i_hLeft = i_hRight;
00098                       i_huLeft = -i_huRight;
00099                       i_bLeft = i_bRight;
00100                   } else if ( i_hRight < dryTol ) {
00101                       i_hRight = i_hLeft;
00102                       i_huRight = -i_huLeft;
00103                       i_bRight = i_bLeft;
00104                   }
00105 
00107                   T uLeft = i_huLeft / i_hLeft;
00109                   T uRight = i_huRight / i_hRight;
00110 
00112                   T waveSpeeds0 = 0., waveSpeeds1 = 0.;
00113 
00114                   //compute the wave speeds
00115                   fWaveComputeWaveSpeeds( i_hLeft, i_hRight,
00116                                           i_huLeft, i_huRight,
00117                                           uLeft, uRight,
00118                                           i_bLeft, i_bRight,
00119 
00120                                           waveSpeeds0, waveSpeeds1 );
00121 
00123                   T fWaves0 = 0., fWaves1 = 0.;
00124 
00125                   //compute the decomposition into f-waves
00126                   fWaveComputeWaveDecomposition( i_hLeft, i_hRight,
00127                                                  i_huLeft, i_huRight,
00128                                                  uLeft, uRight,
00129                                                  i_bLeft, i_bRight,
00130 
00131                                                  waveSpeeds0, waveSpeeds1,
00132                                                  fWaves0, fWaves1);
00133 
00134                   //compute the net-updates
00135                   T hUpdateLeft = 0.;
00136                   T hUpdateRight = 0.;
00137                   T huUpdateLeft = 0.;
00138                   T huUpdateRight = 0.;
00139 
00140                   //1st wave family
00141                   if(waveSpeeds0 < -zeroTol) { //left going
00142                     hUpdateLeft +=  fWaves0;
00143                     huUpdateLeft += fWaves0 * waveSpeeds0;
00144                   }
00145                   else if(waveSpeeds0 > zeroTol) { //right going
00146                     hUpdateRight +=  fWaves0;
00147                     huUpdateRight += fWaves0 * waveSpeeds0;
00148                   }
00149                   else { //split waves
00150                     hUpdateLeft +=   (T).5*fWaves0;
00151                     huUpdateLeft +=  (T).5*fWaves0 * waveSpeeds0;
00152                     hUpdateRight +=  (T).5*fWaves0;
00153                     huUpdateRight += (T).5*fWaves0 * waveSpeeds0;
00154                   }
00155 
00156                   //2nd wave family
00157                   if(waveSpeeds1 > zeroTol) { //right going
00158                     hUpdateRight +=  fWaves1;
00159                     huUpdateRight += fWaves1 * waveSpeeds1;
00160                   }
00161                   else if(waveSpeeds1 < -zeroTol) { //left going
00162                         hUpdateLeft +=  fWaves1;
00163                         huUpdateLeft += fWaves1 * waveSpeeds1;
00164                   }
00165                   else { //split waves
00166                     hUpdateLeft +=   (T).5*fWaves1;
00167                     huUpdateLeft +=  (T).5*fWaves1 * waveSpeeds1;
00168                     hUpdateRight +=  (T).5*fWaves1;
00169                     huUpdateRight += (T).5*fWaves1 * waveSpeeds1;
00170                   }
00171 
00172                   // Set output variables
00173                   o_hUpdateLeft = hUpdateLeft;
00174                   o_hUpdateRight = hUpdateRight;
00175                   o_huUpdateLeft = huUpdateLeft;
00176                   o_huUpdateRight = huUpdateRight;
00177 
00178                   //compute maximum wave speed (-> CFL-condition)
00179                   o_maxWaveSpeed = std::max( std::abs(waveSpeeds0) , std::abs(waveSpeeds1) );
00180     }
00181 
00182         inline
00183         void fWaveComputeWaveSpeeds(
00184             const T i_hLeft,  const T i_hRight,
00185             const T i_huLeft, const T i_huRight,
00186             const T i_uLeft,  const T i_uRight,
00187             const T i_bLeft,  const T i_bRight,
00188 
00189             T &o_waveSpeed0, T &o_waveSpeed1 ) const
00190         {
00191                 //compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i}
00192                 T characteristicSpeed0 = 0., characteristicSpeed1 = 0.;
00193                 characteristicSpeed0 = i_uLeft - std::sqrt(gravity*i_hLeft);
00194                 characteristicSpeed1 = i_uRight + std::sqrt(gravity*i_hRight);
00195 
00196                 //compute "Roe speeds"
00197                 T hRoe = (T).5 * (i_hRight + i_hLeft);
00198                 T uRoe = i_uLeft * std::sqrt(i_hLeft) + i_uRight * std::sqrt(i_hRight);
00199                 uRoe /= std::sqrt(i_hLeft) + std::sqrt(i_hRight);
00200 
00201                 T roeSpeed0 = 0., roeSpeed1 = 0.;
00202                 roeSpeed0 = uRoe - std::sqrt(gravity*hRoe);
00203                 roeSpeed1 = uRoe + std::sqrt(gravity*hRoe);
00204 
00205                 //computer eindfeldt speeds
00206                 o_waveSpeed0 = std::min(characteristicSpeed0, roeSpeed0);
00207                 o_waveSpeed1 = std::max(characteristicSpeed1, roeSpeed1);
00208         }
00209 
00210         inline
00211         void fWaveComputeWaveDecomposition(
00212                         const T i_hLeft,  const T i_hRight,
00213                         const T i_huLeft, const T i_huRight,
00214                         const T i_uLeft,  const T i_uRight,
00215                         const T i_bLeft,  const T i_bRight,
00216                         const T i_waveSpeed0, const T i_waveSpeed1,
00217 
00218                         T &o_fWave0, T &o_fWave1 ) const
00219         {
00220           T lambdaDif = i_waveSpeed1 - i_waveSpeed0;
00221 
00222           //compute the inverse matrix R^{-1}
00223           T Rinv00 = 0., Rinv01 = 0., Rinv10 = 0., Rinv11 = 0.;
00224 
00225           T oneDivLambdaDif = (T)1. /  lambdaDif;
00226           Rinv00 = oneDivLambdaDif *  i_waveSpeed1;
00227           Rinv01 = -oneDivLambdaDif;
00228 
00229           Rinv10 = oneDivLambdaDif * -i_waveSpeed0;
00230           Rinv11 = oneDivLambdaDif;
00231 
00232           //right hand side
00233           T fDif0 = 0., fDif1 = 0.;
00234 
00235           //calculate modified (bathymetry!) flux difference
00236           // f(Q_i) - f(Q_{i-1})
00237           fDif0 = i_huRight - i_huLeft;
00238           fDif1 = i_huRight * i_uRight + (T).5 * gravity * i_hRight * i_hRight
00239                   -(i_huLeft  * i_uLeft  + (T).5 * gravity * i_hLeft  * i_hLeft);
00240 
00241           // \delta x \Psi[2]
00242           T psi = -gravity * (T).5 * (i_hRight + i_hLeft)*(i_bRight - i_bLeft);
00243           fDif1 -= psi;
00244 
00245           //solve linear equations
00246           o_fWave0 = Rinv00 * fDif0 + Rinv01 * fDif1;
00247           o_fWave1 = Rinv10 * fDif0 + Rinv11 * fDif1;
00248         }
00249 };
00250 
00251 }
00252 
00253 #endif // FWAVEVEC_HPP_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends