Add basin.

Makefile

@@ -2,17 +2,18 @@ QT_FLAGS	=-fPIC `pkg-config --cflags --libs Qt5OpenGL`
 QT_LIB		=`pkg-config --libs Qt5OpenGL` -lGL -lGLU
 QT_MOC		= moc
 GPP		= g++
+
 FLAGS 		= -g -O3 -Isrc -Wfatal-errors
 
 QT_INPUT_FILES	= main_window data input physics time geometry initial_state view tank pump pump_tab cloud clouds_tab
 QT_INPUT_MOC_FILES	= main_window physics input view geometry initial_state tank pump_tab pump cloud clouds_tab
 MATH_INPUT_FILES	= poly3 algo spline
-KERNEL_INPUT_FILES = physics time geometry initial_state source input_data
+KERNEL_INPUT_FILES = physics time geometry initial_state source input_data basin
 
 QT_OUTPUT_FILES	= main_window solver output_view output
 QT_OUTPUT_MOC_FILES	= main_window solver output
 MATH_OUTPUT_FILES	= algo
-KERNEL_OUTPUT_FILES = physics time geometry initial_state source input_data solution solver piccard horizontal_problem overland all_vertical_richards richards_evolutive_time richards
+KERNEL_OUTPUT_FILES = physics time geometry initial_state source input_data solution solver piccard horizontal_problem overland all_vertical_richards richards_evolutive_time richards basin
 
 
 QT_INPUT_OBJS	= $(addprefix obj/qt/input/,$(addsuffix .o,$(QT_INPUT_FILES)))

src/debug.hpp

@@ -9,8 +9,8 @@ using namespace std;
 
 class Debug{
 public:
-  static const int level=2;
-  static const size_t ix=36;
+  static const int level=0;
+  static const size_t ix=44;
   static void display(string name,double* u,size_t start,size_t end);
   static void pause();
   static bool debug_Thomas;

src/kernel/basin.cpp

@@ -0,0 +1,219 @@
+#include "basin.hpp"
+
+void
+Basin::compute_leak_direction(){
+  switch(position){
+    case Middle:
+    if(abs(hleft-hright)<delta_both_leak) leak_direction=Both;
+    else leak_direction=(hleft>hright)?Right:Left;
+    break;
+    case BorderLeft:
+    leak_direction=Right;
+    break;
+    case BorderRight:
+    leak_direction=Left;
+    break;
+    default:
+    leak_direction=None;
+    break;
+  };
+}
+
+void Basin::compute_maxs(double* hsoil){
+  switch(leak_direction){
+    case Left:
+    hmax=hleft;
+    xmax_left=xleft;
+    {
+      int x;
+      for(x=xleft+1;x<xright;++x){
+        if(hsoil[x]+Physics::minimal_height_to_runoff>hmax){
+          xmax_right=x;
+          break;
+        }
+      }
+      xmax_right=min(x,(int)xright);
+    }
+    break;
+    case Right:
+    hmax=hright;
+    xmax_right=xright;
+    {
+      int x;
+      for(x=xright-1;x>xleft;--x){
+        if(hsoil[x]+Physics::minimal_height_to_runoff>hmax){
+          xmax_left=x;
+          break;
+        }
+      }
+      xmax_left=max(x,(int)xleft);
+    }
+    break;
+    case Both:
+    hmax=min(hleft,hright);
+    xmax_left=xleft;
+    xmax_right=xright;
+    break;
+    case None:
+    default:
+    hmax=numeric_limits<double>::infinity();
+    xmax_left=xleft;
+    xmax_right=xright;
+    break;
+  };
+}
+
+Basin::Basin(Basin* left_,Basin* right_){
+  left=left_;
+  right=right_;
+  xleft=left->xleft;
+  xright=right->xright;
+  hleft=left->hleft;
+  hright=right->hright;
+  hmin=left->hmax;
+  xacc=left->xacc;
+  if(left->position==BorderLeft){
+    if(right->position==BorderRight) position=BorderBoth;
+    else position=BorderLeft;
+  }
+  else if(right->position==BorderRight) position=BorderRight;
+  else position=Middle;
+  compute_leak_direction();
+
+}
+
+void
+Basin::display(string indent){
+  cout<<indent<<'['<<xleft<<','<<xright<<"] : "<<hmax<<" : ["<<xmax_left<<","<<xmax_right<<"] and "<<hmin<<endl;
+  if(left==nullptr) return;
+  left->display(indent+' ');
+  right->display(indent+' ');
+}
+
+void
+Basin::compute_vflow(double* hsoil,double* hov){
+  if(is_primitive()){
+    vflow=max(0.0,hov[xacc]-hsoil[xacc]-Physics::minimal_height_to_runoff);
+  }
+  else{
+    left->compute_vflow(hsoil,hov);
+    right->compute_vflow(hsoil,hov);
+    vflow=left->vflow+right->vflow;
+  }
+}
+
+void
+Basin::repartition(double* hsoil,double* hov,Direction* rundir){
+
+  overflow=false;
+  vmin=0;
+  if(not is_primitive()){
+    for(size_t ix=xleft;ix<=xright;++ix){
+      vmin+=max(0.0,hmin-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
+    }
+  }
+  if(vflow!=0 and vflow>=vmin){
+    vmax=0;
+    for(size_t ix=xmax_left;ix<=xmax_right;++ix){
+      vmax+=max(0.0,hmax-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
+    }
+    if(vflow-vmax>1e-10){
+      overflow=true;
+      vleak=vflow-vmax;
+      if(not is_primitive()){
+        double delta=hov[right->xacc]-hsoil[right->xacc]-Physics::minimal_height_to_runoff;
+        hov[left->xacc]+=delta;
+        hov[right->xacc]-=delta;
+      }
+      hov[xacc]-=vleak;
+      switch(leak_direction){
+        case Left:
+        hov[xleft-1]+=vleak;
+        runoff(xleft-1,Left,rundir,hsoil,hov);
+        break;
+        case Right:
+        hov[xright+1]+=vleak;
+        runoff(xright+1,Right,rundir,hsoil,hov);
+        break;
+        case Both:
+        hov[xleft-1]+=vleak/2;
+        runoff(xleft-1,Left,rundir,hsoil,hov);
+        hov[xright+1]+=vleak/2;
+        runoff(xright+1,Right,rundir,hsoil,hov);
+        break;
+        default:
+        cerr<<"[Bug] Impossible case ! ... in theory :)"<<endl;
+        exit(-1);
+        break;
+      }
+    }
+  }
+  else{
+    if(not is_primitive()){
+      left->repartition(hsoil,hov,rundir);
+      right->repartition(hsoil,hov,rundir);
+      overflow=left->overflow|right->overflow;
+    }
+  }
+}
+
+void
+Basin::runoff(size_t xstart,Direction dir,Direction* rundir,double* hsoil,double* hov){
+  double r=max(0.0,hov[xstart]-hsoil[xstart]-Physics::minimal_height_to_runoff);
+  int x=xstart;
+  int xstep=(dir==Left)?-1:1;
+  while(r>0 and rundir[x]!=None){
+    hov[x]-=r;
+    x+=xstep;
+    hov[x]+=r;
+    r=max(0.0,hov[x]-hsoil[x]-Physics::minimal_height_to_runoff);
+  }
+
+}
+
+double
+Basin::f(double r,double* hsoil,double* hov,double vext){
+  double res=-vext;
+  for(size_t ix=xleft;ix<=xright;++ix){
+    res+=max(r-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff),0.0);
+  }
+  return res;
+}
+
+double
+Basin::df(double r,double* hsoil,double* hov,double vext){
+  double res=0;
+  for(size_t ix=xleft;ix<=xright;++ix){
+    if(r>min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff)) ++res;
+  }
+  return res;
+}
+
+void
+Basin::flatten(double* hsoil,double* hov){
+  double crit=1;
+  size_t count=0;
+  if(vflow>1e-9){
+    double r=hsoil[xacc]+1;//Physics::minimal_height_to_runoff;
+    while(crit>1e-10 and count<100){
+      ++count;
+      double rold=r;
+      r-=f(r,hsoil,hov,vflow)/df(r,hsoil,hov,vflow);
+      crit=abs(r-rold)+abs(f(r,hsoil,hov,vflow));
+    }
+    if(crit>1e-10){
+      cerr<<" Overland coverge + [TODO] on fait quoi ?"<<endl;
+    }
+    if(r>=hmin or is_primitive()){
+      for(size_t ix=xleft;ix<=xright;++ix){
+        hov[ix]=max(r,min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
+      }
+
+
+    }
+    else{
+      left->flatten(hsoil,hov);
+      right->flatten(hsoil,hov);
+    }
+  }
+}

src/kernel/basin.hpp

@@ -2,100 +2,73 @@
 #define BASIN_HPP
 
 #include <string>
+#include "physics.hpp"
+#include "debug.hpp"
 
 using namespace std;
 
-
-  enum Direction{Left,Right,Both};
-  enum BasinPosition{Middle,BorderLeft,BorderRight};
-
-  struct Summit{
-    size_t ix;
-    double h;
-    bool operator<(const Summit& S) const;
-    Summit(size_t ix,double h);
-  };
-
-
-  class Basin{
-  public:
-    static constexpr double delta_both_leak=0.01;
-    size_t xleft,xright;
-    size_t xbottom; //Used only for primitives bassin
-    Direction leak_direction;
-    BasinPosition position;
-    //Basin* neighbours[2];
-    double hmin,hleft,hright;
-    Basin* left;
-    Basin* right;
-    void compute_leak_direction();
-  public:
-    Basin();
-    Basin(Basin* left,Basin* right);
-    void display(string indent="");
-    bool is_primitive();
-  };
-
-
-
-  inline
-  Summit::Summit(size_t ix_,double h_):ix(ix_),h(h_){
-  }
-
-  inline bool
-  Summit::operator<(const Summit& S) const{
-    return h<S.h;
-  }
-
-  inline Basin::Basin(){
-    left=nullptr;
-    right=nullptr;
-    position=Middle;
-  }
-
-  inline void
-  Basin::compute_leak_direction(){
-    switch(position){
-      case Middle:
-         if(abs(hleft-hright)<delta_both_leak) leak_direction=Both;
-         else leak_direction=(hleft>hright)?Right:Left;
-         break;
-      case BorderLeft:
-         leak_direction=Right;
-         break;
-      case BorderRight:
-         leak_direction=Left;
-         break;
-     default:
-         break;
-    };
-  }
-
-  inline Basin::Basin(Basin* left_,Basin* right_){
-    left=left_;
-    right=right_;
-    xleft=left->xleft;
-    xright=right->xright;
-    hleft=left->hleft;
-    hright=right->hright;
-    hmin=min(left->hmin,right->hmin);
-    if(left->position==BorderLeft) position=BorderLeft;
-    else if(right->position==BorderRight) position=BorderRight;
-    else position=Middle;
-    compute_leak_direction();
-  }
-
-  inline void Basin::display(string indent){
-    cout<<indent<<'['<<xleft<<','<<xright<<']'<<endl;
-    if(left==nullptr) return;
-    left->display(indent+' ');
-    right->display(indent+' ');
-  }
-
-  inline
-  bool Basin::is_primitive(){
-    return left=nullptr;
-  }
+enum Direction{None,Left,Right,Both};
+enum BasinPosition{Middle,BorderLeft,BorderRight,BorderBoth};
+
+struct Summit{
+  size_t ix;
+  double h;
+  bool operator<(const Summit& S) const;
+  Summit(size_t ix,double h);
+};
+
+
+class Basin{
+public:
+  static constexpr double delta_both_leak=0.01;
+  bool overflow;
+  size_t xleft,xright;
+  size_t xmax_left,xmax_right;
+  size_t xacc;
+  Direction leak_direction;
+  BasinPosition position;
+  double hmin,hleft,hright,hmax;
+  Basin* left;
+  Basin* right;
+  double vmin;
+  double vflow;
+  double vmax;
+  double vleak;
+  void compute_leak_direction();
+  void compute_maxs(double* hsoil);
+  void compute_vflow(double* hsoil,double* hov);
+  void repartition(double* hsoil,double* hov,Direction* rundir);
+  void runoff(size_t xstart,Direction dir,Direction* rundir,double* hsoil,double* hov);
+
+  double f(double r,double* hsoil,double* hov,double vext);
+  double df(double r,double* hsoil,double* hov,double vext);
+public:
+  Basin();
+  Basin(Basin* left,Basin* right);
+  void display(string indent="");
+  bool is_primitive();
+    void flatten(double* hsoil,double* hov);
+};
+
+inline
+Summit::Summit(size_t ix_,double h_):ix(ix_),h(h_){
+}
+
+inline bool
+Summit::operator<(const Summit& S) const{
+  return h<S.h;
+}
+
+inline Basin::Basin(){
+  left=nullptr;
+  right=nullptr;
+  position=Middle;
+}
+
+inline bool
+Basin::is_primitive(){
+  return left==nullptr;
+}
 
 
 #endif

src/kernel/geometry.cpp

@@ -1,4 +1,5 @@
 #include "geometry.hpp"
+#include "physics.hpp"
 
 double inf = numeric_limits<double>::infinity();
 
@@ -101,7 +102,6 @@ Geometry::load(fstream& file,bool init) {
   }
   dhsoil[nX-1]=(hsoil[nX-1]-hsoil[nX-2])/dX;
   dhbot[nX-1]=(hbot[nX-1]-hbot[nX-2])/dX;
-  cout<<"Done"<<endl;
 
 }
 
@@ -112,11 +112,9 @@ Geometry::find_basins(const Summit& s,list<Basin*>& basins){
   for(list<Basin*>::iterator it=basins.begin();it!=basins.end();++it){
     Basin* b=*it;
     if(b->xleft==s.ix){
-      cout<<"Right found"<<endl;
       it_right=it;
     }
     if(b->xright==s.ix){
-      cout<<"Left found"<<endl;
       it_left=it;
     }
   }
@@ -142,7 +140,7 @@ void Geometry::compute_basins(){
     Basin* b=new Basin;
     size_t xb=v_min[ib];
     //Find left bound
-    b->xbottom=xb;
+    b->xacc=xb;
     if(xb==0){
       b->xleft=xb;
       b->position=BorderLeft;
@@ -161,37 +159,57 @@ void Geometry::compute_basins(){
     }
     else{
       size_t ix=xb;
-      while(ix<nX and hsoil[ix+1]>hsoil[ix]) ++ix;
+      while(ix<nX-1 and hsoil[ix+1]>hsoil[ix]) ++ix;
       if(ix==nX-1)b->position=BorderRight;
       b->xright=ix;
     }
-    b->hleft=hsoil[b->xleft];
-    b->hright=hsoil[b->xright];
+    b->hleft=hsoil[b->xleft]+Physics::minimal_height_to_runoff;
+    b->hright=hsoil[b->xright]+Physics::minimal_height_to_runoff;
 
     b->compute_leak_direction();
+    b->compute_maxs(hsoil);
 
     // Set hmin
-    b->hmin=hsoil[xb];
+    b->hmin=hsoil[xb]+Physics::minimal_height_to_runoff;
 
     basins_to_merge.push_back(b);
   }
-  //Determine meta basins`
-  sort(summits.begin(),summits.end());
-  for(auto it=summits.begin();it!=summits.end();++it){
-    Summit& s=*it;
-    auto res=find_basins(s,basins_to_merge);
-    list<Basin*>::iterator it_left=res.first;
-    list<Basin*>::iterator it_right=res.second;
-    Basin* b_left=*it_left;
-    Basin* b_right=*it_right;
-
-    basins_to_merge.erase(it_left);
-    basins_to_merge.erase(it_right);
-    Basin* b=new Basin(b_left,b_right);
-    basins_to_merge.push_back(b);
+  //Determine meta basins
+
+  if(summits.empty()){
+    root_basin=new Basin;
+    root_basin->position=BorderBoth;
+    root_basin->left=nullptr;
+    root_basin->right=nullptr;
+    root_basin->xacc=0;
+    root_basin->leak_direction=None;
+    root_basin->xleft=0;
+    root_basin->xright=nX-1;
+    double h=hsoil[0];
+    root_basin->hleft=h;
+    root_basin->hright=h;
+    root_basin->hmin=h;
+
+  }
+  else{
+    sort(summits.begin(),summits.end());
+    for(auto it=summits.begin();it!=summits.end();++it){
+      Summit& s=*it;
+      auto res=find_basins(s,basins_to_merge);
+      list<Basin*>::iterator it_left=res.first;
+      list<Basin*>::iterator it_right=res.second;
+      Basin* b_left=*it_left;
+      Basin* b_right=*it_right;
+
+      basins_to_merge.erase(it_left);
+      basins_to_merge.erase(it_right);
+      Basin* b=new Basin(b_left,b_right);
+      b->compute_maxs(hsoil);
+      basins_to_merge.push_back(b);
+    }
+    root_basin=basins_to_merge.front();
+    root_basin->display();
   }
-  root_basin=basins_to_merge.front();
-  root_basin->display();
 
   // Compute Runoff directions
   runoff_directions=new Direction[nX];
@@ -199,6 +217,16 @@ void Geometry::compute_basins(){
     double dh=(hsoil[ix+1]-hsoil[ix-1])/(2*dX);
     if(abs(dh)<max_slope_both_side_runoff) runoff_directions[ix]=Both;
     else runoff_directions[ix]=(dh>0)?Left:Right;
-    cout<<ix<<" "<<dh<<" "<<runoff_directions[ix]<<endl;
   }
+  for(size_t ib=0;ib<nb;++ib){
+    size_t xb=v_min[ib];
+    runoff_directions[xb]=None;
+  }
+  runoff_directions[0]=Right;
+  runoff_directions[nX-1]=Left;
+  cout<<"Runoff "<<endl;
+  for(size_t ix=1;ix<nX-1;++ix){
+    cout<<ix<<" : "<<runoff_directions[ix]<<endl;
+  }
+  Debug::pause();
 }

src/kernel/geometry.hpp

@@ -21,7 +21,7 @@ extern double inf;
 //! The Geometry class contains all geometric parameters of the domain.
 class Geometry{
 public:
-  static constexpr double max_slope_both_side_runoff=0.01;
+  static constexpr double max_slope_both_side_runoff=0.00;
   //! Horizontal length of the domain
   double lX;
 

src/kernel/overland.cpp

@@ -10,8 +10,11 @@ namespace Kernel{
     compute_flux_soil();
     total_error=0;
     for(size_t ix=0;ix<geometry->nX;++ix){
-      //TODO :: add rain
       hov[ix]=init_hov[ix]-dt*flux_soil[ix];
+      //TODO :: add rain
+      if(ix>89){
+        hov[ix]+=dt*0.0001;
+      }
     }
     apply_flow();
     for(size_t ix=0;ix<geometry->nX;++ix){
@@ -27,7 +30,6 @@ namespace Kernel{
     #ifdef MYDEBUG
     cout<<"   [Overland::run] end"<<endl;
     #endif
-
   }
 
   void Overland::compute_flux_soil(){
@@ -40,13 +42,6 @@ namespace Kernel{
       else if(ix==geometry->nX-1) flux-=Physics::k0*Physics::kr(Pl[ix]+Physics::rho*Physics::g*(l[ix]-geometry->hsoil[ix]))*(hydr[ix]-hydr[ix-1])*geometry->dhsoil[ix]/geometry->dX;
       else flux-=Physics::k0*Physics::kr(Pl[ix]+Physics::rho*Physics::g*(l[ix]-geometry->hsoil[ix]))*(hydr[ix+1]-hydr[ix-1])/2*geometry->dhsoil[ix]/geometry->dX;
       flux_soil[ix]=flux;
-      #ifdef MYDEBUG
-      if(Debug::ix==ix){
-        cout<<"      [Overland::compute_flux_soil] p1 = "<<p1<<endl;
-        cout<<"      [Overland::compute_flux_soil] p2 = "<<p2<<endl;
-        cout<<"      [Overland::compute_flux_soil] flux = "<<flux<<endl;
-      }
-      #endif
 
     }
   }
@@ -56,11 +51,100 @@ namespace Kernel{
     while(overflow){
       overflow=false;
       runoff_to_local_min();
-
+      basin_transfer();
     }
   }
 
   void Overland::runoff_to_local_min(){
+  //  cout<<"      [Overland::runoff_to_local_min]"<<endl;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double dh=hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff;
+      if(dh>0){
+        double dh_left,dh_right;
+        switch(geometry->runoff_directions[ix]){
+          case Left:
+          dh_left=dh;
+          dh_right=0;
+          break;
+          case Right:
+          dh_right=dh;
+          dh_left=0;
+          break;
+          case Both:
+          dh_left=dh/2;
+          dh_right=dh/2;
+          break;
+          default:
+          dh_left=0;
+          dh_right=0;
+          break;
+        }
+        size_t jx=ix;
+        while(dh_left>0){
+          --jx;
+          if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Right){
+            hov[jx]+=dh_left;
+            hov[ix]-=dh_left;
+            break;
+          }
+          else{
+            double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_left);
+            hov[jx]+=temp;
+            hov[ix]-=temp;
+            dh_left-=temp;
+          }
+        }
+        jx=ix;
+        while(dh_right>0){
+          ++jx;
+          if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Left){
+            hov[jx]+=dh_right;
+            hov[ix]-=dh_right;
+            break;
+          }
+          else{
+            double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_right);
+            hov[jx]+=temp;
+            hov[ix]-=temp;
+            dh_right-=temp;
+          }
+        }
 
+      }
+    }
+  /*cout<<"------------------- after runoff ----------------"<<endl;
+    for(size_t ix=0;ix<100;++ix){
+      cout<<ix<<" : "<<max(0.0,hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff)<<endl;;
+    }*/
+  }
+
+  void Overland::basin_transfer(){
+    double vol_tot=0;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      vol_tot+=max(hov[ix]-geometry->hsoil[ix],0.0);
+    }
+    //TODO :: Remove guard
+    int nmax=100;
+    do{
+      --nmax;
+      geometry->root_basin->compute_vflow(geometry->hsoil,hov);
+      geometry->root_basin->repartition(geometry->hsoil,hov,geometry->runoff_directions);
+    }while(geometry->root_basin->overflow and nmax>0);
+    if(nmax==0){
+      cerr<<"Too many loop"<<endl;
+      Debug::pause();//exit(-1);
+    }
+    geometry->root_basin->compute_vflow(geometry->hsoil,hov);
+    double vol_tot_mid=0;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      vol_tot_mid+=max(hov[ix]-geometry->hsoil[ix],0.0);
+    }
+    geometry->root_basin->flatten(geometry->hsoil,hov);
+
+    double vol_tot_new=0;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      vol_tot_new+=max(hov[ix]-geometry->hsoil[ix],0.0);
+    }
+    cout<<vol_tot-vol_tot_new<<endl;
   }
 }

src/kernel/overland.hpp

@@ -12,6 +12,8 @@ namespace Kernel{
     void compute_flux_soil();
     void apply_flow();
     void runoff_to_local_min();
+    void basin_transfer();
+
   public:
     //Input
     double dt;
@@ -22,6 +24,7 @@ namespace Kernel{
     double* hydr;
     double n1_init_hov;
     double* previous_hov;
+
     //Output
     double* Psoil; //Own by Overland object
     double* hov;

src/kernel/physics.hpp

@@ -8,10 +8,13 @@
 
 using namespace std;
 
+
 //! The Physics class contains all physical parameters characterising the soil.
 //! This class contains only static members.
 class Physics{
 public:
+  static constexpr double minimal_height_to_runoff=0.019;
+
   enum Model{BrooksCorey};
 
   //! Set physics model

src/kernel/piccard.cpp

@@ -84,12 +84,12 @@ namespace Kernel{
     new_solution->P=all_vertical_richards.P;
     size_t k=0;
     double previous_error=numeric_limits<double>::infinity();
-    //Debug::pause();
+  //  Debug::pause();
+
     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;
-
       #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard::run] k = "<<k<<endl;
       #endif
@@ -149,12 +149,12 @@ namespace Kernel{
         has_converged=false;
         return;
       }
-      //Debug::pause();
 
     }
     if(error>=tolerence_Piccard){
       #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard]::run not converge"<<endl;
+      exit(0);
       #endif
       has_converged=false;
       //Voir nettoyage memoire

src/kernel/richards.cpp

@@ -94,7 +94,6 @@ namespace Kernel{
       //Computed P is in out_P=temp_P[count%2]
       /*if(Debug::ix==ix) Debug::display("in_P",in_P,0,nZ);
       if(Debug::ix==ix) Debug::display("out_P",out_P,0,nZ);*/
-
       error=error2(in_P,out_P,nZ)/norm_previous_P;
       if(count==1){
         in_P=temp_P[1];

src/kernel/solver.cpp

@@ -101,7 +101,7 @@ namespace Kernel{
 
   bool
   Solver::compute_next_solution(){
-    //cout<<"[solver::next]"<<endl;
+
     compute_sources(number_computed_solutions*Time::dT);
     if(m1>1) m1=m1/2;
     Solution* s=space_solution();
@@ -109,6 +109,7 @@ namespace Kernel{
     solutions[number_computed_solutions++]=s;
     //cout<<"[next] done"<<endl;
     return true;
+
   }
 
   Solution*
@@ -173,7 +174,7 @@ namespace Kernel{
           pumps[ix][iz]+=p;
         }
       }
-      }
+    }
   }
 
 }

src/qt/input/data.cpp

@@ -129,7 +129,7 @@ void
 QtInputData::updateInitialState(){
   for(size_t i=0;i<geometry.nX;++i){
     double x=double(i)/(geometry.nX-1);
-    initial_state->hsat[i]=spline_sat(x)*factor;
+    initial_state->hsat[i]=min(geometry.hsoil[i],spline_sat(x)*factor);
   }
   initial_state->updatePressure();
   updateOverland();

src/qt/output/output_view.cpp

@@ -275,6 +275,7 @@ QtOutputView::drawOverland(){
 void
 QtOutputView::displayInfos(double x,double z){
   cout<<"---------------------------------"<<endl;
+  cout<<"step = "<<step<<endl;
   cout<<"x = "<<x<<endl;
   cout<<"z = "<<z<<endl;
   size_t ix=floor(x/dX+0.5);