L1L2ConstitutiveRelation.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 _L1L2CONSTITUTIVERELATION_H_
#define _L1L2CONSTITUTIVERELATION_H_

#include "ConstitutiveRelation.H"
#include "IceThicknessIBC.H"
#include "NamespaceHeader.H"


#if BISICLES_Z == BISICLES_LAYERED

class L1L2ConstitutiveRelation : public ConstitutiveRelation
{

  GlensFlowRelation glensFlowRelation;

public:
  L1L2ConstitutiveRelation() 
    :m_solverTol(1.0e-6),
     m_additionalVelocitySIAGradSLimit(1.0e-2),
     m_effectiveViscositySIAGradSLimit(1.0e+10),
     m_additionalVelocitySIAOnly(true),
     m_startFromAnalyticMu(false),
     m_layerCoarsening(0)
  {;}
  
  void parseParameters();

  virtual ~L1L2ConstitutiveRelation() {;}

  GlensFlowRelation* getGlensFlowRelationPtr()
  { return &glensFlowRelation;}
  GlensFlowRelation& getGlensFlowRelation()
  { return glensFlowRelation;}

  virtual Real power() const 
  { return glensFlowRelation.power(); }

  virtual ConstitutiveRelation* getNewConstitutiveRelation() const;


  virtual void computeMu(LevelData<FArrayBox>& a_mu,
                         const LevelData<FArrayBox>& a_vel,  const Real& a_scale,
                         const LevelData<FArrayBox>* a_crseVelPtr,
                         int a_nRefCrse,
                         const LevelData<FArrayBox>& a_A,
                         const LevelSigmaCS& a_coordSys,
                         const ProblemDomain& a_domain,
                         const IntVect& a_ghostVect = IntVect::Zero) const;
  
 

  virtual void computeDissipation(LevelData<FArrayBox>& a_dissipation,
                                  const LevelData<FArrayBox>& a_vel, 
                                  const LevelData<FArrayBox>* a_crseVel,
                                  int a_nRefCrse,
                                  const LevelData<FArrayBox>& a_A,
                                  const LevelSigmaCS& a_coordSys,
                                  const ProblemDomain& a_domain,
                                  const IntVect& a_ghostVect = IntVect::Zero) const;
  

  virtual void computeFaceMu(LevelData<FluxBox>& a_mu,
                             LevelData<FArrayBox>& a_vel,  const Real& a_scale,
                             const LevelData<FArrayBox>* a_crseVelPtr,
                             int nRefCrse,
                             const LevelData<FluxBox>& a_A, 
                             const LevelSigmaCS& a_coordSys,
                             const ProblemDomain& a_domain,
                             const IntVect& a_ghostVect = IntVect::Zero) const;




  virtual void modifyTransportCoefficients
  (const LevelData<FArrayBox>& a_cellVel,
   const LevelData<FArrayBox>* a_crseVelPtr,
   const LevelData<FArrayBox>* a_crseDiffusivityPtr,
   int a_nRefCrse,
   const LevelSigmaCS& a_coordSys,
   const DisjointBoxLayout& a_grids,
   const ProblemDomain& a_domain,
   const LevelData<FArrayBox>& a_A,
   const LevelData<FArrayBox>& a_sA,
   const LevelData<FArrayBox>& a_bA,
   LevelData<FluxBox>& a_faceVelAdvection,
   LevelData<FluxBox>& a_faceVelTotal,
   LevelData<FluxBox>& a_faceDiffusivity,
   LevelData<FArrayBox>& a_cellDiffusivity,
   LevelData<FluxBox>& a_layerXYFaceXYVel,
   LevelData<FArrayBox>& a_layerSFaceXYVel) const;
    

  void solverTolerance(Real a_solverTol) 
  { m_solverTol = a_solverTol; }


private:
  virtual void computeMuZ(LevelData<FArrayBox>& a_mu,
                          LevelData<FArrayBox>& a_gradVel,
                          LevelData<FArrayBox>& a_epsSqr,
                          const Vector<Real>& a_sigma,
                          const LevelData<FArrayBox>& a_vel, const Real& a_scale,
                          const LevelData<FArrayBox>* a_crseVelPtr,
                          int a_nRefCrse,
                          const LevelData<FArrayBox>& a_A,
                          const LevelSigmaCS& a_coordSys,
                          const ProblemDomain& a_domain,
                          const IntVect& a_ghostVect = IntVect::Unit) const;
  
  

  //compute mu for given sigma at a_faceDir cell faces
  virtual void computeFaceMuZ(LevelData<FluxBox>& a_mu,
                              LevelData<FluxBox>& a_gradVel,
                              LevelData<FluxBox>& a_epsSqr,
                              const Vector<Real>& a_sigma,
                              LevelData<FArrayBox>& a_vel, const Real& a_scale,
                              const LevelData<FArrayBox>* a_crseVelPtr,
                              int a_nRefCrse,
                              const LevelData<FluxBox>& a_A,
                              const LevelSigmaCS& a_coordSys,
                              const ProblemDomain& a_domain,
                              const IntVect& a_ghostVect) const;
  

  virtual void computeEitherMuZ(FArrayBox& a_mu,
                                const Vector<Real>& a_sigma,
                                const FArrayBox& a_grads,
                                const FArrayBox& a_epsSqr,
                                const FArrayBox& a_A,const Real& a_scale,
                                const FArrayBox& a_H,
                                const Real& a_rhog,
                                const Box& a_box,
                                int a_layerCoarsen) const;


  virtual void computeFaceFluxVelocity(LevelData<FArrayBox>& a_vel,
                                       const LevelData<FArrayBox>* a_crseVelPtr,
                                       int a_nRefCrse,
                                       const LevelData<FArrayBox>& a_A,
                                       const LevelData<FArrayBox>& a_thickness,
                                       const RealVect& a_dx,
                                       LevelData<FluxBox>& a_fluxVel,
                                       LevelData<FluxBox>& a_layerXYFaceXYVel,
                                       LevelData<FArrayBox>& a_layerSFaceXYVel,
                                       const LevelSigmaCS& a_coordSys,
                                       const ProblemDomain& a_domain,
                                       const IntVect& a_cellGhost = IntVect::Zero) const;

 
  // tolerance for the secant solve when computing mu
  Real m_solverTol;
  // if > 0, limit the magnitude of grad(s) when computing the SIA part of u(z) - u(z=base)
  Real m_additionalVelocitySIAGradSLimit; 
  // if > 0, limit the magnitude of grad(s) when computing the effective viscosity
  Real m_effectiveViscositySIAGradSLimit;
  // neglect non SIA contribution to u(z) - u(z=base)  
  bool m_additionalVelocitySIAOnly; 

  // for n = 3, there is an analytic solution to f(mu) = 0; sometimes
  // it is preferable to start from that value
  bool m_startFromAnalyticMu;

  // Set m_layerCoarsening > 0 to skip layers in expensive calculations and interpolate instead
  int m_layerCoarsening;
  

};
#endif 
#include "NamespaceFooter.H"
#endif