Add pumps + gestion of indice_x_richards

src/kernel/all_vertical_richards.cpp

@@ -2,26 +2,36 @@
 
 namespace Kernel{
   void
-  AllVerticalRichards::init(const Geometry* geometry_){
+  AllVerticalRichards::init(const Geometry* geometry_,double** pumps_){
     geometry=geometry_;
+    pumps=pumps_;
     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);
+      richards_evolutive_time[ix].init(ix,geometry,pumps[ix]);
     }
   }
 
+  void
+  AllVerticalRichards::update_indice_x_Richards(){
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      indice_x_Richards[ix]=(error_x[ix]>n1_previous_hydr*max_error_x);
+    }
+  }
+
+
   void
   AllVerticalRichards::init_indice_x_Richards(){
-    for(size_t x=0;x<geometry->nX;++x){
-      indice_x_Richards[x]=true;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      indice_x_Richards[ix]=true;
     }
   }
 
+
   void
   AllVerticalRichards::run(){
     Debug::debug_Thomas=false;
-    
+
     size_t i_left,i_middle,i_right;
 
     #ifdef MYDEBUG

src/kernel/all_vertical_richards.hpp

@@ -7,13 +7,18 @@
 namespace Kernel{
   class AllVerticalRichards{
   private:
+    static constexpr double max_error_x=1e-8;
     const Geometry* geometry;
+    double** pumps;
     bool* indice_x_Richards;
     RichardsEvolutiveTime* richards_evolutive_time;
 
   public:
     //Input
     double dt;
+    double* error_x;
+    double n1_previous_hydr;
+
     double** init_P; //P_0
     double** previous_P; //P_{k-1}
     double* hydr;
@@ -26,8 +31,9 @@ namespace Kernel{
     bool has_converged;
 
     AllVerticalRichards();
-    void init(const Geometry* geometry);
+    void init(const Geometry* geometry,double** pumps);
     void init_indice_x_Richards();
+    void update_indice_x_Richards();
     void compute_hsat();
     void run();
   };

src/kernel/horizontal_problem.cpp

@@ -12,6 +12,16 @@ namespace Kernel{
     F=new double[geometry->nX];
   }
 
+  void
+  HorizontalProblem::compute_error(){
+    total_error=0;
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double e=error_x[ix]=abs(previous_hydr[ix]-hydr[ix]);
+      total_error+=e;
+    }
+  }
+
+
   void
   HorizontalProblem::run(){
     #ifdef MYDEBUG
@@ -21,7 +31,7 @@ namespace Kernel{
       hydr[ix]=previous_hydr[ix];
     }*/
     solve_system();
-    error=error2(previous_hydr,hydr,geometry->nX);
+    compute_error();
     #ifdef MYDEBUG
     cout<<"   [HorizontalProblem::run] end"<<endl;
     #endif

src/kernel/horizontal_problem.hpp

@@ -17,7 +17,7 @@ namespace Kernel{
     double* sub_M;
     double* F;
     void solve_system();
-    double compute_error();
+    void compute_error();
   public:
     //Input
     double* previous_hydr;
@@ -25,20 +25,23 @@ namespace Kernel{
     double* Pl;
     //Output
     double* hydr;
-    double error;
+    double total_error;
+    double* error_x;
+
 
     HorizontalProblem();
 
     void init(const Geometry* geometry);
 
     void run();
+
+    double get_error(size_t ix);
   };
 
   inline HorizontalProblem::HorizontalProblem(){
   }
 
 
-
 }
 
 #endif

src/kernel/piccard.cpp

@@ -9,15 +9,17 @@ namespace Kernel{
     Pl=nullptr;
   }
 
-  void Piccard::init(const Geometry* geometry_){
+  void Piccard::init(const Geometry* geometry_,double** pumps_){
     geometry=geometry_;
+    pumps=pumps_;
     horizontal_problem.init(geometry);
     overland.init(geometry);
-    all_vertical_richards.init(geometry);
+    all_vertical_richards.init(geometry,pumps);
     if(l!=nullptr) delete[] l;
     if(Pl!=nullptr) delete[] Pl;
     l=new double[geometry->nX];
     Pl=new double[geometry->nX];
+    error_x=new double[geometry->nX];
   }
 
   void
@@ -32,15 +34,18 @@ namespace Kernel{
 
     compute_Pl(previous_solution->P);
 
+    n1_previous_hydr=norm1(previous_solution->hydr,geometry->nX);
+
 
     //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.error_x=error_x;
     horizontal_problem.run();
     new_solution->hydr=horizontal_problem.hydr;
-    error+=horizontal_problem.error;
+    error+=horizontal_problem.total_error/n1_previous_hydr;
 
     //Compute Overland
     overland.in_hov=previous_solution->hov;
@@ -58,7 +63,7 @@ namespace Kernel{
     // - P
     // - hsat
     //---------------------------------------------
-    all_vertical_richards.init_indice_x_Richards();
+    all_vertical_richards.n1_previous_hydr=n1_previous_hydr;
     all_vertical_richards.dt=dt;
     all_vertical_richards.init_P=previous_solution->P;
     all_vertical_richards.previous_P=previous_solution->P;
@@ -69,6 +74,8 @@ namespace Kernel{
     //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.error_x=error_x;
+    all_vertical_richards.init_indice_x_Richards();
     all_vertical_richards.run();
     new_solution->hsat=all_vertical_richards.hsat;
     new_solution->P=all_vertical_richards.P;
@@ -102,7 +109,7 @@ namespace Kernel{
       horizontal_problem.run();
       new_solution->hydr=horizontal_problem.hydr;
 
-      error+=horizontal_problem.error;
+      error+=horizontal_problem.total_error/n1_previous_hydr;
 
       #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard::run] error (l+hydr) = \033[31m"<<error<<"\033[0m"<<endl;
@@ -129,6 +136,8 @@ namespace Kernel{
       //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.update_indice_x_Richards();
+
       all_vertical_richards.run();
       new_solution->hsat=all_vertical_richards.hsat;
       new_solution->P=all_vertical_richards.P;

src/kernel/piccard.hpp

@@ -18,6 +18,7 @@ namespace Kernel{
     //! Geometry of the simulation
     const Geometry* geometry;
 
+    double** pumps;
     //! A pointer to an HorizontalProblem object
     HorizontalProblem horizontal_problem;
 
@@ -35,7 +36,11 @@ namespace Kernel{
 
     //! Compute Pl from l and P
     void compute_Pl(double** P);
+
+    double n1_previous_hydr;
+    double* error_x;
   public:
+
     double dt;
     Solution* previous_solution;
     Solution* new_solution;
@@ -45,7 +50,7 @@ namespace Kernel{
     Piccard();
 
     //! Init from a Geometry
-    void init(const Geometry* geometry);
+    void init(const Geometry* geometry,double** pumps);
 
     //! Compute a new Solution form a previous once.
     //! The delta time depends of the number of subdivisions required by Solver::space_solution.

src/kernel/richards.cpp

@@ -22,7 +22,8 @@ namespace Kernel{
   }
 
   void
-  Richards::init(size_t ix_,size_t nZ_,double dZ_){
+  Richards::init(size_t ix_,size_t nZ_,double dZ_,double* pumps_){
+    pumps=pumps_;
     ix=ix_;
     nZ=nZ_;
     dZ=dZ_;
@@ -195,15 +196,18 @@ namespace Kernel{
    }
 
    //Compute BB
-   clear(BB,nZ-1);
-   //TODO : Add BB computation from fpump
+   for(size_t i=1;i<nZ-2;++i){
+     BB[i]=dZ*(pumps[i-1]/6+(2*pumps[i])/3+pumps[i+1]/6);
+   }
+   BB[0]=dZ*(pumps[0]/3+pumps[1]/6);
+   BB[nZ-2]=dZ*(pumps[nZ-3]/6+(2*pumps[nZ-2])/3);
 
    //Compute F
-   F[0]=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];
+   F[0]=R[0]-G[0]-I[0]-S[0]-H[0]+BB[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[i]=div_w[i]*dZ+R[i]-G[i]-I[i]-S[i]-H[i]+BB[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];
+   F[nZ-2]=div_w[nZ-2]*dZ+R[nZ-2]-G[nZ-2]-I[nZ-2]-S[nZ-2]-H[nZ-2]+BB[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

src/kernel/richards.hpp

@@ -17,6 +17,8 @@ namespace Kernel{
     static constexpr size_t max_iterations_Richards=50;
     static constexpr double tolerence_Richards=1e-10;
 
+    double* pumps;
+
     double* sup_A;
     double* diag_A;
     double* sub_A;
@@ -55,7 +57,7 @@ namespace Kernel{
     bool has_converged;
 
     Richards();
-    void init(size_t ix,size_t nZ,double dZ);
+    void init(size_t ix,size_t nZ,double dZ,double* pumps);
 
     void run();
   };

src/kernel/richards_evolutive_time.cpp

@@ -2,13 +2,14 @@
 
 namespace Kernel{
   void
-  RichardsEvolutiveTime::init(size_t ix_,const Geometry* geometry){
+  RichardsEvolutiveTime::init(size_t ix_,const Geometry* geometry,double* pumps_){
     ix=ix_;
     nZ=geometry->nZ[ix];
     Z=geometry->Z[ix];
     dX=geometry->dX;
+    pumps=pumps_;
     div_w=new double[nZ];
-    richards.init(ix,nZ,geometry->dZ[ix]);
+    richards.init(ix,nZ,geometry->dZ[ix],pumps);
   }
 
   void

src/kernel/richards_evolutive_time.hpp

@@ -15,13 +15,15 @@ namespace Kernel{
 
     double flux_bot;
     double* div_w;
-
+    double* pumps;
+    
     void compute_flux_bot();
     void compute_div_w();
     Richards richards;
   public:
     //Input
     double dt;
+
     double* init_P;
     double* previous_P; //relatively to k
 
@@ -43,7 +45,7 @@ namespace Kernel{
     bool has_converged;
 
     RichardsEvolutiveTime();
-    void init(size_t ix,const Geometry* geometry);
+    void init(size_t ix,const Geometry* geometry,double* pumps);
     void run();
   };
 

src/kernel/solver.cpp

@@ -8,8 +8,28 @@ namespace Kernel{
     input_data.load(file,true);
     file.close();
 
+    //HACK for pump
+    for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
+      Pump& P=*(*it);
+      double lX=input_data.geometry.lX;
+      P.left_init*=lX;
+      P.right_init*=lX;
+      P.delta_left_init*=lX;
+      P.delta_right_init*=lX;
+      P.left_final*=lX;
+      P.right_final*=lX;
+      P.delta_left_final*=lX;
+      P.delta_right_final*=lX;
+
+    }
+
+    pumps=new double*[input_data.geometry.nX];
+    for(size_t ix=0;ix<input_data.geometry.nX;++ix){
+      pumps[ix]=new double[input_data.geometry.nZ[ix]];
+    }
+
     //Init Piccard object
-    piccard.init(&input_data.geometry);
+    piccard.init(&input_data.geometry,pumps);
 
     //Create pool of solutions
     for(size_t i=0;i<Time::nT;++i){
@@ -72,6 +92,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();
     if(s==nullptr) return false;
@@ -92,6 +113,7 @@ namespace Kernel{
     piccard.new_solution=get_new_solution();
     piccard.dt=Time::dT/m1*3600;
 
+
     while(m1<=max_time_subdivisions){
       #ifdef MYDEBUG
       if(Debug::level>1) cout<<" [Solver::spaceSolution] n1 = "<<n1<<" and m1 = "<<m1<<endl;
@@ -125,4 +147,18 @@ namespace Kernel{
     #endif
     return piccard.new_solution;
   }
+
+  void Solver::compute_sources(double t){
+    for(size_t ix=0;ix<input_data.geometry.nX;++ix){
+      double x=ix*input_data.geometry.dX;
+      for(size_t iz=0;iz<input_data.geometry.nZ[ix];++iz){
+        double z=input_data.geometry.hbot[ix]+iz*input_data.geometry.dZ[ix];
+        pumps[ix][iz]=0;
+        for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
+          double p=(*it)->value(x,z,t/Time::T);
+          pumps[ix][iz]+=p;
+        }
+      }
+      }
+  }
 }

src/kernel/solver.hpp

@@ -23,6 +23,9 @@ namespace Kernel{
   */
   class Solver{
   private:
+
+    double** pumps;
+
     //! Piccard object used for Piccard algorithm.
     Piccard piccard;
 
@@ -56,6 +59,8 @@ namespace Kernel{
     //! Try to compute the solution at time t+1 from the solution at time t
     Solution* space_solution();
 
+    void compute_sources(double t);
+
   public:
     //! Construct a Solver from the input file input_filename
     Solver(string input_filename);
@@ -63,6 +68,8 @@ namespace Kernel{
     //! Return the geometry of the simulation, which is given in input file
     const Geometry& get_geometry() const;
 
+    const Source& get_source() const;
+
     //! Return the vertical factor (used for drawing)
     double get_vertical_factor() const;
 
@@ -93,6 +100,11 @@ namespace Kernel{
     return input_data.geometry;
   }
 
+  inline const Source&
+  Solver::get_source() const{
+    return input_data.source;
+  }
+
   inline double
   Solver::get_vertical_factor() const{
     return input_data.factor;

src/kernel/source.cpp

@@ -110,6 +110,22 @@ Pump::get_bottom_delta(double t){
   return (1-t)*delta_bottom_init+t*delta_bottom_final;
 }
 
+double
+Pump::value(double x,double z,double t){
+  double amplitude=interpol(amplitude_init,amplitude_final,t);
+  double delta_left=interpol(delta_left_init,delta_left_final,t);
+  double left=interpol(left_init,left_final,t);
+  double delta_right=interpol(delta_right_init,delta_right_final,t);
+  double right=interpol(right_init,right_final,t);
+  double delta_top=interpol(delta_top_init,delta_top_final,t);
+  double top=interpol(top_init,top_final,t);
+  double delta_bottom=interpol(delta_bottom_init,delta_bottom_final,t);
+  double bottom=interpol(bottom_init,bottom_final,t);
+  double fx=bump(x,delta_left,left,delta_right,right);
+  double fz=bump(z,delta_bottom,bottom,delta_top,top);
+  return amplitude*fx*fz;
+}
+
 Cloud::Cloud(){
   amplitude_init=1.e-4;
   left_init=0.45;

src/kernel/source.hpp

@@ -4,8 +4,11 @@
 #include <list>
 #include <fstream>
 #include <iostream>
+#include "math/algo.hpp"
+
 using namespace std;
 
+
 class Pump{
 public:
   double amplitude_init;
@@ -26,6 +29,7 @@ public:
   Pump();
   void save(fstream& file);
   void load(fstream& file);
+  double value(double x,double z,double t);
 };
 
 class Cloud{

src/math/algo.cpp

@@ -24,6 +24,15 @@ norm2(double* u,size_t n){
   return sqrt(r);
 }
 
+double
+norm1(double* u,size_t n){
+  double r=0;
+  for(size_t i=0;i<n;++i){
+    r+=abs(u[i]);
+  }
+  return r;
+}
+
 double
 error2(double* u,double* v,size_t n){
   double r=0;
@@ -49,3 +58,20 @@ display(string name,double* u,size_t n){
   }
   cout<<"--------------------------"<<endl;
 }
+
+
+double bump(double x,double delta_left,double left,double delta_right,double right){
+  double left2=left-delta_left;
+  double right2=right+delta_right;
+  if(x>left){
+    if(x<right) return 1;
+    if(x>right2) return 0;
+    return (right2-x)/delta_right;
+  }
+  else if(x<left2) return 0;
+  return (x-left2)/delta_left;
+}
+
+double interpol(double x0,double x1,double t){
+  return x0*(1-t)+x1*t;
+}

src/math/algo.hpp

@@ -18,8 +18,12 @@ using namespace std;
 */
 
 void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x);
+double norm1(double* u,size_t n);
 double norm2(double* u,size_t n);
 double error2(double* u,double* v,size_t n);
 void clear(double* u,size_t n);
 void display(string name,double* u,size_t n);
+double bump(double x,double delta_left,double left,double delta_right,double right);
+double interpol(double x0,double x1,double t);
+
 #endif

src/qt/output/output_view.cpp

@@ -110,6 +110,10 @@ QtOutputView::paintGL(){
   }
   glEnd();
 
+  //Draw pumps
+  for(auto it=solver->get_source().pumps.begin();it!=solver->get_source().pumps.end();++it){
+    drawPump(*it);
+  }
 
   drawVectors();
 }
@@ -215,6 +219,42 @@ QtOutputView::drawVectors(){
   }
 }
 
+void
+QtOutputView::drawPump(Pump* pump){
+  double time=(step*Time::dT)/Time::T;
+
+  double a=pump->get_amplitude(time);
+
+  if(a==0) return;
+  double l=pump->get_left(time);
+  double r=pump->get_right(time);
+  double b=pump->get_bottom(time);//*factor_inv;
+  double t=pump->get_top(time);//*factor_inv;
+  double dl=pump->get_left_delta(time);
+  double dr=pump->get_right_delta(time);
+  double db=pump->get_bottom_delta(time);//*factor_inv;
+  double dt=pump->get_top_delta(time);//*factor_inv;
+  glLineWidth(2);
+  if(a>0) glColor3f(0,0.6,1);
+  else glColor3f(1,0.2,0);
+  glLineStipple(3, 0xAAAA);
+  glBegin(GL_LINE_LOOP);
+  glVertex3f(l,b,1);
+  glVertex3f(r,b,1);
+  glVertex3f(r,t,1);
+  glVertex3f(l,t,1);
+  glEnd();
+  glLineWidth(1);
+  glEnable(GL_LINE_STIPPLE);
+  glBegin(GL_LINE_LOOP);
+  glVertex3f(l-dl,b-db,1);
+  glVertex3f(r+dr,b-db,1);
+  glVertex3f(r+dr,t+dt,1);
+  glVertex3f(l-dl,t+dt,1);
+  glEnd();
+  glDisable(GL_LINE_STIPPLE);
+}
+
 void
 QtOutputView::displayInfos(double x,double z){
   cout<<"---------------------------------"<<endl;

src/qt/output/output_view.hpp

@@ -39,6 +39,7 @@ public:
   void displayInfos(double x,double y);
   void drawVector(size_t ix,size_t iz,double u,double v);
   void drawVectors();
+  void drawPump(Pump* pump);
 };
 
 inline

src/qt/output/solver.cpp

@@ -8,7 +8,11 @@ QtSolver::run(){
   QTime start=QTime::currentTime();
   cout<<"[QtSolver::run] start"<<endl;
   for(size_t i=1;i<Time::nT;++i){
-    compute_next_solution();
+    bool has_converged=compute_next_solution();
+    if(!has_converged){
+      cout<<"Does not converge ... Game Over"<<endl;
+      break;
+    }
   }
   cout<<"[QtSolver::run] end"<<endl;
   QTime end=QTime::currentTime();

todo.txt

@@ -1,8 +1,10 @@
 [A finir]
-HorizontalProblem::run
+
 Overland::run
 indice_x_Richards
 
 
 [Limitation]
 Ici k0 est constant -> impact sur horizontal problem
+
+Rmove hack for pumps in solver load file