123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180 |
- #include "solver.hpp"
- namespace Kernel{
- Solver::Solver(string filename){
- fstream file;
- file.open(filename.c_str(),fstream::in|fstream::binary);
- input_data.load(file,true);
- input_data.geometry.compute_basins();
- file.close();
- //HACK for pump
- for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
- Pump& P=*(*it);
- double lX=input_data.geometry.lX;
- P.left_init*=lX;
- P.right_init*=lX;
- P.delta_left_init*=lX;
- P.delta_right_init*=lX;
- P.left_final*=lX;
- P.right_final*=lX;
- P.delta_left_final*=lX;
- P.delta_right_final*=lX;
- }
- pumps=new double*[input_data.geometry.nX];
- for(size_t ix=0;ix<input_data.geometry.nX;++ix){
- pumps[ix]=new double[input_data.geometry.nZ[ix]];
- }
- //Init Piccard object
- piccard.init(&input_data.geometry,pumps);
- //Create pool of solutions
- for(size_t i=0;i<Time::nT;++i){
- Solution* s=new Solution(input_data.geometry);
- solutions_stack.push(s);
- }
- //Init solutions
- solutions=new Solution*[Time::nT];
- for(size_t i=1;i<Time::nT;++i){
- solutions[i]=nullptr;
- }
- //Init first solution
- solutions[0]=new Solution(input_data.geometry);
- init_first_solution();
- m1=1;
- }
- void
- Solver::init_first_solution(){
- Solution& s=*solutions[0];
- //Init pressure
- for(size_t x=0;x<input_data.geometry.nX;++x){
- for(size_t z=0;z<input_data.geometry.nZ[x];++z){
- s.P[x][z]=input_data.initial_state->Pinit[x][z];
- }
- }
- //Init hydraulic head
- for(size_t x=0;x<input_data.geometry.nX;++x){
- s.hydr[x]=input_data.initial_state->hsat[x]+Psat/(Physics::rho*Physics::g);
- }
- //[Obsolete]
- ////Init overland level
- //for(size_t x=0;x<input_data.geometry.nX;++x){
- // s.hov[x]=Physics::invPsol(input_data.geometry.hsoil[x],s.P[x][input_data.geometry.nZ[x]-1]); //TODO[Chrisophe] à initiliser
- //}
- for(size_t x=0;x<input_data.geometry.nX;++x){
- s.hov[x]=input_data.initial_state->hov[x];
- }
- //Init l
- for(size_t x=0;x<input_data.geometry.nX;++x){
- double t=input_data.initial_state->hsat[x];
- s.l[x]=t;
- s.hsat[x]=t;
- }
- //Init Pl
- for(size_t x=0;x<input_data.geometry.nX;++x){
- s.Pl[x]=Psat;
- }
- //Initial state of input_data will be no longer used.
- input_data.remove_initial_state();
- number_computed_solutions=1;
- }
- bool
- Solver::compute_next_solution(){
- compute_sources(number_computed_solutions*Time::dT);
- if(m1>1) m1=m1/2;
- Solution* s=space_solution();
- if(s==nullptr) return false;
- solutions[number_computed_solutions++]=s;
- //cout<<"[next] done"<<endl;
- return true;
- }
- Solution*
- Solver::space_solution(){
- #ifdef MYDEBUG
- cout<<" [Solver::spaceSolution] start"<<endl;
- #endif
- size_t n1=0; //we work on time interval [n-1+n1/m1,n-1+(n1+1)/m1]
- piccard.previous_solution=solutions[number_computed_solutions-1];
- piccard.new_solution=get_new_solution();
- piccard.dt=Time::dT/m1*3600;
- while(m1<=max_time_subdivisions){
- #ifdef MYDEBUG
- if(Debug::level>1) cout<<" [Solver::spaceSolution] n1 = "<<n1<<" and m1 = "<<m1<<endl;
- #endif
- piccard.run();
- if(!piccard.has_converged){
- if(n1!=0){//S_in is different form solution[n-1]
- release_solution(piccard.previous_solution);
- }
- m1*=2;
- n1=0;
- piccard.dt/=2;
- piccard.previous_solution=solutions[number_computed_solutions-1];
- }
- else{
- if(n1!=0){//S_in is different form solution[n-1]
- release_solution(piccard.previous_solution);
- }
- ++n1;
- if(n1==m1) break;
- piccard.previous_solution=piccard.new_solution;
- piccard.new_solution=get_new_solution();
- }
- }
- if(m1>max_time_subdivisions){
- #ifdef MYDEBUG
- cout<<" [Solver::spaceSolution] not converge"<<endl;
- #endif
- release_solution(piccard.new_solution);
- return nullptr;
- }
- #ifdef MYDEBUG
- cout<<" [Solver::spaceSolution] converge"<<endl;
- #endif
- return piccard.new_solution;
- }
- void Solver::compute_sources(double t){
- for(size_t ix=0;ix<input_data.geometry.nX;++ix){
- double x=ix*input_data.geometry.dX;
- for(size_t iz=0;iz<input_data.geometry.nZ[ix];++iz){
- double z=input_data.geometry.hbot[ix]+iz*input_data.geometry.dZ[ix];
- pumps[ix][iz]=0;
- for(auto it=input_data.source.pumps.begin();it!=input_data.source.pumps.end();++it){
- double p=(*it)->value(x,z,t/Time::T);
- pumps[ix][iz]+=p;
- }
- }
- }
- }
- }
|