|
SWE
|
00001 00063 /* 00064 * TODO: store nLow/nHigh variables in array[2] 00065 */ 00066 00067 00068 #ifndef AUGRIE_HPP_ 00069 #define AUGRIE_HPP_ 00070 00071 #include "WavePropagation.hpp" 00072 00073 #include <cassert> 00074 #include <cmath> 00075 #include <iostream> 00076 #include <string> 00077 #include <vector> 00078 00088 //#define correctrarefactions 00089 //#define complexsteadystatewave 00090 00091 namespace solver { 00092 template <typename T> class AugRie; 00093 } 00094 00100 template <typename T> class solver::AugRie: public WavePropagation<T> { 00101 private: //explicit for unit tests 00102 //use nondependent names (template base class) 00103 using solver::WavePropagation<T>::zeroTol; 00104 using solver::WavePropagation<T>::g; 00105 using solver::WavePropagation<T>::dryTol; 00106 00107 using solver::WavePropagation<T>::hLeft; 00108 using solver::WavePropagation<T>::hRight; 00109 using solver::WavePropagation<T>::huLeft; 00110 using solver::WavePropagation<T>::huRight; 00111 using solver::WavePropagation<T>::bLeft; 00112 using solver::WavePropagation<T>::bRight; 00113 using solver::WavePropagation<T>::uLeft; 00114 using solver::WavePropagation<T>::uRight; 00115 00116 using solver::WavePropagation<T>::wetDryState; 00117 using solver::WavePropagation<T>::DryDry; 00118 using solver::WavePropagation<T>::WetWet; 00119 using solver::WavePropagation<T>::WetDryInundation; 00120 using solver::WavePropagation<T>::WetDryWall; 00121 using solver::WavePropagation<T>::WetDryWallInundation; 00122 using solver::WavePropagation<T>::DryWetInundation; 00123 using solver::WavePropagation<T>::DryWetWall; 00124 using solver::WavePropagation<T>::DryWetWallInundation; 00125 00126 using solver::WavePropagation<T>::storeParameters; 00127 00128 00130 const T newtonTol; 00132 const int maxNumberOfNewtonIterations; 00133 00135 T hMiddle; 00136 00138 T middleStateSpeeds[2]; 00139 00140 00144 enum RiemannStructure { 00145 DrySingleRarefaction, 00146 SingleRarefactionDry, 00147 ShockShock, 00148 ShockRarefaction, 00149 RarefactionShock, 00150 RarefactionRarefaction 00151 }; 00152 00154 RiemannStructure riemannStructure; 00155 00156 public: 00166 AugRie( T i_dryTolerance = (T) 0.01, 00167 T i_gravity = (T) 9.81, 00168 T i_newtonTolerance = (T) 0.000001, 00169 int i_maxNumberOfNewtonIterations = 10, 00170 T i_zeroTolerance = (T) 0.00001 ): 00171 WavePropagation<T>( i_dryTolerance, i_gravity, i_zeroTolerance ), 00172 newtonTol(i_newtonTolerance), 00173 maxNumberOfNewtonIterations(i_maxNumberOfNewtonIterations) { 00174 00175 #ifndef NDEBUG 00176 #ifndef SUPPRESS_SOLVER_DEBUG_OUTPUT 00177 //print some information about the used solver 00178 std::cout << " *** solver::AugRie created" << std::endl 00179 << " zeroTolerance=" << zeroTol << std::endl 00180 << " gravity=" << g << std::endl 00181 << " dryTolerance=" << dryTol << std::endl 00182 << " newtonTolerance=" << newtonTol << std::endl 00183 << " maxNumberNumberOfNewtonIterations=" << maxNumberOfNewtonIterations << std::endl 00184 << "\n optional features:" 00185 #ifdef correctrarefactions 00186 << " correctrarefactions" 00187 #endif 00188 #ifdef complexsteadystatewave 00189 << " complexsteadystatewave" 00190 #endif 00191 << "\n ***\n\n"; 00192 #endif 00193 #endif 00194 00195 } 00196 //declare variables which are used over and over again 00197 #if 0 00198 T sqrt_g_h[2]; 00199 T sqrt_h[2]; 00200 00201 #define sqrt_g_hLeft (sqrt_g_h[0]) 00202 #define sqrt_g_hRight (sqrt_g_h[1]) 00203 00204 #define sqrt_hLeft (sqrt_h[0]) 00205 #define sqrt_hRight (sqrt_h[1]) 00206 00207 #else 00208 T sqrt_g_hLeft; 00209 T sqrt_g_hRight; 00210 00211 T sqrt_hLeft; 00212 T sqrt_hRight; 00213 #endif 00214 00215 00216 T sqrt_g; 00217 00218 00224 void computeNetUpdates ( const T &i_hLeft, const T &i_hRight, 00225 const T &i_huLeft, const T &i_huRight, 00226 const T &i_bLeft, const T &i_bRight, 00227 00228 T &o_hUpdateLeft, 00229 T &o_hUpdateRight, 00230 T &o_huUpdateLeft, 00231 T &o_huUpdateRight, 00232 T &o_maxWaveSpeed ) { 00233 // store parameters to member variables 00234 storeParameters( i_hLeft, i_hRight, 00235 i_huLeft, i_huRight, 00236 i_bLeft, i_bRight ); 00237 00238 //set speeds to zero (will be determined later) 00239 uLeft = uRight = 0.; 00240 00241 //reset net updates and the maximum wave speed 00242 o_hUpdateLeft = o_hUpdateRight = o_huUpdateLeft = o_huUpdateRight = (T)0.; 00243 o_maxWaveSpeed = (T)0.; 00244 00245 //determine the wet/dry state and compute local variables correspondingly 00246 determineWetDryState(); 00247 00248 if (wetDryState == DryDry) //nothing to do in a dry region 00249 return; 00250 00251 //precompute some terms which are fixed during 00252 //the computation after some specific point 00253 sqrt_g = std::sqrt(g); 00254 sqrt_hLeft = std::sqrt(hLeft); 00255 sqrt_hRight = std::sqrt(hRight); 00256 00257 sqrt_g_hLeft = sqrt_g*sqrt_hLeft; 00258 sqrt_g_hRight = sqrt_g*sqrt_hRight; 00259 00260 00261 //where to store the three waves 00262 T fWaves[3][2]; 00263 //and their speeds 00264 T waveSpeeds[3]; 00265 00266 //compute the augmented decomposition 00267 // (thats the place where the computational work is done..) 00268 computeWaveDecomposition( fWaves, 00269 waveSpeeds ); 00270 00271 00272 //compute the updates from the three propagating waves 00273 //A^-\delta Q = \sum{s[i]<0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i 00274 //A^+\delta Q = \sum{s[i]>0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i 00275 for (int waveNumber = 0; waveNumber < 3; waveNumber++) { 00276 if (waveSpeeds[waveNumber] < -zeroTol) { //left going 00277 o_hUpdateLeft += fWaves[waveNumber][0]; 00278 o_huUpdateLeft += fWaves[waveNumber][1]; 00279 } 00280 00281 else if (waveSpeeds[waveNumber] > zeroTol) { //right going 00282 o_hUpdateRight += fWaves[waveNumber][0]; 00283 o_huUpdateRight += fWaves[waveNumber][1]; 00284 } 00285 else { //TODO: this case should not happen mathematically, but it does. Where is the bug? Machine accuracy only? 00286 o_hUpdateLeft += (T).5*fWaves[waveNumber][0]; 00287 o_huUpdateLeft += (T).5*fWaves[waveNumber][1]; 00288 00289 o_hUpdateRight += (T).5*fWaves[waveNumber][0]; 00290 o_huUpdateRight += (T).5*fWaves[waveNumber][1]; 00291 } 00292 00293 //no wave speeds => zero strength fWaves 00294 // assert( std::fabs(fWaves[waveNumber][0]) < zeroTol ); 00295 // assert( std::fabs(fWaves[waveNumber][1]) < zeroTol ); 00296 } 00297 00298 #ifndef NDEBUG 00299 if(wetDryState == DryWetWall) 00300 assert( std::fabs(o_hUpdateLeft) < zeroTol && std::fabs(o_huUpdateLeft) < zeroTol ); 00301 else if(wetDryState == WetDryWall) 00302 assert( std::fabs(o_hUpdateRight) < zeroTol && std::fabs(o_huUpdateRight) < zeroTol ); 00303 #endif 00304 00305 //compute maximum wave speed (-> CFL-condition) 00306 waveSpeeds[0] = std::fabs(waveSpeeds[0]); 00307 waveSpeeds[1] = std::fabs(waveSpeeds[1]); 00308 waveSpeeds[2] = std::fabs(waveSpeeds[2]); 00309 00310 o_maxWaveSpeed = std::max(waveSpeeds[0], waveSpeeds[1]); 00311 o_maxWaveSpeed = std::max(o_maxWaveSpeed, waveSpeeds[2]); 00312 } 00313 00314 private: 00318 void determineWetDryState() { 00319 //compute speeds or set them to zero (dry cells) 00320 if (hLeft > dryTol) { 00321 uLeft = huLeft / hLeft; 00322 } 00323 else { 00324 bLeft += hLeft; 00325 hLeft = huLeft = uLeft = 0; 00326 } 00327 00328 if (hRight > dryTol) { 00329 uRight = huRight / hRight; 00330 } 00331 else { 00332 bRight += hRight; 00333 hRight = huRight = uRight = 0; 00334 } 00335 00336 //test for simple wet/wet case since this is most probably the 00337 //most frequently executed branch 00338 if (hLeft >= dryTol && hRight >= dryTol) { 00339 wetDryState = WetWet; 00340 } 00341 00342 //check for the dry/dry-case 00343 else if (hLeft < dryTol && hRight < dryTol) 00344 wetDryState = DryDry; 00345 00346 //we have a shoreline: one cell dry, one cell wet 00347 00348 //check for simple inundation problems 00349 // (=> dry cell lies lower than the wet cell) 00350 else if (hLeft < dryTol && hRight + bRight > bLeft) 00351 wetDryState = DryWetInundation; 00352 00353 else if (hRight < dryTol && hLeft + bLeft > bRight) 00354 wetDryState = WetDryInundation; 00355 00356 //dry cell lies higher than the wet cell 00357 else if (hLeft < dryTol) { 00358 //lets check if the momentum is able to overcome the difference in height 00359 // => solve homogeneous Riemann-problem to determine the middle state height 00360 // which would arise if there is a wall (wall-boundary-condition) 00361 // \cite[ch. 6.8.2]{george2006finite}) 00362 // \cite[ch. 5.2]{george2008augmented} 00363 computeMiddleState( hRight, hRight, 00364 -uRight, uRight, 00365 -huRight, huRight, 00366 maxNumberOfNewtonIterations ); 00367 00368 if (hMiddle + bRight > bLeft) { 00369 //momentum is large enough, continue with the original values 00370 // bLeft = hMiddle + bRight; 00371 wetDryState = DryWetWallInundation; 00372 } 00373 else { 00374 //momentum is not large enough, use wall-boundary-values 00375 hLeft = hRight; 00376 uLeft = -uRight; 00377 huLeft = -huRight; 00378 bLeft = bRight = (T)0.; 00379 wetDryState = DryWetWall; 00380 } 00381 } 00382 00383 else if (hRight < dryTol) { 00384 //lets check if the momentum is able to overcome the difference in height 00385 // => solve homogeneous Riemann-problem to determine the middle state height 00386 // which would arise if there is a wall (wall-boundary-condition) 00387 // \cite[ch. 6.8.2]{george2006finite}) 00388 // \cite[ch. 5.2]{george2008augmented} 00389 computeMiddleState( hLeft, hLeft, 00390 uLeft, -uLeft, 00391 huLeft, -huLeft, 00392 maxNumberOfNewtonIterations ); 00393 00394 if (hMiddle + bLeft > bRight) { 00395 //momentum is large enough, continue with the original values 00396 // bRight = hMiddle + bLeft; 00397 wetDryState = WetDryWallInundation; 00398 } 00399 else { 00400 hRight = hLeft; 00401 uRight = -uLeft; 00402 huRight = -huLeft; 00403 bRight = bLeft = (T)0.; 00404 wetDryState = WetDryWall; 00405 } 00406 } 00407 //done with all cases 00408 else assert(false); 00409 00410 //limit the effect of the source term if there is a "wall" 00411 //\cite[end of ch. 5.2?]{george2008augmented} 00412 //\cite[rpn2ez_fast_geo.f][levequeclawpack] 00413 if ( wetDryState == DryWetWallInundation ) 00414 bLeft = hRight + bRight; 00415 00416 else if( wetDryState == WetDryWallInundation ) 00417 bRight = hLeft + bLeft; 00418 } 00419 00426 inline void computeWaveDecomposition( T o_fWaves[3][2], 00427 T o_waveSpeeds[3] ) { 00428 //compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i} (char. speeds) 00429 T characteristicSpeeds[2]; 00430 characteristicSpeeds[0] = uLeft - sqrt_g_hLeft; 00431 characteristicSpeeds[1] = uRight + sqrt_g_hRight; 00432 00433 //compute "Roe speeds" 00434 T hRoe = (T).5 * (hRight + hLeft); 00435 T uRoe = uLeft * sqrt_hLeft + uRight * sqrt_hRight; 00436 uRoe /= sqrt_hLeft + sqrt_hRight; 00437 00438 T roeSpeeds[2]; 00439 //optimization for dumb compilers 00440 T sqrt_g_hRoe = std::sqrt(g*hRoe); 00441 roeSpeeds[0] = uRoe - sqrt_g_hRoe; 00442 roeSpeeds[1] = uRoe + sqrt_g_hRoe; 00443 00444 //compute the middle state of the homogeneous Riemann-Problem 00445 if(wetDryState != WetDryWall && wetDryState != DryWetWall) { //case WDW and DWW was computed in determineWetDryState already 00446 computeMiddleState( hLeft, hRight, 00447 uLeft, uRight, 00448 huLeft, huRight 00449 ); 00450 } 00451 00452 //compute extended eindfeldt speeds (einfeldt speeds + middle state speeds) 00453 // \cite[ch. 5.2]{george2008augmented}, \cite[ch. 6.8]{george2006finite} 00454 T extEinfeldtSpeeds[2]= {(T)0, (T)0}; 00455 if( wetDryState == WetWet || 00456 wetDryState == WetDryWall || 00457 wetDryState == DryWetWall ) { 00458 extEinfeldtSpeeds[0] = std::min(characteristicSpeeds[0], roeSpeeds[0]); 00459 extEinfeldtSpeeds[0] = std::min(extEinfeldtSpeeds[0], middleStateSpeeds[1]); 00460 00461 extEinfeldtSpeeds[1] = std::max(characteristicSpeeds[1], roeSpeeds[1]); 00462 extEinfeldtSpeeds[1] = std::max(extEinfeldtSpeeds[1], middleStateSpeeds[0]); 00463 } 00464 else if (hLeft < dryTol) { //ignore undefined speeds 00465 extEinfeldtSpeeds[0] = std::min(roeSpeeds[0], middleStateSpeeds[1]); 00466 extEinfeldtSpeeds[1] = std::max(characteristicSpeeds[1], roeSpeeds[1]); 00467 00468 assert(middleStateSpeeds[0] < extEinfeldtSpeeds[1]); 00469 } 00470 else if (hRight < dryTol) { //ignore undefined speeds 00471 extEinfeldtSpeeds[0] = std::min(characteristicSpeeds[0], roeSpeeds[0]); 00472 extEinfeldtSpeeds[1] = std::max(roeSpeeds[1], middleStateSpeeds[0]); 00473 00474 assert(middleStateSpeeds[1] > extEinfeldtSpeeds[0]); 00475 } 00476 else { 00477 assert(false); 00478 } 00479 00480 //HLL middle state 00481 // \cite[theorem 3.1]{george2006finite}, \cite[ch. 4.1]{george2008augmented} 00482 T hLLMiddleHeight = huLeft - huRight + extEinfeldtSpeeds[1] * hRight - extEinfeldtSpeeds[0] * hLeft; 00483 hLLMiddleHeight /= extEinfeldtSpeeds[1] - extEinfeldtSpeeds[0]; 00484 hLLMiddleHeight = std::max(hLLMiddleHeight, (T)0.); 00485 00486 //define eigenvalues 00487 T eigenValues[3]; 00488 eigenValues[0] = extEinfeldtSpeeds[0]; 00489 //eigenValues[1] --> corrector wave 00490 eigenValues[2] = extEinfeldtSpeeds[1]; 00491 00492 //define eigenvectors 00493 T eigenVectors[3][3]; 00494 00495 //set first and third eigenvector 00496 eigenVectors[0][0] = (T)1; 00497 eigenVectors[0][2] = (T)1; 00498 00499 eigenVectors[1][0] = eigenValues[0]; 00500 eigenVectors[1][2] = eigenValues[2]; 00501 00502 eigenVectors[2][0] = eigenValues[0]*eigenValues[0]; 00503 eigenVectors[2][2] = eigenValues[2]*eigenValues[2]; 00504 00505 00506 //compute rarefaction corrector wave 00507 // \cite[ch. 6.7.2]{george2006finite}, \cite[ch. 5.1]{george2008augmented} 00508 bool strongRarefaction = false; 00509 #ifdef correctrarefactions 00510 if ( (riemannStructure == ShockRarefaction || 00511 riemannStructure == RarefactionShock || 00512 riemannStructure == RarefactionRarefaction) && 00513 hMiddle > dryTol ) { //limit to ensure non-negative depth 00514 00515 //TODO: GeoClaw, riemann_aug_JCP; No explicit boundaries for "strong rarefaction" in literature? 00516 T rareBarrier[2] = {.5, .9}; 00517 00518 //lower rare barrier in the case of a transsonic rarefaction 00519 if ( (riemannStructure == RarefactionShock || riemannStructure == RarefactionRarefaction) && 00520 eigenValues[0]*middleStateSpeeds[0] < (T).0 ) { //transsonic rarefaction, first wave family 00521 rareBarrier[0] = (T).2; 00522 } 00523 else if ( (riemannStructure == ShockRarefaction || riemannStructure == RarefactionRarefaction) && 00524 eigenValues[2]*middleStateSpeeds[1] < (T)0. ) { //transsonic rarefaction, second wave family 00525 rareBarrier[0] = (T).2; 00526 } 00527 00528 00529 T sqrt_g_hMiddle = std::sqrt(g*hMiddle); 00530 //determine the max. rarefaction size (distance between head and tail) 00531 T rareFactionSize[2]; 00532 rareFactionSize[0] = (T)3. * (sqrt_g_hLeft - sqrt_g_hMiddle); 00533 rareFactionSize[1] = (T)3. * (sqrt_g_hRight - sqrt_g_hMiddle); 00534 00535 T maxRareFactionSize = std::max(rareFactionSize[0], rareFactionSize[1]); 00536 00537 //set the eigenvalue of the corrector wave in the case of a "strong rarefaction" 00538 if ( maxRareFactionSize > rareBarrier[0] * (eigenValues[2] - eigenValues[0]) && 00539 maxRareFactionSize < rareBarrier[1] * (eigenValues[2] - eigenValues[0]) 00540 ) { 00541 strongRarefaction = true; 00542 if (rareFactionSize[0] > rareFactionSize[1]) 00543 eigenValues[1] = middleStateSpeeds[0]; 00544 else 00545 eigenValues[1] = middleStateSpeeds[1]; 00546 } 00547 00548 //TODO: implemented in clawpack, why? 00549 // if ( hMiddle < std::min(hLeft, hRight)/5.) { //middle state in an HLL solve 00550 // strongRarefaction = false; 00551 // } 00552 } 00553 #endif 00554 00555 //set 2nd eigenvector 00556 if(strongRarefaction == false) { 00557 eigenValues[1] = (T).5*(eigenValues[0]+eigenValues[2]); 00558 eigenVectors[0][1] = 0.; 00559 eigenVectors[1][1] = 0.; 00560 eigenVectors[2][1] = 1.; 00561 } 00562 else { 00563 eigenVectors[0][1] = (T)1.; 00564 eigenVectors[1][1] = eigenValues[1]; 00565 eigenVectors[2][1] = eigenValues[1]*eigenValues[1]; 00566 } 00567 00568 //compute the jump in state 00569 T rightHandSide[3]; 00570 rightHandSide[0] = hRight - hLeft; 00571 rightHandSide[1] = huRight - huLeft; 00572 rightHandSide[2] = huRight * uRight + (T).5 * g * hRight * hRight 00573 -( huLeft * uLeft + (T).5 * g * hLeft * hLeft ); 00574 00575 //compute steady state wave 00576 // \cite[ch. 4.2.1 \& app. A]{george2008augmented}, \cite[ch. 6.2 \& ch. 4.4]{george2006finite} 00577 T steadyStateWave[2]; 00578 T hBar = (hLeft + hRight)*(T).5; 00579 #ifdef complexsteadystatewave 00580 T lLambdaBar = (uLeft + uRight)*(uLeft + uRight)*.25 - g*hBar; 00581 00582 T lLambdaTilde = std::max((T)0., uLeft*uRight) - g*hBar; 00583 00584 //near sonic as defined in geoclaw (TODO: literature?) 00585 if( ( std::fabs(lLambdaBar) < zeroTol ) || 00586 ( lLambdaBar*lLambdaTilde < zeroTol ) || 00587 ( lLambdaBar*eigenValues[0]*eigenValues[1] < zeroTol ) || 00588 ( std::min( std::fabs( eigenValues[0]), std::fabs(eigenValues[2]) ) < zeroTol ) || 00589 ( eigenValues[0] < (T)0. && middleStateSpeeds[0] > (T)0. ) || 00590 ( eigenValues[2] > (T)0. && middleStateSpeeds[1] < (T)0. ) || 00591 ( (uLeft + sqrt(g*hLeft)) * (uRight+sqrt(g*hRight)) < (T)0. ) || 00592 ( (uLeft - sqrt(g*hLeft)) * (uRight-sqrt(g*hRight)) < (T)0. ) ) { 00593 #endif 00594 steadyStateWave[0] = -(bRight - bLeft); 00595 steadyStateWave[1] = -g * hBar * (bRight - bLeft); 00596 #ifdef complexsteadystatewave 00597 } 00598 else { 00599 steadyStateWave[0] = g * (hBar / lLambdaBar) * (bRight - bLeft); 00600 T hTilde = hBar * lLambdaTilde / lLambdaBar; 00601 00602 //set bounds for problems far from steady state 00603 hTilde = std::max(std::min(hLeft, hRight), hTilde); 00604 hTilde = std::min(std::max(hLeft, hRight), hTilde); 00605 steadyStateWave[1] = -g*hTilde*(bRight - bLeft); 00606 } 00607 #endif 00608 00609 //preserve depth-positivity 00610 // \cite[ch. 6.5.2]{george2006finite}, \cite[ch. 4.2.3]{george2008augmented} 00611 if(eigenValues[0] < -zeroTol && eigenValues[2] > zeroTol) { //subsonic 00612 steadyStateWave[0] = std::max(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[0]); 00613 steadyStateWave[0] = std::min(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[2]); 00614 } 00615 else if(eigenValues[0] > zeroTol) { //supersonic right TODO: motivation? 00616 steadyStateWave[0] = std::max(steadyStateWave[0], -hLeft); 00617 steadyStateWave[0] = std::min(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[0]); 00618 } 00619 else if(eigenValues[2] < -zeroTol) { //supersonic left TODO: motivation? 00620 steadyStateWave[0] = std::max(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[2]); 00621 steadyStateWave[0] = std::min(steadyStateWave[0], hRight); 00622 } 00623 00624 //Limit the effect of the source term 00625 // \cite[ch. 6.4.2]{george2006finite} 00626 steadyStateWave[1] = std::min(steadyStateWave[1], g*std::max(-hLeft*(bRight-bLeft) , -hRight*(bRight-bLeft))); 00627 steadyStateWave[1] = std::max(steadyStateWave[1], g*std::min(-hLeft*(bRight-bLeft) , -hRight*(bRight-bLeft))); 00628 00629 //no source term in the case of a wall 00630 if( wetDryState == WetDryWall || wetDryState == DryWetWall ) { 00631 assert(std::fabs(steadyStateWave[0])<zeroTol); 00632 assert(std::fabs(steadyStateWave[1])<zeroTol); 00633 } 00634 00635 rightHandSide[0] -= steadyStateWave[0]; 00636 //rightHandSide[1]: no source term 00637 rightHandSide[2] -= steadyStateWave[1]; 00638 00639 //everything is ready, solve the equations! 00640 T beta[3]; 00641 solveLinearEquation(eigenVectors, rightHandSide, beta); 00642 00643 //compute f-waves and wave-speeds 00644 if (wetDryState == WetDryWall) { //zero ghost updates (wall boundary) 00645 //care about the left going wave (0) only 00646 o_fWaves[0][0] = beta[0] * eigenVectors[1][0]; 00647 o_fWaves[0][1] = beta[0] * eigenVectors[2][0]; 00648 00649 //set the rest to zero 00650 o_fWaves[1][0] = o_fWaves[1][1] = (T)0.; 00651 o_fWaves[2][0] = o_fWaves[2][1] = (T)0.; 00652 00653 o_waveSpeeds[0] = eigenValues[0]; 00654 o_waveSpeeds[1] = o_waveSpeeds[2] = (T)0.; 00655 00656 assert(eigenValues[0] < zeroTol); 00657 } 00658 else if (wetDryState == DryWetWall) { //zero ghost updates (wall boundary) 00659 //care about the right going wave (2) only 00660 o_fWaves[2][0] = beta[2] * eigenVectors[1][2]; 00661 o_fWaves[2][1] = beta[2] * eigenVectors[2][2]; 00662 00663 //set the rest to zero 00664 o_fWaves[0][0] = o_fWaves[0][1] = (T)0.; 00665 o_fWaves[1][0] = o_fWaves[1][1] = (T)0.; 00666 00667 o_waveSpeeds[2] = eigenValues[2]; 00668 o_waveSpeeds[0] = o_waveSpeeds[1] = 0.; 00669 00670 assert(eigenValues[2] > -zeroTol); 00671 } 00672 else { 00673 //compute f-waves (default) 00674 for (int waveNumber = 0; waveNumber < 3; waveNumber++) { 00675 o_fWaves[waveNumber][0] = beta[waveNumber] * eigenVectors[1][waveNumber]; //select 2nd and 00676 o_fWaves[waveNumber][1] = beta[waveNumber] * eigenVectors[2][waveNumber]; //3rd component of the augmented decomposition 00677 } 00678 00679 o_waveSpeeds[0] = eigenValues[0]; 00680 o_waveSpeeds[1] = eigenValues[1]; 00681 o_waveSpeeds[2] = eigenValues[2]; 00682 } 00683 } 00684 00697 inline void computeMiddleState( const T &i_hLeft, const T &i_hRight, 00698 const T &i_uLeft, const T &i_uRight, 00699 const T &i_huLeft, const T &i_huRight, 00700 const int &i_maxNumberOfNewtonIterations = 1 ) { 00701 00702 //set everything to zero 00703 hMiddle = (T)0.; 00704 middleStateSpeeds[0] = (T)0.; 00705 middleStateSpeeds[1] = (T)0.; 00706 00707 //compute local square roots 00708 //(not necessarily the same ones as in computeNetUpdates!) 00709 T l_sqrt_g_hRight = std::sqrt(g*i_hRight); 00710 T l_sqrt_g_hLeft = std::sqrt(g*i_hLeft); 00711 00712 //single rarefaction in the case of a wet/dry interface 00713 if (i_hLeft < dryTol) { 00714 middleStateSpeeds[1] = middleStateSpeeds[0] = i_uRight - (T)2*l_sqrt_g_hRight; 00715 riemannStructure = DrySingleRarefaction; 00716 return; 00717 } 00718 else if (i_hRight < dryTol) { 00719 middleStateSpeeds[0] = middleStateSpeeds[1] = i_uLeft + (T)2*l_sqrt_g_hLeft; 00720 riemannStructure = SingleRarefactionDry; 00721 return; 00722 } 00723 00724 //determine the wave structure of the Riemann-problem 00725 riemannStructure = determineRiemannStructure( i_hLeft, i_hRight, 00726 i_uLeft, i_uRight ); 00727 00728 //will be computed later 00729 T sqrt_g_hMiddle = 0.; 00730 00731 if (riemannStructure == ShockShock) { 00732 /*compute All-Shock Riemann Solution 00733 * \cite[ch. 13.7]{leveque2002finite} 00734 * 00735 *compute middle state h_m 00736 * => Solve non-linear scalar equation by Newton's method 00737 * u_l - (h_m - h_l) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_l} \right) } 00738 * = u_r + (h_m - h_r) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_r} \right) } 00739 * 00740 * Therefore determine the root of phi_{ss}: 00741 *\begin{equation} 00742 * \begin{matrix} 00743 * \phi_{ss}(h) =&& u_r - u_l + 00744 * (h-h_l) \sqrt{ 00745 * \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right) 00746 * } \\ 00747 * && +(h-h_r) \sqrt{ 00748 * \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right) 00749 * } 00750 * \end{matrix} 00751 *\end{equation} 00752 * 00753 *\begin{equation} 00754 * \begin{matrix} 00755 * \phi_{ss}'(h) = &&\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right) } 00756 * +\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right) }\\ 00757 * && -\frac{g}{4} 00758 * \frac{h-h_l} 00759 * { 00760 * h^2\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right) } 00761 * }- 00762 * \frac{g}{4} 00763 * \frac{h-h_r} 00764 * { 00765 * h^2\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right) } 00766 * } 00767 * \end{matrix} 00768 *\end{equation} 00769 */ 00770 00771 // hMiddle = (hLeft + hRight)*(T)0.618; //first estimate 00772 hMiddle = std::min(i_hLeft, i_hRight); 00773 00774 T l_sqrtTermH[2] = {0,0}; 00775 00776 for (int i = 0; i < i_maxNumberOfNewtonIterations; i++) { 00777 // sqrtTermHLow = std::sqrt((T).5 * g * ((T)1/hMiddle + (T)1/hLeft)); 00778 l_sqrtTermH[0] = std::sqrt((T).5 * g * ((hMiddle+i_hLeft)/(hMiddle*i_hLeft))); 00779 // sqrtTermHHigh = std::sqrt((T).5 * g * ((T)1/hMiddle + (T)1/hRight)); 00780 l_sqrtTermH[1] = std::sqrt((T).5 * g * ((hMiddle+i_hRight)/(hMiddle*i_hRight))); 00781 00782 T phi = i_uRight - i_uLeft + (hMiddle - i_hLeft) * l_sqrtTermH[0] + (hMiddle - i_hRight) * l_sqrtTermH[1]; 00783 00784 if (std::fabs(phi) < newtonTol) 00785 break; 00786 00787 T derivativePhi = l_sqrtTermH[0] + l_sqrtTermH[1] 00788 - (T).25*g*(hMiddle - i_hLeft) / (l_sqrtTermH[0]*hMiddle*hMiddle) 00789 - (T).25*g*(hMiddle - i_hRight) / (l_sqrtTermH[1]*hMiddle*hMiddle); 00790 00791 hMiddle = hMiddle - phi/derivativePhi; //Newton step 00792 assert(hMiddle >= dryTol); 00793 00794 // if(i == i_maxNumberOfNewtonIterations-1) { 00795 // std::cerr << "Newton-Method did not converge" << std::endl; 00796 // std::cerr << "std::fabs(phi): " << std::fabs(phi) << std::endl; 00797 // std::cerr << "hMiddle: " << hMiddle << std::endl; 00798 // assert(false); 00799 // } 00800 } 00801 00802 sqrt_g_hMiddle = std::sqrt(g*hMiddle); 00803 00804 //compute middle speed u_m 00805 //\begin{equation} 00806 // \label{eq:hmshock1stfamsimp} 00807 // u_m = u_l - (h_m - h_l) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_l} \right) } 00808 //\end{equation} 00809 // 00810 //\begin{equation} 00811 // \label{eq:hmshock2ndfamsimp} 00812 // u_m = u_r + (h_m - h_r) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_r} \right) } 00813 //\end{equation} 00814 00815 // T uMiddleEstimates[2]; 00816 // uMiddleEstimates[0] = uLeft - (hMiddle - hLeft) * l_sqrtTermH[0]; 00817 // uMiddleEstimates[1] = i_uRight + (hMiddle - hRight) * l_sqrtTermH[1]; 00818 // uMiddle = (T).5 * (uMiddleEstimates[0] + uMiddleEstimates[1]); 00819 00820 //middle state speeds as they are implemented in clawpack, TODO: why? 00821 // middleStateSpeeds[0] = uMiddleEstimates[0] - sqrt_g_hMiddle; 00822 // middleStateSpeeds[1] = uMiddleEstimates[1] + sqrt_g_hMiddle; 00823 } 00824 00825 if (riemannStructure == RarefactionRarefaction) { 00826 //Compute All-Rarefaction Riemann Solution 00827 // \cite[ch. 13.8.6]{leveque2002finite} 00828 00829 //compute middle state height h_m 00830 hMiddle = std::max((T)0., i_uLeft - i_uRight + (T)2.*(l_sqrt_g_hLeft + l_sqrt_g_hRight)); //std::max -> Text after formula (13.56), page 279 00831 hMiddle = (T)1/((T)16*g)* hMiddle * hMiddle; 00832 00833 sqrt_g_hMiddle = std::sqrt(g*hMiddle); 00834 00835 //middle state speeds as they are implemented in clawpack, why? 00836 // middleStateSpeeds[0] = uLeft + (T)2.*l_sqrt_g_hLow - (T)3.*sqrt_g_hMiddle; 00837 // middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hHigh + (T)3.*sqrt_g_hMiddle; 00838 } 00839 00840 if (riemannStructure == ShockRarefaction || riemannStructure == RarefactionShock) { 00841 //compute dam-break Riemann-solution 00842 //TODO: reference 00843 T hMin, hMax; 00844 if (riemannStructure == ShockRarefaction) { 00845 hMin = i_hLeft; 00846 hMax = i_hRight; 00847 } 00848 else { 00849 hMin = i_hRight; 00850 hMax = i_hLeft; 00851 } 00852 00853 //hMiddle = (hLeft + hRight)*(T)0.618; //first estimate 00854 hMiddle = hMin; 00855 00856 sqrt_g_hMiddle = std::sqrt(g*hMiddle); 00857 T sqrt_g_hMax = std::sqrt(g*hMax); 00858 for (int i = 0; i < i_maxNumberOfNewtonIterations; i++) { 00859 /* 00860 *compute middle state h_m 00861 * => Solve non-linear scalar equation by Newton's method 00862 * 00863 * Therefore determine the root of phi_{sr}/phi_{rs}: 00864 *\begin{equation} 00865 * \begin{matrix} 00866 * h_{min} = \min(h_l, h_r)\\ 00867 * h_{max} = \max(h_l, h_r) 00868 * \end{matrix} 00869 *\end{equation} 00870 *\begin{equation} 00871 * \begin{matrix} 00872 * \phi_{sr}(h) = \phi_{rs}(h) =&& u_r - u_l + (h - h_{min}) \sqrt{ 00873 * \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_{min}} \right) 00874 * } \\ 00875 * && + 2 (\sqrt{gh} - \sqrt{gh_{max}}) 00876 * \end{matrix} 00877 *\end{equation} 00878 *\begin{equation} 00879 * \phi'_{sr}(h) = \phi'_{rs}(h) = \sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_{min}} \right) } 00880 * -\frac{g}{4} \frac{h-h_{min}} 00881 * { h^2 00882 * \sqrt{ 00883 * \frac{g}{2} \left(\frac{1}{h} + \frac{1}{h_{min}} \right) 00884 * } 00885 * } 00886 * + \sqrt{\frac{g}{h}} 00887 *\end{equation} 00888 */ 00889 00890 T sqrtTermHMin = std::sqrt( (T).5*g* ((hMiddle+hMin)/(hMiddle*hMin)) ); 00891 00892 T phi = i_uRight - i_uLeft + (hMiddle - hMin) * sqrtTermHMin + (T)2*(sqrt_g_hMiddle - sqrt_g_hMax); 00893 00894 if (std::fabs(phi) < newtonTol) 00895 break; 00896 00897 T derivativePhi = sqrtTermHMin - (T).25*g * (hMiddle - hMin) / (hMiddle * hMiddle * sqrtTermHMin) + sqrt_g/sqrt_g_hMiddle; 00898 00899 hMiddle = hMiddle - phi/derivativePhi; //Newton step 00900 00901 #ifndef NDEBUG 00902 if (hMiddle < hMin) { 00903 std::cout << phi << std::endl; 00904 std::cout << derivativePhi << std::endl; 00905 std::cerr << "hMiddle(" << hMiddle << ") < hMin(" << hMin << ")" << std::endl; 00906 assert(false); 00907 } 00908 #endif 00909 00910 // if(i == i_maxNumberOfNewtonIterations-1) { 00911 // std::cerr << "Newton-Method did not converge" << std::endl; 00912 // std::cerr << "std::fabs(phi): " << std::fabs(phi) << std::endl; 00913 // std::cerr << "hMiddle: " << hMiddle << std::endl; 00914 // assert(false); 00915 // } 00916 00917 sqrt_g_hMiddle = std::sqrt(g*hMiddle); 00918 } 00919 00920 // //middle state speeds as they are implemented in clawpack, TODO: why? 00921 // if (riemannStructure == ShockRarefaction) { 00922 // uMiddle = i_uRight - (T)2.*( l_sqrt_g_hHigh - sqrt_g_hMiddle ); 00923 // middleStateSpeeds[0] = i_uRight - (T)2.*l_sqrt_g_hHigh + sqrt_g_hMiddle; 00924 // middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hHigh + (T)3.*sqrt_g_hMiddle; 00925 // } 00926 // else { 00927 // uMiddle = uLeft + (T)2.*( l_sqrt_g_hLow - sqrt_g_hMiddle ); 00928 // middleStateSpeeds[0] = uLeft + (T)2.*l_sqrt_g_hLow - (T)3.*sqrt_g_hMiddle; 00929 // middleStateSpeeds[1] = uLeft + (T)2.*l_sqrt_g_hLow - sqrt_g_hMiddle; 00930 // } 00931 } 00932 00933 middleStateSpeeds[0] = i_uLeft + (T)2.*l_sqrt_g_hLeft - (T)3.*sqrt_g_hMiddle; 00934 middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hRight + (T)3.*sqrt_g_hMiddle; 00935 00936 assert(hMiddle >= 0); 00937 } 00938 00950 inline RiemannStructure determineRiemannStructure( const T &i_hLeft, const T &i_hRight, 00951 const T &i_uLeft, const T &i_uRight ) const { 00952 T hMin = std::min(i_hLeft, i_hRight); 00953 T hMax = std::max(i_hLeft, i_hRight); 00954 00955 T uDif = i_uRight - i_uLeft; 00956 00957 if( 0 <= (T)2.0*(std::sqrt(g*hMin) - std::sqrt(g*hMax)) + uDif) 00958 return RarefactionRarefaction; 00959 00960 if( (hMax - hMin) * std::sqrt(g*(T)0.5*( 1/hMax + 1/hMin )) + uDif <= 0 ) 00961 return ShockShock; 00962 00963 if( i_hLeft < i_hRight ) 00964 return ShockRarefaction; 00965 00966 return RarefactionShock; 00967 } 00968 00977 static inline void solveLinearEquation( const T i_matrix[3][3], 00978 const T i_b[3], 00979 T o_x[3] ) { 00980 //#if 1 00981 // compute inverse of 3x3 matrix 00982 const T m[3][3] = 00983 { { (i_matrix[1][1]*i_matrix[2][2] - i_matrix[1][2]*i_matrix[2][1]), -(i_matrix[0][1]*i_matrix[2][2] - i_matrix[0][2]*i_matrix[2][1]), (i_matrix[0][1]*i_matrix[1][2] - i_matrix[0][2]*i_matrix[1][1]) }, 00984 { -(i_matrix[1][0]*i_matrix[2][2] - i_matrix[1][2]*i_matrix[2][0]), (i_matrix[0][0]*i_matrix[2][2] - i_matrix[0][2]*i_matrix[2][0]), -(i_matrix[0][0]*i_matrix[1][2] - i_matrix[0][2]*i_matrix[1][0]) }, 00985 { (i_matrix[1][0]*i_matrix[2][1] - i_matrix[1][1]*i_matrix[2][0]), -(i_matrix[0][0]*i_matrix[2][1] - i_matrix[0][1]*i_matrix[2][0]), (i_matrix[0][0]*i_matrix[1][1] - i_matrix[0][1]*i_matrix[1][0]) } 00986 }; 00987 T d = (i_matrix[0][0]*m[0][0] + i_matrix[0][1]*m[1][0] + i_matrix[0][2]*m[2][0]); 00988 00989 #ifndef NDEBUG 00990 if (std::fabs(d) < 0.000000001) { 00991 std::cerr << "Division close to zero!" << std::endl; 00992 std::cout << "Matrix: " << std::endl; 00993 std::cout << i_matrix[0][0] << ", " << i_matrix[0][1] << ", " << i_matrix[0][2] << std::endl; 00994 std::cout << i_matrix[1][0] << ", " << i_matrix[1][1] << ", " << i_matrix[1][2] << std::endl; 00995 std::cout << i_matrix[2][0] << ", " << i_matrix[2][1] << ", " << i_matrix[2][2] << std::endl; 00996 std::cout << "b: " << std::endl; 00997 std::cout << i_b[0] << ", " << i_b[1] << ", " << i_b[2] << std::endl; 00998 std::cout << "d: " << d << std::endl; 00999 assert(false); 01000 } 01001 #endif 01002 01003 // m stores not really the inverse matrix, but the inverse multiplied by d 01004 T s = (T)1/d; 01005 01006 // compute m*i_b 01007 o_x[0] = (m[0][0]*i_b[0] + m[0][1]*i_b[1] + m[0][2]*i_b[2])*s; 01008 o_x[1] = (m[1][0]*i_b[0] + m[1][1]*i_b[1] + m[1][2]*i_b[2])*s; 01009 o_x[2] = (m[2][0]*i_b[0] + m[2][1]*i_b[1] + m[2][2]*i_b[2])*s; 01010 01011 // return; 01012 //#else 01013 // T origDet = computeDeterminant(i_matrix); 01014 // 01015 // if (std::fabs(origDet) > zeroTol) { 01016 // T modifiedMatrix[3][3]; 01017 // 01018 // for (int column = 0; column < 3; column++) 01019 // { 01020 // memcpy(modifiedMatrix, i_matrix, sizeof(T)*3*3); 01021 // 01022 // //set a column of the matrix to i_b 01023 // modifiedMatrix[0][column] = i_b[0]; 01024 // modifiedMatrix[1][column] = i_b[1]; 01025 // modifiedMatrix[2][column] = i_b[2]; 01026 // 01027 // o_x[column] = computeDeterminant(modifiedMatrix) / origDet; 01028 // } 01029 // } 01030 // else { 01031 // std::cerr << "Warning: Linear dependent eigenvectors! (using Jacobi Solver)" << std::endl; 01032 // std::cerr << "Determinant: " << origDet << std::endl; 01033 // std::cerr << i_matrix[0][0] << "\t" << i_matrix[0][1] << "\t" << i_matrix[0][2] << std::endl; 01034 // std::cerr << i_matrix[1][0] << "\t" << i_matrix[1][1] << "\t" << i_matrix[1][2] << std::endl; 01035 // std::cerr << i_matrix[2][0] << "\t" << i_matrix[2][1] << "\t" << i_matrix[2][2] << std::endl; 01036 //#if 0 01037 // T xTemp[3]; 01038 // xTemp[0] = xTemp[1] = xTemp[2] = 0.; 01039 // for (int m = 0; m < 20; m++) { 01040 // for (int row = 0; row < 3; row++) { 01041 // o_x[row] = 0.; 01042 // for (int col = 0; col < 3; col++) { 01043 // if (col != row) 01044 // o_x[row] += i_matrix[row][col]*xTemp[col]; 01045 // } 01046 // o_x[row] = (i_b[row] - o_x[row])/i_matrix[row][row]; 01047 // } 01048 // 01049 // if (fabs(o_x[0]-xTemp[0]) + fabs(o_x[1]-xTemp[1]) + fabs(o_x[2]-xTemp[2]) < zeroTol*10.) 01050 // break; 01051 // else { 01052 // std::cout << "Error: " << fabs(o_x[0]-xTemp[0]) + fabs(o_x[1]-xTemp[1]) + fabs(o_x[2]-xTemp[2]) << std::endl; 01053 // xTemp[0] = o_x[0]; 01054 // xTemp[1] = o_x[1]; 01055 // xTemp[2] = o_x[2]; 01056 // } 01057 // } 01058 // std::cout << "Solution:" << std::endl; 01059 // std::cout << "\t" << o_x[0] << std::endl; 01060 // std::cout << "x=\t" << o_x[1] << std::endl; 01061 // std::cout << "\t" << o_x[2] << std::endl; 01062 // std::cout << "*********" << std::endl; 01063 // std::cout << "\t" << i_b[0] << std::endl; 01064 // std::cout << "b=\t" << i_b[1] << std::endl; 01065 // std::cout << "\t" << i_b[2] << std::endl; 01066 // std::cout << "***A*x****" << std::endl; 01067 // for (int row = 0; row < 3; row++) { 01068 // std::cout << "\t" << 01069 // i_matrix[row][0] * o_x[0] + 01070 // i_matrix[row][1] * o_x[1] + 01071 // i_matrix[row][2] * o_x[2] << std::endl; 01072 // } 01073 //#endif 01074 // assert(false); 01075 // } 01076 //#endif 01077 } 01078 //#if 0 01079 // /** 01080 // * Compute the determinant of a 3x3 matrix 01081 // * using formula: |A|=a11 * (a22a33-a32a23) + a12 * (a23a31-a33a21) + a13 * (a21a32-a31a22) 01082 // * |A|=a00 * (a11a22-a21a12) + a01 * (a12a20-a22a10) + a02 * (a10a21-a20a11) 01083 // */ 01084 // T computeDeterminant(T matrix[3][3]) const 01085 // { 01086 // 01087 // T determinant = matrix[0][0] * ( matrix[1][1] * matrix[2][2] - matrix[2][1] * matrix[1][2] ) 01088 // + matrix[0][1] * ( matrix[1][2] * matrix[2][0] - matrix[2][2] * matrix[1][0] ) 01089 // + matrix[0][2] * ( matrix[1][0] * matrix[2][1] - matrix[2][0] * matrix[1][1] ); 01090 // 01091 // return determinant; 01092 // } 01093 //#endif 01094 01095 01096 #undef sqrt_g_hLeft 01097 #undef sqrt_g_hRight 01098 #undef sqrt_hLeft 01099 #undef sqrt_hRight 01100 01101 }; 01102 01103 #undef complexsteadystatewave 01104 #undef correctrarefactions 01105 01106 01107 #endif /* AUGRIE_HPP_ */
1.7.6.1
