WavePropagation2d.h

Go to the documentation of this file.

 
#ifndef TSUNAMI_LAB_PATCHES_WAVE_PROPAGATION_2D
#define TSUNAMI_LAB_PATCHES_WAVE_PROPAGATION_2D
 
#include "WavePropagation.h"
 
namespace tsunami_lab {
  namespace patches {
    template<class T>
    class WavePropagation2d;
  }
}
 
template<class T>
class tsunami_lab::patches::WavePropagation2d : public WavePropagation<T> {
  private:
    unsigned short m_step = 0;
 
    t_idx m_nx = 0;
 
    t_idx m_ny = 0;
 
    t_idx m_stride = 0;
 
    t_idx m_innerOffset = 0;
 
    t_real * m_h[2] = { nullptr, nullptr };
 
    t_real * m_hu[2] = { nullptr, nullptr };
 
    t_real * m_hv[2] = { nullptr, nullptr };
 
    t_real * m_bm = nullptr;
 
    inline t_idx idx( t_idx i_ix, t_idx i_iy ) const {
      return ( i_ix + 1 ) + ( i_iy + 1 ) * m_stride;
    }
 
  public:
    WavePropagation2d( t_idx i_nx, t_idx i_ny ) : WavePropagation<T>() {
      m_nx          = i_nx;
      m_ny          = i_ny;
      m_stride      = i_nx + 2;
      m_innerOffset = m_stride + 1;
 
      const t_idx l_total = m_stride * ( m_ny + 2 );
 
      // allocate state buffers (zero-initialised, including ghosts)
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        m_h [l_st] = new t_real[ l_total ]();
        m_hu[l_st] = new t_real[ l_total ]();
        m_hv[l_st] = new t_real[ l_total ]();
      }
 
      // bathymetry (zero-initialised)
      m_bm = new t_real[ l_total ]();
    }
 
    ~WavePropagation2d() {
      for( unsigned short l_st = 0; l_st < 2; l_st++ ) {
        delete[] m_h [l_st];
        delete[] m_hu[l_st];
        delete[] m_hv[l_st];
      }
      delete[] m_bm;
    }
 
    void timeStep( t_real i_scaling ) {
      // pointers to old (Q^n) and new (Q^{n+1}) state
      t_real * l_hOld  = m_h [m_step];
      t_real * l_huOld = m_hu[m_step];
      t_real * l_hvOld = m_hv[m_step];
 
      m_step = ( m_step + 1 ) % 2;
      t_real * l_hNew  = m_h [m_step];
      t_real * l_huNew = m_hu[m_step];
      t_real * l_hvNew = m_hv[m_step];
 
      const t_idx l_s = m_stride;
 
      // initialise Q^{n+1} = Q^n on the interior
      for( t_idx l_j = 1; l_j <= m_ny; l_j++ ) {
        for( t_idx l_i = 1; l_i <= m_nx; l_i++ ) {
          const t_idx l_id = l_i + l_j * l_s;
          l_hNew [l_id] = l_hOld [l_id];
          l_huNew[l_id] = l_huOld[l_id];
          l_hvNew[l_id] = l_hvOld[l_id];
        }
      }
 
      // x-sweep: vertical edges between cells (i, j) and (i+1, j)
      // for i ∈ [0, m_nx], j ∈ [1, m_ny] — touches the left/right ghosts
      for( t_idx l_j = 1; l_j <= m_ny; l_j++ ) {
        for( t_idx l_i = 0; l_i <= m_nx; l_i++ ) {
          const t_idx l_idL = l_i       + l_j * l_s;
          const t_idx l_idR = ( l_i+1 ) + l_j * l_s;
 
          t_real l_hL  = l_hOld [l_idL];
          t_real l_hR  = l_hOld [l_idR];
          t_real l_huL = l_huOld[l_idL];
          t_real l_huR = l_huOld[l_idR];
          t_real l_bL  = m_bm   [l_idL];
          t_real l_bR  = m_bm   [l_idR];
 
          const bool l_dryL = l_hL <= 0;
          const bool l_dryR = l_hR <= 0;
          if( l_dryL && l_dryR ) continue;
 
          // wet-dry edge: mirror dry side from wet partner (Spec 3.2.1)
          if( l_dryL ) {
            l_hL  = l_hR;  l_huL = -l_huR;  l_bL = l_bR;
          } else if( l_dryR ) {
            l_hR  = l_hL;  l_huR = -l_huL;  l_bR = l_bL;
          }
 
          t_real l_netUpdates[2][2];
          this->solver.netUpdates( l_hL,  l_hR,
                                   l_huL, l_huR,
                                   l_bL,  l_bR,
                                   l_netUpdates[0], l_netUpdates[1] );
 
          if( !l_dryL ) {
            l_hNew [l_idL] -= i_scaling * l_netUpdates[0][0];
            l_huNew[l_idL] -= i_scaling * l_netUpdates[0][1];
          }
          if( !l_dryR ) {
            l_hNew [l_idR] -= i_scaling * l_netUpdates[1][0];
            l_huNew[l_idR] -= i_scaling * l_netUpdates[1][1];
          }
        }
      }
 
      // y-sweep: horizontal edges between cells (i, j) and (i, j+1)
      // for i ∈ [1, m_nx], j ∈ [0, m_ny] — touches the bottom/top ghosts.
      // hv plays the role of the normal momentum here; hu is untouched.
      for( t_idx l_j = 0; l_j <= m_ny; l_j++ ) {
        for( t_idx l_i = 1; l_i <= m_nx; l_i++ ) {
          const t_idx l_idB = l_i + l_j       * l_s;     // below
          const t_idx l_idT = l_i + ( l_j+1 ) * l_s;     // above (top)
 
          t_real l_hB  = l_hOld [l_idB];
          t_real l_hT  = l_hOld [l_idT];
          t_real l_hvB = l_hvOld[l_idB];
          t_real l_hvT = l_hvOld[l_idT];
          t_real l_bB  = m_bm   [l_idB];
          t_real l_bT  = m_bm   [l_idT];
 
          const bool l_dryB = l_hB <= 0;
          const bool l_dryT = l_hT <= 0;
          if( l_dryB && l_dryT ) continue;
 
          if( l_dryB ) {
            l_hB  = l_hT;  l_hvB = -l_hvT;  l_bB = l_bT;
          } else if( l_dryT ) {
            l_hT  = l_hB;  l_hvT = -l_hvB;  l_bT = l_bB;
          }
 
          t_real l_netUpdates[2][2];
          this->solver.netUpdates( l_hB,  l_hT,
                                   l_hvB, l_hvT,
                                   l_bB,  l_bT,
                                   l_netUpdates[0], l_netUpdates[1] );
 
          if( !l_dryB ) {
            l_hNew [l_idB] -= i_scaling * l_netUpdates[0][0];
            l_hvNew[l_idB] -= i_scaling * l_netUpdates[0][1];
          }
          if( !l_dryT ) {
            l_hNew [l_idT] -= i_scaling * l_netUpdates[1][0];
            l_hvNew[l_idT] -= 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];
      t_real * l_hv = m_hv[m_step];
      const t_idx l_s = m_stride;
 
      // Note: the four corner ghost cells (i ∈ {0, m_nx+1} ∧ j ∈ {0, m_ny+1})
      // are intentionally not written — see class docstring. A future
      // CTU/transverse extension would have to fill them in addition.
 
      // left edge: i = 0 mirrors i = 1
      for( t_idx l_j = 1; l_j <= m_ny; l_j++ ) {
        l_h [          l_j*l_s] = l_h [1 + l_j*l_s];
        l_hu[          l_j*l_s] = ( i_left == BoundaryCondition::REFLECTING )
                                  ? -l_hu[1 + l_j*l_s] : l_hu[1 + l_j*l_s];
        l_hv[          l_j*l_s] = l_hv[1 + l_j*l_s];
        m_bm[          l_j*l_s] = m_bm[1 + l_j*l_s];
      }
 
      // right edge: i = m_nx + 1 mirrors i = m_nx
      for( t_idx l_j = 1; l_j <= m_ny; l_j++ ) {
        l_h [(m_nx+1) + l_j*l_s] = l_h [m_nx + l_j*l_s];
        l_hu[(m_nx+1) + l_j*l_s] = ( i_right == BoundaryCondition::REFLECTING )
                                   ? -l_hu[m_nx + l_j*l_s] : l_hu[m_nx + l_j*l_s];
        l_hv[(m_nx+1) + l_j*l_s] = l_hv[m_nx + l_j*l_s];
        m_bm[(m_nx+1) + l_j*l_s] = m_bm[m_nx + l_j*l_s];
      }
 
      // bottom edge: j = 0 mirrors j = 1
      for( t_idx l_i = 1; l_i <= m_nx; l_i++ ) {
        l_h [l_i              ] = l_h [l_i + l_s];
        l_hu[l_i              ] = l_hu[l_i + l_s];
        l_hv[l_i              ] = ( i_bottom == BoundaryCondition::REFLECTING )
                                  ? -l_hv[l_i + l_s] : l_hv[l_i + l_s];
        m_bm[l_i              ] = m_bm[l_i + l_s];
      }
 
      // top edge: j = m_ny + 1 mirrors j = m_ny
      for( t_idx l_i = 1; l_i <= m_nx; l_i++ ) {
        l_h [l_i + (m_ny+1)*l_s] = l_h [l_i + m_ny*l_s];
        l_hu[l_i + (m_ny+1)*l_s] = l_hu[l_i + m_ny*l_s];
        l_hv[l_i + (m_ny+1)*l_s] = ( i_top == BoundaryCondition::REFLECTING )
                                   ? -l_hv[l_i + m_ny*l_s] : l_hv[l_i + m_ny*l_s];
        m_bm[l_i + (m_ny+1)*l_s] = m_bm[l_i + m_ny*l_s];
      }
    }
 
    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_stride;
    }
 
    t_real const * getHeight() {
      return m_h[m_step] + m_innerOffset;
    }
 
    t_real const * getMomentumX() {
      return m_hu[m_step] + m_innerOffset;
    }
 
    t_real const * getMomentumY() {
      return m_hv[m_step] + m_innerOffset;
    }
 
    t_real const * getBathymetry() {
      return m_bm + m_innerOffset;
    }
 
    void setHeight( t_idx  i_ix,
                    t_idx  i_iy,
                    t_real i_h ) {
      m_h[m_step][ idx( i_ix, i_iy ) ] = i_h;
    }
 
    void setMomentumX( t_idx  i_ix,
                       t_idx  i_iy,
                       t_real i_hu ) {
      m_hu[m_step][ idx( i_ix, i_iy ) ] = i_hu;
    }
 
    void setMomentumY( t_idx  i_ix,
                       t_idx  i_iy,
                       t_real i_hv ) {
      m_hv[m_step][ idx( i_ix, i_iy ) ] = i_hv;
    }
 
    void setBathymetry( t_idx  i_ix,
                        t_idx  i_iy,
                        t_real i_b ) {
      m_bm[ idx( i_ix, i_iy ) ] = i_b;
    }
};
 
#endif