123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219 |
- #include "basin.hpp"
- void
- Basin::compute_leak_direction(){
- switch(position){
- case Middle:
- if(abs(hleft-hright)<delta_both_leak) leak_direction=Both;
- else leak_direction=(hleft>hright)?Right:Left;
- break;
- case BorderLeft:
- leak_direction=Right;
- break;
- case BorderRight:
- leak_direction=Left;
- break;
- default:
- leak_direction=None;
- break;
- };
- }
- void Basin::compute_maxs(double* hsoil){
- switch(leak_direction){
- case Left:
- hmax=hleft;
- xmax_left=xleft;
- {
- int x;
- for(x=xleft+1;x<xright;++x){
- if(hsoil[x]+Physics::minimal_height_to_runoff>hmax){
- xmax_right=x;
- break;
- }
- }
- xmax_right=min(x,(int)xright);
- }
- break;
- case Right:
- hmax=hright;
- xmax_right=xright;
- {
- int x;
- for(x=xright-1;x>xleft;--x){
- if(hsoil[x]+Physics::minimal_height_to_runoff>hmax){
- xmax_left=x;
- break;
- }
- }
- xmax_left=max(x,(int)xleft);
- }
- break;
- case Both:
- hmax=min(hleft,hright);
- xmax_left=xleft;
- xmax_right=xright;
- break;
- case None:
- default:
- hmax=numeric_limits<double>::infinity();
- xmax_left=xleft;
- xmax_right=xright;
- break;
- };
- }
- Basin::Basin(Basin* left_,Basin* right_){
- left=left_;
- right=right_;
- xleft=left->xleft;
- xright=right->xright;
- hleft=left->hleft;
- hright=right->hright;
- hmin=left->hmax;
- xacc=left->xacc;
- if(left->position==BorderLeft){
- if(right->position==BorderRight) position=BorderBoth;
- else position=BorderLeft;
- }
- else if(right->position==BorderRight) position=BorderRight;
- else position=Middle;
- compute_leak_direction();
- }
- void
- Basin::display(string indent){
- cout<<indent<<'['<<xleft<<','<<xright<<"] : "<<hmax<<" : ["<<xmax_left<<","<<xmax_right<<"] and "<<hmin<<endl;
- if(left==nullptr) return;
- left->display(indent+' ');
- right->display(indent+' ');
- }
- void
- Basin::compute_vflow(double* hsoil,double* hov){
- if(is_primitive()){
- vflow=max(0.0,hov[xacc]-hsoil[xacc]-Physics::minimal_height_to_runoff);
- }
- else{
- left->compute_vflow(hsoil,hov);
- right->compute_vflow(hsoil,hov);
- vflow=left->vflow+right->vflow;
- }
- }
- void
- Basin::repartition(double* hsoil,double* hov,Direction* rundir){
- overflow=false;
- vmin=0;
- if(not is_primitive()){
- for(size_t ix=xleft;ix<=xright;++ix){
- vmin+=max(0.0,hmin-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
- }
- }
- if(vflow!=0 and vflow>=vmin){
- vmax=0;
- for(size_t ix=xmax_left;ix<=xmax_right;++ix){
- vmax+=max(0.0,hmax-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
- }
- if(vflow-vmax>1e-10){
- overflow=true;
- vleak=vflow-vmax;
- if(not is_primitive()){
- double delta=hov[right->xacc]-hsoil[right->xacc]-Physics::minimal_height_to_runoff;
- hov[left->xacc]+=delta;
- hov[right->xacc]-=delta;
- }
- hov[xacc]-=vleak;
- switch(leak_direction){
- case Left:
- hov[xleft-1]+=vleak;
- runoff(xleft-1,Left,rundir,hsoil,hov);
- break;
- case Right:
- hov[xright+1]+=vleak;
- runoff(xright+1,Right,rundir,hsoil,hov);
- break;
- case Both:
- hov[xleft-1]+=vleak/2;
- runoff(xleft-1,Left,rundir,hsoil,hov);
- hov[xright+1]+=vleak/2;
- runoff(xright+1,Right,rundir,hsoil,hov);
- break;
- default:
- cerr<<"[Bug] Impossible case ! ... in theory :)"<<endl;
- exit(-1);
- break;
- }
- }
- }
- else{
- if(not is_primitive()){
- left->repartition(hsoil,hov,rundir);
- right->repartition(hsoil,hov,rundir);
- overflow=left->overflow|right->overflow;
- }
- }
- }
- void
- Basin::runoff(size_t xstart,Direction dir,Direction* rundir,double* hsoil,double* hov){
- double r=max(0.0,hov[xstart]-hsoil[xstart]-Physics::minimal_height_to_runoff);
- int x=xstart;
- int xstep=(dir==Left)?-1:1;
- while(r>0 and rundir[x]!=None){
- hov[x]-=r;
- x+=xstep;
- hov[x]+=r;
- r=max(0.0,hov[x]-hsoil[x]-Physics::minimal_height_to_runoff);
- }
- }
- double
- Basin::f(double r,double* hsoil,double* hov,double vext){
- double res=-vext;
- for(size_t ix=xleft;ix<=xright;++ix){
- res+=max(r-min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff),0.0);
- }
- return res;
- }
- double
- Basin::df(double r,double* hsoil,double* hov,double vext){
- double res=0;
- for(size_t ix=xleft;ix<=xright;++ix){
- if(r>min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff)) ++res;
- }
- return res;
- }
- void
- Basin::flatten(double* hsoil,double* hov){
- double crit=1;
- size_t count=0;
- if(vflow>1e-9){
- double r=hsoil[xacc]+1;//Physics::minimal_height_to_runoff;
- while(crit>1e-10 and count<100){
- ++count;
- double rold=r;
- r-=f(r,hsoil,hov,vflow)/df(r,hsoil,hov,vflow);
- crit=abs(r-rold)+abs(f(r,hsoil,hov,vflow));
- }
- if(crit>1e-10){
- cerr<<" Overland coverge + [TODO] on fait quoi ?"<<endl;
- }
- if(r>=hmin or is_primitive()){
- for(size_t ix=xleft;ix<=xright;++ix){
- hov[ix]=max(r,min(hov[ix],hsoil[ix]+Physics::minimal_height_to_runoff));
- }
- }
- else{
- left->flatten(hsoil,hov);
- right->flatten(hsoil,hov);
- }
- }
- }
|