Start to add overland

Makefile

@@ -32,7 +32,7 @@ EXE_INPUT = RichardsFastSlowInput
 EXE_OUTPUT = RichardsFastSlowOutput
 EXE_SOLVER = RichardsFastSlowSolver
 
-all: $(EXE_OUTPUT)
+all: $(EXE_OUTPUT) $(EXE_INPUT)
 input: $(EXE_INPUT)
 output: $(EXE_OUTPUT)
 solver: $(EXE_SOLVER)

src/debug.hpp

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

src/kernel/all_vertical_richards.cpp

@@ -15,7 +15,7 @@ namespace Kernel{
   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);
+      indice_x_Richards[ix]=(error_x[ix]>max_error_x);
     }
   }
 
@@ -39,12 +39,23 @@ namespace Kernel{
     #endif
 
     for(size_t ix=0;ix<geometry->nX;++ix){
+
+      P[ix][geometry->nZ[ix]-1]=Psoil[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];
+
+
+        //Update Psoil
+  /*      r.P[geometry->nZ[ix]-1]=Psoil[ix];
+        r.init_P[geometry->nZ[ix]-1]=Psoil[ix];
+        r.previous_P[geometry->nZ[ix]-1]=Psoil[ix];*/
+
+
+
         if(ix==0){
           i_left=0;
           i_middle=1;

src/kernel/all_vertical_richards.hpp

@@ -17,8 +17,7 @@ namespace Kernel{
     //Input
     double dt;
     double* error_x;
-    double n1_previous_hydr;
-
+  
     double** init_P; //P_0
     double** previous_P; //P_{k-1}
     double* hydr;

src/kernel/horizontal_problem.cpp

@@ -10,13 +10,16 @@ namespace Kernel{
     diag_M=new double[geometry->nX];
     sup_M=new double[geometry->nX-1];
     F=new double[geometry->nX];
+    hydr=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]);
+      double e=abs(previous_hydr[ix]-hydr[ix])/n1_init_hydr;
+      error_x[ix]=e;
       total_error+=e;
     }
   }

src/kernel/horizontal_problem.hpp

@@ -23,6 +23,7 @@ namespace Kernel{
     double* previous_hydr;
     double* l;
     double* Pl;
+    double n1_init_hydr;
     //Output
     double* hydr;
     double total_error;

src/kernel/overland.cpp

@@ -7,13 +7,42 @@ namespace Kernel{
     #ifdef MYDEBUG
     cout<<"   [Overland::run] start"<<endl;
     #endif
+    compute_flux_soil();
+    total_error=0;
     for(size_t ix=0;ix<geometry->nX;++ix){
-      hov[ix]=in_hov[ix];
-      Psoil[ix]=previous_P[ix][geometry->nZ[ix]-1];
+      //TODO :: add rain
+      hov[ix]=init_hov[ix]-dt*flux_soil[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]);
     }
-    error=0;
+    #ifdef MYDEBUG
+    cout<<"   [Overland::run] error : "<<total_error<<endl;
+    #endif
     #ifdef MYDEBUG
     cout<<"   [Overland::run] end"<<endl;
     #endif
   }
+
+  void Overland::compute_flux_soil(){
+    for(size_t ix=0;ix<geometry->nX;++ix){
+      double p1=P[ix][geometry->nZ[ix]-1];
+      double p2=P[ix][geometry->nZ[ix]-2];
+
+      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;
+      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
+
+    }
+  }
 }

src/kernel/overland.hpp

@@ -2,22 +2,29 @@
 #define OVERLAND_HPP
 
 #include "geometry.hpp"
+#include "physics.hpp"
 
 namespace Kernel{
   class Overland{
   private:
     const Geometry* geometry;
+    double* flux_soil;
+    void compute_flux_soil();
   public:
     //Input
-    double* in_hov;
-    double** previous_P;
+    double dt;
+    double* init_hov;
+    double** P;
     double* l;
     double* Pl;
     double* hydr;
+    double n1_init_hov;
+    double* previous_hov;
     //Output
     double* Psoil; //Own by Overland object
     double* hov;
-    double error;
+    double* error_x;
+    double total_error;
 
     Overland();
     void init(const Geometry* geometry);
@@ -32,7 +39,10 @@ namespace Kernel{
   Overland::init(const Geometry* geometry_){
     geometry=geometry_;
     Psoil=new double[geometry->nX];
+    flux_soil=new double[geometry->nX];
+    hov=new double[geometry->nX];
   }
 
+
 }
 #endif

src/kernel/physics.cpp

@@ -138,20 +138,20 @@ Physics::load(fstream& file){
 }
 
 double
-Physics::Psol(double hsoil,double hov){
+Physics::Psoil(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){
+Physics::dPsoil(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){
+Physics::invPsoil(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

@@ -86,9 +86,9 @@ public:
   //! 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 Psoil(double hsoil,double hov);
+  static double dPsoil(double hsoil,double hov);
+  static double invPsoil(double hsoil,double Psol);
 
   static double tilde_a(double l,double hbot);
   static double tilde_c (double l,double hbot);

src/kernel/piccard.cpp

@@ -34,40 +34,43 @@ namespace Kernel{
 
     compute_Pl(previous_solution->P);
 
-    n1_previous_hydr=norm1(previous_solution->hydr,geometry->nX);
-
+    double n1_init_hydr=norm1(previous_solution->hydr,geometry->nX);
+    double n1_init_hov=norm1(previous_solution->hov,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.n1_init_hydr=n1_init_hydr;
     horizontal_problem.error_x=error_x;
     horizontal_problem.run();
-    new_solution->hydr=horizontal_problem.hydr;
-    error+=horizontal_problem.total_error/n1_previous_hydr;
+    swap(new_solution->hydr,horizontal_problem.hydr);
+    error+=horizontal_problem.total_error;
 
     //Compute Overland
-    overland.in_hov=previous_solution->hov;
-    overland.previous_P=previous_solution->P;
+    overland.dt=dt;
+    overland.init_hov=previous_solution->hov;
+    overland.previous_hov=previous_solution->hov;
+    overland.P=previous_solution->P;
     overland.l=previous_solution->hsat;
     overland.Pl=Pl;
-    overland.hydr=horizontal_problem.hydr;
-    overland.hov=new_solution->hov;
+    overland.hydr=new_solution->hydr;
+    overland.n1_init_hov=n1_init_hov;
+    overland.error_x=error_x;
     overland.run();
-    new_solution->hov=overland.hov;
-    error+=overland.error;
+    swap(new_solution->hov,overland.hov);
+    error+=overland.total_error;
 
     //----------------------------------------------
     // Apply AllVerticalRichads algorithm to obtain
     // - P
     // - hsat
     //---------------------------------------------
-    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;
-    all_vertical_richards.hydr=horizontal_problem.hydr;
+    all_vertical_richards.hydr=new_solution->hydr;
     all_vertical_richards.l=previous_solution->hsat;
     all_vertical_richards.Pl=Pl;
     all_vertical_richards.Psoil=overland.Psoil;
@@ -79,12 +82,9 @@ namespace Kernel{
     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=numeric_limits<double>::infinity();
-    #ifdef MYDEBUG
-    if(Debug::level>1) cout<<"  [Piccard::run] error (horizontal) = \033[31m"<<error<<"\033[0m"<<endl;
-    #endif
+    //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;
@@ -96,40 +96,43 @@ namespace Kernel{
 
       compute_l(l,all_vertical_richards.hsat,error);
 
-      #ifdef MYMYDEBUG
+      #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard::run] error (l) = \033[31m"<<error<<"\033[0m"<<endl;
       #endif
 
       compute_Pl(all_vertical_richards.P);
-
-      horizontal_problem.previous_hydr=horizontal_problem.hydr;
+      horizontal_problem.previous_hydr=new_solution->hydr;
       horizontal_problem.l=l;
       horizontal_problem.Pl=Pl;
-      horizontal_problem.hydr=new_solution->hydr;
       horizontal_problem.run();
-      new_solution->hydr=horizontal_problem.hydr;
+      swap(new_solution->hydr,horizontal_problem.hydr);
+
+      //new_solution->hydr=horizontal_problem.hydr;
 
-      error+=horizontal_problem.total_error/n1_previous_hydr;
+      error+=horizontal_problem.total_error;
 
       #ifdef MYDEBUG
-      if(Debug::level>1) cout<<"  [Piccard::run] error (l+hydr) = \033[31m"<<error<<"\033[0m"<<endl;
+      if(Debug::level>1) cout<<"  [Piccard::run] error (l+hori) = \033[31m"<<error<<"\033[0m"<<endl;
       #endif
 
       //Compute Overland
-      overland.in_hov=previous_solution->hov;
-      overland.previous_P=all_vertical_richards.P;
+      overland.previous_hov=new_solution->hov;
+      overland.P=new_solution->P;
       overland.l=l;
       overland.Pl=Pl;
-      overland.hydr=horizontal_problem.hydr;
-      overland.hov=new_solution->hov;
+      overland.hydr=new_solution->hydr;
       overland.run();
-      new_solution->hov=overland.hov;
-      error+=overland.error;
+      swap(new_solution->hov,overland.hov);
+      error+=overland.total_error;
+
+      #ifdef MYDEBUG
+      if(Debug::level>1) cout<<"  [Piccard::run] error (l+hori+overland) = \033[31m"<<error<<"\033[0m"<<endl;
+      #endif
 
       //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.hydr=new_solution->hydr;
       all_vertical_richards.l=l;
       all_vertical_richards.Pl=Pl;
       all_vertical_richards.Psoil=overland.Psoil;
@@ -146,6 +149,8 @@ namespace Kernel{
         has_converged=false;
         return;
       }
+      //Debug::pause();
+
     }
     if(error>=tolerence_Piccard){
       #ifdef MYDEBUG

src/kernel/piccard.hpp

@@ -37,7 +37,6 @@ namespace Kernel{
     //! Compute Pl from l and P
     void compute_Pl(double** P);
 
-    double n1_previous_hydr;
     double* error_x;
   public:
 

src/kernel/richards.cpp

@@ -220,7 +220,7 @@ namespace Kernel{
    diag_A[nZ-2]+=(diag_B[nZ-2]+diag_C[nZ-2]);
    //Debug::display("sup_A",sup_A,0,nZ-2);
   Thomas(nZ-1,sub_A,diag_A,sup_A,F,out_P);
-  out_P[nZ-1]=in_P[nZ-1];
+  out_P[nZ-1]=P[nZ-1];
   #ifdef MYDEBUG
   if(Debug::ix==ix){
     cout<<"       [Richards::solve_system] out = "<<out_P<<endl;

src/kernel/solver.cpp

@@ -66,11 +66,16 @@ namespace Kernel{
       s.hydr[x]=input_data.initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
     }
 
-    //Init overland level
+    //[Obsolete]
+    ////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
+    //}
     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
+      s.hov[x]=input_data.initial_state->hov[x];
     }
 
+
     //Init l
     for(size_t x=0;x<input_data.geometry.nX;++x){
       double t=input_data.initial_state->hsat[x];
@@ -120,6 +125,10 @@ namespace Kernel{
       #endif
       piccard.run();
       if(!piccard.has_converged){
+        if(n1!=0){//S_in is different form solution[n-1]
+          release_solution(piccard.previous_solution);
+        }
+
         m1*=2;
         n1=0;
         piccard.dt/=2;

src/kernel/solver.hpp

@@ -85,6 +85,9 @@ namespace Kernel{
 
   inline Solution*
   Solver::get_new_solution(){
+    if(solutions_stack.empty()){
+      cerr<<"[Error] Solution stack is empty"<<endl;
+    }
     Solution* s=solutions_stack.top();
     solutions_stack.pop();
     return s;

src/qt/input/data.cpp

@@ -133,6 +133,7 @@ QtInputData::updateInitialState(){
   }
   initial_state->updatePressure();
   updateOverland();
+  updateOverlandAndPressure();
 }
 
 void
@@ -168,7 +169,6 @@ QtInputData::load(fstream& file){
   updateSplineSoil();
   updateGeometry();
   updateInitialState();
-  updateOverland();
 }
 
 QtInputData::~QtInputData(){
@@ -202,13 +202,37 @@ QtInputData::updateOverland(){
     if(geometry.hsoil[x_ov]<y){
       size_t x_left=findSoilOver(y,x_ov,-1);
       size_t x_right=findSoilOver(y,x_ov,1);
-        for(size_t x=x_left;x<=x_right;++x){
+      for(size_t x=x_left;x<=x_right;++x){
         initial_state->hov[x]=max(initial_state->hov[x],y);
       }
     }
   }
 }
 
+void
+QtInputData::updateOverlandAndPressure(){
+  for(size_t x=0;x<geometry.nX;++x){
+    double hsoil=geometry.hsoil[x];
+    double hov=initial_state->hov[x];
+    if(hov==hsoil){
+      //No overland water and set hov to coincide with Psoil
+      double Psoil=initial_state->Pinit[x][geometry.nZ[x]-1];
+      double d=Physics::invPsoil(hsoil,Psoil);
+      initial_state->hov[x]=Physics::invPsoil(hsoil,Psoil);
+    }
+    else{
+    double Psoil=Physics::Psoil(hsoil,hov);
+    double* Pinit_x=initial_state->Pinit[x];
+    for(size_t z=0;z<geometry.nZ[x];++z){
+      Pinit_x[z]=max(Pinit_x[z],20000*(geometry.Z[x][z]-geometry.hsoil[x])+Psoil);
+    }
+      //Overland water and set pressure in consequence
+
+
+    }
+  }
+}
+
 size_t
 QtInputData::findSoilOver(double y,int x_ov,int x_step){
   int x=x_ov;

src/qt/input/data.hpp

@@ -29,6 +29,7 @@ public:
   void updateGeometry();
   void updateInitialState();
   void updateOverland();
+  void updateOverlandAndPressure();
   double evalSplineBot(double x) const;
   double evalSplineSat(double x) const;
   double evalSplineSoil(double x) const;

src/qt/input/view.cpp

@@ -317,5 +317,5 @@ QtInputView::moveSelected(double x,double y){
 
 void
 QtInputView::mouseReleaseEvent(QMouseEvent* event){
-  if(selected<3*QtInputData::np) emit geometryChanged();
+  if(selected<3*QtInputData::np+QtInputData::np_ov) emit geometryChanged();
 }

src/qt/output/main_window.cpp

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

src/qt/output/output_view.cpp

@@ -115,6 +115,7 @@ QtOutputView::paintGL(){
     drawPump(*it);
   }
 
+drawOverland();
   drawVectors();
 }
 
@@ -255,6 +256,22 @@ QtOutputView::drawPump(Pump* pump){
   glDisable(GL_LINE_STIPPLE);
 }
 
+
+void
+QtOutputView::drawOverland(){
+  const Geometry& G=solver->get_geometry();
+  double dX=G.dX;
+  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],0);
+    glVertex3f(i*dX,solver->get_solution(step)->hov[i],0);
+    glVertex3f(i*dX+dX,solver->get_solution(step)->hov[i+1],0);
+    glVertex3f(i*dX+dX,G.hsoil[i+1],0);
+  }
+  glEnd();
+}
+
 void
 QtOutputView::displayInfos(double x,double z){
   cout<<"---------------------------------"<<endl;

src/qt/output/output_view.hpp

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

todo.txt

@@ -8,3 +8,12 @@ indice_x_Richards
 Ici k0 est constant -> impact sur horizontal problem
 
 Rmove hack for pumps in solver load file
+
+[Interface]
+Add flux bot comme nuage mais sans temps (etenacheite sol)
+Voir problem locale . ,
+
+[Optimisation]
+
+[Overland]
+Add rain