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 };
 
  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;
        }
      }
    }
 
    ~WavePropagation1d() {
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        delete[] m_h[l_st];
        delete[] m_hu[l_st];
      }
    }
 
    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;
 
        // compute net-updates
        t_real l_netUpdates[2][2];
 
        this->solver.netUpdates( l_hOld[l_ceL],
                                  l_hOld[l_ceR],
                                  l_huOld[l_ceL],
                                  l_huOld[l_ceR],
                                  l_netUpdates[0],
                                  l_netUpdates[1] );
 
        // update the cells' quantities
        l_hNew[l_ceL]  -= i_scaling * l_netUpdates[0][0];
        l_huNew[l_ceL] -= i_scaling * l_netUpdates[0][1];
 
        l_hNew[l_ceR]  -= i_scaling * l_netUpdates[1][0];
        l_huNew[l_ceR] -= i_scaling * l_netUpdates[1][1];
      }
    }
 
    void setGhostOutflow() {
      t_real * l_h = m_h[m_step];
      t_real * l_hu = m_hu[m_step];
 
      // set left boundary
      l_h[0] = l_h[1];
      l_hu[0] = l_hu[1];
 
      // set right boundary
      l_h[m_nCells+1] = l_h[m_nCells];
      l_hu[m_nCells+1] = l_hu[m_nCells];
    }
 
    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;
    }
 
    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 ) {};
};
 
#endif