123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158 |
- #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);
- }
|