tracefacefluxwithgravity.cpp

Go to the documentation of this file.

#include "tracefacefluxwithgravity.h"


 TraceFaceFluxWithGravity::TraceFaceFluxWithGravity(unsigned numFaces,Function1D &fMobW,Function1D &fMobGrav,VecDouble &cK,double gc) 
  :m_fMobW(fMobW),m_fMobGrav(fMobGrav),m_cK(cK),m_Vel(numFaces,3)
{
  m_gc = gc;
  m_maxK = cK.linfty_norm();
}


double TraceFaceFluxWithGravity::fluxAtFace(OrthoMesh::Face_It &faceIt,int face,int cell1,int cell2,double &u1,double &u2) 
{
  assert( (unsigned) face < m_Vel.m());
  assert( (unsigned) cell1 < m_cK.size());
  assert( (unsigned) cell2 < m_cK.size());
  assert(m_cK.size() == getSw().size());
  VecDouble &cVSw = getSw();
  double S1 = cVSw(cell1);
  double S2 = cVSw(cell2);
  
  double dMobW1 = 1-m_fMobW(S1);
  double dMobW2 = 1-m_fMobW(S2);
  double *vel = &m_Vel(face,0);
  double flux = (dMobW1+dMobW2)/2.0*faceIt->normal_multiply(vel);
  
  //If the face has normal in the Z direction, then add the gravity term.
  if (faceIt->getNormalNonZeroComponent() == Z && !faceIt->at_boundary() )
  {
    VecDouble &cK = getPermeability();
    double K1 = cK(cell1);
    double K2 = cK(cell2);
    flux += (m_fMobGrav(S1) + m_fMobGrav(S2)) * K1*K2/(K1+K2) * m_gc; 
  }
  return (u1+u2)/(S1+S2)*flux/(S;
}



void  TraceFaceFluxWithGravity::maxGlobalCharVelocity(double vel[3]) 
{
  VecDouble minVel(3);
  VecDouble maxVel(3);

  NumericMethods::getMinMaxValueByColumn(m_Vel,minVel,maxVel);
  vel[X]= NumericMethods::maxMod(minVel(X),maxVel(X))*m_fMobW.getMaxNormGrad(0,1);
  vel[Y]= NumericMethods::maxMod(minVel(Y),maxVel(Y))*m_fMobW.getMaxNormGrad(0,1);
  vel[Z]= NumericMethods::maxMod(minVel(Z),maxVel(Z))*m_fMobW.getMaxNormGrad(0,1) + m_maxK*fabs(m_fMobGrav.getMaxNormGrad(0,1)*m_gc);
}



double TraceFaceFluxWithGravity::maxLocalCharVelocity(OrthoMesh::Face_It &faceIt,int face,int cell1,int cell2,double &u1,double &u2)
{
  assert(u1<=u2);
  double *vel = &(m_Vel(face,0));
  
  if (faceIt->getNormalNonZeroComponent() != Z)
    return  faceIt->normal_multiply(vel)*m_fMobW.getMaxNormGrad(u1,u2);
  else
  {
    double term1 = faceIt->normal_multiply(vel)*m_fMobW.getMaxNormGrad(u1,u2);

    //Get the max permeability
    double K = fmax(getPermeability()(cell1), getPermeability()(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_fMobGrav.getMinMaxGradValues(u1,u2,dMin,dMax); 

    double result1 = term1 + K*dMin*m_gc;
    double result2 = term1 + K*dMax*m_gc;

    if (fabs(result1) < fabs(result2))
      return result2;
    else
      return result1;
  }
}


void TraceFaceFluxWithGravity::setVelocity(DynamicBase &dynMod) 
{
  dynMod.getVelocitiesAtFaces(m_Vel);
}

---