LevelSigmaCS.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef _LEVELSIGMACS_H_
#define _LEVELSIGMACS_H_

#include "DisjointBoxLayout.H"
#include "LevelData.H"
#include "FArrayBox.H"
#include "NodeFArrayBox.H"
#include "FluxBox.H"
#include "RealVect.H"
#include "Interval.H"
#include "AMRIO.H"


#include "NamespaceHeader.H"

#define BISICLES_LAYERED 0
#define BISICLES_FULLZ 1
#if CH_SPACEDIM == 2
// for now, BISICLES_LAYERED is the same as 2D, but
// if we one day look at flow line models there could a 1D variant
#define BISICLES_Z 0
#elif CH_SPACEDIM == 3
// for now, BISICLES_FULLZ is the same as 3D, but
// if we one day look at flow line models there could a 2D variant
#define BISICLES_Z 1
#elif CH_SPACEDIM == 1
// for now, same as 2D
// for now, BISICLES_LAYERED is the same as 2D, but
// if we one day look at flow line models there could a 1D variant
#define BISICLES_Z 0
#endif


class LevelSigmaCS
{
public:
  LevelSigmaCS();

  LevelSigmaCS(const DisjointBoxLayout& a_grids,
               const RealVect& a_dx, 
               const IntVect& a_ghostVect = IntVect::Unit);

  LevelSigmaCS( const DisjointBoxLayout& a_grids,
                const RealVect& a_dx,
                const LevelSigmaCS& a_fineCS, 
                int a_nRef
                );
  virtual ~LevelSigmaCS();

  void define(const DisjointBoxLayout& a_grids,
              const RealVect& a_dx, 
              const IntVect& a_ghostVect = IntVect::Unit);

  void define(const LevelSigmaCS& a_fineCS, 
              int a_nRef);


  RealVect realCoord(const RealVect& a_Xi, const DataIndex& a_index) const;
  

  RealVect mappedCoord(const RealVect& a_x, const DataIndex& a_index) const;


  void getNodeRealCoordinates(LevelData<NodeFArrayBox>& a_nodeCoords) const;


  // return cell spacing dx
  RealVect dx() const {return m_dx;}

  LevelData<FArrayBox>& getH(); 

  const LevelData<FArrayBox>& getH() const;
  
  LevelData<FluxBox>& getFaceH(); 

  const LevelData<FluxBox>& getFaceH() const;

  const LevelData<FArrayBox>& getTopography() const;

  LevelData<FArrayBox>& getTopography();


  void setTopography(const LevelData<FArrayBox>& a_topography);

  //void setUnitShift(const RealVect& a_unitShift)
  //{
  //  m_unitShift = a_unitShift;
  //}



  void setBackgroundSlope(const RealVect& a_backgroundSlope)
  {
    m_backgroundSlope = a_backgroundSlope;
  }

  const RealVect& getBackgroundSlope() const
  {
    return m_backgroundSlope;
  }


  //sets the surface elevation (at cell centres)
  void setSurfaceHeight(const LevelData<FArrayBox>& a_surface);
  //sets the surface gradient (at cell centres)
  void setGradSurface(const LevelData<FArrayBox>& a_gradSurface);
  //sets the surface gradient (at cell faces)
  void setGradSurfaceFace(const LevelData<FluxBox>& a_gradSurface);
  
  
  //returns the surface elevation (at cell centres)
  const LevelData<FArrayBox>& getSurfaceHeight() const ;
  //returns the surface gradient (at cell centres)
  const LevelData<FArrayBox>& getGradSurface() const;
  //returns the surface gradient (at cell faces)
  const LevelData<FluxBox>& getGradSurfaceFace() const;
  //returns the thickness over flotation (at cell centres)
  const LevelData<FArrayBox>& getThicknessOverFlotation() const;


  void getSurfaceHeight(LevelData<FArrayBox>& a_zSurface) const;


  const LevelData<BaseFab<int> >& getFloatingMask() const {return m_floatingMask;}


  void recomputeGeometry( const LevelSigmaCS* a_crseCoords, 
                          const int a_refRatio) ;
 

  void recomputeGeometryFace( const LevelSigmaCS* a_crseCoords, 
                              const int a_refRatio) ;


  const LevelData<FArrayBox>& deltaFactors() const;


  const LevelData<FluxBox>& faceDeltaFactors() const;


  Real seaLevel() const {return m_seaLevel;}

  void setSeaLevel(const Real& a_seaLevel) {m_seaLevel = a_seaLevel;}

  Real iceDensity() const {return m_iceDensity;}
  Real waterDensity() const {return m_waterDensity;}
  Real gravity() const {return m_gravity;}

  void setIceDensity(const Real& a_iceDensity) 
  {m_iceDensity = a_iceDensity;}
  void setWaterDensity(const Real& a_waterDensity) 
  {m_waterDensity = a_waterDensity;}
  void setGravity(const Real& a_gravity) 
  {m_gravity = a_gravity;}
  

  bool isDefined() const {return m_isDefined;}
  
  LevelSigmaCS* makeCoarser(int a_coarsenFactor) const;
                   
  const DisjointBoxLayout& grids() const {return m_grids;}

  const LayoutData<bool>&  anyFloating() const { return m_anyFloating;}


  enum ThicknessInterpolationMethod{
    STD_THICKNESS_INTERPOLATION_METHOD = 0,
    SMOOTH_SURFACE_THICKNESS_INTERPOLATION_METHOD = 1,
    MAX_THICKNESS_INTERPOLATION_METHOD = 2};

  /*** 
       Interpolating thickness and topography independently tends to
       move the grounding line, and results in steep surfaces there
       which are not seen in interpolation of the surface itself. 
       interpFromCoarse interpolates topography, the upper and lower surfaces
       then adjusts base (and so thickness) to maintain the grounding line.
       
       optionally, topography and thickness can be set elsewhere,
       and the mask can be recomputed, in either the valid region, or in
       both valid region and ghost region. 
 */
  void interpFromCoarse(const LevelSigmaCS& a_crseCoords,
                        const int a_refinementRatio, 
                        const bool a_interpolateTopography = true, 
                        const bool a_interpolateThickness = true, 
                        const bool a_preserveMask = true,
                        const bool a_interpolateTopographyGhost = true, 
                        const bool a_interpolateThicknessGhost = true, 
                        const bool a_preserveMaskGhost = true, 
                        int a_thicknessInterpolationMethod = STD_THICKNESS_INTERPOLATION_METHOD);





  // exchange topography accounting for the unit shift.
  void exchangeTopography();
  // exchange a_level accounting for the unit shift.
  void unitShiftExchange(LevelData<FArrayBox>& a_level);

  IntVect ghostVect() const {return m_ghostVect;}

protected:

  DisjointBoxLayout m_grids;

  DisjointBoxLayout m_Hgrids;

  // cell spacing
  RealVect m_dx;

  IntVect m_ghostVect;

  // for periodic problems :  a vector with components
  // (topography(x + Lx , y) - topography(x , y), 
  // (topography(x , y + Ly) - topography(x , y)
  RealVect m_unitShift;

  RealVect m_backgroundSlope;

  LevelData<FArrayBox> m_topography;

  LevelData<FArrayBox> m_H;

  LevelData<FluxBox> m_faceH;

  LevelData<FArrayBox> m_surface;

  LevelData<FArrayBox> m_gradSurface;

  LevelData<FluxBox> m_gradSurfaceFace;

  LevelData<FArrayBox> m_thicknessOverFlotation;

  LevelData<BaseFab<int> > m_floatingMask;

  LayoutData<bool> m_anyFloating; 
  LevelData<FArrayBox> m_deltaFactors;

  LevelData<FluxBox> m_faceDeltaFactors;

  bool m_isDefined;


  // sea level
  Real m_seaLevel;

  Real m_iceDensity;

  Real  m_waterDensity;

  Real  m_gravity;

  // set a reasonable set of default values
  void setDefaultValues();
  
  void computeDeltaFactors();

  void computeSurface(const LevelSigmaCS* a_crseCoords, 
                      const int a_refRatio);

  void computeFloatingMask(const LevelData<FArrayBox>& a_surface);


#if BISICLES_Z == BISICLES_LAYERED

public:
  // some poor man's multidim stuff that is only useful when we
  // are treating the problem as a stack of layers with constant
  // sigma at the boundaries which never get refined

  // vector of sigma values at the layer interfaces  
  const Vector<Real>& getFaceSigma() const
  { return m_faceSigma; }

  //vector of delta sigma's (at the layer midpoints)
  const Vector<Real>& getDSigma() const
  { return m_dSigma;}

  // vector of sigma values at the layer midpoints
  const Vector<Real>& getSigma() const
  { return m_sigma; }

  void setFaceSigma(const Vector<Real>& a_faceSigma)
  {
    m_faceSigma = a_faceSigma;
    m_sigma.resize(m_faceSigma.size()-1);
    m_dSigma.resize(m_faceSigma.size()-1);
    for (unsigned int i = 0; i <  m_sigma.size(); ++i)
      {
        //layer midpoint
        m_sigma[i] = 0.5 * (a_faceSigma[i] + a_faceSigma[i+1]);
        //layer delta sigma
        m_dSigma[i] = (a_faceSigma[i+1] - a_faceSigma[i]);
      }
  }

protected:

  Vector<Real> m_faceSigma;
  Vector<Real> m_sigma;
  Vector<Real> m_dSigma;

#endif


};


// specializations of IO routines for LayoutData<SigmaCS>

void WriteSigmaMappedUGHDF5(const string&               a_fileRoot,
                            const DisjointBoxLayout&    a_grids,
                            const LevelData<FArrayBox>& a_data,
                            const LevelSigmaCS&  a_CoordSys,
                            const Box& a_domainBox);

void WriteSigmaMappedUGHDF5(const string&               a_fileRoot,
                            const DisjointBoxLayout&    a_grids,
                            const LevelData<FArrayBox>& a_data,
                            const LevelSigmaCS&  a_CoordSys,
                            const Vector<string>& a_vectNames,
                            const Box& a_domainBox);

void
WriteSigmaMappedAMRHierarchyHDF5(const string& a_fileRoot,
                                 const Vector<DisjointBoxLayout>& a_vectGrids,
                                 const Vector<LevelData<FArrayBox>* > & a_vectData,
                                 const Vector<string>& a_vectNames,
                                 const Vector<const LevelSigmaCS* >& a_vectCoordSys,
                                 const Box& a_baseDomain,
                                 const Real& a_dt,
                                 const Real& a_time,
                                 const Vector<int>& a_vectRatio,
                                 const int& a_numLevels);

#include "NamespaceFooter.H"

#endif