RichardsFastSlow/src/kernel/overland.cpp

#include "overland.hpp"
#include "debug.hpp"

namespace Kernel{

  void Overland::run(){
    #ifdef MYDEBUG
    cout<<"   [Overland::run] start"<<endl;
    #endif
    compute_flux_soil();
    total_error=0;
    for(size_t ix=0;ix<geometry->nX;++ix){
      //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]);
    }
    #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

    }
  }
}