facefluxwithgravity.cpp

Go to the documentation of this file.

#include "facefluxwithgravity.h"
#include "numericmethods.h"

FaceFluxWithGravity::FaceFluxWithGravity(OrthoMesh &mesh,Function1D &fMobW,Function1D &DfMobW,Function1D &fMobGrav,Function1D &DfMobGrav,const VecDouble  &cK,double gc)
  :FaceFluxFunction(1),m_fMobW(fMobW),m_DfMobW(DfMobW),m_fMobGrav(fMobGrav),m_DfMobGrav(DfMobGrav),m_Vn(mesh.numFaces()),m_mesh(mesh),m_cK(cK)
{

  m_gc = gc;
  m_maxK = m_cK.linfty_norm();
}





void FaceFluxWithGravity::exactFluxAtFace(VecDouble &vFlux, const FaceInfo &face,const VecDouble &Q1,const VecDouble &Q2)
{
  assert( (unsigned) face.index < m_Vn.size());
  assert( (unsigned) face.cell1 < m_cK.size());
  assert( (unsigned) face.cell2 < m_cK.size());

  
     
  double v = m_Vn(face.index);
  double fMin,fMax;
  double ul = Q1(0);
  double ur = Q2(0);
  

  if (v > 0)
  {
    if (ul > ur)
    {
      m_fMobW.getMinMaxValues(ur,ul,fMin,fMax);
      vFlux(0)=fMax*v;
    }
    else 
    {
      m_fMobW.getMinMaxValues(ul,ur,fMin,fMax);
      vFlux(0)=fMin*v;
    }

  }
  else
  {
    if (ul > ur)
    {
      m_fMobW.getMinMaxValues(ur,ul,fMin,fMax);
      vFlux(0)=fMin*v;
    }
    else 
    {
      m_fMobW.getMinMaxValues(ul,ur,fMin,fMax);
      vFlux(0)=fMax*v;
    }
  }

  
  //If the face has normal in the Z direction, then add the gravity term.
  if (face.normal == Z && !face.at_boundary )
  {
    const VecDouble &cK = getPermeability();
    double K1 = cK(face.cell1);
    double K2 = cK(face.cell2);
    vFlux(0) += (m_fMobGrav(ul) + m_fMobGrav(ur)) * K1*K2/(K1+K2) * m_gc; 
  }
}



 void FaceFluxWithGravity::fluxAtFace(VecDouble &vFlux, const FaceInfo &face,const VecDouble &Q1,const VecDouble &Q2) 
{
  assert( (unsigned) face.index < m_Vn.size());
  assert( (unsigned) face.cell1 < m_cK.size());
  assert( (unsigned) face.cell2 < m_cK.size());

  const double& u1 = Q1(0);
  const double& u2 = Q2(0);
     
  double dMobW1 = m_fMobW(u1);
  double dMobW2 = m_fMobW(u2);

  double v = m_Vn(face.index);
  vFlux(0) = (dMobW1+dMobW2)/2.0*v;
  
  //If the face has normal in the Z direction, then add the gravity term.
  if (face.normal == Z && !(face.at_boundary) )
  {
    const VecDouble &cK = getPermeability();
    double K1 = cK(face.cell1);
    double K2 = cK(face.cell2);
    vFlux(0) += (m_fMobGrav(u1) + m_fMobGrav(u2)) * K1*K2/(K1+K2) * m_gc;
    //flux += (K1*m_fMobGrav(u1) + K2*m_fMobGrav(u2))/2.0 * m_gc; 
  }
}









void  FaceFluxWithGravity::maxGlobalCharVelocity(double vel[3]) 
{
  static VecDouble normalVel(3);
  m_mesh.getMaxAbsVelocitiesInFacesByComponents(normalVel,m_Vn);
  vel[X]= normalVel(X)*m_DfMobW.getMaxNorm(0,1);
  vel[Y]= normalVel(Y)*m_DfMobW.getMaxNorm(0,1);
  vel[Z]= normalVel(Z)*m_DfMobW.getMaxNorm(0,1) + m_maxK*fabs(m_DfMobGrav.getMaxNorm(0,1)*m_gc);
  printf("Max Vel %g, GradMob  %g\n",std::max(normalVel(X), std::max(normalVel(Y),normalVel(Z))),m_DfMobW.getMaxNorm(0,1));
}



double FaceFluxWithGravity::maxLocalCharVelocity(const FaceInfo &face,const VecDouble &vQ1,const VecDouble &vQ2)
{
  double u1=vQ1(0);
  double u2=vQ2(0);
  if (u1>u2)
    std::swap(u1,u2);

  double v=m_Vn(face.index);


  
  if (face.normal != Z)
    return  v*m_DfMobW.getMaxNorm(u1,u2);
  else
  {
    double term1 = v*m_DfMobW.getMaxNorm(u1,u2);

    //Get the max permeability
    double K = fmax(getPermeability()(face.cell1), getPermeability()(face.cell2));

    //The m_fMobGrav function can have negative derivative.
    //An the gravity term of the flux can be negative or positive.
    double dMax,dMin;
    m_DfMobGrav.getMinMaxValues(u1,u2,dMin,dMax); 
    double dMaxAbs = std::max(fabs(dMin),fabs(dMax));
    
    return  fabs(term1) + fabs(K*dMaxAbs*m_gc);
  }
}





 void FaceFluxWithGravity::updateDynamicData(DynamicBase &dynMod) 
{
  dynMod.getNormalVelocityAtFaces(m_Vn);
}

---