Commit 7a2eac2d authored by Coon, Ethan's avatar Coon, Ethan
Browse files

starts to namespace, rename files to match some convention. imports utils...

starts to namespace, rename files to match some convention.  imports utils from utils branch.  Builds first cut c-array based demo around BareGroundFluxes.  Adds to CMake to include new drivers
parent 9bb3b803
......@@ -4,10 +4,11 @@ set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
enable_testing()
project(ELM_physics CXX)
project(ELM CXX)
# need netcdf
list(APPEND CMAKE_MODULE_PATH "${ELM_physics_SOURCE_DIR}/config/cmake")
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/config/cmake")
include(AddImportedLibrary)
if (BUILD_SHARED_LIBS)
......
add_subdirectory(utils)
add_subdirectory(cc)
include_directories(${NetCDF_INCLUDE_DIR})
add_executable (ELM_driver driver.cc)
target_link_libraries (ELM_driver LINK_PUBLIC physics ${NetCDF_LIBRARIES})
install(TARGETS ELM_driver)
add_test(NAME ELM_driver COMMAND ELM_driver)
# include_directories(${NetCDF_INCLUDE_DIR})
# add_executable (ELM_driver driver.cc)
# target_link_libraries (ELM_driver LINK_PUBLIC physics ${NetCDF_LIBRARIES})
# install(TARGETS ELM_driver)
# add_test(NAME ELM_driver COMMAND ELM_driver)
include_directories(${ELM_UTILS_SOURCE_DIR})
add_executable (ELM_cc driver.cc)
target_link_libraries (ELM_cc LINK_PUBLIC elm_physics )
install(TARGETS ELM_cc)
add_test(NAME ELM_cc COMMAND ELM_cc)
#include "clm_constants.h"
#include "landtype.h"
#include "BareGroundFluxes.h"
#include "array.hh"
int main(int argc, char** argv)
{
const int ncells = 3;
const int ntimes = 2;
// instantiate data
ELM::LandType land;
ELM::Array<int,1> frac_vec_nosno_in(ncells);
ELM::Array<double,1> forc_u(ncells);
ELM::Array<double,1> forc_v(ncells);
ELM::Array<double,1> forc_q_in(ncells);
ELM::Array<double,1> forc_th_in(ncells);
ELM::Array<double,1> forc_hgt_u_patch_in(ncells);
ELM::Array<double,1> thm_in(ncells);
ELM::Array<double,1> thv_in(ncells);
ELM::Array<double,1> t_grnd(ncells);
ELM::Array<double,1> qg(ncells);
ELM::Array<double,1> z0mg_in(ncells);
ELM::Array<double,1> dlrad(ncells);
ELM::Array<double,1> ulrad(ncells);
ELM::Array<double,1> forc_hgt_t_patch(ncells);
ELM::Array<double,1> forc_hgt_q_patch(ncells);
ELM::Array<double,1> z0mg(ncells);
ELM::Array<double,1> zii(ncells);
ELM::Array<double,1> beta(ncells);
ELM::Array<double,1> z0hg(ncells);
ELM::Array<double,1> z0qg(ncells);
ELM::Array<int,1> snl(ncells, 1); // fix me!
ELM::Array<double,1> forc_rho(ncells);
ELM::Array<double,1> soilbeta(ncells);
ELM::Array<double,1> dqgdT(ncells);
ELM::Array<double,1> htvp(ncells);
ELM::Array<double,1> t_h2osfc(ncells);
ELM::Array<double,1> qg_snow(ncells);
ELM::Array<double,1> qg_soil(ncells);
ELM::Array<double,1> qg_h2osfc(ncells);
ELM::Array<double,2> t_soisno(ncells, ELM::nlevsno+1); // is this correct? This is required by BareGroundFluxes_impl.hh:108
ELM::Array<double,1> forc_pbot(ncells);
ELM::Array<double,1> cgrnds(ncells);
ELM::Array<double,1> cgrndl(ncells);
ELM::Array<double,1> cgrnd(ncells);
ELM::Array<double,1> eflx_sh_grnd(ncells);
ELM::Array<double,1> eflx_sh_tot(ncells);
ELM::Array<double,1> eflx_sh_snow(ncells);
ELM::Array<double,1> eflx_sh_soil(ncells);
ELM::Array<double,1> eflx_sh_h2osfc(ncells);
ELM::Array<double,1> qflx_evap_soi(ncells);
ELM::Array<double,1> qflx_evap_tot(ncells);
ELM::Array<double,1> qflx_ev_snow(ncells);
ELM::Array<double,1> qflx_ev_soil(ncells);
ELM::Array<double,1> qflx_ev_h2osfc(ncells);
ELM::Array<double,1> t_ref2m(ncells);
ELM::Array<double,1> t_ref2m_r(ncells);
ELM::Array<double,1> q_ref2m(ncells);
ELM::Array<double,1> rh_ref2m(ncells);
ELM::Array<double,1> rh_ref2m_r(ncells);
// initialize
ELM::Array<ELM::BareGroundFluxes, 1> bg_fluxes(ncells);
for (int c=0; c!=ncells; ++c) {
bg_fluxes[c].InitializeFlux(land,
frac_vec_nosno_in[c],
forc_u[c],
forc_v[c],
forc_q_in[c],
forc_th_in[c],
forc_hgt_u_patch_in[c],
thm_in[c],
thv_in[c],
t_grnd[c],
qg[c],
z0mg_in[c],
dlrad[c],
ulrad[c]);
}
// iterate in time
for (int time=0; time!=ntimes; ++time) {
for (int c=0; c!=ncells; ++c) {
bg_fluxes[c].StabilityIteration(land,
forc_hgt_t_patch[c],
forc_hgt_q_patch[c],
z0mg[c],
zii[c],
beta[c],
z0hg[c],
z0qg[c]);
bg_fluxes[c].ComputeFlux(land,
snl[c],
forc_rho[c],
soilbeta[c],
dqgdT[c],
htvp[c],
t_h2osfc[c],
qg_snow[c],
qg_soil[c],
qg_h2osfc[c],
t_soisno[c],
forc_pbot[c],
cgrnds[c],
cgrndl[c],
cgrnd[c],
eflx_sh_grnd[c],
eflx_sh_tot[c],
eflx_sh_snow[c],
eflx_sh_soil[c],
eflx_sh_h2osfc[c],
qflx_evap_soi[c],
qflx_evap_tot[c],
qflx_ev_snow[c],
qflx_ev_soil[c],
qflx_ev_h2osfc[c],
t_ref2m[c],
t_ref2m_r[c],
q_ref2m[c],
rh_ref2m[c],
rh_ref2m_r[c]);
}
}
return 0;
}
project(ELM_UTILS)
#target_include_directories (elm_utils PUBLIC ${ELM_UTILS_SOURCE_DIR})
//! A set of utilities for testing ELM kernels in C++
#ifndef ELM_KERNEL_TEST_ARRAY_HH_
#define ELM_KERNEL_TEST_ARRAY_HH_
#include <algorithm>
#include <array>
#include <numeric>
#include <assert.h>
#include <memory>
namespace ELM {
namespace Impl {
template<typename T>
class Data_ {
public:
// assigment operator
void operator=(T t) {
for (auto& d : *this) {
d = t;
}
}
// iterators
T const * begin() const { return &(d_[0]); }
T const * end() const { return &(d_[len_]); }
T * begin() { return &(d_[0]); }
T * end() { return &(d_[len_]); }
size_t size() const { return len_; }
T const * data() const { return d_; }
T * data() { return d_; }
protected:
//
// Constructors are all protected -- do not use this class directly!
//
// construct from a total length
Data_(int len)
: len_(len),
do_(std::shared_ptr<T>(new T[len], std::default_delete<T[]>()))
{
d_ = do_.get();
}
// construct and initialize
Data_(int len, T d)
: Data_(len)
{
*this = d;
}
// construct non-owning view
Data_(int len, T* d)
: len_(len),
d_(d) {}
// copy construct
Data_(const Data_& other) = default;
protected:
// global length
const int len_;
// owning data
std::shared_ptr<T> do_;
// non-owning data
T* d_;
};
} // namespace Impl
// templated class for multi-dimensional array
template<typename T, size_t D>
class Array : public Impl::Data_<T> {};
// 1D specialization
template<typename T>
class Array<T,1> : public Impl::Data_<T> {
static const size_t dim = 1;
public:
// forward construction
Array(int N)
: Impl::Data_<T>(N) {}
Array(std::array<int,1> N)
: Impl::Data_<T>(std::get<0>(N)) {}
// forward construction
Array(int N, T t)
: Impl::Data_<T>(N, t) {}
Array(std::array<int,1> N, T t)
: Impl::Data_<T>(std::get<0>(N), t) {}
// forward construction
Array(int N, T* d)
: Impl::Data_<T>(N, d) {}
Array(std::array<int,1> N, T* d)
: Impl::Data_<T>(std::get<0>(N), d) {}
// forward construction
Array(const Array<T,1>& other) = default;
// accessors
T& operator()(const int i) { assert(0 <= i && i < len_); return d_[i]; }
const T& operator()(const int i) const { assert(0 <= i && i < len_); return d_[i]; }
T& operator[](const int i) { assert(0 <= i && i < len_); return d_[i]; }
const T& operator[](int i) const { assert(0 <= i && i < len_); return d_[i]; }
// shape
int extent(int d) const {
assert(d < 1 && "Array::extent requested for dimension greater than is stored.");
return shape()[d];
}
std::array<int,1> shape() const {
return { len_ };
}
protected:
using Impl::Data_<T>::len_;
using Impl::Data_<T>::d_;
};
// 2D specialization
template<typename T>
class Array<T,2> : public Impl::Data_<T> {
static const size_t dim = 2;
public:
// forward construction
Array(int M, int N)
: Impl::Data_<T>(N*M),
M_(M),
N_(N)
{}
Array(std::array<int, 2> N)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N)),
M_(std::get<0>(N)),
N_(std::get<1>(N))
{}
// forward construction
Array(int M, int N, T t)
: Impl::Data_<T>(N*M, t),
M_(M),
N_(N)
{}
Array(std::array<int, 2> N, T t)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N), t),
M_(std::get<0>(N)),
N_(std::get<1>(N))
{}
// forward construction
Array(int M, int N, T* d)
: Impl::Data_<T>(N*M, d),
M_(M),
N_(N)
{}
Array(std::array<int, 2> N, T* d)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N), d),
M_(std::get<0>(N)),
N_(std::get<1>(N))
{}
// forward construction
Array(const Array<T,2>& other) = default;
T& operator()(int i, int j) { assert(0 <= i && i < M_ && 0 <= j && j < N_); return d_[j+i*N_]; }
const T& operator()(int i, int j) const { assert(0 <= i && i < M_ && 0 <= j && j < N_); return d_[j+i*N_]; }
Array<T,1> operator[](int i) { assert(0 <= i && i < M_); return Array<T,1>(N_, &d_[i*N_]); }
const Array<T,1> operator[](int i) const { assert(0 <= i && i < M_); return Array<T,1>(N_, &d_[i*N_]); }
// shape
int extent(size_t d) const {
assert(d < 2 && "Array::extent requested for dimension greater than is stored.");
return shape()[d];
}
std::array<int,2> shape() const {
return { M_, N_ };
}
protected:
int M_, N_;
using Impl::Data_<T>::d_;
};
// 3D specialization
template<typename T>
class Array<T,3> : public Impl::Data_<T> {
static const size_t dim = 3;
public:
// forward construction
Array(int M, int N, int P)
: Impl::Data_<T>(N*M*P),
M_(M),
N_(N),
P_(P)
{}
Array(std::array<int, 3> N)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N)),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N))
{}
// forward construction
Array(int M, int N, int P, T t)
: Impl::Data_<T>(N*M*P, t),
M_(M),
N_(N),
P_(P)
{}
Array(std::array<int, 3> N, T t)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N), t),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N))
{}
// forward construction
Array(int M, int N, int P, T* d)
: Impl::Data_<T>(N*M*P, d),
M_(M),
N_(N),
P_(P)
{}
Array(std::array<int, 3> N, T* d)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N), d),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N))
{}
// forward construction
Array(const Array<T,3>& other) = default;
T& operator()(int i, int j, int k) {
assert(0 <= i && i < M_ &&
0 <= j && j < N_ &&
0 <= k && k < P_);
return d_[k + P_*(j+i*N_)];
}
const T& operator()(int i, int j, int k) const {
assert(0 <= i && i < M_ &&
0 <= j && j < N_ &&
0 <= k && k < P_);
return d_[k + P_*(j+i*N_)];
}
Array<T,2> operator[](int i) {
assert(0 <= i && i < M_);
return Array<T,2>(N_, P_, &d_[i*N_*P_]);
}
const Array<T,2> operator[](int i) const {
assert(0 <= i && i < M_);
return Array<T,2>(N_, P_, &d_[i*N_*P_]);
}
// shape
int extent(int d) const {
assert(d < 3 && "Array::extent requested for dimension greater than is stored.");
return shape()[d];
}
std::array<int,3> shape() const {
return { M_, N_, P_ };
}
protected:
int M_, N_, P_;
using Impl::Data_<T>::d_;
};
// 4D specialization
template<typename T>
class Array<T,4> : public Impl::Data_<T> {
static const size_t dim = 4;
public:
// forward construction
Array(int M, int N, int P, int Q)
: Impl::Data_<T>(N*M*P*Q),
M_(M),
N_(N),
P_(P),
Q_(Q)
{}
Array(std::array<int, 4> N)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N) * std::get<3>(N)),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N)),
Q_(std::get<3>(N))
{}
// forward construction
Array(int M, int N, int P, int Q, T t)
: Impl::Data_<T>(N*M*P*Q, t),
M_(M),
N_(N),
P_(P),
Q_(Q)
{}
Array(std::array<int, 4> N, T t)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N) * std::get<3>(N), t),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N)),
Q_(std::get<3>(N))
{}
// forward construction
Array(int M, int N, int P, int Q, T* d)
: Impl::Data_<T>(N*M*P*Q, d),
M_(M),
N_(N),
P_(P),
Q_(Q)
{}
Array(std::array<int, 4> N, T* d)
: Impl::Data_<T>(std::get<0>(N) * std::get<1>(N) * std::get<2>(N) * std::get<3>(N), d),
M_(std::get<0>(N)),
N_(std::get<1>(N)),
P_(std::get<2>(N)),
Q_(std::get<3>(N))
{}
// forward construction
Array(const Array<T,4>& other) = default;
T& operator()(int i, int j, int k, int l) {
assert(0 <= i && i < M_ &&
0 <= j && j < N_ &&
0 <= k && k < P_ &&
0 <= l && l < Q_);
return d_[l + Q_*(k + P_*(j+i*N_))];
}
const T& operator()(int i, int j, int k, int l) const {
assert(0 <= i && i < M_ &&
0 <= j && j < N_ &&
0 <= k && k < P_ &&
0 <= l && l < Q_);
return d_[l + Q_*(k + P_*(j+i*N_))];
}
Array<T,3> operator[](int i) {
assert(0 <= i && i < M_);
return Array<T,3>(N_, P_, Q_, &d_[i*N_*P_*Q_]);
}
const Array<T,3> operator[](int i) const {
assert(0 <= i && i < M_);
return Array<T,3>(N_, P_, Q_, &d_[i*N_*P_*Q_]);
}
// shape
int extent(int d) const {
assert(d < 4 && "Array::extent requested for dimension greater than is stored.");
return shape()[d];
}
std::array<int,4> shape() const {
return { M_, N_, P_, Q_ };
}
protected:
int M_, N_, P_, Q_;
using Impl::Data_<T>::d_;
};
//
// Construct a non-owning Array of a different shape but with the same data as
// an old shape.
//
// NOTE, this does no transposing!
//
template<size_t D1, size_t D2, typename T>
Array<T,D2> reshape(Array<T,D1>& arr_in, const std::array<int, D2>& new_shape) {
int new_length = std::accumulate(new_shape.begin(), new_shape.end(), 1, std::multiplies<int>());
assert(new_length == arr_in.size() && "Invalid Array reshape");
// construct and return
return Array<T,D2>(new_shape, (T*) arr_in.begin());
}
//
// Copies from an Array into an Array-like object.
//
template<typename T, typename Array_type>
void
deep_copy(Array_type& arr, const Array<T,1>& arr_in)
{
assert(arr.extent(0) == arr_in.extent(0));
for (int i=0; i!=arr_in.extent(0); ++i) {
arr(i) = arr_in(i);
}
}
template<typename T, typename Array_type>
void