Commit 04d1abce authored by Coon, Ethan's avatar Coon, Ethan
Browse files

work toward a dynamically allocated but fixed size array

parent e33ffbe9
......@@ -39,11 +39,6 @@ static const int NUM_PROCS = 6;
static const int NUM_PROCS_X = 2;
static const int NUM_PROCS_Y = 3;
using MatrixStatePFT = MatrixStatic<n_grid_cells, n_pfts>;
using MatrixStateSoilColumn = MatrixStatic<n_grid_cells, n_levels_snow>;
using MatrixForc = MatrixStatic<n_max_times,n_grid_cells>;
using VectorColumn = VectorStatic<n_grid_cells>;
using VectorColumnInt = VectorStatic<n_grid_cells,int>;
} // namespace
} // namespace
......@@ -86,17 +81,17 @@ int main(int argc, char ** argv)
const double n_melt = 0.7;
// phenology input
ELM::Utils::MatrixStatePFT elai;
ELM::Utils::MatrixStatePFT esai;
ELM::Utils::Matrix<double> elai(n_grid_cells, n_pfts);
ELM::Utils::Matrix<double> esai(n_grid_cells, n_pfts);
ELM::Utils::read_phenology("../links/surfacedataWBW.nc", n_months, n_pfts, 0, elai, esai);
ELM::Utils::read_phenology("../links/surfacedataBRW.nc", n_months, n_pfts, n_months, elai, esai);
// forcing input
ELM::Utils::MatrixForc forc_rain;
ELM::Utils::MatrixForc forc_snow;
ELM::Utils::MatrixForc forc_air_temp;
ELM::Utils::Matrix<double> forc_rain(n_max_times, n_grid_cells);
ELM::Utils::Matrix<double> forc_snow(n_max_times, n_grid_cells);
ELM::Utils::Matrix<double> forc_air_temp(n_max_times, n_grid_cells);
const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 0, n_grid_cells, forc_rain, forc_snow, forc_air_temp);
ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
ELM::Utils::Matrix<double> forc_irrig(n_max_times, n_grid_cells, 0.);
double qflx_floodg = 0.0;
......@@ -104,47 +99,47 @@ int main(int argc, char ** argv)
//
// NOTE: in a real case, these would be populated, but we don't actually
// need them to be for these kernels. --etc
auto z = ELM::Utils::MatrixStateSoilColumn(0.);
auto zi = ELM::Utils::MatrixStateSoilColumn(0.);
auto dz = ELM::Utils::MatrixStateSoilColumn(0.);
auto z = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto zi = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto dz = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
// state variables that require ICs and evolve (in/out)
auto h2ocan = ELM::Utils::MatrixStatePFT(); h2ocan = 0.;
auto swe_old = ELM::Utils::MatrixStateSoilColumn(0.);
auto h2osoi_liq = ELM::Utils::MatrixStateSoilColumn(0.);
auto h2osoi_ice = ELM::Utils::MatrixStateSoilColumn(0.);
auto t_soisno = ELM::Utils::MatrixStateSoilColumn(0.);
auto frac_iceold = ELM::Utils::MatrixStateSoilColumn(0.);
auto t_grnd = ELM::Utils::VectorColumn(0.);
auto h2osno = ELM::Utils::VectorColumn(0.); h2osno = 0.;
auto snow_depth = ELM::Utils::VectorColumn(0.);
auto snow_level = ELM::Utils::VectorColumnInt(0.); // note this tracks the snow_depth
auto h2osfc = ELM::Utils::VectorColumn(0.);
auto frac_h2osfc = ELM::Utils::VectorColumn(0.); frac_h2osfc = 0.;
auto h2ocan = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto swe_old = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto h2osoi_liq = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto h2osoi_ice = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto t_soisno = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto frac_iceold = ELM::Utils::Matrix<double>(n_grid_cells, n_levels_snow, 0.);
auto t_grnd = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto h2osno = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto snow_depth = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto snow_level = ELM::Utils::Vector<int>(n_grid_cells, 0); // note this tracks the snow_depth
auto h2osfc = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto frac_h2osfc = ELM::Utils::Vector<double>(n_grid_cells, 0.);
// output fluxes by pft
auto qflx_prec_intr = ELM::Utils::MatrixStatePFT();
auto qflx_irrig = ELM::Utils::MatrixStatePFT();
auto qflx_prec_grnd = ELM::Utils::MatrixStatePFT();
auto qflx_snwcp_liq = ELM::Utils::MatrixStatePFT();
auto qflx_snwcp_ice = ELM::Utils::MatrixStatePFT();
auto qflx_snow_grnd_patch = ELM::Utils::MatrixStatePFT();
auto qflx_rain_grnd = ELM::Utils::MatrixStatePFT();
auto qflx_prec_intr = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_irrig = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_prec_grnd = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_snwcp_liq = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_snwcp_ice = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_snow_grnd_patch = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
auto qflx_rain_grnd = ELM::Utils::Matrix<double>(n_grid_cells, n_pfts, 0.);
// FIXME: I have no clue what this is... it is inout on WaterSnow. For now I
// am guessing the data structure. Ask Scott. --etc
auto integrated_snow = ELM::Utils::VectorColumn(0.);
auto integrated_snow = ELM::Utils::Vector<double>(n_grid_cells, 0.);
// output fluxes, state by the column
auto qflx_snow_grnd_col = ELM::Utils::VectorColumn();
auto qflx_snow_h2osfc = ELM::Utils::VectorColumn();
auto qflx_h2osfc2topsoi = ELM::Utils::VectorColumn();
auto qflx_floodc = ELM::Utils::VectorColumn();
auto qflx_snow_grnd_col = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto qflx_snow_h2osfc = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto qflx_h2osfc2topsoi = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto qflx_floodc = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto frac_sno_eff = ELM::Utils::VectorColumn();
auto frac_sno = ELM::Utils::VectorColumn();
auto frac_sno_eff = ELM::Utils::Vector<double>(n_grid_cells, 0.);
auto frac_sno = ELM::Utils::Vector<double>(n_grid_cells, 0.);
#ifdef TRACE
std::ofstream soln_file;
soln_file.open("test_CanopyHydrology_module.soln");
......
......@@ -11,28 +11,26 @@ namespace ELM {
namespace Utils {
template<size_t N, typename T=double>
template<typename T=double>
class VectorStatic {
public:
VectorStatic() :
do_(std::shared_ptr<std::array<T,N>>(new std::array<T,N>())) {
VectorStatic(std::size_t N) :
do_(std::shared_ptr<T>(new T[N])),
N_(N) {
d_ = do_.get();
}
VectorStatic(T t) :
VectorStatic()
VectorStatic(std::size_t N, T t) :
VectorStatic(N)
{ *this = t; }
VectorStatic(const std::shared_ptr<std::array<T,N> >& d) :
do_(d) {
d_ = do_.get();
}
VectorStatic(std::array<T,N>* d) :
VectorStatic(std::size_t N, T* d) :
do_(nullptr),
d_(d)
d_(d),
N_(N)
{}
VectorStatic(const VectorStatic<N,T>& other) = default;
VectorStatic(const VectorStatic& other) = default;
T& operator()(size_t i) { assert(0 <= i && i < N); return (*d_)[i]; }
const T& operator()(size_t i) const { assert(0 <= i && i < N); return (*d_)[i]; }
......@@ -41,251 +39,58 @@ class VectorStatic {
const T& operator[](size_t i) const { assert(0 <= i && i < N); return (*d_)[i]; }
void operator=(T t) {
for (size_t i=0; i!=N; ++i) {
for (size_t i=0; i!=N_; ++i) {
(*d_)[i] = t;
}
}
double const * begin() const { return &((*d_)[0]); }
double const * end() const { return &((*d_)[N]); }
std::size_t size() const { return N_; }
private:
std::array<T,N>* d_;
std::shared_ptr< std::array<T,N> > do_;
const std::size_t N_;
T[] d_;
std::shared_ptr<T[]> do_;
};
template<size_t ROW, size_t COL, typename T=double>
template<typename T=double>
class MatrixStatic {
using Data_t = std::array<std::array<T,COL>, ROW>;
public:
MatrixStatic() :
d_(std::shared_ptr<Data_t>(new Data_t())) {}
MatrixStatic(std::size_t M, std::size_t N) :
d_(std::shared_ptr<T>(new T[N*M])),
M_(M),
N_(N),
NM_(N*M)
{}
MatrixStatic(T t) :
MatrixStatic()
MatrixStatic(std::size_t M, std::size_t N, T t) :
MatrixStatic(M,N)
{ *this = t; }
MatrixStatic(const MatrixStatic<ROW,COL,T>& other) = default;
MatrixStatic(const MatrixStatic& other) = default;
T& operator()(size_t i, size_t j) { assert(0 <= i && i < ROW && 0 <= j && j < COL); return (*d_)[i][j]; }
const T& operator()(size_t i, size_t j) const { assert(0 <= i && i < ROW && 0 <= j && j < COL); return (*d_)[i][j]; }
T& operator()(size_t i, size_t j) { assert(0 <= i && i < M_ && 0 <= j && j < N_); return (*d_)[j+i*N_]; }
const T& operator()(size_t i, size_t j) const { assert(0 <= i && i < M_ && 0 <= j && j < N_); return (*d_)[j+i*N_]; }
VectorStatic<COL,T> operator[](size_t i) { assert(0 <= i && i < ROW); return &((*d_)[i]); }
const VectorStatic<COL,T> operator[](size_t i) const { assert(0 <= i && i < ROW); return &((*d_)[i]); }
VectorStatic<T> operator[](size_t i) { assert(0 <= i && i < M_); return VectorStatic<T>(N_, &(*d_)[i*N_]); }
const VectorStatic<COL,T> operator[](size_t i) const { assert(0 <= i && i < M_); return VectorStatic<T>(N_, &(*d_)[i*N_]); }
void operator=(T t) {
for (size_t i=0; i!=ROW; ++i) {
for (size_t j=0; j!=COL; ++j) {
(*d_)[i][j] = t;
}
}
}
double const * begin() const { return &(*d_)[0][0]; }
double const * end() const { return &(*d_)[ROW-1][COL-1] +1; }
private:
std::shared_ptr<Data_t> d_;
};
enum struct Ordering { C, FORTRAN };
template <typename T=double, Ordering O=Ordering::C>
class Matrix {
public:
Matrix(size_t nrows, size_t ncols) :
nrows_(nrows),
ncols_(ncols),
d_(nrows*ncols) {}
inline
const T& operator()(size_t i, size_t j) const {
return d_[ncols_*i + j];
}
inline
T& operator()(size_t i, size_t j) {
return d_[ncols_*i + j];
}
void operator=(T val) {
for (size_t lcv=0; lcv!=nrows_*ncols_; ++lcv) {
d_[lcv] = val;
}
}
typename std::vector<T>::const_iterator begin() const {
return d_.begin();
}
typename std::vector<T>::const_iterator end() const {
return d_.end();
}
/*
inline
T* operator[](size_t i) {
return &d_[i*ncols];
}
*/
private:
size_t nrows_, ncols_;
std::vector<T> d_;
};
template <typename T>
class Matrix<T,Ordering::FORTRAN> {
public:
Matrix(size_t nrows, size_t ncols) :
nrows_(nrows),
ncols_(ncols),
d_(nrows*ncols) {}
inline
const T& operator()(size_t i, size_t j) const {
return d_[nrows_*j + i];
}
inline
T& operator()(size_t i, size_t j) {
return d_[nrows_*j + i];
}
/*
inline
T* operator[](size_t i) {
return &d_[i*ncols];
}
*/
void operator=(T val) {
for (size_t lcv=0; lcv!=nrows_*ncols_; ++lcv) {
d_[lcv] = val;
}
}
typename std::vector<T>::const_iterator begin() const {
return d_.begin();
}
typename std::vector<T>::const_iterator end() const {
return d_.end();
}
private:
size_t nrows_, ncols_;
std::vector<T> d_;
};
template <typename T=double, Ordering O=Ordering::C>
class TensorRank3 {
public:
TensorRank3(size_t n_first, size_t n_second, size_t n_third) :
n_first_(n_first),
n_second_(n_second),
n_third_(n_third),
d_(n_first*n_second*n_third) {}
inline
const T& operator()(size_t i, size_t j, size_t k) const {
return d_[n_third_ * (n_second_*i + j) + k];
}
inline
T& operator()(size_t i, size_t j, size_t k) {
return d_[n_third_ * (n_second_*i + j) + k];
}
/*
inline
T* operator[](size_t i) {
return &d_[i*n_second_*n_third_];
}
*/
void operator=(T val) {
for (size_t lcv=0; lcv!=n_first_*n_second_*n_third_; ++lcv) {
d_[lcv] = val;
for (std::size_t i=0; i!=NM_; ++i) {
(*d_)[i] = t;
}
}
typename std::vector<T>::const_iterator begin() const {
return d_.begin();
}
typename std::vector<T>::const_iterator end() const {
return d_.end();
}
double const * begin() const { return &(*d_)[0]; }
double const * end() const { return &(*d_)[M_-1][N_-1] +1; }
private:
size_t n_first_, n_second_, n_third_;
std::vector<T> d_;
std::shared_ptr<T[]> d_;
};
template <typename T>
class TensorRank3<T,Ordering::FORTRAN> {
public:
TensorRank3(size_t n_first, size_t n_second, size_t n_third) :
n_first_(n_first),
n_second_(n_second),
n_third_(n_third),
d_(n_first*n_second*n_third) {}
inline
const T& operator()(size_t i, size_t j, size_t k) const {
return d_[n_first_*(n_second_*k + j) + i];
}
inline
T& operator()(size_t i, size_t j, size_t k) {
return d_[n_first_*(n_second_*k + j) + i];
}
/*
inline
T* operator[](size_t k) {
return &d_[k*n_first_*n_second_];
}
*/
void operator=(T val) {
for (size_t lcv=0; lcv!=n_first_*n_second_*n_third_; ++lcv) {
d_[lcv] = val;
}
}
typename std::vector<T>::const_iterator begin() const {
return d_.begin();
}
typename std::vector<T>::const_iterator end() const {
return d_.end();
}
private:
size_t n_first_, n_second_, n_third_;
std::vector<T> d_;
};
} // namespace Utils
} // namespace ELM
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment