RichardsFastSlow/src/kernel/solver.cpp

#include "solver.hpp"

namespace Kernel{

  Solver::Solver(string filename){
    fstream file;
    file.open(filename.c_str(),fstream::in|fstream::binary);
    input_data.load(file,true);
    file.close();

    //HACK for pump
    for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
      Pump& P=*(*it);
      double lX=input_data.geometry.lX;
      P.left_init*=lX;
      P.right_init*=lX;
      P.delta_left_init*=lX;
      P.delta_right_init*=lX;
      P.left_final*=lX;
      P.right_final*=lX;
      P.delta_left_final*=lX;
      P.delta_right_final*=lX;

    }

    pumps=new double*[input_data.geometry.nX];
    for(size_t ix=0;ix<input_data.geometry.nX;++ix){
      pumps[ix]=new double[input_data.geometry.nZ[ix]];
    }

    //Init Piccard object
    piccard.init(&input_data.geometry,pumps);

    //Create pool of solutions
    for(size_t i=0;i<Time::nT;++i){
      Solution* s=new Solution(input_data.geometry);
      solutions_stack.push(s);
    }

    //Init solutions
    solutions=new Solution*[Time::nT];
    for(size_t i=1;i<Time::nT;++i){
      solutions[i]=nullptr;
    }

    //Init first solution
    solutions[0]=new Solution(input_data.geometry);
    init_first_solution();

    m1=1;
  }

  void
  Solver::init_first_solution(){
    Solution& s=*solutions[0];

    //Init pressure
    for(size_t x=0;x<input_data.geometry.nX;++x){
      for(size_t z=0;z<input_data.geometry.nZ[x];++z){
        s.P[x][z]=input_data.initial_state->Pinit[x][z];
      }
    }

    //Init hydraulic head
    for(size_t x=0;x<input_data.geometry.nX;++x){
      s.hydr[x]=input_data.initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
    }

    //Init overland level
    for(size_t x=0;x<input_data.geometry.nX;++x){
      s.hov[x]=Physics::invPsol(input_data.geometry.hsoil[x],s.P[x][input_data.geometry.nZ[x]-1]); //TODO[Chrisophe] à initiliser
    }

    //Init l
    for(size_t x=0;x<input_data.geometry.nX;++x){
      double t=input_data.initial_state->hsat[x];
      s.l[x]=t;
      s.hsat[x]=t;
    }

    //Init Pl
    for(size_t x=0;x<input_data.geometry.nX;++x){
      s.Pl[x]=Psat;
    }

    //Initial state of input_data will be no longer used.
    input_data.remove_initial_state();

    number_computed_solutions=1;
  }

  bool
  Solver::compute_next_solution(){
    //cout<<"[solver::next]"<<endl;
    compute_sources(number_computed_solutions*Time::dT);
    if(m1>1) m1=m1/2;
    Solution* s=space_solution();
    if(s==nullptr) return false;
    solutions[number_computed_solutions++]=s;
    //cout<<"[next] done"<<endl;
    return true;
  }

  Solution*
  Solver::space_solution(){
    #ifdef MYDEBUG
    cout<<" [Solver::spaceSolution] start"<<endl;
    #endif

    size_t n1=0; //we work on time interval [n-1+n1/m1,n-1+(n1+1)/m1]

    piccard.previous_solution=solutions[number_computed_solutions-1];
    piccard.new_solution=get_new_solution();
    piccard.dt=Time::dT/m1*3600;


    while(m1<=max_time_subdivisions){
      #ifdef MYDEBUG
      if(Debug::level>1) cout<<" [Solver::spaceSolution] n1 = "<<n1<<" and m1 = "<<m1<<endl;
      #endif
      piccard.run();
      if(!piccard.has_converged){
        m1*=2;
        n1=0;
        piccard.dt/=2;
        piccard.previous_solution=solutions[number_computed_solutions-1];
      }
      else{
        if(n1!=0){//S_in is different form solution[n-1]
          release_solution(piccard.previous_solution);
        }
        ++n1;
        if(n1==m1) break;
        piccard.previous_solution=piccard.new_solution;
        piccard.new_solution=get_new_solution();
      }
    }
    if(m1>max_time_subdivisions){
      #ifdef MYDEBUG
      cout<<" [Solver::spaceSolution] not converge"<<endl;
      #endif
      release_solution(piccard.new_solution);
      return nullptr;
    }
    #ifdef MYDEBUG
    cout<<" [Solver::spaceSolution] converge"<<endl;
    #endif
    return piccard.new_solution;
  }

  void Solver::compute_sources(double t){
    for(size_t ix=0;ix<input_data.geometry.nX;++ix){
      double x=ix*input_data.geometry.dX;
      for(size_t iz=0;iz<input_data.geometry.nZ[ix];++iz){
        double z=input_data.geometry.hbot[ix]+iz*input_data.geometry.dZ[ix];
        pumps[ix][iz]=0;
        for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
          double p=(*it)->value(x,z,t/Time::T);
          pumps[ix][iz]+=p;
        }
      }
      }
  }
}