ConstitutiveRelation.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 _CONSTITUTIVERELATION_H_
#define _CONSTITUTIVERELATION_H_

#include "LevelData.H"
#include "FArrayBox.H"
#include "FluxBox.H"
#include "LevelSigmaCS.H"
#include "CellToEdge.H"
#include "NamespaceHeader.H"




class 
ConstitutiveRelation
{

public:

  static ConstitutiveRelation* parse(const char* a_prefix);

  ConstitutiveRelation() {;}

  virtual ~ConstitutiveRelation() {;}
 

  virtual void computeMu(LevelData<FArrayBox>& a_mu,
                         const LevelData<FArrayBox>& a_vel, const Real& a_scale,
                         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 = 0;

  

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


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


  virtual ConstitutiveRelation* getNewConstitutiveRelation() const = 0;

  void computeStrainRateInvariant(LevelData<FArrayBox>& a_epsilonSquared,
                                  const LevelData<FArrayBox>& a_velocity,
                                  const LevelData<FArrayBox>* a_crseVel,
                                  int nRefCrse,
                                  const LevelSigmaCS& a_coordSys,
                                  const IntVect& a_ghostVect = IntVect::Zero) const;

  void computeStrainRateInvariantFace(LevelData<FluxBox>& a_epsilonSquared,
                                      LevelData<FArrayBox>& a_velocity,
                                      const LevelData<FArrayBox>* a_crseVel,
                                      int a_nRefCrse,
                                      const LevelSigmaCS& a_coordSys,
                                      const IntVect& a_ghostVect = IntVect::Zero) const;  

  void computeStrainRateInvariant(LevelData<FArrayBox>& a_epsilonSquared,
                                  LevelData<FArrayBox>& a_gradVelocity,
                                  const LevelData<FArrayBox>& a_velocity,
                                  const LevelData<FArrayBox>* a_crseVel,
                                  int a_nRefCrse,
                                  const LevelSigmaCS& a_coordSys,
                                  const IntVect& a_ghostVect = IntVect::Zero) const;

  void computeStrainRateInvariantFace(LevelData<FluxBox>& a_epsilonSquared,
                                      LevelData<FluxBox>& a_gradVelocity,
                                      LevelData<FArrayBox>& a_velocity,
                                      const LevelData<FArrayBox>* a_crseVel,
                                      int a_nRefCrse,
                                      const LevelSigmaCS& a_coordSys,
                                      const IntVect& a_ghostVect = IntVect::Zero) const;  

protected:

};
 
class RateFactor 
{

public:
  
  virtual ~RateFactor(){;}

  virtual void computeA(FArrayBox& a_A, 
                        const FArrayBox& a_thetaStar,
                        const FArrayBox& a_pressure,
                        const Box& a_box) const = 0;

  virtual RateFactor* getNewRateFactor() const = 0;


} ;


class ConstantRateFactor : public RateFactor  
{
  Real m_A; // the constant rate factor
public:
  
  // ConstantRateFactor(Real a_seconds_per_unit_time)
  //   : m_A(9.2e-18 / a_seconds_per_unit_time ) // A measured in Pa^{n} a^{-1}
  // {
  // }
 
  ConstantRateFactor(Real a_A) : m_A(a_A)
  {
  }


  void computeA(FArrayBox& a_A, 
                const FArrayBox& a_theta, 
                const FArrayBox& a_pressure,
                const Box& a_box) const
  {
    a_A.setVal(m_A, a_box, 0);
  }

  RateFactor* getNewRateFactor() const 
  {
    ConstantRateFactor* newPtr = new ConstantRateFactor(m_A);
    return static_cast<RateFactor*>(newPtr);
  }


};


class ArrheniusRateFactor : public RateFactor  
{
  Real m_n;
  Real m_enhance;
  Real m_B0;
  Real m_theta_r;
  Real m_K;
  Real m_C;
  Real m_R;
  // (insert funny comment about lazy creeps here...)
  Real m_Q;

public:
  ArrheniusRateFactor(Real a_seconds_per_unit_time);
 
  void setDefaultParameters(Real a_seconds_per_unit_time);
 
  
  void setParameters(Real a_n,
                     Real a_enhance,
                     Real a_B0,
                     Real a_theta_r,
                     Real a_K,
                     Real a_C,
                     Real a_R,
                     Real a_Q);

  //compute the (temperature-dependent) rate factor A
  void computeA(FArrayBox& a_A, 
                const FArrayBox& a_thetaStar,
                const FArrayBox& a_pressure,
                const Box& a_box) const;

  RateFactor* getNewRateFactor() const; 
  

};


class PatersonRateFactor : public RateFactor  
{

public:
  Real m_E;
  Real m_A0;
  Real m_T0;
  Real m_R;
  // (insert funny comment about two lazy creeps here...)
  Real m_Qm, m_Qp;

public:

  PatersonRateFactor(Real a_seconds_per_unit_time);
 
  void setDefaultParameters(Real a_seconds_per_unit_time);
 
  void setA0(Real a_A0)
  {
    m_A0 = a_A0;
  }

  void setParameters(Real a_E, Real a_A0, Real a_T0, Real a_R, 
                     Real a_Qm, Real a_Qp);

  //compute the (pressure-corrected temperature- 
  //and pressure- dependent) rate factor A
  void computeA(FArrayBox& a_A, 
                const FArrayBox& a_thetaPC, 
                const FArrayBox& a_pressure,
                const Box& a_box) const;

  RateFactor* getNewRateFactor() const; 
  

};



class ZwingerRateFactor : public RateFactor  
{
  
  Real m_E;
  Real m_A0;
  Real m_T0;
  Real m_R;
  // (insert funny comment about two lazy creeps here...)
  Real m_Qm, m_Qp;

public:

  ZwingerRateFactor(Real a_seconds_per_unit_time);
 
  void setDefaultParameters(Real a_seconds_per_unit_time);
 

  void setParameters(Real a_E, Real a_A0, Real a_T0, Real a_R, 
                     Real a_Qm, Real a_Qp);

  //compute the (pressure-corrected temperature- 
  //and pressure- dependent) rate factor A
  void computeA(FArrayBox& a_A, 
                const FArrayBox& a_thetaPC, 
                const FArrayBox& a_pressure,
                const Box& a_box) const;

  RateFactor* getNewRateFactor() const; 
  

};

// $\f \tau^2 \tau _ij = 2 A \dot{\epsilon}_{ij} \f$ 



class GlensFlowRelation: public ConstitutiveRelation
{
public:

  GlensFlowRelation();

  virtual ~GlensFlowRelation();



  virtual void computeMu(LevelData<FArrayBox>& a_mu,
                         const LevelData<FArrayBox>& a_vel,  const Real& a_scale,
                         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 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_crseVel,
                             int a_nRefCrse,
                             const LevelData<FluxBox>& a_A, 
                             const LevelSigmaCS& a_coordSys,
                             const ProblemDomain& a_domain,
                             const IntVect& a_ghostVect = IntVect::Zero) const;


  virtual void setDefaultParameters();

  virtual void setParameters(Real a_n,
                             Real a_epsSqr0, 
                             Real a_delta);

  virtual ConstitutiveRelation* getNewConstitutiveRelation() const;
  
  Real m_n;

  virtual Real power() const {return m_n;}

  //RateFactor* m_rateFactor;

  Real m_epsSqr0, m_delta;

protected:
  void computeMu0(FArrayBox& a_mu0,
                  const FArrayBox& a_A,const Real& a_scale,
                  const Box& a_box) const;
  

};



class constMuRelation: public ConstitutiveRelation
{
public:
  constMuRelation();

  virtual ~constMuRelation();


  virtual void computeMu(LevelData<FArrayBox>& a_mu,
                         const LevelData<FArrayBox>& a_vel, const Real& a_scale,
                         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 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_crseVel,
                             int a_nRefCrse, 
                             const LevelData<FluxBox>& a_A, 
                             const LevelSigmaCS& a_coordSys,
                             const ProblemDomain& a_domain,
                             const IntVect& a_ghostVect = IntVect::Zero) const;



  virtual void setDefaultParameters();

  virtual void setConstVal(Real a_mu) {m_mu = a_mu;}

  virtual ConstitutiveRelation* getNewConstitutiveRelation() const;

  // constant value of mu
  Real m_mu;

  virtual Real power() const{return 1.0;}

protected:

};



#include "NamespaceFooter.H"

#endif