RichardsFastSlow/src/initial_state.cpp

#include "initial_state.hpp"

Tank::Tank(){
  saturation=0.8;
  left=0.4;right=0.6;
  bottom=0.6;top=0.7;
  delta_left=0.05;
  delta_right=0.05;
  delta_bottom=0.05;
  delta_top=0.05;
}

void
Tank::save(fstream& file){
  file.write((char*)&saturation,sizeof(double));
  file.write((char*)&left,sizeof(double));
  file.write((char*)&right,sizeof(double));
  file.write((char*)&bottom,sizeof(double));
  file.write((char*)&top,sizeof(double));
  file.write((char*)&delta_left,sizeof(double));
  file.write((char*)&delta_right,sizeof(double));
  file.write((char*)&delta_bottom,sizeof(double));
  file.write((char*)&delta_top,sizeof(double));
}

void
Tank::load(fstream& file){
  file.read((char*)&saturation,sizeof(double));
  file.read((char*)&left,sizeof(double));
  file.read((char*)&right,sizeof(double));
  file.read((char*)&bottom,sizeof(double));
  file.read((char*)&top,sizeof(double));
  file.read((char*)&delta_left,sizeof(double));
  file.read((char*)&delta_right,sizeof(double));
  file.read((char*)&delta_bottom,sizeof(double));
  file.read((char*)&delta_top,sizeof(double));
}

InitialState::InitialState(){
  geometry=nullptr;
  hsat=nullptr;
  Pinit=nullptr;
}

InitialState::InitialState(Geometry* _geometry){
  geometry=_geometry;
  hsat=nullptr;
  Pinit=nullptr;
}

Tank*
InitialState::addTank(){
  Tank* tank=new Tank;
  tanks.push_back(tank);
  //updatePressure();
  return tank;
}

void
InitialState::removeTank(Tank* tank){
  for(auto it=tanks.begin();it!=tanks.end();++it){
    if(*it==tank){
      delete *it;
      tanks.erase(it);
      return;
    }
  }
  //updatePressure();
}

InitialState::~InitialState(){
  if(hsat!=nullptr){
    delete[] hsat;
    for(size_t i=0;i<geometry->nX;++i){
      delete[] Pinit[i];
    }
    delete[] Pinit;
    for(auto it=tanks.begin();it!=tanks.end();++it){
      delete(*it);
    }
  }
}

void
InitialState::updatePressure(){
  for(size_t i=0;i<geometry->nX;++i){
    double x=double(i)/(geometry->nX-1);
    size_t nZ=geometry->nZ[i];
    double temp=hsat[i]+Psat/(Physics::g*Physics::rho); //TODO : replace Physics::model_data[0] by Psat in future
    //cout<<i<<" -> "<<nZ<<endl;
    for(size_t j=0;j<nZ;++j){
      double z=geometry->Z[i][j];
      Pinit[i][j]=Physics::rho*Physics::g*(temp-z);
      for(auto it=tanks.begin();it!=tanks.end();++it){
        Tank* tank=*it;
        double S=tank->saturation;
        double l=tank->left;
        double r=tank->right;
        double t=tank->top;
        double b=tank->bottom;
        double dl=tank->delta_left;
        double dr=tank->delta_right;
        double dt=tank->delta_top;
        double db=tank->delta_bottom;

        double Ps=Physics::s_inv(S);
        double Px,Pz;
        if(x<=l){
          Px=(Ps-Physics::Psec)*(x-l)/dl+Ps;
        }
        else if(x>=r){
          Px=(Physics::Psec-Ps)*(x-r)/dr+Ps;
        }
        else{
          Px=Ps;
        }
        if(z<=b){
          Pz=(Ps-Physics::Psec)*(z-b)/db+Ps;
        }
        else if(z>=t){
          Pz=(Physics::Psec-Ps)*(z-t)/dt+Ps;
        }
        else{
          Pz=Ps;
        }
        Pinit[i][j]=max(max(min(Px,Pz),Pinit[i][j]),Physics::Psec);
      }
      //TODO : voir considrer max avec psec
    }
  }
}

void
InitialState::save(fstream& file){
  file.write((char*)hsat,geometry->nX*sizeof(double));
  file.write((char*)hov,geometry->nX*sizeof(double));
  for(size_t i=0;i<geometry->nX;++i){
    file.write((char*)Pinit[i],geometry->nZ[i]*sizeof(double));
  }
  size_t nt=tanks.size();
  file.write((char*)&nt,sizeof(double));
  for(auto it=tanks.begin();it!=tanks.end();++it){
    (*it)->save(file);
  }
}

void
InitialState::load(fstream& file,bool init){
  if(init){
    hsat=new double[geometry->nX];
    hov=new double[geometry->nX];
    Pinit=new double*[geometry->nX];
  }
  file.read((char*)hsat,geometry->nX*sizeof(double));
  file.read((char*)hov,geometry->nX*sizeof(double));
  for(size_t i=0;i<geometry->nX;++i){
    if(init) Pinit[i]=new double[geometry->nZ[i]];
    file.read((char*)Pinit[i],geometry->nZ[i]*sizeof(double));
  }
  size_t nt;
  file.read((char*)&nt,sizeof(size_t));
  for(size_t i=0;i<nt;++i){
    Tank* tank=addTank();
    tank->load(file);
  }
}