RichardsFastSlow/src/kernel/richards_evolutive_time.cpp

#include "richards_evolutive_time.hpp"

namespace Kernel{
  void
  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],pumps);
  }

  void
  RichardsEvolutiveTime::compute_flux_bot(){
    flux_bot=0;
  }

  void
  RichardsEvolutiveTime::compute_div_w(){

    double cst=Physics::rho*Physics::g;
    double G_left=Pl_left/cst+l_left;
    double G_middle=Pl_middle/cst+l_middle;
    double G_right=Pl_right/cst+l_right;
    double m_left=max(G_left,G_middle);
    double m_right=max(G_middle,G_right);

    for(size_t i=0;i<nZ;++i){
      div_w[i]=(Physics::kr(cst*(m_right-Z[i]))*Physics::k0*(hydr_right-hydr_middle)-Physics::kr(cst*(m_left-Z[i]))*Physics::k0*(hydr_middle-hydr_left))/(dX*dX);
    }


  }

  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;
    size_t n2=0;

    double dt_sub=dt;

    richards.div_w=div_w;
    richards.dt=dt_sub;
    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;
      }
      else{
        richards.previous_P=richards.P;
        richards.near_P=richards.P;
      }
      richards.P=P;

      richards.run();
      P=richards.P;

      if(!richards.has_converged){
        m2*=2;
        dt_sub/=2;
        richards.dt=dt_sub;
        n2=0;
      }
      else{
        ++n2;
      }
    }
    has_converged=(m2<=max_Richards_time_subdivisions);
    #ifdef MYDEBUG
    if(ix==Debug::ix){
      if(Debug::level>1){
        cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] out = "<<P<<endl;
        cout<<"    [RichardsEvolutiveTime ix="<<ix<<"]";
        if(not has_converged) cout<<" not";
        cout<<" converge"<<endl;
      }
      cout<<"    [RichardsEvolutiveTime ix="<<ix<<"] stop"<<endl;
    }
    #endif
  }
}