RichardsFastSlow/src/kernel/solver.hpp

#ifndef KERNEL_HPP
#define KERNEL_HPP
#include <cstdio>
#include <cstring>
#include <fstream>
#include <stack>
#include "physics.hpp"
#include "geometry.hpp"
#include "time.hpp"
#include "input_data.hpp"
#include "solution.hpp"
#include "piccard.hpp"

using namespace std;

namespace Kernel{

  //! The Solver class.
  /*!
  This is the entry class to compute solutions of the problem.
  From an input file, it computes a Solution to each time step
  of the simulation.
  */
  class Solver{
  private:

    double** pumps;

    //! Piccard object used for Piccard algorithm.
    Piccard piccard;

    //! Maximal number of time subintervals allowed to compute a new Solution form a previous one.
    static constexpr size_t max_time_subdivisions=32;

    //! Datas stored in the input file will be loaded in input_data
    InputData input_data;

    //! A stack of pointers to valid Solution object.
    stack<Solution*> solutions_stack;

    //! Array of pointer to Solutions computed for each time step
    Solution** solutions;

    //! Numbers of Solutions computed so far
    size_t number_computed_solutions;

    //! Number of time subintervals used to compute a new Solution from the previous one
    size_t m1;

    //! Get a newSolution object from the stack #solutions_stacks
    Solution* get_new_solution();

    //! Returns a solution that is no longer used on the stack
    void release_solution(Solution*);

    //! Init the solution at time t=0 from the input file
    void init_first_solution();

    //! Try to compute the solution at time t+1 from the solution at time t
    Solution* space_solution();

    void compute_sources(double t);

  public:
    //! Construct a Solver from the input file input_filename
    Solver(string input_filename);

    //! Return the geometry of the simulation, which is given in input file
    const Geometry& get_geometry() const;

    const Source& get_source() const;

    //! Return the vertical factor (used for drawing)
    double get_vertical_factor() const;

    //! Return tjhe number of solution computed until now, relatively to the time
    size_t get_computed_solutions_number() const;

    //! Return the solution computed at time t. Return nullptr if the solution is not yet computed.
    const Solution* get_solution(size_t t) const;

    //! Compute the next Solution. Return false if the schema can't obtain such a soluition.
    bool compute_next_solution();
  };

  inline Solution*
  Solver::get_new_solution(){
    if(solutions_stack.empty()){
      cerr<<"[Error] Solution stack is empty"<<endl;
    }
    Solution* s=solutions_stack.top();
    solutions_stack.pop();
    return s;
  }

  inline void
  Solver::release_solution(Solution* s){
    solutions_stack.push(s);
  }

  inline const Geometry&
  Solver::get_geometry() const{
    return input_data.geometry;
  }

  inline const Source&
  Solver::get_source() const{
    return input_data.source;
  }

  inline double
  Solver::get_vertical_factor() const{
    return input_data.factor;
  }

  inline size_t
  Solver::get_computed_solutions_number() const{
    return number_computed_solutions;
  }

  inline const Solution*
  Solver::get_solution(size_t n) const{
    return solutions[n];
  }
}
#endif