RichardsFastSlow/src/kernel/horizontal_problem.cpp

#include "horizontal_problem.hpp"
#include "debug.hpp"

namespace Kernel{

  void
  HorizontalProblem::init(const Geometry* geometry_){
    geometry=geometry_;
    sub_M=new double[geometry->nX-1];
    diag_M=new double[geometry->nX];
    sup_M=new double[geometry->nX-1];
    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
    cout<<"   [HorizontalProblem::run] start"<<endl;
    #endif
    /*for(size_t ix=0;ix<geometry->nX;++ix){
      hydr[ix]=previous_hydr[ix];
    }*/
    solve_system();
    compute_error();
    #ifdef MYDEBUG
    cout<<"   [HorizontalProblem::run] end"<<endl;
    #endif
  }

  void
  HorizontalProblem::solve_system(){
    for(size_t ix=1;ix<geometry->nX-1;++ix){
      double dx=geometry->dX;
      double hbot=geometry->hbot[ix];
      double lx=l[ix];
      double a=-Physics::tilde_a(lx,hbot)*Physics::k0/(dx*dx);
      double c=Physics::tilde_c(lx,hbot)*Physics::k0/dx;
      double dhbot=geometry->dhbot[ix];
      sub_M[ix-1]=a-(dhbot*c)/2;
      diag_M[ix]=1-2*a;
      sup_M[ix]=a+(dhbot*c)/2;
      F[ix]=Pl[ix]/(Physics::rho*Physics::g)+lx;
    }
    diag_M[0]=1;
    sup_M[0]=-1;
    F[0]=0;
    diag_M[geometry->nX-1]=1;
    sub_M[geometry->nX-2]=-1;
    F[geometry->nX-1]=0;
    Thomas(geometry->nX,sub_M,diag_M,sup_M,F,hydr);
  }

}