RichardsFastSlow/src/physics.cpp

#include "physics.hpp"

double Physics::model_datas[6];

//-------------
// Constructor
//-------------

Physics::Physics(PhysicModel model){
  switch(model){
  case BrooksCorey:
    s=&s_BC;
    ds=&ds_BC;
    s_ds=&s_ds_BC;
    kr=&kr_BC;
    dkr=&dkr_BC;
    kr_dkr=&kr_dkr_BC;
    break;
  default:
    assert(false);
  };

}
//------------------------
// Brooks and Corey model
//------------------------

#define Psat model_datas[0]
#define sres model_datas[1]
#define lambda model_datas[2]
#define alpha model_datas[3]

double
Physics::s_BC(double P){
  if(P>=Psat) return 1;
  return sres+(1-sres)*pow(Psat/P,lambda);
}

double
Physics::ds_BC(double P){
  if(P>=Psat) return 0;
  return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
}

void
Physics::s_ds_BC(double P,double& v,double& dv){
  if(P>=Psat){
    v=1;
    dv=0;
  }
  else{
    double t=(1-sres)*pow(Psat/P,lambda);
    v=sres+t;
    dv=-(lambda*t)/P;
  }
}

double
Physics::kr_BC(double P){
  if(P>=Psat) return 1;
  return pow(Psat/P,alpha);
}

double
Physics::dkr_BC(double P){
  if(P>=Psat) return 0;
  return -alpha*pow(Psat/P,alpha)/P;
}

void
Physics::kr_dkr_BC(double P,double& v,double& dv){
  if(P>=Psat){
    v=1;
    dv=0;
  }
  else{
    double t=pow(Psat/P,alpha);
    v=t;
    dv=-(alpha*t)/P;
  }
}