123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149 |
- #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){
- hov[ix]=init_hov[ix]-dt*flux_soil[ix];
- //TODO :: add rain
- if(ix>89){
- hov[ix]+=dt*0.0001;
- }
- }
- apply_flow();
- for(size_t ix=0;ix<geometry->nX;++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;
- }
- }
- void Overland::apply_flow(){
- bool overflow=true;
- while(overflow){
- overflow=false;
- runoff_to_local_min();
- basin_transfer();
- }
- }
- void Overland::runoff_to_local_min(){
- // cout<<" [Overland::runoff_to_local_min]"<<endl;
- for(size_t ix=0;ix<geometry->nX;++ix){
- double dh=hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff;
- if(dh>0){
- double dh_left,dh_right;
- switch(geometry->runoff_directions[ix]){
- case Left:
- dh_left=dh;
- dh_right=0;
- break;
- case Right:
- dh_right=dh;
- dh_left=0;
- break;
- case Both:
- dh_left=dh/2;
- dh_right=dh/2;
- break;
- default:
- dh_left=0;
- dh_right=0;
- break;
- }
- size_t jx=ix;
- while(dh_left>0){
- --jx;
- if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Right){
- hov[jx]+=dh_left;
- hov[ix]-=dh_left;
- break;
- }
- else{
- double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_left);
- hov[jx]+=temp;
- hov[ix]-=temp;
- dh_left-=temp;
- }
- }
- jx=ix;
- while(dh_right>0){
- ++jx;
- if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Left){
- hov[jx]+=dh_right;
- hov[ix]-=dh_right;
- break;
- }
- else{
- double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_right);
- hov[jx]+=temp;
- hov[ix]-=temp;
- dh_right-=temp;
- }
- }
- }
- }
- /*cout<<"------------------- after runoff ----------------"<<endl;
- for(size_t ix=0;ix<100;++ix){
- cout<<ix<<" : "<<max(0.0,hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff)<<endl;;
- }*/
- }
- void Overland::basin_transfer(){
- double vol_tot=0;
- for(size_t ix=0;ix<geometry->nX;++ix){
- vol_tot+=max(hov[ix]-geometry->hsoil[ix],0.0);
- }
- //TODO :: Remove guard
- int nmax=100;
- do{
- --nmax;
- geometry->root_basin->compute_vflow(geometry->hsoil,hov);
- geometry->root_basin->repartition(geometry->hsoil,hov,geometry->runoff_directions);
- }while(geometry->root_basin->overflow and nmax>0);
- if(nmax==0){
- cerr<<"Too many loop"<<endl;
- Debug::pause();//exit(-1);
- }
- geometry->root_basin->compute_vflow(geometry->hsoil,hov);
- double vol_tot_mid=0;
- for(size_t ix=0;ix<geometry->nX;++ix){
- vol_tot_mid+=max(hov[ix]-geometry->hsoil[ix],0.0);
- }
- geometry->root_basin->flatten(geometry->hsoil,hov);
- double vol_tot_new=0;
- for(size_t ix=0;ix<geometry->nX;++ix){
- vol_tot_new+=max(hov[ix]-geometry->hsoil[ix],0.0);
- }
- }
- }
|