RichardsFastSlow/src/qt/input/data.cpp

#include "qt/input/data.hpp"

QtInputData::QtInputData(){
  depth=1;
  nZ_max=100;
  initPhysics();
  initSplines();
  initGeometry();
  initInitialState();
}

void QtInputData::initPhysics(){
  Physics::g=9.81;
  Physics::rho=1000;
  Physics::phi=0.3;
  Physics::k0=3e-5;
  Physics::nivrivsat=0.02;
  Physics::Psec=-10000;
  Physics::model_data[0]=-2000;
  Physics::model_data[1]=0;
  Physics::model_data[2]=3;
  Physics::model_data[3]=11;
  Physics::model_data[4]=0;
  Physics::model_data[5]=0;
  Physics::setModel(Physics::BrooksCorey);
}

void
QtInputData::initSplines() {
  // Control points of splines
  for(size_t i=0;i<np;++i){
    double x=double(i)/(np-1);
    point[i].x=x;
    point[i].y=0.2;
    point[np+i].x=x;
    point[np+i].y=0.8;
    point[2*np+i].x=x;
    point[2*np+i].y=0.5;
  }
  // Control points of oveland points
  for(size_t i=0;i<np_ov;++i){
    double x=double(i+0.5)/np_ov;
    point[3*np+i].x=x;
    point[3*np+i].y=0.75;
  }

  // Splines
  spline_bot.setPoints(point,np);
  spline_bot.compute();
  spline_soil.setPoints(&point[np],np);
  spline_soil.compute();
  spline_sat.setPoints(&point[2*np],np);
  spline_sat.compute();
}

void
QtInputData::initGeometry(){
  geometry.lX=10;
  geometry.nX=100;
  geometry.dX=geometry.lX/(geometry.nX-1);
  geometry.hsoil=new double[nmax_Qt];
  geometry.dhsoil=new double[nmax_Qt];
  geometry.hbot=new double[nmax_Qt];
  geometry.dhbot=new double[nmax_Qt];
  geometry.nZ=new size_t[nmax_Qt];
  geometry.dZ=new double[nmax_Qt];
  geometry.Z=new double*[nmax_Qt];
  for(size_t i=0;i<nmax_Qt;++i){
    geometry.Z[i]=new double[nmax_Qt];
  }
  updateGeometry();
}

void
QtInputData::initInitialState(){
  initial_state->hsat=new double[nmax_Qt];
  initial_state->hov=new double[nmax_Qt];
  initial_state->Pinit=new double*[nmax_Qt];
  for(size_t i=0;i<nmax_Qt;++i){
    initial_state->Pinit[i]=new double[nmax_Qt];
  }
  updateInitialState();
}

void
QtInputData::updateGeometry(){
  double prev_factor=factor;
  double hs_max=-inf,hb_min=inf;
  double v;
  geometry.dX=geometry.lX/(geometry.nX-1);
  for(size_t k=0;k<geometry.nX;++k){
    double x=k*geometry.dX/geometry.lX;
    v=geometry.hsoil[k]=spline_soil(x);
    hs_max=max(v,hs_max);
    geometry.dhsoil[k]=spline_soil.derivate(x);
    v=geometry.hbot[k]=spline_bot(x);
    hb_min=min(v,hb_min);
    geometry.dhbot[k]=spline_bot.derivate(x);
  }
  double dZ_avg=depth/nZ_max;
  factor=depth/(hs_max-hb_min);
  for(size_t k=0;k<geometry.nX;++k){
    geometry.hsoil[k]*=factor;
    geometry.hbot[k]*=factor;
  }
  geometry.initZ(dZ_avg,false);
  double r=factor/prev_factor;
  for(auto it=initial_state->tanks.begin();it!=initial_state->tanks.end();++it){
    Tank* tank=*it;
    tank->top*=r;
    tank->bottom*=r;
    tank->delta_top*=r;
    tank->delta_bottom*=r;
  }
  for(auto it=source.pumps.begin();it!=source.pumps.end();++it){
    Pump* pump=*it;
    pump->top_init*=r;
    pump->bottom_init*=r;
    pump->delta_top_init*=r;
    pump->delta_bottom_init*=r;
    pump->top_final*=r;
    pump->bottom_final*=r;
    pump->delta_top_final*=r;
    pump->delta_bottom_final*=r;
  }
}

void
QtInputData::updateInitialState(){
  for(size_t i=0;i<geometry.nX;++i){
    double x=double(i)/(geometry.nX-1);
    initial_state->hsat[i]=spline_sat(x)*factor;
  }
  initial_state->updatePressure();
  updateOverland();
}

void
QtInputData::save(fstream& file){
  InputData::save(file);
  file.write((char*)&depth,sizeof(double));
  file.write((char*)&nZ_max,sizeof(size_t));
  size_t temp=np;
  file.write((char*)&temp,sizeof(size_t));
  for(size_t i=0;i<3*np+np_ov;++i){
    file.write((char*)&point[i].x,sizeof(double));
    file.write((char*)&point[i].y,sizeof(double));
  }
}

void
QtInputData::load(fstream& file){
  InputData::load(file);
  file.read((char*)&depth,sizeof(double));
  file.read((char*)&nZ_max,sizeof(size_t));
  size_t temp;
  file.read((char*)&temp,sizeof(size_t));
  if(temp!=np){
    cerr<<"Incorrect number of points"<<endl;
    return;
  }
  for(size_t i=0;i<3*np+np_ov;++i){
    file.read((char*)&point[i].x,sizeof(double));
    file.read((char*)&point[i].y,sizeof(double));
  }
  updateSplineBot();
  updateSplineSat();
  updateSplineSoil();
  updateGeometry();
  updateInitialState();
  updateOverland();
}

QtInputData::~QtInputData(){
  delete[] geometry.hsoil;
  delete[] geometry.dhsoil;
  delete[] geometry.hbot;
  delete[] geometry.dhbot;
  delete[] geometry.nZ;
  for(size_t i=0;i<nmax_Qt;++i){
    delete[] geometry.Z[i];
  }
  delete[] geometry.Z;
  geometry.hsoil=nullptr; //Bypass geometry destructor
  delete[] initial_state->hsat;
  for(size_t i=0;i<nmax_Qt;++i){
    delete[] initial_state->Pinit[i];
  }
  delete[] initial_state->Pinit;
  initial_state->hsat=nullptr; //Bypassfirst initial_state destructor
}

void
QtInputData::updateOverland(){
  for(size_t x=0;x<geometry.nX;++x){
    initial_state->hov[x]=geometry.hsoil[x];
  }
  for(size_t i=0;i<np_ov;++i){
    Point& P=point[3*np+i];
    size_t x_ov=P.x*(geometry.nX-1);
    double y=P.y*factor;
    if(geometry.hsoil[x_ov]<y){
      size_t x_left=findSoilOver(y,x_ov,-1);
      size_t x_right=findSoilOver(y,x_ov,1);
        for(size_t x=x_left;x<=x_right;++x){
        initial_state->hov[x]=max(initial_state->hov[x],y);
      }
    }
  }
}

size_t
QtInputData::findSoilOver(double y,int x_ov,int x_step){
  int x=x_ov;
  while(x>=0 and x<geometry.nX and geometry.hsoil[x]<y){
    x+=x_step;
  }
  x=min(x,(int)geometry.nX-1);
  x=max(x,0);
  return x;
}