123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107 |
- #include "all_vertical_richards.hpp"
- namespace Kernel{
- void
- AllVerticalRichards::init(const Geometry* geometry_,double** pumps_){
- geometry=geometry_;
- pumps=pumps_;
- indice_x_Richards=new bool[geometry->nX];
- richards_evolutive_time=new RichardsEvolutiveTime[geometry->nX];
- for(size_t ix=0;ix<geometry->nX;++ix){
- richards_evolutive_time[ix].init(ix,geometry,pumps[ix]);
- }
- }
- void
- AllVerticalRichards::update_indice_x_Richards(){
- for(size_t ix=0;ix<geometry->nX;++ix){
- indice_x_Richards[ix]=(error_x[ix]>n1_previous_hydr*max_error_x);
- }
- }
- void
- AllVerticalRichards::init_indice_x_Richards(){
- for(size_t ix=0;ix<geometry->nX;++ix){
- indice_x_Richards[ix]=true;
- }
- }
- void
- AllVerticalRichards::run(){
- Debug::debug_Thomas=false;
- size_t i_left,i_middle,i_right;
- #ifdef MYDEBUG
- cout<<" [AllVerticalRichards::run] start"<<endl;
- #endif
- for(size_t ix=0;ix<geometry->nX;++ix){
- if(indice_x_Richards[ix]){
- RichardsEvolutiveTime& r=richards_evolutive_time[ix];
- r.dt=dt;
- r.init_P=init_P[ix];
- r.previous_P=previous_P[ix];
- r.P=P[ix];
- if(ix==0){
- i_left=0;
- i_middle=1;
- i_right=2;
- }
- else if(ix==geometry->nX-1){
- i_left=geometry->nX-3;
- i_middle=geometry->nX-2;
- i_right=geometry->nX-1;
- }
- else{
- i_left=ix-1;
- i_middle=ix;
- i_right=ix+1;
- }
- r.hydr_left=hydr[i_left];
- r.l_left=l[i_left];
- r.Pl_left=Pl[i_left];
- r.hydr_middle=hydr[i_middle];
- r.l_middle=l[i_middle];
- r.Pl_middle=Pl[i_middle];
- r.hydr_right=hydr[i_right];
- r.l_right=l[i_right];
- r.Pl_right=Pl[i_right];
- r.run();
- P[ix]=r.P;
- }
- }
- //[If Paraller wait for sync]
- has_converged=true;
- for(size_t ix=0;ix<geometry->nX;++ix){
- if(!richards_evolutive_time[ix].has_converged){
- has_converged=false;
- }
- }
- if(has_converged) compute_hsat();
- #ifdef MYDEBUG
- cout<<" [AllVerticalRichards::run] ";
- if(not has_converged) cout<<"not";
- cout<<"converge"<<endl<<" [AllVerticalRichards::run] end"<<endl;
- #endif
- }
- void
- AllVerticalRichards::compute_hsat(){
- for(size_t ix=0;ix<geometry->nX;++ix){
- double* Px=previous_P[ix];
- size_t iz=0;
- for(;iz<geometry->nZ[ix] and Px[iz]>Psat;++iz);
- if(iz>0 and Px[iz]<=Psat){
- hsat[ix]=(Psat-Px[iz-1])*geometry->dZ[ix]/(Px[iz]-Px[iz-1])+geometry->Z[ix][iz-1];
- }
- else{
- hsat[ix]=(iz==0)?geometry->hbot[ix]:geometry->hsoil[ix];
- }
- }
- }
- }
|