Work on flow, correct overland problem artefact comming from bad derivatives.

src/kernel/basin.hpp

@@ -0,0 +1,101 @@
+#ifndef BASIN_HPP
+#define BASIN_HPP
+
+#include <string>
+
+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;
+  }
+
+
+#endif

src/kernel/geometry.cpp

@@ -91,4 +91,114 @@ Geometry::load(fstream& file,bool init) {
     if(init) Z[i]=new double[nZ[i]];
     file.read((char*)Z[i],nZ[i]*sizeof(double));
   }
+  //Comoute dhsoil and dhbot
+  dhsoil[0]=(hsoil[1]-hsoil[0])/dX;
+  dhbot[0]=(hbot[1]-hbot[0])/dX;
+
+  for(size_t i=1;i<nX-1;++i){
+    dhsoil[i]=(hsoil[i+1]-hsoil[i-1])/(2*dX);
+    dhbot[i]=(hbot[i+1]-hbot[i-1])/(2*dX);
+  }
+  dhsoil[nX-1]=(hsoil[nX-1]-hsoil[nX-2])/dX;
+  dhbot[nX-1]=(hbot[nX-1]-hbot[nX-2])/dX;
+  cout<<"Done"<<endl;
+
+}
+
+pair<list<Basin*>::iterator,list<Basin*>::iterator>
+Geometry::find_basins(const Summit& s,list<Basin*>& basins){
+  list<Basin*>::iterator it_left;
+  list<Basin*>::iterator it_right;
+  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;
+    }
+  }
+  return pair<list<Basin*>::iterator,list<Basin*>::iterator>(it_left,it_right);
+}
+
+void Geometry::compute_basins(){
+  //Find local min point of hsoil
+  vector<size_t> v_min;
+  if(hsoil[0]<hsoil[1]) v_min.push_back(0);
+  for(size_t ix=1;ix<nX-1;++ix){
+    double h=hsoil[ix];
+    if(hsoil[ix-1]>h and h<hsoil[ix+1]) v_min.push_back(ix);
+  }
+  if(hsoil[nX-2]>hsoil[nX-1]) v_min.push_back(nX-1);
+  size_t nb=v_min.size();
+
+  //Primitive Basins
+  vector<Summit> summits;
+  list<Basin*> basins_to_merge;
+
+  for(size_t ib=0;ib<nb;++ib){
+    Basin* b=new Basin;
+    size_t xb=v_min[ib];
+    //Find left bound
+    b->xbottom=xb;
+    if(xb==0){
+      b->xleft=xb;
+      b->position=BorderLeft;
+    }
+    else{
+      size_t ix=xb;
+      while(ix>0 and hsoil[ix-1]>hsoil[ix]) --ix;
+      b->xleft=ix;
+      if(ix==0) b->position=BorderLeft;
+      if(ix>0) summits.emplace_back(ix,hsoil[ix]);
+    }
+    //Find right bound
+    if(xb==nX-1){
+      b->xright=xb;
+      b->position=BorderRight;
+    }
+    else{
+      size_t ix=xb;
+      while(ix<nX 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->compute_leak_direction();
+
+    // Set hmin
+    b->hmin=hsoil[xb];
+
+    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);
+  }
+  root_basin=basins_to_merge.front();
+  root_basin->display();
+
+  // Compute Runoff directions
+  runoff_directions=new Direction[nX];
+  for(size_t ix=1;ix<nX-1;++ix){
+    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;
+  }
 }

src/kernel/geometry.hpp

@@ -4,18 +4,24 @@
 #include <iostream>
 #include <fstream>
 #include <limits>
+#include <list>
+#include <vector>
+#include <algorithm>
 
 #include "math/spline.hpp"
-
+#include "basin.hpp"
 using namespace std;
 
 using Func = double (*)(double);
 
 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;
   //! Horizontal length of the domain
   double lX;
 
@@ -48,6 +54,9 @@ public:
   //! 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;
 
+  Direction* runoff_directions;
+
+  Basin* root_basin;
   //!
   Geometry();
 
@@ -58,6 +67,10 @@ public:
   void save(fstream& file);
 
   void load(fstream& file,bool init);
+
+  pair<list<Basin*>::iterator,list<Basin*>::iterator> find_basins(const Summit& s,list<Basin*>& basins);
+
+  void compute_basins();
 };
 
 

src/kernel/overland.cpp

@@ -12,9 +12,13 @@ namespace Kernel{
     for(size_t ix=0;ix<geometry->nX;++ix){
       //TODO :: add rain
       hov[ix]=init_hov[ix]-dt*flux_soil[ix];
+    }
+    apply_flow();
+    for(size_t ix=0;ix<geometry->nX;++ix){
       double e=abs(hov[ix]-previous_hov[ix])/n1_init_hov;
       error_x[ix]+=e;
       total_error+=e;
+
       Psoil[ix]=Physics::Psoil(geometry->hsoil[ix],hov[ix]);
     }
     #ifdef MYDEBUG
@@ -23,6 +27,7 @@ namespace Kernel{
     #ifdef MYDEBUG
     cout<<"   [Overland::run] end"<<endl;
     #endif
+
   }
 
   void Overland::compute_flux_soil(){
@@ -33,7 +38,7 @@ namespace Kernel{
       double flux=Physics::k0/2*(Physics::kr(p1)+Physics::kr(p2))*((p1-p2)/(geometry->dZ[ix]*Physics::rho*Physics::g)+1);
       if(ix==0) flux-=Physics::k0*Physics::kr(Pl[0]+Physics::rho*Physics::g*(l[0]-geometry->hsoil[0]))*(hydr[1]-hydr[0])*geometry->dhsoil[0]/geometry->dX;
       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;
+      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){
@@ -45,4 +50,17 @@ namespace Kernel{
 
     }
   }
+
+  void Overland::apply_flow(){
+    bool overflow=true;
+    while(overflow){
+      overflow=false;
+      runoff_to_local_min();
+
+    }
+  }
+
+  void Overland::runoff_to_local_min(){
+
+  }
 }

src/kernel/overland.hpp

@@ -10,6 +10,8 @@ namespace Kernel{
     const Geometry* geometry;
     double* flux_soil;
     void compute_flux_soil();
+    void apply_flow();
+    void runoff_to_local_min();
   public:
     //Input
     double dt;

src/kernel/solver.cpp

@@ -6,6 +6,7 @@ namespace Kernel{
     fstream file;
     file.open(filename.c_str(),fstream::in|fstream::binary);
     input_data.load(file,true);
+    input_data.geometry.compute_basins();
     file.close();
 
     //HACK for pump
@@ -92,6 +93,10 @@ namespace Kernel{
     input_data.remove_initial_state();
 
     number_computed_solutions=1;
+
+
+
+
   }
 
   bool
@@ -170,4 +175,5 @@ namespace Kernel{
       }
       }
   }
+
 }

src/kernel/solver.hpp

@@ -10,7 +10,7 @@
 #include "input_data.hpp"
 #include "solution.hpp"
 #include "piccard.hpp"
-
+#include "basin.hpp"
 using namespace std;
 
 namespace Kernel{
@@ -61,6 +61,7 @@ namespace Kernel{
 
     void compute_sources(double t);
 
+
   public:
     //! Construct a Solver from the input file input_filename
     Solver(string input_filename);

src/qt/output/main_window.cpp

@@ -1,7 +1,7 @@
 #include "qt/output/main_window.hpp"
 
 QtMainWindow::QtMainWindow():QMainWindow(){
-  string filename="inputs/test.input";
+  string filename="inputs/ov1.input";
   output=new QtOutput(filename);
   setCentralWidget(output);
 }

todo.txt

@@ -17,3 +17,6 @@ Voir problem locale . ,
 
 [Overland]
 Add rain
+
+[Save and load]
+Remove dhbot and dhsoil