123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245 |
- #include "qt/input/data.hpp"
- QtInputData::QtInputData(){
- depth=1;
- nZ_max=100;
- initPhysics();
- initSplines();
- initGeometry();
- initInitialState();
- }
- void QtInputData::initPhysics(){
- Physics::g=9.81;
- Physics::rho=1000;
- Physics::phi=0.3;
- Physics::k0=3e-5;
- Physics::nivrivsat=0.02;
- Physics::Psec=-10000;
- Physics::model_data[0]=-2000;
- Physics::model_data[1]=0;
- Physics::model_data[2]=3;
- Physics::model_data[3]=11;
- Physics::model_data[4]=0;
- Physics::model_data[5]=0;
- Physics::setModel(Physics::BrooksCorey);
- }
- void
- QtInputData::initSplines() {
- // Control points of splines
- for(size_t i=0;i<np;++i){
- double x=double(i)/(np-1);
- point[i].x=x;
- point[i].y=0.2;
- point[np+i].x=x;
- point[np+i].y=0.8;
- point[2*np+i].x=x;
- point[2*np+i].y=0.5;
- }
- // Control points of oveland points
- for(size_t i=0;i<np_ov;++i){
- double x=double(i+0.5)/np_ov;
- point[3*np+i].x=x;
- point[3*np+i].y=0.75;
- }
- // Splines
- spline_bot.setPoints(point,np);
- spline_bot.compute();
- spline_soil.setPoints(&point[np],np);
- spline_soil.compute();
- spline_sat.setPoints(&point[2*np],np);
- spline_sat.compute();
- }
- void
- QtInputData::initGeometry(){
- geometry.lX=10;
- geometry.nX=100;
- geometry.dX=geometry.lX/(geometry.nX-1);
- geometry.hsoil=new double[nmax_Qt];
- geometry.dhsoil=new double[nmax_Qt];
- geometry.hbot=new double[nmax_Qt];
- geometry.dhbot=new double[nmax_Qt];
- geometry.nZ=new size_t[nmax_Qt];
- geometry.dZ=new double[nmax_Qt];
- geometry.Z=new double*[nmax_Qt];
- for(size_t i=0;i<nmax_Qt;++i){
- geometry.Z[i]=new double[nmax_Qt];
- }
- updateGeometry();
- }
- void
- QtInputData::initInitialState(){
- initial_state->hsat=new double[nmax_Qt];
- initial_state->hov=new double[nmax_Qt];
- initial_state->Pinit=new double*[nmax_Qt];
- for(size_t i=0;i<nmax_Qt;++i){
- initial_state->Pinit[i]=new double[nmax_Qt];
- }
- updateInitialState();
- }
- void
- QtInputData::updateGeometry(){
- double prev_factor=factor;
- double hs_max=-inf,hb_min=inf;
- double v;
- geometry.dX=geometry.lX/(geometry.nX-1);
- for(size_t k=0;k<geometry.nX;++k){
- double x=k*geometry.dX/geometry.lX;
- v=geometry.hsoil[k]=spline_soil(x);
- hs_max=max(v,hs_max);
- geometry.dhsoil[k]=spline_soil.derivate(x);
- v=geometry.hbot[k]=spline_bot(x);
- hb_min=min(v,hb_min);
- geometry.dhbot[k]=spline_bot.derivate(x);
- }
- double dZ_avg=depth/nZ_max;
- factor=depth/(hs_max-hb_min);
- for(size_t k=0;k<geometry.nX;++k){
- geometry.hsoil[k]*=factor;
- geometry.hbot[k]*=factor;
- }
- geometry.initZ(dZ_avg,false);
- double r=factor/prev_factor;
- for(auto it=initial_state->tanks.begin();it!=initial_state->tanks.end();++it){
- Tank* tank=*it;
- tank->top*=r;
- tank->bottom*=r;
- tank->delta_top*=r;
- tank->delta_bottom*=r;
- }
- for(auto it=source.pumps.begin();it!=source.pumps.end();++it){
- Pump* pump=*it;
- pump->top_init*=r;
- pump->bottom_init*=r;
- pump->delta_top_init*=r;
- pump->delta_bottom_init*=r;
- pump->top_final*=r;
- pump->bottom_final*=r;
- pump->delta_top_final*=r;
- pump->delta_bottom_final*=r;
- }
- }
- void
- QtInputData::updateInitialState(){
- for(size_t i=0;i<geometry.nX;++i){
- double x=double(i)/(geometry.nX-1);
- initial_state->hsat[i]=min(geometry.hsoil[i],spline_sat(x)*factor);
- }
- initial_state->updatePressure();
- updateOverland();
- updateOverlandAndPressure();
- }
- void
- QtInputData::save(fstream& file){
- InputData::save(file);
- file.write((char*)&depth,sizeof(double));
- file.write((char*)&nZ_max,sizeof(size_t));
- size_t temp=np;
- file.write((char*)&temp,sizeof(size_t));
- for(size_t i=0;i<3*np+np_ov;++i){
- file.write((char*)&point[i].x,sizeof(double));
- file.write((char*)&point[i].y,sizeof(double));
- }
- }
- void
- QtInputData::load(fstream& file){
- InputData::load(file);
- file.read((char*)&depth,sizeof(double));
- file.read((char*)&nZ_max,sizeof(size_t));
- size_t temp;
- file.read((char*)&temp,sizeof(size_t));
- if(temp!=np){
- cerr<<"Incorrect number of points"<<endl;
- return;
- }
- for(size_t i=0;i<3*np+np_ov;++i){
- file.read((char*)&point[i].x,sizeof(double));
- file.read((char*)&point[i].y,sizeof(double));
- }
- updateSplineBot();
- updateSplineSat();
- updateSplineSoil();
- updateGeometry();
- updateInitialState();
- }
- QtInputData::~QtInputData(){
- delete[] geometry.hsoil;
- delete[] geometry.dhsoil;
- delete[] geometry.hbot;
- delete[] geometry.dhbot;
- delete[] geometry.nZ;
- for(size_t i=0;i<nmax_Qt;++i){
- delete[] geometry.Z[i];
- }
- delete[] geometry.Z;
- geometry.hsoil=nullptr; //Bypass geometry destructor
- delete[] initial_state->hsat;
- for(size_t i=0;i<nmax_Qt;++i){
- delete[] initial_state->Pinit[i];
- }
- delete[] initial_state->Pinit;
- initial_state->hsat=nullptr; //Bypassfirst initial_state destructor
- }
- void
- QtInputData::updateOverland(){
- for(size_t x=0;x<geometry.nX;++x){
- initial_state->hov[x]=geometry.hsoil[x];
- }
- for(size_t i=0;i<np_ov;++i){
- Point& P=point[3*np+i];
- size_t x_ov=P.x*(geometry.nX-1);
- double y=P.y*factor;
- if(geometry.hsoil[x_ov]<y){
- size_t x_left=findSoilOver(y,x_ov,-1);
- size_t x_right=findSoilOver(y,x_ov,1);
- for(size_t x=x_left;x<=x_right;++x){
- initial_state->hov[x]=max(initial_state->hov[x],y);
- }
- }
- }
- }
- void
- QtInputData::updateOverlandAndPressure(){
- for(size_t x=0;x<geometry.nX;++x){
- double hsoil=geometry.hsoil[x];
- double hov=initial_state->hov[x];
- if(hov==hsoil){
- //No overland water and set hov to coincide with Psoil
- double Psoil=initial_state->Pinit[x][geometry.nZ[x]-1];
- double d=Physics::invPsoil(hsoil,Psoil);
- initial_state->hov[x]=Physics::invPsoil(hsoil,Psoil);
- }
- else{
- double Psoil=Physics::Psoil(hsoil,hov);
- double* Pinit_x=initial_state->Pinit[x];
- for(size_t z=0;z<geometry.nZ[x];++z){
- Pinit_x[z]=max(Pinit_x[z],20000*(geometry.Z[x][z]-geometry.hsoil[x])+Psoil);
- }
- //Overland water and set pressure in consequence
- }
- }
- }
- size_t
- QtInputData::findSoilOver(double y,int x_ov,int x_step){
- int x=x_ov;
- while(x>=0 and x<geometry.nX and geometry.hsoil[x]<y){
- x+=x_step;
- }
- x=min(x,(int)geometry.nX-1);
- x=max(x,0);
- return x;
- }
|