WavePropagation1d.h

Go to the documentation of this file.

 
#ifndef TSUNAMI_LAB_PATCHES_WAVE_PROPAGATION_1D
#define TSUNAMI_LAB_PATCHES_WAVE_PROPAGATION_1D
 
#include "WavePropagation.h"
 
namespace tsunami_lab {
  namespace patches {
    template<class T>
    class WavePropagation1d;
  }
}
 
template<class T>
class tsunami_lab::patches::WavePropagation1d: public WavePropagation<T> {
  private:
    unsigned short m_step = 0;
 
    t_idx m_nCells = 0;
 
    t_real * m_h[2] = { nullptr, nullptr };
 
    t_real * m_hu[2] = { nullptr, nullptr };
 
    // bathymetry array
    t_real * m_bm = nullptr;
 
  public:
    WavePropagation1d( t_idx i_nCells ) : WavePropagation<T>() {
      m_nCells = i_nCells;
 
      // allocate memory including a single ghost cell on each side
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        m_h[l_st] = new t_real[  m_nCells + 2 ];
        m_hu[l_st] = new t_real[ m_nCells + 2 ];
      }
 
      // init to zero (including both ghost cells at indices 0 and m_nCells+1)
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        for( t_idx l_ce = 0; l_ce < m_nCells + 2; l_ce++ ) {
          m_h[l_st][l_ce] = 0;
          m_hu[l_st][l_ce] = 0;
        }
      }
 
      // init bathymetry array to 0 (value-initialised)
      m_bm = new t_real[m_nCells + 2]();
    }
 
    ~WavePropagation1d() {
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        delete[] m_h[l_st];
        delete[] m_hu[l_st];
      }
      delete[] m_bm;
    }
 
    void timeStep( t_real i_scaling ) {
      // pointers to old and new data
      t_real * l_hOld = m_h[m_step];
      t_real * l_huOld = m_hu[m_step];
 
      m_step = (m_step+1) % 2;
      t_real * l_hNew =  m_h[m_step];
      t_real * l_huNew = m_hu[m_step];
 
      // init new cell quantities
      for( t_idx l_ce = 1; l_ce < m_nCells+1; l_ce++ ) {
        l_hNew[l_ce] = l_hOld[l_ce];
        l_huNew[l_ce] = l_huOld[l_ce];
      }
 
      // iterate over edges and update with Riemann solutions
      for( t_idx l_ed = 0; l_ed < m_nCells+1; l_ed++ ) {
        // determine left and right cell-id
        t_idx l_ceL = l_ed;
        t_idx l_ceR = l_ed+1;
 
        t_real l_hL_ref = l_hOld[l_ceL];
        t_real l_hR_ref = l_hOld[l_ceR];
        t_real l_huL_ref = l_huOld[l_ceL];
        t_real l_huR_ref = l_huOld[l_ceR];
        t_real l_bmL_ref = m_bm[l_ceL];
        t_real l_bmR_ref = m_bm[l_ceR];
 
        // wet/dry classification by water column: h > 0 is wet, h == 0 is dry
        const bool l_dryL = l_hL_ref <= 0;
        const bool l_dryR = l_hR_ref <= 0;
 
        // dry-dry edge: Riemann problem is undefined, no update needed
        if( l_dryL && l_dryR ) continue;
 
        // wet-dry interface: replace the dry side with a mirror of the wet side
        // (Spec 3.2.1) so the edge acts as a solid wall
        if( l_dryL ) {
          l_hL_ref  = l_hR_ref;
          l_huL_ref = -l_huR_ref;
          l_bmL_ref = l_bmR_ref;
        } else if( l_dryR ) {
          l_hR_ref  = l_hL_ref;
          l_huR_ref = -l_huL_ref;
          l_bmR_ref = l_bmL_ref;
        }
 
        // compute net-updates
        t_real l_netUpdates[2][2];
 
        this->solver.netUpdates( l_hL_ref,
                                  l_hR_ref,
                                  l_huL_ref,
                                  l_huR_ref,
                                  l_bmL_ref,
                                  l_bmR_ref,
                                  l_netUpdates[0],
                                  l_netUpdates[1] );
 
        // only apply updates to wet cells; Spec 3.2 requires dry cells to stay dry
        if( !l_dryL ) {
          l_hNew[l_ceL]  -= i_scaling * l_netUpdates[0][0];
          l_huNew[l_ceL] -= i_scaling * l_netUpdates[0][1];
        }
        if( !l_dryR ) {
          l_hNew[l_ceR]  -= i_scaling * l_netUpdates[1][0];
          l_huNew[l_ceR] -= i_scaling * l_netUpdates[1][1];
        }
      }
    }
 
    void setGhostBoundary( BoundaryCondition i_left,
                           BoundaryCondition i_right,
                           BoundaryCondition /*i_bottom*/,
                           BoundaryCondition /*i_top*/ ) {
      t_real * l_h  = m_h[m_step];
      t_real * l_hu = m_hu[m_step];
 
      // left ghost cell mirrors cell 1
      l_h[0]  = l_h[1];
      m_bm[0] = m_bm[1];
      l_hu[0] = (i_left == BoundaryCondition::REFLECTING) ? -l_hu[1] : l_hu[1];
 
      // right ghost cell mirrors cell m_nCells
      l_h[m_nCells + 1]  = l_h[m_nCells];
      m_bm[m_nCells + 1] = m_bm[m_nCells];
      l_hu[m_nCells + 1] = (i_right == BoundaryCondition::REFLECTING)
        ? -l_hu[m_nCells]
        : l_hu[m_nCells];
    }
 
    void setGhostOutflow() {
      setGhostBoundary( BoundaryCondition::OUTFLOW,  BoundaryCondition::OUTFLOW,
                        BoundaryCondition::OUTFLOW,  BoundaryCondition::OUTFLOW );
    }
 
    void setGhostReflection() {
      setGhostBoundary( BoundaryCondition::REFLECTING, BoundaryCondition::REFLECTING,
                        BoundaryCondition::REFLECTING, BoundaryCondition::REFLECTING );
    }
 
    t_idx getStride(){
      return m_nCells+2;
    }
 
    t_real const * getHeight(){
      return m_h[m_step]+1;
    }
 
    t_real const * getMomentumX(){
      return m_hu[m_step]+1;
    }
 
    t_real const * getMomentumY(){
      return nullptr;
    }
 
    t_real const *getBathymetry() {
        return m_bm + 1;
    }
 
    void setHeight( t_idx  i_ix,
                    t_idx,
                    t_real i_h ) {
      m_h[m_step][i_ix+1] = i_h;
    }
 
    void setMomentumX( t_idx  i_ix,
                       t_idx,
                       t_real i_hu ) {
      m_hu[m_step][i_ix+1] = i_hu;
    }
 
    void setMomentumY( t_idx,
                       t_idx,
                       t_real ) {};
 
    
    void setBathymetry(t_idx i_ix, t_idx /*i_iy*/, t_real i_b) {
      m_bm[i_ix + 1] = i_b;
    };
};
 
#endif