1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 |
- #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
- }
- }
|