RichardsFastSlow/src/kernel/physics.cpp

#include "physics.hpp"
//------------------------
// Brooks and Corey model
//------------------------

double Physics::g;
double Physics::rho;
double Physics::phi;
double Physics::k0;
double Physics::nivrivsat;
double Physics::Psec;
Physics::Model Physics::model;
double (*Physics::s)(double);
double (*Physics::ds)(double);
void (*Physics::s_ds)(double,double&,double&);
double (*Physics::s_inv)(double);
double (*Physics::kr)(double);
double (*Physics::dkr)(double);
void (*Physics::kr_dkr)(double,double&,double&);
double Physics::model_data[6];


#define sres Physics::model_data[1]
#define lambda Physics::model_data[2]
#define alpha Physics::model_data[3]

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

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::s_inv_BC(double S){
  return pow((S-sres)/(1-sres),-1/lambda)*Psat;
}

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;
  }
}

double
Physics::tilde_a(double l,double hbot){
  double t=(l-hbot);
  return t*t/(3*k0);
}

double
Physics::tilde_c(double l,double hbot){
  return (l-hbot)/(2*k0);
}

void
Physics::save(fstream& file){
  file.write((char*)&g,sizeof(double));
  file.write((char*)&rho,sizeof(double));
  file.write((char*)&phi,sizeof(double));
  file.write((char*)&k0,sizeof(double));
  file.write((char*)&Psec,sizeof(double));
  file.write((char*)&nivrivsat,sizeof(double));
  file.write((char*)&model,sizeof(Model));
  for(size_t i=0;i<max_model_parameters;++i){
    file.write((char*)&model_data[i],sizeof(double));
  }
}

void
Physics::load(fstream& file){
  file.read((char*)&g,sizeof(double));
  file.read((char*)&rho,sizeof(double));
  file.read((char*)&phi,sizeof(double));
  file.read((char*)&k0,sizeof(double));
  file.read((char*)&Psec,sizeof(double));
  file.read((char*)&nivrivsat,sizeof(double));
  Model m;
  file.read((char*)&m,sizeof(Model));
  setModel(m);
  for(size_t i=0;i<max_model_parameters;++i){
    file.read((char*)&model_data[i],sizeof(double));
  }
}

double
Physics::Psol(double hsoil,double hov){
  double a=Psat-rho*g*nivrivsat;
  return rho*g*(hov-hsoil)+a*nivrivsat/(hov-hsoil);
}

double
Physics::dPsol(double hsoil,double hov){
  double a=Psat-rho*g*nivrivsat;
  double b=hov-hsoil;
  return rho*g-a*nivrivsat/(b*b);
}

double
Physics::invPsol(double hsoil,double Psol){
  double a=Psat-rho*g*nivrivsat;
  double b=a*nivrivsat;
  return hsoil+(Psol+sqrt(Psol*Psol-4*rho*g*b))/(2*rho*g);
}