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();

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

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