Hydr (with missing files)

backup/obsolete/geometry_spline.cpp

@@ -0,0 +1,72 @@
+#include "geometry_spline.hpp"
+
+void
+GeometrySpline::initDefault(){
+  lX=10;
+  nX=100;
+  nX_max=nmax_Qt;
+  dX=lX/(nX-1);
+  hsoil=new double[nmax_Qt];
+  dhsoil=new double[nmax_Qt];
+  hbot=new double[nmax_Qt];
+  dhbot=new double[nmax_Qt];
+  nZ=new size_t[nmax_Qt];
+  dZ=new double[nmax_Qt];
+  Z=new double*[nmax_Qt];
+  for(size_t i=0;i<nmax_Qt;++i){
+    Z[i]=new double[nmax_Qt];
+  }
+  for(size_t i=0;i<nX;++i){
+    hsoil[i]=0.8;
+    hbot[i]=0.2;
+    dhsoil[i]=0;
+    dhbot[i]=0;
+  }
+  depth=5;
+  nZ_max=100;
+  factor=5/0.6;
+  initZ(depth/nZ_max,false);
+}
+
+void
+GeometrySpline::update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& Ssoil,Spline& Sbot){
+  lX=_lX;
+  depth=_depth;
+  nX=_nX;
+  nZ_max=_nZ_max;
+  dX=lX/(nX-1);
+  double hs_max=-inf,hb_min=inf;
+  double v;
+  for(size_t k=0;k<nX;++k){
+    double x=k*dX/lX;
+    v=hsoil[k]=Ssoil(x);
+    hs_max=max(v,hs_max);
+    dhsoil[k]=Ssoil.derivate(x);
+    v=hbot[k]=Sbot(x);
+    hb_min=min(v,hb_min);
+    dhbot[k]=Sbot.derivate(x);
+  }
+  double dZ_avg=depth/_nZ_max;
+  factor=depth/(hs_max-hb_min);
+  for(size_t k=0;k<nX;++k){
+    hsoil[k]*=factor;
+    hbot[k]*=factor;
+  }
+  initZ(dZ_avg,false);
+}
+
+GeometrySpline::~GeometrySpline(){
+  if(hsoil!=nullptr){
+    delete[] hsoil;
+    delete[] dhsoil;
+    delete[] hbot;
+    delete[] dhbot;
+    delete[] dZ;
+    delete[] nZ;
+    for(size_t i=0;i<nX_max;++i){
+      delete[] Z[i];
+    }
+    delete[] Z;
+  }
+  hsoil=nullptr;
+}

backup/obsolete/geometry_spline.hpp

@@ -0,0 +1,23 @@
+#ifndef GEOMETRY_SPLINE_HPP
+#define GEOMETRY_SPLINE_HPP
+
+#include "geometry.hpp"
+
+class GeometrySpline:public Geometry{
+public:
+  size_t nX_max;
+  double factor;
+  static constexpr size_t const nmax_Qt=400;
+  GeometrySpline();
+  ~GeometrySpline();
+  void initDefault(Geometry* ge);
+  void update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& hbot,Spline& hsoil);
+};
+
+inline
+GeometrySpline::GeometrySpline():Geometry(){
+
+}
+
+
+#endif

backup/obsolete/initial_state_spline.cpp

@@ -0,0 +1,26 @@
+#include "initial_state_spline.hpp"
+
+void
+InitialStateSpline::update(Spline* _Ssat){
+  Ssat=_Ssat;
+  if(hsat!=nullptr){
+    delete[] hsat;
+    for(size_t i=0;i<geometry_spline->nX;++i){
+      delete[] Pinit[i];
+    }
+    delete[] Pinit;
+  }
+  size_t nX=geometry_spline->nX;
+  hsat=new double[nX];
+  Pinit=new double*[nX];
+
+  for(size_t i=0;i<nX;++i){
+    size_t nZ=geometry_spline->nZ[i];
+    Pinit[i]=new double[nZ];
+  }
+  for(size_t i=0;i<geometry_spline->nX;++i){
+    double x=double(i)/(geometry_spline->nX-1);
+    hsat[i]=(*Ssat)(x)*geometry_spline->factor;
+  }
+  updatePressure();
+}

backup/obsolete/initial_state_spline.hpp

@@ -0,0 +1,21 @@
+#ifndef INITIAL_STATE_SPLINE_HPP
+#define INITIAL_STATE_SPLINE_HPP
+
+#include "initial_state.hpp"
+#include "geometry_spline.hpp"
+
+class InitialStateSpline:public InitialState{
+private:
+  Spline* Ssat;
+  GeometrySpline* geometry_spline;
+public:
+  InitialStateSpline(GeometrySpline* geometry);
+  void update(Spline* Ssat);
+};
+
+inline
+InitialStateSpline::InitialStateSpline(GeometrySpline* g):InitialState((Geometry*)g){
+  geometry_spline=g;
+}
+
+#endif

backup/obsolete/input_geometry_curves.cpp

@@ -0,0 +1,163 @@
+#include "input_geometry_curves.hpp"
+
+QtInputGeometryCurves::QtInputGeometryCurves(){
+  margin=10;
+  radius=8;
+  min_d=0.01;
+  initPoints();
+  hbot.setPoints(point,np);
+  hbot.compute();
+  hsoil.setPoints(&point[np],np);
+  hsoil.compute();
+}
+
+void
+QtInputGeometryCurves::initPoints(){
+  for(size_t i=0;i<np;++i){
+    double x=double(i)/(np-1);
+    point[i].x=x;
+    point[i].y=0.7;
+    point[i+np].x=x;
+    point[i+np].y=0.3;
+  }
+}
+
+int
+QtInputGeometryCurves::get_x(float x,int w){
+  return margin+x*w;
+}
+
+int
+QtInputGeometryCurves::get_y(float y,int h){
+  return y*h;
+}
+
+void
+QtInputGeometryCurves::paintEvent(QPaintEvent* event){
+  QPainter painter(this);
+  painter.setRenderHint(QPainter::Antialiasing);
+  int w=width()-2*margin;
+  int h=height();
+  // Hbot
+  painter.setBrush(Qt::SolidPattern);
+  painter.setPen(QColor(102,102,102));
+  drawSpline(hbot,painter);
+  painter.setPen(QColor(154,77,0));
+  drawSpline(hsoil,painter);
+  painter.setPen(Qt::black);
+
+  //Control points
+  for(size_t i=0;i<2*np;++i){
+    Point& P=point[i];
+    painter.drawEllipse(get_x(P.x,w)-radius,get_y(P.y,h)-radius,2*radius,2*radius);
+  }
+}
+
+void
+QtInputGeometryCurves::mousePressEvent(QMouseEvent* event){
+  int mx=event->x();
+  int my=event->y();
+  selected=findPoint(mx,my);
+}
+
+void
+QtInputGeometryCurves::mouseMoveEvent(QMouseEvent* event){
+  if(selected<2*np){
+    int w=width()-2*margin;
+    int h=height();
+    float mx=event->x();
+    float my=event->y();
+    float x=(mx-margin)/w;
+    float y=my/h;
+    moveSelected(x,y);
+    update();
+  }
+}
+
+void
+QtInputGeometryCurves::mouseReleaseEvent(QMouseEvent* event){
+  selected=2*np;
+}
+
+size_t
+QtInputGeometryCurves::findPoint(int x,int y){
+  int w=width()-2*margin;
+  int h=height();
+  for(size_t i=0;i<2*np;++i){
+    Point& P=point[i];
+    int px=get_x(P.x,w);
+    int py=get_y(P.y,h);
+    int dx=x-px;
+    int dy=y-py;
+    int d=dx*dx+dy*dy;
+    if(d<=radius*radius) return i;
+  }
+  return 2*np;
+}
+
+void
+QtInputGeometryCurves::moveSelected(float x,float y){
+  double xp=point[selected].x;
+  double yp=point[selected].y;
+  if(selected!=0 and selected!=np-1){
+    float xmin=point[selected-1].x+min_d;
+    float xmax=point[selected+1].x-min_d;
+    if(xmin<=x and x<=xmax){
+      point[selected].x=x;
+    }
+  }
+  if(0<=y and y<=1){
+    point[selected].y=y;
+  }
+  
+  if(selected<np) hbot.compute();
+  else hsoil.compute();
+  double w=width()-2*margin;
+  bool ok=true;
+  for(int i=1;i<=w;++i){
+    double x=i/w;
+    double ybot=hbot(x);
+    double ysoil=hsoil(x);
+    if(ybot<min_d or ybot>1-min_d or ysoil<min_d or ysoil>1-min_d or ybot<ysoil+min_d){
+      point[selected].x=xp;
+      point[selected].y=yp;
+      if(selected<np) hbot.compute();
+      else hsoil.compute();
+      return;
+    }
+  }
+}
+
+void
+QtInputGeometryCurves::drawSpline(Spline& S,QPainter& painter){
+  double w=width()-2*margin;
+  double h=height();
+  double xp=0;
+  double yp=S(0);
+  for(int i=1;i<=w;i+=2){
+    double x=i/w;
+    double y=S(x);
+    painter.drawLine(margin+xp*w,yp*h,margin+x*w,y*h);
+    xp=x;
+    yp=y;
+  }
+}
+
+void
+QtInputGeometryCurves::save(fstream& file){
+  for(size_t i=0;i<2*np;++i){
+    file.write((char*)&point[i].x,sizeof(double));
+    file.write((char*)&point[i].y,sizeof(double));
+  }
+}
+
+void
+QtInputGeometryCurves::load(fstream& file){
+  for(size_t i=0;i<2*np;++i){
+    file.read((char*)&point[i].x,sizeof(double));
+    file.read((char*)&point[i].y,sizeof(double));
+  }
+  hbot.compute();
+  hsoil.compute();
+  update();
+}

backup/obsolete/input_geometry_curves.hpp

@@ -0,0 +1,61 @@
+#ifndef QT_INPUT_GEOMETRY_CURVES_HPP
+#define QT_INPUT_GEOMETRY_CURVES_HPP
+
+#include <iostream>
+#include <fstream>
+#include <QWidget>
+#include <QMouseEvent>
+#include <QPainter>
+
+#include "math/point.hpp"
+#include "math/spline.hpp"
+
+using namespace std;
+
+static const size_t np = 10;
+
+class QtInputGeometryCurves:public QWidget{
+private:
+  int margin;
+  int radius;
+  float min_d;
+  size_t selected;
+  Point point[2*np];
+  Spline hsoil;
+  Spline hbot;
+  
+  size_t findPoint(int x,int y);
+  int get_x(float x,int w);
+  int get_y(float y,int h);
+  void moveSelected(float x,float y);
+  void drawSpline(Spline& S,QPainter& painter);
+public:
+  QtInputGeometryCurves();
+  ~QtInputGeometryCurves();
+  void save(fstream& file);
+  void load(fstream& file);
+  Spline& getSsoil();
+  Spline& getSbot();
+protected:
+  void initPoints();
+  void paintEvent(QPaintEvent* event);
+  void mousePressEvent(QMouseEvent* event) override;
+  void mouseMoveEvent(QMouseEvent* event) override;
+  void mouseReleaseEvent(QMouseEvent* event) override;
+};
+
+inline
+QtInputGeometryCurves::~QtInputGeometryCurves(){
+}
+
+inline Spline&
+QtInputGeometryCurves::getSsoil(){
+  return hsoil;
+}
+
+inline Spline&
+QtInputGeometryCurves::getSbot(){
+  return hbot;
+}
+
+#endif

backup/obsolete/input_overland.hpp

@@ -0,0 +1 @@
+

backup/obsolete/input_sources.cpp

@@ -0,0 +1,7 @@
+#include "input_sources.hpp"
+
+QtInputSources::QtInputSources(){
+}
+
+QtInputSources::~QtInputSources(){
+}

backup/obsolete/input_sources.hpp

@@ -0,0 +1,18 @@
+#ifndef QT_INPUT_SOURCES
+#define QT_INPUT_SOURCES
+
+#include <iostream>
+#include <fstream>
+
+#include <QWidget>
+
+#include "source.hpp"
+
+class QtInputSources:public QWidget{
+private:
+public:
+  QtInputSources();
+  ~QtInputSources();
+};
+
+#endif

backup/obsolete/spline.hpp

@@ -0,0 +1,11 @@
+#ifndef QSPLINE_HPP
+#define QSPLINE_HPP
+
+#include <GL/glut.h>
+
+class QSpline:public Spline{
+public:
+  void QSpline(float z);
+  void draw_points_GL()
+};
+#endif

backup/obsolete/view.cpp

@@ -0,0 +1,111 @@
+#include "view.hpp"
+
+void
+QtView::initializeGL(){
+  glClearColor(1,1,1,1);
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_LIGHTING);
+  glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
+  glEnable(GL_COLOR_MATERIAL);
+}
+
+void
+QtView::resizeGL(int w,int h){
+  glViewport(0,0,w,h);
+  glMatrixMode(GL_PROJECTION);
+  glLoadIdentity();
+  gluOrtho2D(0,1,0,1);
+  glMatrixMode(GL_MODELVIEW);
+  glLoadIdentity();
+}
+
+void
+QtView::paintGL(){
+  Geometry& G=data->geometry;
+  glBegin(GL_TRIANGLES);
+  for(size_t i=0;i<G.nX-1;++i){
+    size_t max_left=G.nZ[i];
+    size_t max_right=G.nZ[i+1];
+    size_t left=0;
+    size_t right=0;
+      while(left<max_left-1 and right<max_right-1){
+      if(G.Z[i][left+1]<G.Z[i+1][right+1]){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+      else{
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(left==max_left-1){
+      while(right<max_right-1){
+        drawTriangle(i,left,i+1,right,i+1,right+1);
+        ++right;
+      }
+    }
+    if(right==max_right-1){
+      while(left<max_left-1){
+        drawTriangle(i,left,i+1,right,i,left+1);
+        ++left;
+      }
+    }
+  }
+  glEnd();
+  drawOverland();
+}
+
+void
+QtView::setColor(size_t ix,size_t iz){
+  double p=getP(ix,iz);
+  double s=Physics::s(p);
+  double r=(1-s);
+  double g=0.7*(1-s);
+  double b=0.8*s+0.4*(1-s);
+  glColor3f(r,g,b);
+}
+
+void
+QtView::drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3){
+  Geometry& G=data->geometry;
+  double factor_inv=1/data->factor;
+  double dX=1/double(G.nX-1);
+  double x1=ix1*dX;
+  double x2=ix2*dX;
+  double x3=ix3*dX;
+  double y1=G.Z[ix1][iz1]*factor_inv;
+  double y2=G.Z[ix2][iz2]*factor_inv;
+  double y3=G.Z[ix3][iz3]*factor_inv;
+  setColor(ix1,iz1);
+  glVertex3f(x1,y1,0);
+  setColor(ix2,iz2);
+  glVertex3f(x2,y2,0);
+  setColor(ix3,iz3);
+  glVertex3f(x3,y3,0);
+}
+
+void
+QtView::drawOverland(){
+  Geometry& G=data->geometry;
+  double dX=1/double(G.nX-1);
+  double factor_inv=1/data->factor;
+  glColor3f(0.5,0.5,1);
+  glBegin(GL_QUADS);
+  for(size_t i=0;i<G.nX-1;++i){
+    glVertex3f(i*dX,G.hsoil[i]*factor_inv,0);
+    glVertex3f(i*dX,data->initial_state->hov[i]*factor_inv,0);
+    glVertex3f(i*dX+dX,data->initial_state->hov[i+1]*factor_inv,0);
+    glVertex3f(i*dX+dX,G.hsoil[i+1]*factor_inv,0);
+  }
+  glEnd();
+/*  glColor3f(0.2,0.2,1);
+
+  glLineWidth(2);
+  glBegin(GL_LINE_STRIP);
+  for(size_t i=0;i<G.nX;++i){
+    glVertex2f(i*dX,data->initial_state.hov[i]*factor_inv);
+  }
+  glEnd();*/
+
+}

backup/obsolete/view.hpp

@@ -0,0 +1,39 @@
+#ifndef QT_VIEW_HPP
+#define QT_VIEW_HPP
+
+#include <iostream>
+
+#include <QWidget>
+#include <QPainter>
+#include <QOpenGLWidget>
+#include <QMouseEvent>
+
+#include <GL/glut.h>
+
+#include "input_data.hpp"
+
+using namespace std;
+
+class QtView:public QOpenGLWidget{
+protected:
+  InputData* data;
+  void drawOverland();
+public:
+  QtView(InputData* data);
+  void setGeometry(Geometry* geometry);
+  void initializeGL();
+  void paintGL();
+  void resizeGL(int x,int h);
+  void drawTriangle(size_t ix1,size_t iz1,size_t ix2,size_t iz2,size_t ix3,size_t iz3);
+  virtual double getP(size_t ix,size_t iz)=0;
+  void setColor(size_t ix,size_t iz);
+
+};
+
+inline
+QtView::QtView(InputData* d):QOpenGLWidget(){
+  data=d;
+}
+
+
+#endif

backup/old/kernel.cpp

@@ -0,0 +1,489 @@
+#include "kernel.hpp"
+
+Kernel::Kernel(string filename){
+  fstream file;
+  file.open(filename.c_str(),fstream::in|fstream::binary);
+  InputData::load(file,true);
+
+  initSolutions();
+
+  Solution& S0=*solution[0];
+
+  //Init pressure
+  for(size_t x=0;x<geometry.nX;++x){
+    for(size_t z=0;z<geometry.nZ[x];++z){
+      S0.P[x][z]=initial_state->Pinit[x][z];
+    }
+  }
+  //Init hydraulic head
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hydr[x]=initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
+  }
+
+  //Init overland level
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.hov[x]=Physics::invPsol(geometry.hsoil[x],S0.P[x][geometry.nZ[x]-1]); //TODO[Chrisophe] à initiliser
+  }
+
+  //Init l
+  for(size_t x=0;x<geometry.nX;++x){
+    double t=initial_state->hsat[x];
+    S0.l[x]=t;
+    S0.hsat[x]=t;
+  }
+
+  //Init Pl
+  for(size_t x=0;x<geometry.nX;++x){
+    S0.Pl[x]=Psat;
+  }
+
+
+  indice_x_Richards=new bool[geometry.nX];
+
+  div_w=new double*[geometry.nX];
+  P_temp=new double*[geometry.nX];
+  for(size_t i=0;i<geometry.nX;++i){
+    div_w[i]=new double[geometry.nZ[i]];
+    P_temp[i]=new double[geometry.nZ[i]];
+  }
+  t_sub_A=new double*[geometry.nX];
+  t_sub_B=new double*[geometry.nX];
+  t_sub_C=new double*[geometry.nX];
+  t_diag_A=new double*[geometry.nX];
+  t_diag_B=new double*[geometry.nX];
+  t_diag_C=new double*[geometry.nX];
+  t_sup_A=new double*[geometry.nX];
+  t_sup_B=new double*[geometry.nX];
+  t_sup_C=new double*[geometry.nX];
+  t_F=new double*[geometry.nX];
+  t_G=new double*[geometry.nX];
+  t_S=new double*[geometry.nX];
+  t_H=new double*[geometry.nX];
+  t_R=new double*[geometry.nX];
+  t_I=new double*[geometry.nX];
+  t_BB=new double*[geometry.nX];
+
+  for(size_t i=0;i<geometry.nX;++i){
+    size_t nZ=geometry.nZ[i];
+    t_sub_A[i]=new double[nZ-2];
+    t_sub_B[i]=new double[nZ-2];
+    t_sub_C[i]=new double[nZ-2];
+    t_diag_A[i]=new double[nZ-1];
+    t_diag_B[i]=new double[nZ-1];
+    t_diag_C[i]=new double[nZ-1];
+    t_sup_A[i]=new double[nZ-2];
+    t_sup_B[i]=new double[nZ-2];
+    t_sup_C[i]=new double[nZ-2];
+    t_F[i]=new double[nZ-1];
+    t_G[i]=new double[nZ-1];
+    t_S[i]=new double[nZ-1];
+    t_H[i]=new double[nZ-1];
+    t_R[i]=new double[nZ-1];
+    t_I[i]=new double[nZ-1];
+    t_BB[i]=new double[nZ-1];
+  }
+
+  delete initial_state;
+  initial_state=nullptr;
+
+  file.close();
+  step=0;
+  m=1;
+  scheme_error=0;
+}
+
+void
+Kernel::initSolutions(){
+  Solution* S0=new Solution();
+  S0->init(geometry);
+  solution[0]=S0;
+  for(size_t i=1;i<=Time::nT;++i) solution[i]=nullptr;
+  for(size_t i=0;i<Time::nT+2;++i){
+    Solution* S=new Solution();
+    S->init(geometry);
+    pool_solutions.push(S);
+  }
+}
+
+bool
+Kernel::next(){
+  //cout<<"[Kernel::next]"<<endl;
+  if(m>1) m=m/2;
+  spaceSolution();
+  if(scheme_error>0) return false;
+  ++step;
+  //cout<<"[next] done"<<endl;
+  return true;
+}
+
+void
+Kernel::spaceSolution(){
+  //cout<<"[Kernel::spaceSolution]"<<endl;
+  size_t n=step+1; //We want solution at time n
+  size_t k=0; //we work on time interval [n-1+k/m,n-1+(k+1)/m]
+  Solution *S=solution[step];
+  while(m<=max_time_subdivisions){
+    //cout<<" k,m = "<<k<<','<<m<<endl;
+    S=Piccard(S);
+    if(S==nullptr){
+      m*=2;
+      k=0;
+      S=solution[step];
+    }
+    else{
+      ++k;
+      if(k==m) break;
+    }
+  }
+  if(m>max_time_subdivisions){
+    cerr<<"Max time subdivision reached"<<endl;
+    scheme_error=1;
+  }
+  //cout<<"[spaceSolution] done"<<endl;
+}
+
+Solution*
+Kernel::Piccard(Solution* Sin){
+  //Return a valid solution or nullptr if not converge
+  //cout<<"[Kernel::Piccard]"<<endl;
+  memcpy(Sin.l,Sin.hsat,sizeof(double)*geometry.nX);
+  //Compute Pl with l=hsat of input solution and P=pressure of input_solution
+
+  Solution* Spl=pool_solutions.pop();
+  compute_Pl(*Spl,Sin.hsat,Sin.P);
+  double dt=Time::dT/m; // util ????
+
+  size_t count=0;
+
+  //Specicy to compute Richards on all columns
+  for(size_t i=0;i<geometry.nX;++i) indice_x_Richards[i]=true;
+
+  double error=0;
+  Solution* Snew=allVerticalRichards(Sin,Sin,Spl);
+  if(Snew==nullptr) return false; //allVerti... must compute sol_out.P and sol_out.hsat
+  //cout<<"[Piccard] done"<<endl;
+  //Quitte de facon anticipe
+  return true;
+  compute_l(Snew,error);
+  compute_Pl(Snew,Snew.l,Snew.P);
+  if(!horizontalProblem(Sin,Snew,error)) return false; //Set hydr of Snew voir mettre hydr_in
+  if(!overlandProblem(Sin,Sin,Snew,error)) return false; //Set hov of Snew
+
+  //Comput error
+  double error_prev=error;
+
+  while(error>=tolerence_Piccard and count<max_iterations_Piccard and (abs(error_prev-error)>tolerence_Piccard/100 or error<oscilation_Piccard)){
+    Solution* Sold=Snew;
+    error_prev=error;
+    ++count;
+    error=0;
+    Snew=allVerticalRichards(Sin,Sold,Sold);
+    if(Snew==nullptr) return false;
+    compute_l(Snew,error);
+    compute_Pl(Snew,Snew.l,Snew.P);
+    if(!horizontalProblem(Sold,Snew,error)) return false;
+    if(!overlandProblem(Sin,Sold,Snew,error)) return false;
+  }
+
+  compute_mass(Snew);
+  Snew.nstep_Piccard=count;
+  //cout<<"[Piccard] done"<<endl;
+  if(error>=tolerence_Piccard){
+    return nullptr;
+    //Faire attention à la pile -> rempiler des solutions
+  }
+  return Snew;
+}
+
+
+void
+Kernel::compute_l(Solution& S,double& error) {
+  //Update S.l using S.hsat
+  bool e=0;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double a=S.l[ix];
+    double b=max(S.hsat[ix],a);
+    S.l[ix]=b;
+    double t=b-a;
+    e+=t*t;
+  }
+  error+=sqrt(e);
+}
+
+void
+Kernel::compute_Pl(Solution& S,const double* h,double** P){
+  //cout<<"[Kernel::compute_Pl]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double* Px=P[ix];
+    if(h[ix]==geometry.hsoil[ix]){
+      S.Pl[ix]=Px[geometry.nZ[ix]-1];
+      //if(S.Pl[ix]!=-2000) cout<<"error "<<ix<<endl;
+    }
+    else{
+      size_t a=(h[ix]-geometry.hbot[ix])/geometry.dZ[ix];
+      double p1=Px[a];
+      double p2=Px[a+1];
+      S.Pl[ix]=p1+(p2-p1)/geometry.dZ[ix]*(h[ix]-geometry.Z[ix][a]);
+      //cout<<ix<<" : "<<S.Pl[ix]<<endl;
+
+    }
+  }
+}
+
+bool
+Kernel::allVerticalRichards(Solution& S_0,Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::allVerticalRichards]"<<endl;
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    if(indice_x_Richards[ix]){
+      if(!Richards1dEvolutiveTime(ix,S_0,S_s,S_sp)) return false;
+    }
+  }
+  //cout<<"[Kernel::allVerticalRichards] done"<<endl;
+  return true;
+}
+
+bool
+Kernel::Richards1dEvolutiveTime(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp){
+  //cout<<"[Kernel::Richards1dEvolutiveTime]"<<endl;
+  //cout<<" ix = "<<ix<<endl;
+  //Warning flux_bot and div_w are always zero
+  double flux_bot=bottomFluxRichards(ix,S_s);
+  compute_div_w(ix,S_s);
+  size_t q=1;
+  size_t l=0;
+  bool conv;
+  //cout<<"Time::dT "<<Time::dT<<endl;
+  double dt=Time::dT/m*3600; //in second
+  while(l<q and q<=max_Richards_time_subdivisions){
+    //cout<<" l,q = "<<l<<','<<q<<endl;
+    if(l=0){
+      conv=Richards1d(ix,S_0,S_s,S_sp,dt,flux_bot);
+    }
+    else{
+      conv=Richards1d(ix,S_0,S_sp,S_sp,dt,flux_bot);
+    }
+    if(!conv){
+      q*=2;
+      dt/=2;
+      l=0;
+    }
+    else{
+      ++l;
+    }
+  }
+  //{cout<<"Continue ?"<<endl;char a;cin>>a;}
+
+  if(q>max_Richards_time_subdivisions) return false;
+  return true;
+}
+
+void
+Kernel::compute_mass(Solution& S){
+
+}
+
+bool
+Kernel::horizontalProblem(Solution& S,double& error){
+  //Compute hydr from P, l, Pl, sources, time
+  //return error otherwise
+  return true;
+}
+
+bool
+Kernel::overlandProblem(const Solution& S_in,Solution& S_out,double& error){
+  //Compute hov
+  //return error otherwise
+  return true;
+}
+
+void
+Kernel::saturationInterface(const Solution& S_in,Solution& S_out){
+/*  Solution& S_in=*solution[sol_in];
+  Solution& S_out=*solution[sol_out];
+  for(size_t ix=0;ix<geometry.nX;++ix){
+    double* Px_in=S_in.P[ix];
+    size_t iz=0;
+    for(;iz<geometry.nZ[ix] and Px_in[iz]>Psat;++iz);
+    if(iz>0 and Px_in[iz]<=Psat){
+      S_out.hsat[ix]=(Psat-Px_in[iz-1])*geometry.dZ[ix]/(Px_in[iz]-Px_in[iz-1])+geometry.Z[ix][iz-1];
+    }
+    else{
+      S_out.hsat[ix]=(iz==0)?geometry.hbot[ix]:geometry.hsoil[ix];
+    }
+  }*/
+}
+
+double
+Kernel::bottomFluxRichards(size_t ix,const Solution& S){
+  return 0;
+}
+
+void
+Kernel::compute_div_w(size_t ix,const Solution& S){
+  double* div_w_ix=div_w[ix];
+  for(size_t iz=0;iz<geometry.nZ[ix];++iz){
+    div_w_ix[iz]=0;
+  }
+}
+
+bool
+Kernel::Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot){
+  //cout<<"[Richards1d]"<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double* P0=S_0.P[ix];
+  double* Ps=S_s.P[ix];
+  double* Psp=S_sp.P[ix];
+  double n_P0=norm2(P0,nZ);
+  if(weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,1,n_P0)) return true;
+  return weightedRichards1d(ix,P0,Ps,Psp,dt,flux_bot,0.5,n_P0);
+}
+
+bool
+Kernel::weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0){
+//  cout<<"[weightedRichards1d]"<<endl;
+  //cout<<"r = "<<r<<endl;
+  //cout<<"ix = "<<ix<<endl;
+  size_t nZ=geometry.nZ[ix];
+  double n0=norm2(P0,nZ);
+  double* Pi=Ps;
+  double* Po=P_temp[ix];
+  size_t count=0;
+  double e=+inf;
+  while(e>=tolerence_Richards and count<max_iterations_Richards){
+    ++count;
+    solveSystemRichards(ix,P0,Pi,Po,flux_bot,dt,count);
+    e=error2(Pi,Po,nZ)/n0;
+    //cout<<count<<" : "<<e<<" vs "<<tolerence_Richards<<endl;
+    if(count==1){
+      Pi=Po;
+      Po=Psp;
+    }
+    else{
+      swap(Pi,Po);
+    }
+  }
+  if(e<tolerence_Richards){
+    //cout<<"[converge]"<<endl;
+    memcpy(Psp,Pi,nZ*sizeof(double));
+    return true;
+  }
+  //cout<<"[don't converge]"<<endl;
+  return false;
+}
+
+void
+Kernel::solveSystemRichards(size_t ix,double* P0,double* P,double* Pout,double flux_bot,double dt,size_t count){
+  //TODO treat fpump
+  /*cout<<"[solveSystemRichards]"<<endl;
+  cout<<"P0 = "<<P0<<endl;
+  cout<<"Pin = "<<P<<endl;
+  cout<<"Pout = "<<Pout<<endl;*/
+
+  size_t nZ=geometry.nZ[ix];
+  assert(nZ>=3);
+  double dZ=geometry.dZ[ix];
+
+
+
+  double* sup_A=t_sup_A[ix];
+  double* diag_A=t_diag_A[ix];
+  double* sub_A=t_sub_A[ix];
+  double* sup_B=t_sup_B[ix];
+  double* diag_B=t_diag_B[ix];
+  double* sub_B=t_sub_B[ix];
+  double* sup_C=t_sup_C[ix];
+  double* diag_C=t_diag_C[ix];
+  double* sub_C=t_sub_C[ix];
+  double* F=t_F[ix];
+  double* G=t_G[ix];
+  double* S=t_S[ix];
+  double* H=t_H[ix];
+  double* R=t_R[ix];
+  double* I=t_I[ix];
+  double* BB=t_BB[ix];
+
+  //Compute A
+  diag_A[0]=(Physics::phi*dZ*Physics::ds(P[0]))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    diag_A[i]=(Physics::phi*dZ*Physics::ds(P[i]))/dt;
+  }
+
+
+  //Compute B
+
+  //TODO : A optimiser boucle par rapport aux appels de Kr
+  double alpha=Physics::k0/(2*Physics::rho*Physics::g*dZ);
+  diag_B[0]=alpha*(Physics::kr(P[0])+Physics::kr(P[1]));
+  sup_B[0]=-diag_B[0];
+  diag_B[nZ-2]=alpha*(Physics::kr(P[nZ-3])+2*Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]));
+  sub_B[nZ-3]=-alpha*(Physics::kr(P[nZ-3])+Physics::kr(P[nZ-2]));
+  for(size_t i=1;i<nZ-2;++i){
+    diag_B[i]=alpha*(Physics::kr(P[i-1])+2*Physics::kr(P[i])+Physics::kr(P[i+1]));
+    sub_B[i-1]=-alpha*(Physics::kr(P[i-1])+Physics::kr(P[i]));
+    sup_B[i]=-alpha*(Physics::kr(P[i])+Physics::kr(P[i+1]));
+  }
+
+
+  //Compute C
+  double hk=Physics::k0/2;
+  double temp=1/(dZ*Physics::rho*Physics::g);
+  diag_C[0]=-hk*Physics::dkr(P[0])*(temp*(P[1]-P[0])+1);
+  sup_C[0]=-hk*Physics::dkr(P[1])*(temp*(P[1]-P[0])+1);
+  diag_C[nZ-2]=hk*Physics::dkr(P[nZ-3])*temp*(-P[nZ-3]+2*P[nZ-2]-P[nZ-1]);
+  sub_C[nZ-3]=hk*Physics::dkr(P[nZ-3])*(temp*(P[nZ-2]-P[nZ-3])+1);
+  for(size_t i=1;i<nZ-2;++i){
+    diag_C[i]=hk*Physics::dkr(P[i])*temp*(-P[i-1]+2*P[i]-P[i+1]);
+    sub_C[i-1]=hk*Physics::dkr(P[i-1])*(temp*(P[i]-P[i-1])+1);
+    sup_C[i]=-hk*Physics::dkr(P[i+1])*(temp*(P[i+1]-P[i])+1);
+  }
+
+  //Compute G
+  clear(G,nZ-1);
+  G[nZ-2]=-alpha*(Physics::kr(P[nZ-2])+Physics::kr(P[nZ-1]))*P[nZ-1];
+
+  //Compute S
+  S[0]=-hk*(Physics::kr(P[0])+Physics::kr(P[1]));
+  for(size_t i=1;i<nZ-1;++i){
+    S[i]=hk*(Physics::kr(P[i-1])-Physics::kr(P[i+1]));
+  }
+
+  //Compute H
+  clear(H,nZ-1);
+  H[nZ-2]=-hk*Physics::dkr(P[nZ-1])*(temp*(P[nZ-1]-P[nZ-2])+1);
+
+  //Compute R
+  clear(R,nZ-1);
+  R[0]=-flux_bot;
+
+  //Compute I
+  I[0]=(Physics::phi*dZ*(Physics::s(P[0])-Physics::s(P0[0])))/(2*dt);
+  for(size_t i=1;i<nZ-1;++i){
+    I[i]=(Physics::phi*dZ*(Physics::s(P[i])-Physics::s(P0[i])))/dt;
+  }
+
+  //Compute BB
+  clear(BB,nZ-1);
+  //TODO : Add BB computation from fpump
+
+  //Compute F
+  F[0]=div_w[ix][0]*dZ+R[0]-G[0]-I[0]-S[0]-H[0]+(diag_A[0]+diag_C[0])*P[0]+sup_C[0]*P[1];
+  for(size_t i=1;i<nZ-2;++i){
+    F[i]=div_w[ix][i]*dZ+R[i]-G[i]-I[i]-S[i]-H[i]+(diag_A[i]+diag_C[i])*P[i]+sub_C[i-1]*P[i-1]+sup_C[i]*P[i+1];
+  }
+  F[nZ-2]=div_w[ix][nZ-2]*dZ+R[nZ-2]-G[nZ-2]-I[nZ-2]-S[nZ-2]-H[nZ-2]+(diag_A[nZ-2]+diag_C[nZ-2])*P[nZ-2]+sub_C[nZ-3]*P[nZ-3];
+
+// /  if(ix==39)  display("F",F,nZ-1);
+  //TODO Add contribution of BB in F
+
+  for(size_t i=0;i<nZ-2;++i){
+    sub_A[i]=(sub_B[i]+sub_C[i]);
+    diag_A[i]+=(diag_B[i]+diag_C[i]);
+    sup_A[i]=(sup_B[i]+sup_C[i]);
+  }
+  diag_A[nZ-2]+=(diag_B[nZ-2]+diag_C[nZ-2]);
+
+
+
+ Thomas(nZ-1,sub_A,diag_A,sup_A,F,Pout);
+ Pout[nZ-1]=P[nZ-1];
+}

backup/old/kernel.hpp

@@ -0,0 +1,94 @@
+#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;
+
+class Kernel:public InputData{
+private:
+  static constexpr double tolerence_Piccard=1e-6;
+  static constexpr double oscilation_Piccard=1e-4;
+  static constexpr size_t max_iterations_Piccard=50;
+  static constexpr size_t max_time_subdivisions=32;
+
+  static constexpr size_t max_iterations_Richards=50;
+  static constexpr double tolerence_Richards=1e-10;
+
+  //For which colum we have to compute Richards
+  bool* indice_x_Richards;
+  int scheme_error;
+
+
+
+  void initSolutions();
+
+  void spaceSolution();
+
+  double bottomFluxRichards(size_t ix,const Solution& S_k);
+  Solution* Piccard(Solution* Sin);
+  bool Richards1dEvolutiveTime(size_t ix,const Solution& S_in,const Solution& S_k,Solution& S_kp);
+  bool allVerticalRichards(Solution& S_in,Solution& S_k,Solution& S_kp);
+  void compute_l(Solution& S,double& error);
+  void compute_Pl(Solution& S,const double* h,double** P);
+  void compute_mass(Solution& S);
+  void compute_div_w(size_t ix,const Solution& S);
+  bool horizontalProblem(Solution& S,double& error);
+  bool overlandProblem(const Solution& S_in,Solution& S_out,double& error);
+  void saturationInterface(const Solution& S_in,Solution& S_out);
+  bool Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot);
+  bool weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0);
+  void solveSystemRichards(size_t ix,double* P0,double* Pin,double* Pout,double flux_bot,double dt,size_t count);
+//  double computeError(const Solution& S_in,const Solution& S_out,size_t sol_temp);
+public:
+  size_t step;
+  //Number of sub time intervals of [step,step+1]
+  size_t m;
+  Solution* solution;
+  double** div_w;
+  double** P_temp;
+
+  double** t_sub_A;
+  double** t_sub_B;
+  double** t_sub_C;
+  double** t_diag_A;
+  double** t_diag_B;
+  double** t_diag_C;
+  double** t_sup_A;
+  double** t_sup_B;
+  double** t_sup_C;
+  double** t_F;
+  double** t_G;
+  double** t_S;
+  double** t_H;
+  double** t_R;
+  double** t_I;
+  double** t_BB;
+  size_t nZ(size_t ix);
+  double Z(size_t ix,size_t iz);
+
+  Kernel(string filename);
+  void init(fstream& file);
+  bool next();
+
+};
+
+inline size_t
+Kernel::nZ(size_t ix){
+  return geometry.nZ[ix];
+}
+
+inline double
+Kernel::Z(size_t ix,size_t iz){
+  return geometry.Z[ix][iz];
+}
+
+#endif

src/debug.hpp

@@ -0,0 +1,10 @@
+#ifndef DEBUG_HPP
+#define DEBUG_HPP
+
+class Debug{
+public:
+  static const int level=1;
+  static const size_t ix=0;
+};
+
+#endif

src/kernel/all_vertical_richards.cpp

@@ -0,0 +1,89 @@
+#include "all_vertical_richards.hpp"
+
+namespace Kernel{
+  void
+  AllVerticalRichards::init(const Geometry* geometry_){
+    geometry=geometry_;
+    indice_x_Richards=new bool[geometry->nX];
+    richards_evolutive_time=new RichardsEvolutiveTime[geometry->nX];
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      richards_evolutive_time[ix].init(ix,geometry);
+    }
+  }
+
+  void
+  AllVerticalRichards::init_indice_x_Richards(){
+    for(size_t x=0;x<geometry->nX;++x){
+      indice_x_Richards[x]=true;
+    }
+  }
+
+  void
+  AllVerticalRichards::run(){
+    size_t i_left,i_middle,i_right;
+    if(Debug::level>0) cout<<"   [AllVerticalRichards::run] start"<<endl;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      if(indice_x_Richards[ix]){
+        RichardsEvolutiveTime& r=richards_evolutive_time[ix];
+        r.dt=dt;
+        r.init_P=init_P[ix];
+        r.previous_P=previous_P[ix];
+        r.P=P[ix];
+        if(ix==0){
+          i_left=0;
+          i_middle=1;
+          i_right=2;
+        }
+        else if(ix==geometry->nX-1){
+          i_left=geometry->nX-3;
+          i_middle=geometry->nX-2;
+          i_right=geometry->nX-1;
+        }
+        else{
+          i_left=ix-1;
+          i_middle=ix;
+          i_right=ix+1;
+        }
+        r.hydr_left=hydr[i_left];
+        r.l_left=l[i_left];
+        r.Pl_left=Pl[i_left];
+        r.hydr_middle=hydr[i_middle];
+        r.l_middle=l[i_middle];
+        r.Pl_middle=Pl[i_middle];
+        r.hydr_right=hydr[i_right];
+        r.l_right=l[i_right];
+        r.Pl_right=Pl[i_right];
+        r.run();
+        P[ix]=r.P;
+      }
+    }
+    //[If Paraller wait for sync]
+    has_converged=true;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      if(!richards_evolutive_time[ix].has_converged){
+        has_converged=false;
+      }
+    }
+    if(has_converged) compute_hsat();
+    if(Debug::level>0){
+      cout<<"   [AllVerticalRichards::run] ";
+      if(not has_converged) cout<<"not";
+      cout<<"converge"<<endl<<"   [AllVerticalRichards::run] end"<<endl;
+    }
+  }
+
+  void
+  AllVerticalRichards::compute_hsat(){
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double* Px=previous_P[ix];
+      size_t iz=0;
+      for(;iz<geometry->nZ[ix] and Px[iz]>Psat;++iz);
+      if(iz>0 and Px[iz]<=Psat){
+        hsat[ix]=(Psat-Px[iz-1])*geometry->dZ[ix]/(Px[iz]-Px[iz-1])+geometry->Z[ix][iz-1];
+      }
+      else{
+        hsat[ix]=(iz==0)?geometry->hbot[ix]:geometry->hsoil[ix];
+      }
+    }
+  }
+}

src/kernel/all_vertical_richards.hpp

@@ -0,0 +1,42 @@
+#ifndef ALL_VERTICAL_RICHARDS_HPP
+#define ALL_VERTICAL_RICHARDS_HPP
+
+#include "geometry.hpp"
+#include "richards_evolutive_time.hpp"
+
+namespace Kernel{
+  class AllVerticalRichards{
+  private:
+    const Geometry* geometry;
+    bool* indice_x_Richards;
+    RichardsEvolutiveTime* richards_evolutive_time;
+
+  public:
+    //Input
+    double dt;
+    double** init_P; //P_0
+    double** previous_P; //P_{k-1}
+    double* hydr;
+    double* l;
+    double* Pl;
+    double* Psoil;
+    //Output
+    double** P; //P_k
+    double* hsat;
+    bool has_converged;
+
+    AllVerticalRichards();
+    void init(const Geometry* geometry);
+    void init_indice_x_Richards();
+    void compute_hsat();
+    void run();
+  };
+
+  inline
+  AllVerticalRichards::AllVerticalRichards(){
+    richards_evolutive_time=nullptr;
+  }
+
+
+}
+#endif

src/kernel/geometry.cpp

@@ -0,0 +1,94 @@
+#include "geometry.hpp"
+
+double inf = numeric_limits<double>::infinity();
+
+Geometry::Geometry(){
+  lX=0;
+  nX=0;
+  dX=0;
+  hsoil=nullptr;
+  dhsoil=nullptr;
+  hbot=nullptr;
+  dhbot=nullptr;
+  nZ=nullptr;
+  dZ=nullptr;
+  Z=nullptr;
+}
+
+
+Geometry::~Geometry(){
+  if(hsoil!=nullptr){
+    delete[] hsoil;
+    delete[] dhsoil;
+    delete[] hbot;
+    delete[] dhbot;
+    delete[] dZ;
+    delete[] nZ;
+    for(size_t i=0;i<nX;++i){
+      delete[] Z[i];
+    }
+    delete[] Z;
+  }
+}
+
+void
+Geometry::initZ(double dZ_avg,bool init){
+  if(init){
+    nZ=new size_t[nX];
+    dZ=new double[nX];
+    Z=new double*[nX];
+  }
+  for(size_t k=0;k<nX;++k){
+    double d=hsoil[k]-hbot[k];
+    size_t n=d/dZ_avg;
+    double dz=d/n;
+    dZ[k]=dz;
+    nZ[k]=n+1;
+    if(init) Z[k]=new double[n+1];
+    for(size_t j=0;j<n+1;++j){
+      Z[k][j]=hbot[k]+j*dz;
+    }
+  }
+}
+
+void
+Geometry::save(fstream& file) {
+  file.write((char*)&lX,sizeof(double));
+  file.write((char*)&nX,sizeof(size_t));
+  file.write((char*)&dX,sizeof(double));
+  file.write((char*)hsoil,nX*sizeof(double));
+  file.write((char*)dhsoil,nX*sizeof(double));
+  file.write((char*)hbot,nX*sizeof(double));
+  file.write((char*)dhbot,nX*sizeof(double));
+  file.write((char*)nZ,nX*sizeof(size_t));
+  file.write((char*)dZ,nX*sizeof(double));
+  for(size_t i=0;i<nX;++i){
+    file.write((char*)Z[i],nZ[i]*sizeof(double));
+  }
+}
+
+void
+Geometry::load(fstream& file,bool init) {
+  file.read((char*)&lX,sizeof(double));
+  file.read((char*)&nX,sizeof(size_t));
+  file.read((char*)&dX,sizeof(double));
+  if(init){
+    hsoil=new double[nX];
+    dhsoil=new double[nX];
+    hbot=new double[nX];
+    dhbot=new double[nX];
+    nZ=new size_t[nX];
+    dZ=new double[nX];
+    Z=new double*[nX];
+  }
+  file.read((char*)hsoil,nX*sizeof(double));
+  file.read((char*)dhsoil,nX*sizeof(double));
+  file.read((char*)hbot,nX*sizeof(double));
+  file.read((char*)dhbot,nX*sizeof(double));
+  file.read((char*)nZ,nX*sizeof(size_t));
+  file.read((char*)dZ,nX*sizeof(double));
+  for(size_t i=0;i<nX;++i){
+    if(init) Z[i]=new double[nZ[i]];
+    file.read((char*)Z[i],nZ[i]*sizeof(double));
+  }
+}

src/kernel/geometry.hpp

@@ -0,0 +1,64 @@
+#ifndef GEOMETRY_HPP
+#define GEOMETRY_HPP
+
+#include <iostream>
+#include <fstream>
+#include <limits>
+
+#include "math/spline.hpp"
+
+using namespace std;
+
+using Func = double (*)(double);
+
+extern double inf;
+
+//! The Geometry class contains all geometric parameters of the domain.
+class Geometry{
+public:
+  //! Horizontal length of the domain
+  double lX;
+
+  //! Number of horizontal steps
+  size_t nX;
+
+  //! Horizontal step
+  double dX;
+
+  //! Level of the soil depending on X, e.g, hsoil[k]=level of the soil at X=k*dX.
+  //! Vector of size nX.
+  double* hsoil;
+
+  //! Derivative of the soil depending on X, vector of size nX.
+  double* dhsoil;
+
+  //! Level of the bottom depending on X, e.g, hbot[k]=level of the bottom at X=k*dX.
+  //! Vector of size nX.
+  double* hbot;
+
+  //! Derivative of the bottom depending on X, vector of size nX.
+  double* dhbot;
+
+
+  //! Number of vertical step at a given X, vector of size nX.
+  size_t* nZ;
+  //! Vertical step at a given X, vector of size nX.
+  double* dZ;
+
+  //! Vertical considered positions at a given X, vector of vectors of size nX. For each k, Z[k] is a vector of size nZ[k]
+  double** Z;
+
+  //!
+  Geometry();
+
+  ~Geometry();
+
+  void initZ(double dZ_avg,bool init=true);
+
+  void save(fstream& file);
+
+  void load(fstream& file,bool init);
+};
+
+
+#endif

src/kernel/horizontal_problem.cpp

@@ -0,0 +1,50 @@
+#include "horizontal_problem.hpp"
+#include "debug.hpp"
+
+namespace Kernel{
+
+  void
+  HorizontalProblem::init(const Geometry* geometry_){
+    geometry=geometry_;
+    sub_M=new double[geometry->nX-1];
+    diag_M=new double[geometry->nX];
+    sup_M=new double[geometry->nX-1];
+    F=new double[geometry->nX];
+  }
+
+  void
+  HorizontalProblem::run(){
+    if(Debug::level>0) cout<<"   [HorizontalProblem::run] start"<<endl;
+    /*for(size_t ix=0;ix<geometry->nX;++ix){
+      hydr[ix]=previous_hydr[ix];
+    }*/
+    solve_system();
+    error=error2(previous_hydr,hydr,geometry->nX);
+    if(Debug::level>0) cout<<"   [HorizontalProblem::run] end"<<endl;
+
+  }
+
+  void
+  HorizontalProblem::solve_system(){
+    for(size_t ix=1;ix<geometry->nX-1;++ix){
+      double dx=geometry->dX;
+      double hbot=geometry->hbot[ix];
+      double lx=l[ix];
+      double a=-Physics::tilde_a(lx,hbot)*Physics::k0/(dx*dx);
+      double c=Physics::tilde_c(lx,hbot)*Physics::k0/dx;
+      double dhbot=geometry->dhbot[ix];
+      sub_M[ix-1]=a-(dhbot*c)/2;
+      diag_M[ix]=1-2*a;
+      sup_M[ix]=a+(dhbot*c)/2;
+      F[ix]=Pl[ix]/(Physics::rho*Physics::g)+lx;
+    }
+    diag_M[0]=1;
+    sup_M[0]=-1;
+    F[0]=0;
+    diag_M[geometry->nX-1]=1;
+    sub_M[geometry->nX-2]=-1;
+    F[geometry->nX-1]=0;
+    Thomas(geometry->nX,sub_M,diag_M,sup_M,F,hydr);
+  }
+
+}

src/kernel/horizontal_problem.hpp

@@ -0,0 +1,44 @@
+#ifndef HORIZONTAL_PROBLEM_HPP
+#define HORIZONTAL_PROBLEM_HPP
+
+#include "geometry.hpp"
+#include "physics.hpp"
+
+
+namespace Kernel{
+  /*! Class for Horzontal problem.
+   *! Here we code only homogeneus Neumann conditions.
+  */
+  class HorizontalProblem{
+  private:
+    const Geometry* geometry;
+    double* sup_M;
+    double* diag_M;
+    double* sub_M;
+    double* F;
+    void solve_system();
+    double compute_error();
+  public:
+    //Input
+    double* previous_hydr;
+    double* l;
+    double* Pl;
+    //Output
+    double* hydr;
+    double error;
+
+    HorizontalProblem();
+
+    void init(const Geometry* geometry);
+
+    void run();
+  };
+
+  inline HorizontalProblem::HorizontalProblem(){
+  }
+
+
+
+}
+
+#endif

src/kernel/initial_state.cpp

@@ -0,0 +1,166 @@
+#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);
+  }
+}

src/kernel/initial_state.hpp

@@ -0,0 +1,45 @@
+#ifndef INITIAL_STATE_HPP
+#define INITIAL_STATE_HPP
+
+#include <list>
+#include <fstream>
+#include "physics.hpp"
+#include "geometry.hpp"
+#include "math/spline.hpp"
+
+using namespace std;
+
+class Tank{
+public:
+  double saturation;
+  double left,right,bottom,top;
+  double delta_left,delta_right,delta_bottom,delta_top;
+  Tank();
+  void save(fstream& file);
+  void load(fstream& file);
+};
+
+class InitialState{
+private:
+  Geometry* geometry;
+public:
+  //! Vector of size nX
+  double* hsat;
+
+
+  double* hov;
+
+  //! Vector of vector of size nX. For each k, Pinit[k] is a vector of size geometry->nZ[k]
+  double** Pinit;
+  list<Tank*> tanks;
+  InitialState();
+  InitialState(Geometry* geometry);
+  ~InitialState();
+  Tank* addTank();
+  void removeTank(Tank* tank);
+  void updatePressure();
+  void save(fstream& file);
+  void load(fstream& file,bool init);
+};
+
+#endif

src/kernel/input_data.cpp

@@ -0,0 +1,21 @@
+#include "input_data.hpp"
+
+void
+InputData::load(fstream& file,bool init){
+  Physics::load(file);
+  Time::load(file);
+  geometry.load(file,init);
+  initial_state->load(file,init);
+  source.load(file);
+  file.read((char*)&factor,sizeof(double));
+}
+
+void
+InputData::save(fstream& file){
+  Physics::save(file);
+  Time::save(file);
+  geometry.save(file);
+  initial_state->save(file);
+  source.save(file);
+  file.write((char*)&factor,sizeof(double));
+}

src/kernel/input_data.hpp

@@ -0,0 +1,43 @@
+#ifndef INPUT_DATA_HPP
+#define INPUT_DATA_HPP
+
+#include "physics.hpp"
+#include "time.hpp"
+#include "geometry.hpp"
+#include "initial_state.hpp"
+#include "source.hpp"
+
+
+class InputData{
+public:
+  Geometry geometry;
+  InitialState* initial_state;
+  Source source;
+  double factor;
+  InputData();
+  ~InputData();
+  void load(fstream& file,bool init=false);
+  void save(fstream& file);
+  void remove_initial_state();
+};
+
+
+
+inline
+InputData::InputData(){
+  initial_state=new InitialState(&geometry);
+}
+
+inline InputData::~InputData(){
+  remove_initial_state();
+}
+
+inline void
+InputData::remove_initial_state() {
+  if(initial_state!=nullptr){
+    delete initial_state;
+  }
+}
+
+
+#endif

src/kernel/overland.cpp

@@ -0,0 +1,15 @@
+#include "overland.hpp"
+#include "debug.hpp"
+
+namespace Kernel{
+
+  void Overland::run(){
+    if(Debug::level>0) cout<<"   [Overland::run] start"<<endl;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      hov[ix]=in_hov[ix];
+      Psoil[ix]=previous_P[ix][geometry->nZ[ix]-1];
+    }
+    error=0;
+    if(Debug::level>0) cout<<"   [Overland::run] end"<<endl;
+  }
+}

src/kernel/overland.hpp

@@ -0,0 +1,38 @@
+#ifndef OVERLAND_HPP
+#define OVERLAND_HPP
+
+#include "geometry.hpp"
+
+namespace Kernel{
+  class Overland{
+  private:
+    const Geometry* geometry;
+  public:
+    //Input
+    double* in_hov;
+    double** previous_P;
+    double* l;
+    double* Pl;
+    double* hydr;
+    //Output
+    double* Psoil; //Own by Overland object
+    double* hov;
+    double error;
+
+    Overland();
+    void init(const Geometry* geometry);
+    void run();
+  };
+
+  inline
+  Overland::Overland(){
+  }
+
+  inline void
+  Overland::init(const Geometry* geometry_){
+    geometry=geometry_;
+    Psoil=new double[geometry->nX];
+  }
+
+}
+#endif

src/kernel/physics.cpp

@@ -0,0 +1,158 @@
+#include "physics.hpp"
+//------------------------
+// Brooks and Corey model
+//------------------------
+
+double Physics::g;
+double Physics::rho;
+double Physics::phi;
+double Physics::k0;
+double Physics::nivrivsat;
+double Physics::Psec;
+Physics::Model Physics::model;
+double (*Physics::s)(double);
+double (*Physics::ds)(double);
+void (*Physics::s_ds)(double,double&,double&);
+double (*Physics::s_inv)(double);
+double (*Physics::kr)(double);
+double (*Physics::dkr)(double);
+void (*Physics::kr_dkr)(double,double&,double&);
+double Physics::model_data[6];
+
+
+#define sres Physics::model_data[1]
+#define lambda Physics::model_data[2]
+#define alpha Physics::model_data[3]
+
+void Physics::setModel(Model m){
+  model=m;
+  switch(model){
+    case BrooksCorey:
+    s=&s_BC;
+    ds=&ds_BC;
+    s_ds=&s_ds_BC;
+    s_inv=&s_inv_BC;
+    kr=&kr_BC;
+    dkr=&dkr_BC;
+    kr_dkr=&kr_dkr_BC;
+    break;
+    default:
+    assert(false);
+  };
+}
+
+double
+Physics::s_BC(double P){
+  if(P>=Psat) return 1;
+  return sres+(1-sres)*pow(Psat/P,lambda);
+}
+
+double
+Physics::ds_BC(double P){
+  if(P>=Psat) return 0;
+  return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
+}
+
+void
+Physics::s_ds_BC(double P,double& v,double& dv){
+  if(P>=Psat){
+    v=1;
+    dv=0;
+  }
+  else{
+    double t=(1-sres)*pow(Psat/P,lambda);
+    v=sres+t;
+    dv=-(lambda*t)/P;
+  }
+}
+
+double
+Physics::s_inv_BC(double S){
+  return pow((S-sres)/(1-sres),-1/lambda)*Psat;
+}
+
+double
+Physics::kr_BC(double P){
+  if(P>=Psat) return 1;
+  return pow(Psat/P,alpha);
+}
+
+double
+Physics::dkr_BC(double P){
+  if(P>=Psat) return 0;
+  return -alpha*pow(Psat/P,alpha)/P;
+}
+
+void
+Physics::kr_dkr_BC(double P,double& v,double& dv){
+  if(P>=Psat){
+    v=1;
+    dv=0;
+  }
+  else{
+    double t=pow(Psat/P,alpha);
+    v=t;
+    dv=-(alpha*t)/P;
+  }
+}
+
+double
+Physics::tilde_a(double l,double hbot){
+  double t=(l-hbot);
+  return t*t/(3*k0);
+}
+
+double
+Physics::tilde_c(double l,double hbot){
+  return (l-hbot)/(2*k0);
+}
+
+void
+Physics::save(fstream& file){
+  file.write((char*)&g,sizeof(double));
+  file.write((char*)&rho,sizeof(double));
+  file.write((char*)&phi,sizeof(double));
+  file.write((char*)&k0,sizeof(double));
+  file.write((char*)&Psec,sizeof(double));
+  file.write((char*)&nivrivsat,sizeof(double));
+  file.write((char*)&model,sizeof(Model));
+  for(size_t i=0;i<max_model_parameters;++i){
+    file.write((char*)&model_data[i],sizeof(double));
+  }
+}
+
+void
+Physics::load(fstream& file){
+  file.read((char*)&g,sizeof(double));
+  file.read((char*)&rho,sizeof(double));
+  file.read((char*)&phi,sizeof(double));
+  file.read((char*)&k0,sizeof(double));
+  file.read((char*)&Psec,sizeof(double));
+  file.read((char*)&nivrivsat,sizeof(double));
+  Model m;
+  file.read((char*)&m,sizeof(Model));
+  setModel(m);
+  for(size_t i=0;i<max_model_parameters;++i){
+    file.read((char*)&model_data[i],sizeof(double));
+  }
+}
+
+double
+Physics::Psol(double hsoil,double hov){
+  double a=Psat-rho*g*nivrivsat;
+  return rho*g*(hov-hsoil)+a*nivrivsat/(hov-hsoil);
+}
+
+double
+Physics::dPsol(double hsoil,double hov){
+  double a=Psat-rho*g*nivrivsat;
+  double b=hov-hsoil;
+  return rho*g-a*nivrivsat/(b*b);
+}
+
+double
+Physics::invPsol(double hsoil,double Psol){
+  double a=Psat-rho*g*nivrivsat;
+  double b=a*nivrivsat;
+  return hsoil+(Psol+sqrt(Psol*Psol-4*rho*g*b))/(2*rho*g);
+}

src/kernel/physics.hpp

@@ -0,0 +1,106 @@
+#ifndef PHYSICS_HPP
+#define PHYSICS_HPP
+
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+//! The Physics class contains all physical parameters characterising the soil.
+//! This class contains only static members.
+class Physics{
+public:
+  enum Model{BrooksCorey};
+
+  //! Set physics model
+  static void setModel(Model model);
+
+  //! Gravity acceleration (m.s^-2)
+  static double g;
+  //! Fluid density (g.l^(-1))
+  static double rho;
+
+  //! Porosity of the soil
+  static double phi;
+  //! Conductivity of the saturated soil
+  static double k0;
+
+  static double Psec;
+  //! Characterise the water pressure at the bottom of the overland water
+  static double nivrivsat;
+
+  //! Model used
+  static Model model;
+
+  //! Return the saturation in function of the pressure
+  static double (*s)(double);
+  //! Return the derivtive of the saturation in function of the pressure
+  static double (*ds)(double);
+  //! Set the saturation and its derivative in function of the pressure
+  static void (*s_ds)(double,double&,double&);
+
+  static double (*s_inv)(double);
+
+  //! Return the relative conductivity in function of the pressure
+  static double (*kr)(double);
+  //! Return the derivtive of the relative conductivity in function of the pressure
+  static double (*dkr)(double);
+  //! Set the relative conductivity and its derivative in function of the pressure
+  static void (*kr_dkr)(double,double&,double&);
+
+  //---------------------
+  // Models descriptions
+  //---------------------
+
+  //! Datas used to define the model
+  static double model_data[6];
+
+  //! Max model Datas
+  static const size_t max_model_parameters=6;
+
+  //------------------------
+  // Brooks and Corey model
+  //------------------------
+
+  //model_data[0] -> psat : minimal pressure such that s(psat)=1
+  //model_data[1] -> sres : residual pressure
+
+  //model_data[2] -> lambda
+  //model_data[3] -> alpha
+
+  //! Brooks and Corey saturation map
+  static double s_BC(double P);
+  //! Brooks and Corey derivative of the saturation map
+  static double ds_BC(double P);
+  //! Brooks and Corey saturation and its derivative setter
+  static void (s_ds_BC)(double P,double& v,double& dv);
+
+  static double s_inv_BC(double s);
+
+  //! Brooks and Corey relative conductivity map
+  static double kr_BC(double P);
+  //! Brooks and Corey derivative of the relative conductivity map
+  static double dkr_BC(double P);
+  //! Brooks and Corey relative conductivity and its derivative setter
+  static void (kr_dkr_BC)(double P,double& v,double& dv);
+
+  static double Psol(double hsoil,double hov);
+  static double dPsol(double hsoil,double hov);
+  static double invPsol(double hsoil,double Psol);
+
+  static double tilde_a(double l,double hbot);
+  static double tilde_c (double l,double hbot);
+
+  //--------------------
+  // File input/output
+  //--------------------
+
+  static void save(fstream& file);
+  static void load(fstream& file);
+
+};
+
+#define Psat Physics::model_data[0]
+#endif

src/kernel/piccard.cpp

@@ -0,0 +1,168 @@
+#include "piccard.hpp"
+
+namespace Kernel{
+
+  Piccard::Piccard(){
+    previous_solution=nullptr;
+    new_solution=nullptr;
+    l=nullptr;
+    Pl=nullptr;
+  }
+
+  void Piccard::init(const Geometry* geometry_){
+    geometry=geometry_;
+    horizontal_problem.init(geometry);
+    overland.init(geometry);
+    all_vertical_richards.init(geometry);
+    if(l!=nullptr) delete[] l;
+    if(Pl!=nullptr) delete[] Pl;
+    l=new double[geometry->nX];
+    Pl=new double[geometry->nX];
+  }
+
+  void
+  Piccard::run(){
+    if(Debug::level>0) cout<<"  [Piccard::run] start"<<endl;
+    double error=0;
+    //Copy s_in.hsat in l
+    memcpy(l,previous_solution->hsat,sizeof(double)*geometry->nX);
+    //Compute Pl from s_in.hsat and s_in.P
+
+    compute_Pl(previous_solution->P);
+
+
+    //Compute hydr from s_in.hydr, s_in.hsat and Pl
+    horizontal_problem.previous_hydr=previous_solution->hydr;
+    horizontal_problem.l=previous_solution->hsat;
+    horizontal_problem.Pl=Pl;
+    horizontal_problem.hydr=new_solution->hydr;
+    horizontal_problem.run();
+    new_solution->hydr=horizontal_problem.hydr;
+    error+=horizontal_problem.error;
+
+    //Compute Overland
+    overland.in_hov=previous_solution->hov;
+    overland.previous_P=previous_solution->P;
+    overland.l=previous_solution->hsat;
+    overland.Pl=Pl;
+    overland.hydr=horizontal_problem.hydr;
+    overland.hov=new_solution->hov;
+    overland.run();
+    new_solution->hov=overland.hov;
+    error+=overland.error;
+
+    //----------------------------------------------
+    // Apply AllVerticalRichads algorithm to obtain
+    // - P
+    // - hsat
+    //---------------------------------------------
+    all_vertical_richards.init_indice_x_Richards();
+    all_vertical_richards.dt=dt;
+    all_vertical_richards.init_P=previous_solution->P;
+    all_vertical_richards.previous_P=previous_solution->P;
+    all_vertical_richards.hydr=horizontal_problem.hydr;
+    all_vertical_richards.l=previous_solution->hsat;
+    all_vertical_richards.Pl=Pl;
+    all_vertical_richards.Psoil=overland.Psoil;
+    //Specify output (may can change during the run of all_vertical_richards)
+    all_vertical_richards.hsat=new_solution->hsat;
+    all_vertical_richards.P=new_solution->P;
+    all_vertical_richards.run();
+    new_solution->hsat=all_vertical_richards.hsat;
+    new_solution->P=all_vertical_richards.P;
+
+    size_t k=0;
+    double previous_error=error;
+
+    while(error>=tolerence_Piccard and k<max_iterations_Piccard and (abs(previous_error-error)>tolerence_Piccard/100 or error<oscilation_Piccard)){
+      previous_error=error;
+      error=0;
+      ++k;
+      if(Debug::level>1) cout<<"  [Piccard::run] k = "<<k<<endl;
+
+      compute_l(l,all_vertical_richards.hsat,error);
+
+      compute_Pl(all_vertical_richards.P);
+
+      horizontal_problem.previous_hydr=horizontal_problem.hydr;
+      horizontal_problem.l=l;
+      horizontal_problem.Pl=Pl;
+      horizontal_problem.hydr=new_solution->hydr;
+      horizontal_problem.run();
+      new_solution->hydr=horizontal_problem.hydr;
+
+      error+=horizontal_problem.error;
+
+      //Compute Overland
+      overland.in_hov=previous_solution->hov;
+      overland.previous_P=all_vertical_richards.P;
+      overland.l=l;
+      overland.Pl=Pl;
+      overland.hydr=horizontal_problem.hydr;
+      overland.hov=new_solution->hov;
+      overland.run();
+      new_solution->hov=overland.hov;
+      error+=overland.error;
+
+      //Voir calcul indice_x_Richards
+      all_vertical_richards.init_P=previous_solution->P; //P_0
+      all_vertical_richards.previous_P=all_vertical_richards.P; //P_{k-1}
+      all_vertical_richards.hydr=horizontal_problem.hydr;
+      all_vertical_richards.l=l;
+      all_vertical_richards.Pl=Pl;
+      all_vertical_richards.Psoil=overland.Psoil;
+      //Specify output (may can change during the run of all_vertical_richards)
+      all_vertical_richards.hsat=new_solution->hsat;
+      all_vertical_richards.P=new_solution->P;
+      all_vertical_richards.run();
+      new_solution->hsat=all_vertical_richards.hsat;
+      new_solution->P=all_vertical_richards.P;
+      //we get P_k : all_vertical_richards.P
+      if(!all_vertical_richards.has_converged){
+        has_converged=false;
+        return;
+      }
+    }
+    if(error>=tolerence_Piccard){
+      if(Debug::level>1) cout<<"  [Piccard]::run not converge"<<endl;
+      has_converged=false;
+      //Voir nettoyage memoire
+      return;
+    }
+    if(Debug::level>1) cout<<"  [Piccard]::run converge"<<endl;
+    swap(l,new_solution->l);
+    swap(Pl,new_solution->Pl);
+    has_converged=true;
+    //Voir nettoyage memoire
+    return;
+  }
+
+  void
+  Piccard::compute_l(double* h,double* hsat,double& error) {
+    bool e=0;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double a=h[ix];
+      double b=max(hsat[ix],a);
+      l[ix]=b;
+      double t=b-a;
+      e+=t*t;
+    }
+    error+=sqrt(e);
+  }
+
+  void
+  Piccard::compute_Pl(double** P){
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double* Px=P[ix];
+      if(l[ix]==geometry->hsoil[ix]){
+        Pl[ix]=Px[geometry->nZ[ix]-1];
+      }
+      else{
+        size_t a=(l[ix]-geometry->hbot[ix])/geometry->dZ[ix];
+        double p1=Px[a];
+        double p2=Px[a+1];
+        Pl[ix]=p1+(p2-p1)/geometry->dZ[ix]*(l[ix]-geometry->Z[ix][a]);
+      }
+    }
+  }
+}

src/kernel/piccard.hpp

@@ -0,0 +1,59 @@
+#ifndef PICCARD_HPP
+#define PICCARD_HPP
+
+#include <cstring>
+#include "solution.hpp"
+#include "horizontal_problem.hpp"
+#include "all_vertical_richards.hpp"
+#include "overland.hpp"
+
+namespace Kernel{
+  //! A class implementing Piccard algorithm
+  class Piccard{
+  private:
+    static constexpr double tolerence_Piccard=1e-6;
+    static constexpr double oscilation_Piccard=1e-4;
+    static constexpr size_t max_iterations_Piccard=50;
+
+    //! Geometry of the simulation
+    const Geometry* geometry;
+
+    //! A pointer to an HorizontalProblem object
+    HorizontalProblem horizontal_problem;
+
+    //! A pointer to an Overland object
+    Overland overland;
+
+    //! A pointer to an AllVerticalRichards object
+    AllVerticalRichards all_vertical_richards;
+
+    double *l;
+    double *Pl;
+
+    //! Compute l from h and hsat
+    void compute_l(double* h,double* hsat,double& error);
+
+    //! Compute Pl from l and P
+    void compute_Pl(double** P);
+  public:
+    double dt;
+    Solution* previous_solution;
+    Solution* new_solution;
+    bool has_converged;
+
+    //! Construct an empty Piccard object
+    Piccard();
+
+    //! Init from a Geometry
+    void init(const Geometry* geometry);
+
+    //! Compute a new Solution form a previous once.
+    //! The delta time depends of the number of subdivisions required by Solver::space_solution.
+    void run();
+
+  };
+
+}
+
+
+#endif

src/kernel/richards.cpp

@@ -0,0 +1,194 @@
+#include "richards.hpp"
+
+namespace Kernel{
+
+  Richards::Richards(){
+    sup_A=nullptr;
+    diag_A=nullptr;
+    sub_A=nullptr;
+    sup_B=nullptr;
+    diag_B=nullptr;
+    sub_B=nullptr;
+    sup_C=nullptr;
+    diag_C=nullptr;
+    sub_C=nullptr;
+    F=nullptr;
+    G=nullptr;
+    S=nullptr;
+    H=nullptr;
+    R=nullptr;
+    I=nullptr;
+    BB=nullptr;
+  }
+
+  void
+  Richards::init(size_t ix_,size_t nZ_,double dZ_){
+    ix=ix_;
+    nZ=nZ_;
+    dZ=dZ_;
+    temp_P[0]=new double[nZ];
+    temp_P[1]=new double[nZ];
+
+    sub_A=new double[nZ-2];
+    sub_B=new double[nZ-2];
+    sub_C=new double[nZ-2];
+    diag_A=new double[nZ-1];
+    diag_B=new double[nZ-1];
+    diag_C=new double[nZ-1];
+    sup_A=new double[nZ-2];
+    sup_B=new double[nZ-2];
+    sup_C=new double[nZ-2];
+    F=new double[nZ-1];
+    G=new double[nZ-1];
+    S=new double[nZ-1];
+    H=new double[nZ-1];
+    R=new double[nZ-1];
+    I=new double[nZ-1];
+    BB=new double[nZ-1];
+
+  }
+
+  void
+  Richards::run(){
+    if(Debug::level>0 and Debug::ix==ix) cout<<"     [Richards::run] start"<<endl;
+    norm_previous_P=norm2(previous_P,nZ);
+    has_converged=weighted_run(1);
+    if(!has_converged) has_converged=weighted_run(0.5);
+    if(Debug::level>0 and Debug::ix==ix) cout<<"     [Richards::run] stop"<<endl;
+
+  }
+
+  bool
+  Richards::weighted_run(double w){
+    if(Debug::level>0 and Debug::ix==ix) cout<<"      [Richards::weighted_run] start"<<endl;
+
+    double error=numeric_limits<double>::infinity();
+    size_t count=0;
+    in_P=near_P;
+    out_P=temp_P[1];
+    while(error>=tolerence_Richards and count<max_iterations_Richards){
+      ++count;
+      if(Debug::level>1 and Debug::ix==ix) cout<<"      [Richards::weighted_run] count = "<<count<<endl;
+      solve_system();
+
+      //Computed P is in out_P=temp_P[count%2]
+      error=error2(in_P,out_P,nZ)/norm_previous_P;
+      if(count==1){
+        in_P=temp_P[1];
+        out_P=temp_P[0];
+      }
+      else{
+        swap(in_P,out_P);
+      }
+   }
+   if(error<tolerence_Richards){
+     //Last computed P is in temp_P[count%2] which is now also in_P
+     assert(in_P==temp_P[count%2]);
+     swap(P,temp_P[count%2]);
+     if(Debug::ix==ix){
+       if(Debug::level>1) cout<<"      [Richards::weighted_run] converge"<<endl;
+       if(Debug::level>0) cout<<"      [Richards::weighted_run] stop"<<endl;
+     }
+     return true;
+   }
+   if(Debug::level>0 and Debug::ix==ix){
+     cout<<"      [Richards::weighted_run] not converge"<<endl;
+     cout<<"      [Richards::weighted_run] stop"<<endl;
+   }
+   return false;
+ }
+
+ void
+ Richards::solve_system(){
+   if(Debug::level>0 and Debug::ix==ix) cout<<"       [Richards::solve_system] start"<<endl;
+   assert(nZ>=3);
+
+   //Compute A
+   diag_A[0]=(Physics::phi*dZ*Physics::ds(in_P[0]))/(2*dt);
+   for(size_t i=1;i<nZ-1;++i){
+     diag_A[i]=(Physics::phi*dZ*Physics::ds(in_P[i]))/dt;
+   }
+
+
+   //Compute B
+
+   //TODO : A optimiser boucle par rapport aux appels de Kr
+   double alpha=Physics::k0/(2*Physics::rho*Physics::g*dZ);
+   diag_B[0]=alpha*(Physics::kr(in_P[0])+Physics::kr(in_P[1]));
+   sup_B[0]=-diag_B[0];
+   diag_B[nZ-2]=alpha*(Physics::kr(in_P[nZ-3])+2*Physics::kr(in_P[nZ-2])+Physics::kr(in_P[nZ-1]));
+   sub_B[nZ-3]=-alpha*(Physics::kr(in_P[nZ-3])+Physics::kr(in_P[nZ-2]));
+   for(size_t i=1;i<nZ-2;++i){
+     diag_B[i]=alpha*(Physics::kr(in_P[i-1])+2*Physics::kr(in_P[i])+Physics::kr(in_P[i+1]));
+     sub_B[i-1]=-alpha*(Physics::kr(in_P[i-1])+Physics::kr(in_P[i]));
+     sup_B[i]=-alpha*(Physics::kr(in_P[i])+Physics::kr(in_P[i+1]));
+   }
+
+
+   //Compute C
+   double hk=Physics::k0/2;
+   double temp=1/(dZ*Physics::rho*Physics::g);
+   diag_C[0]=-hk*Physics::dkr(in_P[0])*(temp*(in_P[1]-in_P[0])+1);
+   sup_C[0]=-hk*Physics::dkr(in_P[1])*(temp*(in_P[1]-in_P[0])+1);
+   diag_C[nZ-2]=hk*Physics::dkr(in_P[nZ-3])*temp*(-in_P[nZ-3]+2*in_P[nZ-2]-in_P[nZ-1]);
+   sub_C[nZ-3]=hk*Physics::dkr(in_P[nZ-3])*(temp*(in_P[nZ-2]-in_P[nZ-3])+1);
+   for(size_t i=1;i<nZ-2;++i){
+     diag_C[i]=hk*Physics::dkr(in_P[i])*temp*(-in_P[i-1]+2*in_P[i]-in_P[i+1]);
+     sub_C[i-1]=hk*Physics::dkr(in_P[i-1])*(temp*(in_P[i]-in_P[i-1])+1);
+     sup_C[i]=-hk*Physics::dkr(in_P[i+1])*(temp*(in_P[i+1]-in_P[i])+1);
+   }
+
+   //Compute G
+   clear(G,nZ-1);
+   G[nZ-2]=-alpha*(Physics::kr(in_P[nZ-2])+Physics::kr(in_P[nZ-1]))*in_P[nZ-1];
+
+   //Compute S
+   S[0]=-hk*(Physics::kr(in_P[0])+Physics::kr(in_P[1]));
+   for(size_t i=1;i<nZ-1;++i){
+     S[i]=hk*(Physics::kr(in_P[i-1])-Physics::kr(in_P[i+1]));
+   }
+
+   //Compute H
+   clear(H,nZ-1);
+   H[nZ-2]=-hk*Physics::dkr(in_P[nZ-1])*(temp*(in_P[nZ-1]-in_P[nZ-2])+1);
+
+   //Compute R
+   clear(R,nZ-1);
+   R[0]=-flux_bot;
+
+   //Compute I
+   I[0]=(Physics::phi*dZ*(Physics::s(in_P[0])-Physics::s(previous_P[0])))/(2*dt);
+   for(size_t i=1;i<nZ-1;++i){
+     I[i]=(Physics::phi*dZ*(Physics::s(in_P[i])-Physics::s(previous_P[i])))/dt;
+   }
+
+   //Compute BB
+   clear(BB,nZ-1);
+   //TODO : Add BB computation from fpump
+
+   //Compute F
+   F[0]=div_w[0]*dZ+R[0]-G[0]-I[0]-S[0]-H[0]+(diag_A[0]+diag_C[0])*in_P[0]+sup_C[0]*in_P[1];
+   for(size_t i=1;i<nZ-2;++i){
+     F[i]=div_w[i]*dZ+R[i]-G[i]-I[i]-S[i]-H[i]+(diag_A[i]+diag_C[i])*in_P[i]+sub_C[i-1]*in_P[i-1]+sup_C[i]*in_P[i+1];
+   }
+   F[nZ-2]=div_w[nZ-2]*dZ+R[nZ-2]-G[nZ-2]-I[nZ-2]-S[nZ-2]-H[nZ-2]+(diag_A[nZ-2]+diag_C[nZ-2])*in_P[nZ-2]+sub_C[nZ-3]*in_P[nZ-3];
+
+ // /  if(ix==39)  display("F",F,nZ-1);
+   //TODO Add contribution of BB in F
+
+   for(size_t i=0;i<nZ-2;++i){
+     sub_A[i]=(sub_B[i]+sub_C[i]);
+     diag_A[i]+=(diag_B[i]+diag_C[i]);
+     sup_A[i]=(sup_B[i]+sup_C[i]);
+   }
+   diag_A[nZ-2]+=(diag_B[nZ-2]+diag_C[nZ-2]);
+
+  Thomas(nZ-1,sub_A,diag_A,sup_A,F,out_P);
+  out_P[nZ-1]=in_P[nZ-1];
+  if(Debug::level>0 and Debug::ix==ix){
+    cout<<"       [Richards::solve_system] out = "<<out_P<<endl;
+    cout<<"       [Richards::solve_system] stop"<<endl;
+  }
+
+ }
+}

src/kernel/richards.hpp

@@ -0,0 +1,65 @@
+#ifndef RICHARDS_HPP
+#define RICHARDS_HPP
+
+#include <iostream>
+#include <limits>
+#include <cassert>
+#include "physics.hpp"
+#include "math/algo.hpp"
+#include "debug.hpp"
+
+using namespace std;
+
+namespace Kernel{
+  class Richards{
+  private:
+
+    static constexpr size_t max_iterations_Richards=50;
+    static constexpr double tolerence_Richards=1e-10;
+
+    double* sup_A;
+    double* diag_A;
+    double* sub_A;
+    double* sup_B;
+    double* diag_B;
+    double* sub_B;
+    double* sup_C;
+    double* diag_C;
+    double* sub_C;
+    double* F;
+    double* G;
+    double* S;
+    double* H;
+    double* R;
+    double* I;
+    double* BB;
+
+    size_t ix;
+    size_t nZ;
+    double dZ;
+
+    double norm_previous_P;
+    bool weighted_run(double w);
+    void solve_system();
+    double* temp_P[2];
+    double* in_P;
+    double* out_P;
+  public:
+    double dt;
+    double flux_bot;
+    double* div_w;
+
+    double* near_P; //for Newton algorithm
+    double* previous_P; //previous relatively to n2
+    double* P;
+    bool has_converged;
+
+    Richards();
+    void init(size_t ix,size_t nZ,double dZ);
+
+    void run();
+  };
+
+
+}
+#endif

src/kernel/richards_evolutive_time.cpp

@@ -0,0 +1,86 @@
+#include "richards_evolutive_time.hpp"
+
+namespace Kernel{
+  void
+  RichardsEvolutiveTime::init(size_t ix_,const Geometry* geometry){
+    ix=ix_;
+    nZ=geometry->nZ[ix];
+    Z=geometry->Z[ix];
+    dX=geometry->dX;
+    div_w=new double[nZ];
+    richards.init(ix,nZ,geometry->dZ[ix]);
+  }
+
+  void
+  RichardsEvolutiveTime::compute_flux_bot(){
+    flux_bot=0;
+  }
+
+  void
+  RichardsEvolutiveTime::compute_div_w(){
+
+    double cst=Physics::rho*Physics::g;
+    double G_left=Pl_left/cst+l_left;
+    double G_middle=Pl_middle/cst+l_middle;
+    double G_right=Pl_right/cst+l_right;
+    double m_left=max(G_left,G_middle);
+    double m_right=max(G_middle,G_right);
+
+    for(size_t i=0;i<nZ;++i){
+      div_w[i]=(Physics::kr(cst*(m_right-Z[i]))*Physics::k0*(hydr_right-hydr_middle)-Physics::kr(cst*(m_left-Z[i]))*Physics::k0*(hydr_middle-hydr_left))/(dX*dX);
+    }
+
+
+  }
+
+  void
+  RichardsEvolutiveTime::run(){
+    if(ix==Debug::ix and Debug::level>0) cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] start"<<endl;
+    compute_flux_bot();
+    compute_div_w();
+    size_t m2=1;
+    size_t n2=0;
+
+    double dt_sub=dt;
+
+    richards.div_w=div_w;
+    richards.dt=dt_sub;
+    richards.flux_bot=flux_bot;
+    //You must specify richards.Px and also create it
+    while(n2<m2 and m2<=max_Richards_time_subdivisions){
+      if(ix==Debug::ix and Debug::level>1) cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] n2 = "<<n2<<" and m2 = "<<m2<<endl;
+      if(n2==0){
+        richards.previous_P=init_P;
+        richards.near_P=previous_P;
+      }
+      else{
+        richards.previous_P=richards.P;
+        richards.near_P=richards.P;
+      }
+      richards.P=P;
+
+      richards.run();
+      P=richards.P;
+
+      if(!richards.has_converged){
+        m2*=2;
+        dt_sub/=2;
+        richards.dt=dt_sub;
+        n2=0;
+      }
+      else{
+        ++n2;
+      }
+    }
+    has_converged=(m2<=max_Richards_time_subdivisions);
+    if(ix==Debug::ix and Debug::level>0){
+      if(Debug::level>1){
+        cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] out = "<<P<<endl;
+        cout<<"    [RichardsEvolutiveTime ix="<<ix<<"]";
+        if(not has_converged) cout<<"not";
+        cout<<" converge"<<endl;
+      }
+      cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] stop"<<endl;
+    }
+  }
+}

src/kernel/richards_evolutive_time.hpp

@@ -0,0 +1,55 @@
+#ifndef RICHARDS_EVOLUTIVE_TIME_HPP
+#define RICHARDS_EVOLUTIVE_TIME_HPP
+
+#include "geometry.hpp"
+#include "richards.hpp"
+
+namespace Kernel{
+  class RichardsEvolutiveTime{
+  private:
+    static constexpr size_t max_Richards_time_subdivisions=256;
+    size_t ix;
+    size_t nZ;
+    double* Z;
+    double dX;
+
+    double flux_bot;
+    double* div_w;
+
+    void compute_flux_bot();
+    void compute_div_w();
+    Richards richards;
+  public:
+    //Input
+    double dt;
+    double* init_P;
+    double* previous_P; //relatively to k
+
+    double hydr_left;
+    double l_left;
+    double Pl_left;
+
+    double hydr_middle;
+    double l_middle;
+    double Pl_middle;
+
+    double hydr_right;
+    double l_right;
+    double Pl_right;
+
+    //Output
+    double* P;
+    double hsat;
+    bool has_converged;
+
+    RichardsEvolutiveTime();
+    void init(size_t ix,const Geometry* geometry);
+    void run();
+  };
+
+  inline
+  RichardsEvolutiveTime::RichardsEvolutiveTime(){
+
+  }
+}
+#endif

src/kernel/solution.cpp

@@ -0,0 +1,14 @@
+#include "solution.hpp"
+
+Solution::Solution(const Geometry& geometry){
+  size_t nX=geometry.nX;
+  P=new double*[nX];
+  for(size_t ix=0;ix<nX;++ix){
+    P[ix]=new double[geometry.nZ[ix]];
+  }
+  hydr=new double[nX];
+  hov=new double[nX];
+  hsat=new double[nX];
+  l=new double[nX];
+  Pl=new double[nX];
+}

src/kernel/solution.hpp

@@ -0,0 +1,23 @@
+#ifndef SOLUTION_HPP
+#define SOLUTION_HPP
+
+#include "geometry.hpp"
+
+class Solution{
+public:
+  double** P;
+  double* hydr;
+  double* hov;
+  double* hsat;
+  double* l;
+  double* Pl;
+  double mass;
+  size_t nstep_Piccard;
+
+  Solution(const Geometry& geometry);
+
+};
+
+
+
+#endif

src/kernel/solver.cpp

@@ -0,0 +1,120 @@
+#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(){
+    if(Debug::level>0) cout<<" [Solver::spaceSolution] start"<<endl;
+    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){
+
+      if(Debug::level>1) cout<<" [Solver::spaceSolution] n1 = "<<n1<<" and m1 = "<<m1<<endl;
+      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){
+      if(Debug::level>0) cout<<" [Solver::spaceSolution] not converge"<<endl;
+      release_solution(piccard.new_solution);
+      return nullptr;
+    }
+    if(Debug::level>0) cout<<" [Solver::spaceSolution] converge"<<endl;
+    return piccard.new_solution;
+  }
+}

src/kernel/solver.hpp

@@ -0,0 +1,111 @@
+#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:
+    //! 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();
+
+  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;
+
+    //! 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(){
+    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 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

src/kernel/source.cpp

@@ -0,0 +1,257 @@
+#include "source.hpp"
+
+Pump::Pump(){
+  amplitude_init=1.e-4;
+  left_init=0.45;
+  right_init=0.55;
+  bottom_init=0.45;
+  top_init=0.55;
+  delta_left_init=0.05;
+  delta_right_init=0.05;
+  delta_top_init=0.05;
+  delta_bottom_init=0.05;
+  amplitude_final=1.e-4;
+  left_final=0.45;
+  right_final=0.55;
+  bottom_final=0.45;
+  top_final=0.55;
+  delta_left_final=0.05;
+  delta_right_final=0.05;
+  delta_top_final=0.05;
+  delta_bottom_final=0.05;
+}
+
+void
+Pump::save(fstream& file){
+  file.write((char*)&amplitude_init,sizeof(double));
+  file.write((char*)&left_init,sizeof(double));
+  file.write((char*)&right_init,sizeof(double));
+  file.write((char*)&bottom_init,sizeof(double));
+  file.write((char*)&top_init,sizeof(double));
+  file.write((char*)&delta_left_init,sizeof(double));
+  file.write((char*)&delta_right_init,sizeof(double));
+  file.write((char*)&delta_bottom_init,sizeof(double));
+  file.write((char*)&delta_top_init,sizeof(double));
+  file.write((char*)&amplitude_final,sizeof(double));
+  file.write((char*)&left_final,sizeof(double));
+  file.write((char*)&right_final,sizeof(double));
+  file.write((char*)&bottom_final,sizeof(double));
+  file.write((char*)&top_final,sizeof(double));
+  file.write((char*)&delta_left_final,sizeof(double));
+  file.write((char*)&delta_right_final,sizeof(double));
+  file.write((char*)&delta_bottom_final,sizeof(double));
+  file.write((char*)&delta_top_final,sizeof(double));
+}
+
+void
+Pump::load(fstream& file){
+  file.read((char*)&amplitude_init,sizeof(double));
+  file.read((char*)&left_init,sizeof(double));
+  file.read((char*)&right_init,sizeof(double));
+  file.read((char*)&bottom_init,sizeof(double));
+  file.read((char*)&top_init,sizeof(double));
+  file.read((char*)&delta_left_init,sizeof(double));
+  file.read((char*)&delta_right_init,sizeof(double));
+  file.read((char*)&delta_bottom_init,sizeof(double));
+  file.read((char*)&delta_top_init,sizeof(double));
+  file.read((char*)&amplitude_final,sizeof(double));
+  file.read((char*)&left_final,sizeof(double));
+  file.read((char*)&right_final,sizeof(double));
+  file.read((char*)&bottom_final,sizeof(double));
+  file.read((char*)&top_final,sizeof(double));
+  file.read((char*)&delta_left_final,sizeof(double));
+  file.read((char*)&delta_right_final,sizeof(double));
+  file.read((char*)&delta_bottom_final,sizeof(double));
+  file.read((char*)&delta_top_final,sizeof(double));
+}
+
+double
+Pump::get_amplitude(double t){
+  return (1-t)*amplitude_init+t*amplitude_final;
+}
+
+double
+Pump::get_left(double t){
+  return (1-t)*left_init+t*left_final;
+}
+
+double
+Pump::get_right(double t){
+  return (1-t)*right_init+t*right_final;
+}
+
+double
+Pump::get_top(double t){
+  return (1-t)*top_init+t*top_final;
+}
+
+double
+Pump::get_bottom(double t){
+  return (1-t)*bottom_init+t*bottom_final;
+}
+
+double
+Pump::get_left_delta(double t){
+  return (1-t)*delta_left_init+t*delta_left_final;
+}
+
+double
+Pump::get_right_delta(double t){
+  return (1-t)*delta_right_init+t*delta_right_final;
+}
+
+double
+Pump::get_top_delta(double t){
+  return (1-t)*delta_top_init+t*delta_top_final;
+}
+
+double
+Pump::get_bottom_delta(double t){
+  return (1-t)*delta_bottom_init+t*delta_bottom_final;
+}
+
+Cloud::Cloud(){
+  amplitude_init=1.e-4;
+  left_init=0.45;
+  right_init=0.55;
+  delta_left_init=0.05;
+  delta_right_init=0.05;
+  amplitude_final=1.e-4;
+  left_final=0.45;
+  right_final=0.55;
+  delta_left_final=0.05;
+  delta_right_final=0.05;
+}
+
+
+void
+Cloud::save(fstream& file){
+  file.write((char*)&amplitude_init,sizeof(double));
+  file.write((char*)&left_init,sizeof(double));
+  file.write((char*)&right_init,sizeof(double));
+  file.write((char*)&delta_left_init,sizeof(double));
+  file.write((char*)&delta_right_init,sizeof(double));
+  file.write((char*)&amplitude_final,sizeof(double));
+  file.write((char*)&left_final,sizeof(double));
+  file.write((char*)&right_final,sizeof(double));
+  file.write((char*)&delta_left_final,sizeof(double));
+  file.write((char*)&delta_right_final,sizeof(double));
+}
+
+void
+Cloud::load(fstream& file){
+  file.read((char*)&amplitude_init,sizeof(double));
+  file.read((char*)&left_init,sizeof(double));
+  file.read((char*)&right_init,sizeof(double));
+  file.read((char*)&delta_left_init,sizeof(double));
+  file.read((char*)&delta_right_init,sizeof(double));
+  file.read((char*)&amplitude_final,sizeof(double));
+  file.read((char*)&left_final,sizeof(double));
+  file.read((char*)&right_final,sizeof(double));
+  file.read((char*)&delta_left_final,sizeof(double));
+  file.read((char*)&delta_right_final,sizeof(double));
+}
+
+double
+Cloud::get_amplitude(double t){
+  return (1-t)*amplitude_init+t*amplitude_final;
+}
+
+double
+Cloud::get_amplitude_max(){
+  return max(amplitude_init,amplitude_final);
+}
+
+double
+Cloud::get_left(double t){
+  return (1-t)*left_init+t*left_final;
+}
+
+double
+Cloud::get_right(double t){
+  return (1-t)*right_init+t*right_final;
+}
+
+double
+Cloud::get_left_delta(double t){
+  return (1-t)*delta_left_init+t*delta_left_final;
+}
+
+double
+Cloud::get_right_delta(double t){
+  return (1-t)*delta_right_init+t*delta_right_final;
+}
+
+Source::~Source(){
+  for(auto it=pumps.begin();it!=pumps.end();++it){
+    delete *it;
+  }
+  for(auto it=clouds.begin();it!=clouds.end();++it){
+    delete *it;
+  }
+}
+
+Pump*
+Source::addPump(){
+  Pump* pump=new Pump;
+  pumps.push_back(pump);
+  return pump;
+}
+
+void
+Source::removePump(Pump* pump){
+  for(auto it=pumps.begin();it!=pumps.end();++it){
+    if(*it==pump){
+      delete *it;
+      pumps.erase(it);
+      return;
+    }
+  }
+}
+
+Cloud*
+Source::addCloud(){
+  Cloud* cloud=new Cloud;
+  clouds.push_back(cloud);
+  return cloud;
+}
+
+void
+Source::removeCloud(Cloud* cloud){
+  for(auto it=clouds.begin();it!=clouds.end();++it){
+    if(*it==cloud){
+      delete *it;
+      clouds.erase(it);
+      return;
+    }
+  }
+}
+
+void
+Source::save(fstream& file){
+  size_t n=pumps.size();
+  file.write((char*)&n,sizeof(double));
+  for(auto it=pumps.begin();it!=pumps.end();++it){
+    (*it)->save(file);
+  }
+  n=clouds.size();
+  file.write((char*)&n,sizeof(double));
+  for(auto it=clouds.begin();it!=clouds.end();++it){
+    (*it)->save(file);
+  }
+}
+
+void
+Source::load(fstream& file){
+  size_t n;
+  file.read((char*)&n,sizeof(size_t));
+  for(size_t i=0;i<n;++i){
+    Pump* pump=addPump();
+    pump->load(file);
+  }
+  file.read((char*)&n,sizeof(size_t));
+  for(size_t i=0;i<n;++i){
+    Cloud* cloud=addCloud();
+    cloud->load(file);
+  }
+}

src/kernel/source.hpp

@@ -0,0 +1,68 @@
+#ifndef SOURCE_HPP
+#define SOURCE_HPP
+
+#include <list>
+#include <fstream>
+#include <iostream>
+using namespace std;
+
+class Pump{
+public:
+  double amplitude_init;
+  double left_init,right_init,bottom_init,top_init;
+  double delta_left_init,delta_right_init,delta_bottom_init,delta_top_init;
+  double amplitude_final;
+  double left_final,right_final,bottom_final,top_final;
+  double delta_left_final,delta_right_final,delta_bottom_final,delta_top_final;
+  double get_amplitude(double t);
+  double get_left(double t);
+  double get_right(double t);
+  double get_top(double t);
+  double get_bottom(double t);
+  double get_left_delta(double t);
+  double get_right_delta(double t);
+  double get_top_delta(double t);
+  double get_bottom_delta(double t);
+  Pump();
+  void save(fstream& file);
+  void load(fstream& file);
+};
+
+class Cloud{
+public:
+  double amplitude_init;
+  double left_init,right_init;
+  double delta_left_init,delta_right_init;
+  double amplitude_final;
+  double left_final,right_final;
+  double delta_left_final,delta_right_final;
+  double get_amplitude(double t);
+  double get_amplitude_max();
+  double get_left(double t);
+  double get_right(double t);
+  double get_left_delta(double t);
+  double get_right_delta(double t);
+  Cloud();
+  void save(fstream& file);
+  void load(fstream& file);
+
+};
+
+class Source{
+public:
+  list<Pump*> pumps;
+  list<Cloud*> clouds;
+  Source();
+  ~Source();
+  Pump* addPump();
+  Cloud* addCloud();
+  void removePump(Pump* pump);
+  void removeCloud(Cloud* cloud);
+  void save(fstream& file);
+  void load(fstream& file);
+};
+
+inline Source::Source(){
+}
+
+#endif

src/kernel/time.cpp

@@ -0,0 +1,18 @@
+#include "time.hpp"
+
+double Time::T=10;
+size_t Time::nT=1000;
+double Time::dT=Time::T/(Time::nT-1);
+
+void
+Time::save(fstream& file){
+  file.write((char*)&T,sizeof(double));
+  file.write((char*)&nT,sizeof(size_t));
+}
+
+void
+Time::load(fstream& file){
+  file.read((char*)&T,sizeof(double));
+  file.read((char*)&nT,sizeof(size_t));
+  dT=T/(nT-1);
+}

src/kernel/time.hpp

@@ -0,0 +1,26 @@
+#ifndef TIME_HPP
+#define TIME_HPP
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+
+class Time{
+public:
+  //! Total duration of the simulation
+  static double T;
+
+  //! Number of time steps
+  static size_t nT;
+
+  //! Time step
+  static double dT;
+
+  //! Save time in file
+  static void save(fstream& file);
+
+  //! Load time from file
+  static void load(fstream& file);
+};
+#endif

src/qt/output/solver.cpp

@@ -0,0 +1,19 @@
+#include "qt/output/solver.hpp"
+
+QtSolver::QtSolver(string filename):Kernel::Solver(filename){
+}
+
+void
+QtSolver::run(){
+  QTime start=QTime::currentTime();
+  cout<<"[QtSolver::run] start"<<endl;
+  //  for(size_t i=1;i<Time::nT;++i){
+  for(size_t i=1;i<Time::nT;++i){
+    compute_next_solution();
+  }
+  cout<<"[QtSolver::run] end"<<endl;
+  QTime end=QTime::currentTime();
+  float ms=start.msecsTo(end);
+  cout<<"[QtSolver::run] duration : "<<ms/1000<<endl;
+
+}

src/qt/output/solver.hpp

@@ -0,0 +1,24 @@
+#ifndef QT_KERNEL_HPP
+#define QT_KERNEL_HPP
+
+#include <QThread>
+#include <QTime>
+#include "../../kernel/solver.hpp"
+
+class QtSolver:public QObject,public Kernel::Solver{
+  Q_OBJECT
+public:
+  QtSolver(string filename);
+  ~QtSolver();
+public slots:
+  void run();
+signals:
+  void finished(int);
+};
+
+inline
+QtSolver::~QtSolver(){
+}
+
+
+#endif