12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394 |
- #ifndef KERNEL_HPP
- #define KERNEL_HPP
- #include <cstdio>
- #include <cstring>
- #include <fstream>
- #include <stack>
- #include "physics.hpp"
- #include "geometry.hpp"
- #include "time.hpp"
- #include "input_data.hpp"
- #include "solution.hpp"
- #include "piccard.hpp"
- using namespace std;
- class Kernel:public InputData{
- private:
- static constexpr double tolerence_Piccard=1e-6;
- static constexpr double oscilation_Piccard=1e-4;
- static constexpr size_t max_iterations_Piccard=50;
- static constexpr size_t max_time_subdivisions=32;
- static constexpr size_t max_iterations_Richards=50;
- static constexpr double tolerence_Richards=1e-10;
- //For which colum we have to compute Richards
- bool* indice_x_Richards;
- int scheme_error;
- void initSolutions();
- void spaceSolution();
- double bottomFluxRichards(size_t ix,const Solution& S_k);
- Solution* Piccard(Solution* Sin);
- bool Richards1dEvolutiveTime(size_t ix,const Solution& S_in,const Solution& S_k,Solution& S_kp);
- bool allVerticalRichards(Solution& S_in,Solution& S_k,Solution& S_kp);
- void compute_l(Solution& S,double& error);
- void compute_Pl(Solution& S,const double* h,double** P);
- void compute_mass(Solution& S);
- void compute_div_w(size_t ix,const Solution& S);
- bool horizontalProblem(Solution& S,double& error);
- bool overlandProblem(const Solution& S_in,Solution& S_out,double& error);
- void saturationInterface(const Solution& S_in,Solution& S_out);
- bool Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot);
- bool weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0);
- void solveSystemRichards(size_t ix,double* P0,double* Pin,double* Pout,double flux_bot,double dt,size_t count);
- // double computeError(const Solution& S_in,const Solution& S_out,size_t sol_temp);
- public:
- size_t step;
- //Number of sub time intervals of [step,step+1]
- size_t m;
- Solution* solution;
- double** div_w;
- double** P_temp;
- double** t_sub_A;
- double** t_sub_B;
- double** t_sub_C;
- double** t_diag_A;
- double** t_diag_B;
- double** t_diag_C;
- double** t_sup_A;
- double** t_sup_B;
- double** t_sup_C;
- double** t_F;
- double** t_G;
- double** t_S;
- double** t_H;
- double** t_R;
- double** t_I;
- double** t_BB;
- size_t nZ(size_t ix);
- double Z(size_t ix,size_t iz);
- Kernel(string filename);
- void init(fstream& file);
- bool next();
- };
- inline size_t
- Kernel::nZ(size_t ix){
- return geometry.nZ[ix];
- }
- inline double
- Kernel::Z(size_t ix,size_t iz){
- return geometry.Z[ix][iz];
- }
- #endif
|