Debug F, simple version of Thomas, true weight Richards

Makefile

@@ -32,7 +32,7 @@ EXE_INPUT = RichardsFastSlowInput
 EXE_OUTPUT = RichardsFastSlowOutput
 EXE_SOLVER = RichardsFastSlowSolver
 
-all: $(EXE_INPUT) $(EXE_OUTPUT) 
+all: $(EXE_OUTPUT)
 input: $(EXE_INPUT)
 output: $(EXE_OUTPUT)
 solver: $(EXE_SOLVER)
@@ -40,7 +40,7 @@ solver: $(EXE_SOLVER)
 $(EXE_INPUT) : src/main_input.cpp $(KERNEL_INPUT_OBJS) $(MATH_INPUT_OBJS) $(QT_INPUT_MOCS) $(QT_INPUT_OBJS)
 			$(GPP) $(FLAGS) $(QT_FLAGS) $^ $(QT_LIB) -o $@
 
-$(EXE_OUTPUT) : src/main_output.cpp $(KERNEL_OUTPUT_OBJS) $(MATH_OUTPUT_OBJS) $(QT_OUTPUT_MOCS) $(QT_OUTPUT_OBJS)
+$(EXE_OUTPUT) : src/main_output.cpp src/debug.cpp $(KERNEL_OUTPUT_OBJS) $(MATH_OUTPUT_OBJS) $(QT_OUTPUT_MOCS) $(QT_OUTPUT_OBJS)
 			$(GPP) $(FLAGS) $(QT_FLAGS) $^ $(QT_LIB) -o $@
 
 $(EXE_SOLVER) : src/main_solver.cpp $(KERNEL_SOLVER_OBJS) $(MATH_SOLVER_OBJS)

src/debug.cpp

@@ -0,0 +1,3 @@
+#include "debug.hpp"
+
+bool Debug::debug_Thomas;

src/debug.hpp

@@ -1,10 +1,34 @@
 #ifndef DEBUG_HPP
 #define DEBUG_HPP
 
+//#define MYDEBUG
+
+#include <iostream>
+
+using namespace std;
+
 class Debug{
 public:
   static const int level=2;
-  static const size_t ix=51;
+  static const size_t ix=55;
+  static void display(string name,double* u,size_t start,size_t end);
+  static void pause();
+  static bool debug_Thomas;
 };
 
+
+
+
+inline void Debug::display(string name,double* u,size_t start,size_t end){
+  cout<<"===== "<<name<<" ====="<<endl;
+  for(size_t i=start;i<end;++i){
+    cout<<i<<" : "<<u[i]<<endl;
+  }
+  cout<<"============================"<<endl;
+}
+
+inline void Debug::pause(){
+  char a;
+  cin>>a;
+}
 #endif

src/kernel/all_vertical_richards.cpp

@@ -20,8 +20,14 @@ namespace Kernel{
 
   void
   AllVerticalRichards::run(){
+    Debug::debug_Thomas=false;
+    
     size_t i_left,i_middle,i_right;
-    if(Debug::level>0) cout<<"   [AllVerticalRichards::run] start"<<endl;
+
+    #ifdef MYDEBUG
+    cout<<"   [AllVerticalRichards::run] start"<<endl;
+    #endif
+
     for(size_t ix=0;ix<geometry->nX;++ix){
       if(indice_x_Richards[ix]){
         RichardsEvolutiveTime& r=richards_evolutive_time[ix];
@@ -65,11 +71,13 @@ namespace Kernel{
       }
     }
     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;
-    }
+
+    #ifdef MYDEBUG
+    cout<<"   [AllVerticalRichards::run] ";
+    if(not has_converged) cout<<"not";
+    cout<<"converge"<<endl<<"   [AllVerticalRichards::run] end"<<endl;
+    #endif
+
   }
 
   void

src/kernel/horizontal_problem.cpp

@@ -14,14 +14,17 @@ namespace Kernel{
 
   void
   HorizontalProblem::run(){
-    if(Debug::level>0) cout<<"   [HorizontalProblem::run] start"<<endl;
+    #ifdef MYDEBUG
+    cout<<"   [HorizontalProblem::run] start"<<endl;
+    #endif
     /*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;
-
+    #ifdef MYDEBUG
+    cout<<"   [HorizontalProblem::run] end"<<endl;
+    #endif
   }
 
   void

src/kernel/overland.cpp

@@ -4,12 +4,16 @@
 namespace Kernel{
 
   void Overland::run(){
-    if(Debug::level>0) cout<<"   [Overland::run] start"<<endl;
+    #ifdef MYDEBUG
+    cout<<"   [Overland::run] start"<<endl;
+    #endif
     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;
+    #ifdef MYDEBUG
+    cout<<"   [Overland::run] end"<<endl;
+    #endif
   }
 }

src/kernel/piccard.cpp

@@ -22,7 +22,9 @@ namespace Kernel{
 
   void
   Piccard::run(){
-    if(Debug::level>0) cout<<"  [Piccard::run] start"<<endl;
+    #ifdef MYDEBUG
+    cout<<"  [Piccard::run] start"<<endl;
+    #endif
     double error=0;
     //Copy s_in.hsat in l
     memcpy(l,previous_solution->hsat,sizeof(double)*geometry->nX);
@@ -73,15 +75,24 @@ namespace Kernel{
 
     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
     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
 
       compute_l(l,all_vertical_richards.hsat,error);
-  if(Debug::level>1) cout<<"  [Piccard::run] error (l) = \033[31m"<<error<<"\033[0m"<<endl;
+
+      #ifdef MYMYDEBUG
+      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;
@@ -92,7 +103,10 @@ namespace Kernel{
       new_solution->hydr=horizontal_problem.hydr;
 
       error+=horizontal_problem.error;
+
+      #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard::run] error (l+hydr) = \033[31m"<<error<<"\033[0m"<<endl;
+      #endif
 
       //Compute Overland
       overland.in_hov=previous_solution->hov;
@@ -125,12 +139,16 @@ namespace Kernel{
       }
     }
     if(error>=tolerence_Piccard){
+      #ifdef MYDEBUG
       if(Debug::level>1) cout<<"  [Piccard]::run not converge"<<endl;
+      #endif
       has_converged=false;
       //Voir nettoyage memoire
       return;
     }
+    #ifdef MYDEBUG
     if(Debug::level>1) cout<<"  [Piccard]::run converge"<<endl;
+    #endif
     swap(l,new_solution->l);
     swap(Pl,new_solution->Pl);
     has_converged=true;

src/kernel/richards.cpp

@@ -50,28 +50,50 @@ namespace Kernel{
 
   void
   Richards::run(){
-    if(Debug::level>0 and Debug::ix==ix) cout<<"     [Richards::run] start"<<endl;
+    #ifdef MYDEBUG
+    if(Debug::ix==ix) cout<<"     [Richards::run] start"<<endl;
+    #endif
     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;
+
+    #ifdef MYDEBUG
+    if(Debug::ix==ix) cout<<"     [Richards::run] stop"<<endl;
+    #endif
 
   }
 
   bool
   Richards::weighted_run(double w){
-    if(Debug::level>0 and Debug::ix==ix) cout<<"      [Richards::weighted_run] start"<<endl;
-    if(Debug::level>0 and Debug::ix==ix) cout<<"      [Richards::weighted_run] w = "<<w<<endl;
+    Debug::debug_Thomas=(Debug::ix==ix);
+    #ifdef MYDEBUG
+    if(Debug::ix==ix){
+      cout<<"      [Richards::weighted_run] start"<<endl;
+      cout<<"      [Richards::weighted_run] w = "<<w<<endl;
+    }
+    #endif
     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;
+
+      #ifdef MYDEBUG
       if(Debug::level>1 and Debug::ix==ix) cout<<"      [Richards::weighted_run] count = "<<count<<endl;
+      #endif
+
       solve_system();
 
+      if(w<1){
+        for(size_t i=0;i<nZ;++i){
+          out_P[i]=w*out_P[i]+(1-w)*in_P[i];
+        }
+      }
       //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];
@@ -80,29 +102,39 @@ namespace Kernel{
       else{
         swap(in_P,out_P);
       }
-   }
+      #ifdef MYDEBUG
+      if(Debug::ix==ix){
+        cout<<"      [Richards::weighted_run] error = "<<error<<endl;;
+      }
+      #endif
+    }
    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]);
+     #ifdef MYDEBUG
      if(Debug::ix==ix){
        if(Debug::level>1) cout<<"      [Richards::weighted_run] converge"<<endl;
-       if(Debug::level>0) cout<<"      [Richards::weighted_run] stop"<<endl;
+       cout<<"      [Richards::weighted_run] stop"<<endl;
      }
+     #endif
      return true;
    }
-   if(Debug::level>0 and Debug::ix==ix){
+   #ifdef MYDEBUG
+   if(Debug::ix==ix){
      cout<<"      [Richards::weighted_run] not converge"<<endl;
      cout<<"      [Richards::weighted_run] stop"<<endl;
    }
+   #endif
    return false;
  }
 
  void
  Richards::solve_system(){
-   if(Debug::level>0 and Debug::ix==ix) cout<<"       [Richards::solve_system] start"<<endl;
+   #ifdef MYDEBUG
+   if(Debug::ix==ix) cout<<"       [Richards::solve_system] start"<<endl;
    assert(nZ>=3);
-
+   #endif
    //Compute A
    diag_A[0]=(Physics::phi*dZ*Physics::ds(in_P[0]))/(2*dt);
    for(size_t i=1;i<nZ-1;++i){
@@ -167,7 +199,7 @@ namespace Kernel{
    //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];
+   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];
    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];
    }
@@ -182,13 +214,15 @@ namespace Kernel{
      sup_A[i]=(sup_B[i]+sup_C[i]);
    }
    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];
-  if(Debug::level>0 and Debug::ix==ix){
+  #ifdef MYDEBUG
+  if(Debug::ix==ix){
     cout<<"       [Richards::solve_system] out = "<<out_P<<endl;
     cout<<"       [Richards::solve_system] stop"<<endl;
   }
+  #endif
 
  }
 }

src/kernel/richards_evolutive_time.cpp

@@ -35,7 +35,10 @@ namespace Kernel{
 
   void
   RichardsEvolutiveTime::run(){
+    #ifdef MYDEBUG
     if(ix==Debug::ix and Debug::level>0) cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] start"<<endl;
+    #endif
+
     compute_flux_bot();
     compute_div_w();
     size_t m2=1;
@@ -48,7 +51,10 @@ namespace Kernel{
     richards.flux_bot=flux_bot;
     //You must specify richards.Px and also create it
     while(n2<m2 and m2<=max_Richards_time_subdivisions){
+      #ifdef MYDEBUG
       if(ix==Debug::ix and Debug::level>1) cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] n2 = "<<n2<<" and m2 = "<<m2<<endl;
+      #endif
+
       if(n2==0){
         richards.previous_P=init_P;
         richards.near_P=previous_P;
@@ -73,7 +79,8 @@ namespace Kernel{
       }
     }
     has_converged=(m2<=max_Richards_time_subdivisions);
-    if(ix==Debug::ix and Debug::level>0){
+    #ifdef MYDEBUG
+    if(ix==Debug::ix){
       if(Debug::level>1){
         cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] out = "<<P<<endl;
         cout<<"    [RichardsEvolutiveTime ix="<<ix<<"]";
@@ -82,5 +89,6 @@ namespace Kernel{
       }
       cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] stop"<<endl;
     }
+    #endif
   }
 }

src/kernel/richards_evolutive_time.hpp

@@ -7,7 +7,7 @@
 namespace Kernel{
   class RichardsEvolutiveTime{
   private:
-    static constexpr size_t max_Richards_time_subdivisions=16;//256;
+    static constexpr size_t max_Richards_time_subdivisions=256;
     size_t ix;
     size_t nZ;
     double* Z;

src/kernel/solver.cpp

@@ -82,7 +82,10 @@ namespace Kernel{
 
   Solution*
   Solver::space_solution(){
-    if(Debug::level>0) cout<<" [Solver::spaceSolution] start"<<endl;
+    #ifdef MYDEBUG
+    cout<<" [Solver::spaceSolution] start"<<endl;
+    #endif
+
     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];
@@ -90,8 +93,9 @@ namespace Kernel{
     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;
+      #endif
       piccard.run();
       if(!piccard.has_converged){
         m1*=2;
@@ -110,11 +114,15 @@ namespace Kernel{
       }
     }
     if(m1>max_time_subdivisions){
-      if(Debug::level>0) cout<<" [Solver::spaceSolution] not converge"<<endl;
+      #ifdef MYDEBUG
+      cout<<" [Solver::spaceSolution] not converge"<<endl;
+      #endif
       release_solution(piccard.new_solution);
       return nullptr;
     }
-    if(Debug::level>0) cout<<" [Solver::spaceSolution] converge"<<endl;
+    #ifdef MYDEBUG
+    cout<<" [Solver::spaceSolution] converge"<<endl;
+    #endif
     return piccard.new_solution;
   }
 }

src/math/algo.cpp

@@ -1,17 +1,16 @@
 #include "algo.hpp"
 
 void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x){
-  double last=c[0]=c[0]/b[0];
-  for(size_t i=1;i<n-1;++i){
-    last=c[i]=c[i]/(b[i]-a[i-1]*last);
-  }
-  last=d[0]=d[0]/b[0];
+  //d[0]=3e-5;
   for(size_t i=1;i<n;++i){
-    last=d[i]=(d[i]-a[i-1]*last)/(b[i]-a[i-1]*c[i-1]);
+    double w=a[i-1]/b[i-1];
+    b[i]=b[i]-w*c[i-1];
+    d[i]=d[i]-w*d[i-1];
+
   }
-  last=x[n-1]=d[n-1];
-  for(int i=n-1;i>0;--i){
-    last=x[i-1]=d[i-1]-c[i-1]*last;
+  x[n-1]=d[n-1]/b[n-1];
+  for(int i=n-2;i>=0;--i){
+    x[i]=(d[i]-c[i]*x[i+1])/b[i];
   }
 }
 

src/math/algo.hpp

@@ -1,6 +1,7 @@
 #ifndef ALGO_HPP
 #define ALGO_HPP
 
+#include "debug.hpp"
 
 #include <cmath>
 #include <iostream>

src/qt/output/output_view.cpp

@@ -99,6 +99,18 @@ QtOutputView::paintGL(){
   }
   glEnd();
 
+  //Draw hsat
+  glColor3f(0,1,1);
+  glLineWidth(1);
+  glBegin(GL_LINE_STRIP);
+  for(size_t ix=0;ix<nX;ix+=1){
+    double x=ix*dX;
+    double y=solver->get_solution(step)->hsat[ix];
+    glVertex3f(x,y,1);
+  }
+  glEnd();
+
+
   drawVectors();
 }
 

src/qt/output/solver.cpp

@@ -7,9 +7,9 @@ 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);