SWE
/import/home/rettenbs/src/SWE/submodules/swe_solvers/src/solver/FWaveCuda.h
00001 
00053 #include <cmath>
00054 #include <cstdio>
00055 
00056 typedef float T;
00057 
00059 const T zeroTol = (T) 0.0000001;
00060 
00062 const T dryTol = (T) 100.;
00063 
00064 /*
00065  * Visual C++ compiler doesn't allow the definitions above.
00066  * Therefore the pre-compiler should be used.
00067  */
00068 //#define zeroTol (T)0.0000001
00069 //#define dryTol (T)100.
00070 
00071 __device__
00072 void fWaveComputeNetUpdates( const T i_gravity,
00073                              T i_hLeft,  T i_hRight,
00074                              T i_huLeft, T i_huRight,
00075                              T i_bLeft,  T i_bRight,
00076 
00077                              T o_netUpdates[5]);
00078 
00079 __device__
00080 inline void fWaveComputeWaveSpeeds( const T i_gravity,
00081                                     const T i_hLeft,  const T i_hRight,
00082                                     const T i_huLeft, const T i_huRight,
00083                                     const T i_uLeft,  const T i_uRight,
00084                                     const T i_bLeft,  const T i_bRight,
00085 
00086                                     T o_waveSpeeds[2] );
00087 
00088 __device__
00089 inline void fWaveComputeWaveDecomposition( const T i_gravity,
00090                                            const T i_hLeft,  const T i_hRight,
00091                                            const T i_huLeft, const T i_huRight,
00092                                            const T i_uLeft,  const T i_uRight,
00093                                            const T i_bLeft,  const T i_bRight,
00094                                            const T i_waveSpeeds[2],
00095 
00096                                            T o_fWaves[2][2] );
00097 
00121 __device__
00122 void fWaveComputeNetUpdates( const T i_gravity,
00123                              T i_hLeft,  T i_hRight,
00124                              T i_huLeft, T i_huRight,
00125                              T i_bLeft,  T i_bRight,
00126 
00127                              T o_netUpdates[5] ) {
00128   //reset net updates
00129   o_netUpdates[0] = (T)0.; //hUpdateLeft
00130   o_netUpdates[1] = (T)0.; //hUpdateRight
00131   o_netUpdates[2] = (T)0.; //huUpdateLeft
00132   o_netUpdates[3] = (T)0.; //huUpdateRight
00133 
00134   //reset the maximum wave speed
00135   o_netUpdates[4] = (T)0.; //maxWaveSpeed
00136 
00137   // determine the wet dry state and corr. values, if necessary.
00138   if( !(i_hLeft > dryTol && i_hRight > dryTol) ) {
00139     if( i_hLeft < dryTol && i_hRight < dryTol )
00140       return;
00141     else if ( i_hLeft < dryTol ) {
00142       i_hLeft = i_hRight;
00143       i_huLeft = -i_huRight;
00144       i_bLeft = i_bRight;
00145     }
00146     else { //( i_hRight < dryTol )
00147       i_hRight = i_hLeft;
00148       i_huRight = -i_huLeft;
00149       i_bRight = i_bLeft;
00150     }
00151   }
00152 
00154   T uLeft = i_huLeft / i_hLeft;
00156   T uRight = i_huRight / i_hRight;
00157 
00159   T waveSpeeds[2];
00160 
00161   //compute the wave speeds
00162   fWaveComputeWaveSpeeds( i_gravity,
00163                           i_hLeft, i_hRight,
00164                           i_huLeft, i_huRight,
00165                           uLeft, uRight,
00166                           i_bLeft, i_bRight,
00167 
00168                           waveSpeeds );
00169 
00171   T fWaves[2][2];
00172 
00173   //compute the decomposition into f-waves
00174   fWaveComputeWaveDecomposition( i_gravity,
00175                                  i_hLeft, i_hRight,
00176                                  i_huLeft, i_huRight,
00177                                  uLeft, uRight,
00178                                  i_bLeft, i_bRight,
00179 
00180                                  waveSpeeds,
00181                                  fWaves );
00182 
00183   //compute the net-updates
00184   //1st wave family
00185   if(waveSpeeds[0] < -zeroTol) { //left going
00186     o_netUpdates[0] +=  fWaves[0][0]; //hUpdateLeft
00187     o_netUpdates[2] += fWaves[0][1]; //huUpdateLeft
00188   }
00189   else if(waveSpeeds[0] > zeroTol) { //right going
00190     o_netUpdates[1] +=  fWaves[0][0]; //hUpdateRight
00191     o_netUpdates[3] += fWaves[0][1]; //huUpdateRight
00192   }
00193   else { //split waves
00194     o_netUpdates[0] +=   (T).5*fWaves[0][0]; //hUpdateLeft
00195     o_netUpdates[2] +=  (T).5*fWaves[0][1]; //huUpdateLeft
00196     o_netUpdates[1] +=  (T).5*fWaves[0][0]; //hUpdateRight
00197     o_netUpdates[3] += (T).5*fWaves[0][1]; //huUpdateRight
00198   }
00199 
00200   //2nd wave family
00201   if(waveSpeeds[1] < -zeroTol) { //left going
00202       o_netUpdates[0] +=  fWaves[1][0]; //hUpdateLeft
00203       o_netUpdates[2] += fWaves[1][1]; //huUpdateLeft
00204   }
00205   else if(waveSpeeds[1] > zeroTol) { //right going
00206       o_netUpdates[1] += fWaves[1][0]; //hUpdateRight
00207       o_netUpdates[3] += fWaves[1][1]; //huUpdateRight
00208   }
00209   else { //split waves
00210     o_netUpdates[0] +=   (T).5*fWaves[1][0]; //hUpdateLeft
00211     o_netUpdates[2] +=  (T).5*fWaves[1][1]; //huUpdateLeft
00212     o_netUpdates[1] +=  (T).5*fWaves[1][0]; //hUpdateRight
00213     o_netUpdates[3] += (T).5*fWaves[1][1]; //huUpdateRight
00214   }
00215 
00216   //compute maximum wave speed (-> CFL-condition)
00217   o_netUpdates[4] = fmax( fabs(waveSpeeds[0]) , fabs(waveSpeeds[1]) );
00218 }
00219 
00235 __device__
00236 inline void fWaveComputeWaveSpeeds( const T i_gravity,
00237                                     const T i_hLeft,  const T i_hRight,
00238                                     const T i_huLeft, const T i_huRight,
00239                                     const T i_uLeft,  const T i_uRight,
00240                                     const T i_bLeft,  const T i_bRight,
00241 
00242                                     T o_waveSpeeds[2] ) {
00243   //compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i}
00244   T characteristicSpeeds[2];
00245   characteristicSpeeds[0] = i_uLeft - sqrt(i_gravity*i_hLeft);
00246   characteristicSpeeds[1] = i_uRight + sqrt(i_gravity*i_hRight);
00247 
00248   //compute "Roe speeds"
00249   T hRoe = (T).5 * (i_hRight + i_hLeft);
00250   T uRoe = i_uLeft * sqrt(i_hLeft) + i_uRight * sqrt(i_hRight);
00251   uRoe /= sqrt(i_hLeft) + sqrt(i_hRight);
00252 
00253   T roeSpeeds[2];
00254   roeSpeeds[0] = uRoe - sqrt(i_gravity*hRoe);
00255   roeSpeeds[1] = uRoe + sqrt(i_gravity*hRoe);
00256 
00257   //computer eindfeldt speeds
00258   T einfeldtSpeeds[2];
00259   einfeldtSpeeds[0] = fmin(characteristicSpeeds[0], roeSpeeds[0]);
00260   einfeldtSpeeds[1] = fmax(characteristicSpeeds[1], roeSpeeds[1]);
00261 
00262   //set wave speeds
00263   o_waveSpeeds[0] = einfeldtSpeeds[0];
00264   o_waveSpeeds[1] = einfeldtSpeeds[1];
00265 }
00266 
00282 __device__
00283 inline void fWaveComputeWaveDecomposition( const T i_gravity,
00284                                            const T i_hLeft,  const T i_hRight,
00285                                            const T i_huLeft, const T i_huRight,
00286                                            const T i_uLeft,  const T i_uRight,
00287                                            const T i_bLeft,  const T i_bRight,
00288                                            const T i_waveSpeeds[2],
00289 
00290                                            T o_fWaves[2][2] ) {
00291   T lambdaDif = i_waveSpeeds[1] - i_waveSpeeds[0];
00292 
00293   //compute the inverse matrix R^{-1}
00294   T Rinv[2][2];
00295 
00296   T oneDivLambdaDif = (T)1. /  lambdaDif;
00297   Rinv[0][0] = oneDivLambdaDif *  i_waveSpeeds[1];
00298   Rinv[0][1] = -oneDivLambdaDif;
00299 
00300   Rinv[1][0] = oneDivLambdaDif * -i_waveSpeeds[0];
00301   Rinv[1][1] = oneDivLambdaDif;
00302 
00303   //right hand side
00304   T fDif[2];
00305 
00306   //calculate modified (bathymetry!) flux difference
00307   // f(Q_i) - f(Q_{i-1})
00308   fDif[0] = i_huRight - i_huLeft;
00309   fDif[1] = i_huRight * i_uRight + (T).5 * i_gravity * i_hRight * i_hRight
00310           -(i_huLeft  * i_uLeft  + (T).5 * i_gravity * i_hLeft  * i_hLeft);
00311 
00312   // \delta x \Psi[2]
00313   T psi = -i_gravity * (T).5 * (i_hRight + i_hLeft)*(i_bRight - i_bLeft);
00314   fDif[1] -= psi;
00315 
00316   //solve linear equations
00317   T beta[2];
00318   beta[0] = Rinv[0][0] * fDif[0] + Rinv[0][1] * fDif[1];
00319   beta[1] = Rinv[1][0] * fDif[0] + Rinv[1][1] * fDif[1];
00320 
00321   //return f-waves
00322   o_fWaves[0][0] = beta[0];
00323   o_fWaves[0][1] = beta[0] * i_waveSpeeds[0];
00324 
00325   o_fWaves[1][0] = beta[1];
00326   o_fWaves[1][1] = beta[1] * i_waveSpeeds[1];
00327 
00328   /*
00329    * in the case of a "zero strength" wave set the wave speeds to zero
00330    *   => non present waves don't affect the CFL-condition.
00331    * TODO: Seems not to affect the time step size
00332    */
00333 //  if( fmax(fabs(beta[0]), fabs(beta[1])) < zeroTol ) {
00334 //    io_waveSpeeds[0] = (T) 0.;
00335 //    io_waveSpeeds[1] = (T) 0.;
00336 //  }
00337 }
00338 
00339 // undefine the constants (see above)
00340 #undef zeroTol
00341 #undef dryTol
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends