Skip to content

Commit cc0f137

Browse files
pdziekanpdziekan
committed
pbl with anisotropic smagorinsky test
1 parent bbeb5c5 commit cc0f137

File tree

6 files changed

+86
-35
lines changed

6 files changed

+86
-35
lines changed

‎libmpdata++-config.cmake‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ set(libmpdataxx_INCLUDE_DIRS "${CMAKE_CURRENT_LIST_DIR}/../../include/")
3232

3333
############################################################################################
3434
# debug mode compiler flags
35-
set(libmpdataxx_CXX_FLAGS_DEBUG "${libmpdataxx_CXX_FLAGS_DEBUG} -std=c++14 -DBZ_DEBUG -g -Wno-enum-compare -Wfatal-errors") #TODO: -Og if compiler supports it?
35+
set(libmpdataxx_CXX_FLAGS_DEBUG "${libmpdataxx_CXX_FLAGS_DEBUG} -std=c++17 -DBZ_DEBUG -g -Wno-enum-compare -Wfatal-errors") #TODO: -Og if compiler supports it?
3636

3737

3838
############################################################################################
@@ -42,7 +42,7 @@ if(
4242
CMAKE_CXX_COMPILER_ID STREQUAL "Clang" OR
4343
CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang"
4444
)
45-
set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=c++14 -DNDEBUG -Ofast -march=native -Wno-enum-compare")
45+
set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=c++17 -DNDEBUG -Ofast -march=native -Wno-enum-compare")
4646

4747
# preventing Kahan summation from being optimised out
4848
if (
@@ -58,23 +58,23 @@ if(
5858
CMAKE_CXX_COMPILER_ID STREQUAL "Intel"
5959
)
6060
# flags taken from -fast but without -static
61-
set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=gnu++14 -DNDEBUG -xHOST -O3 -ipo -no-prec-div -fp-model fast=2")
61+
set(libmpdataxx_CXX_FLAGS_RELEASE "${libmpdataxx_CXX_FLAGS_RELEASE} -std=gnu++17 -DNDEBUG -xHOST -O3 -ipo -no-prec-div -fp-model fast=2")
6262
endif()
6363

6464

6565
############################################################################################
66-
# C++14
66+
# C++17
6767
include(CheckCXXSourceCompiles)
68-
set(CMAKE_REQUIRED_FLAGS "-std=c++14")
68+
set(CMAKE_REQUIRED_FLAGS "-std=c++17")
6969
check_cxx_source_compiles("
7070
#include <type_traits>
7171
auto f() { return 1;}
7272
template <bool a, class b> using ei=std::enable_if<a,b>;
7373
struct a {a(int){}};struct b:a {using a::a;};
74-
int main(){b i(1);}
75-
" CXX14_SUPPORTED)
76-
if (NOT CXX14_SUPPORTED)
77-
message(FATAL_ERROR "C++14 compatibility test failed - please update your compiler or point CMake to another one with -DCMAKE_CXX_COMPILER=...")
74+
int main(){b i(1); if constexpr(true){;}}
75+
" CXX17_SUPPORTED)
76+
if (NOT CXX17_SUPPORTED)
77+
message(FATAL_ERROR "C++17 compatibility test failed - please update your compiler or point CMake to another one with -DCMAKE_CXX_COMPILER=...")
7878
endif()
7979
unset(CMAKE_REQUIRED_FLAGS)
8080

‎libmpdata++/solvers/detail/boussinesq_sgs_common.hpp‎

Lines changed: 43 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -22,9 +22,10 @@ namespace libmpdataxx
2222
using real_t = typename ct_params_t::real_t;
2323

2424
protected:
25+
2526
// member fields
2627
real_t prandtl_num;
27-
typename parent_t::arr_t &rcdsn_num, &full_tht, &tdef_sq, &mix_len, &hflux_srfc;
28+
typename parent_t::arr_t &rcdsn_num, &full_tht, &tdef_sq, &mix_len, &mix_len_h, &mix_len_v, &hflux_srfc;
2829
arrvec_t<typename parent_t::arr_t> &grad_tht;
2930

3031
template <int nd = ct_params_t::n_dims>
@@ -74,20 +75,48 @@ namespace libmpdataxx
7475

7576
calc_rcdsn_num();
7677

77-
this->k_m(this->ijk) = where(
78+
if constexpr (static_cast<sgs_scheme_t>(ct_params_t::sgs_scheme) == smg)
79+
this->k_m(this->ijk) = where(
7880
rcdsn_num(this->ijk) / prandtl_num < 1,
7981
pow(this->smg_c * mix_len(this->ijk), 2)
8082
* sqrt(tdef_sq(this->ijk) * (1 - rcdsn_num(this->ijk) / prandtl_num)),
8183
0
8284
);
85+
else // smgani
86+
{
87+
this->k_m[0](this->ijk) = where(
88+
rcdsn_num(this->ijk) / prandtl_num < 1,
89+
pow(this->smg_c * mix_len_h(this->ijk), 2)
90+
* sqrt(tdef_sq(this->ijk) * (1 - rcdsn_num(this->ijk) / prandtl_num)),
91+
0
92+
);
93+
this->k_m[1](this->ijk) = where(
94+
rcdsn_num(this->ijk) / prandtl_num < 1,
95+
pow(this->smg_c * mix_len_v(this->ijk), 2)
96+
* sqrt(tdef_sq(this->ijk) * (1 - rcdsn_num(this->ijk) / prandtl_num)),
97+
0
98+
);
99+
}
83100
// one level above surface
84101
auto ijp1 = this->ijk;
85102
ijp1.lbound(ct_params_t::n_dims - 1) = 1;
86103
ijp1.ubound(ct_params_t::n_dims - 1) = 1;
87104

88-
this->k_m(ij) = this->k_m(ijp1);
105+
if constexpr (static_cast<sgs_scheme_t>(ct_params_t::sgs_scheme) == smg)
106+
this->k_m(ij) = this->k_m(ijp1);
107+
else
108+
{
109+
this->k_m[0](ij) = this->k_m[0](ijp1);
110+
this->k_m[1](ij) = this->k_m[1](ijp1);
111+
}
89112

90-
this->xchng_sclr(this->k_m, this->ijk, 1);
113+
if constexpr (static_cast<sgs_scheme_t>(ct_params_t::sgs_scheme) == smg)
114+
this->xchng_sclr(this->k_m, this->ijk, 1);
115+
else
116+
{
117+
this->xchng_sclr(this->k_m[0], this->ijk, 1);
118+
this->xchng_sclr(this->k_m[1], this->ijk, 1);
119+
}
91120

92121
// havo to use modified ijkm due to shared-memory parallelisation, otherwise overlapping ranges
93122
// would lead to double multiplications
@@ -129,10 +158,12 @@ namespace libmpdataxx
129158
prandtl_num(p.prandtl_num),
130159
rcdsn_num(args.mem->tmp[__FILE__][0][0]),
131160
full_tht(args.mem->tmp[__FILE__][1][0]),
132-
tdef_sq(args.mem->tmp[__FILE__][2][0]),
133-
mix_len(args.mem->tmp[__FILE__][3][0]),
134-
grad_tht(args.mem->tmp[__FILE__][4]),
135-
hflux_srfc(args.mem->tmp[__FILE__][5][0])
161+
mix_len(args.mem->tmp[__FILE__][2][0]),
162+
mix_len_h(args.mem->tmp[__FILE__][3][0]),
163+
mix_len_v(args.mem->tmp[__FILE__][4][0]),
164+
grad_tht(args.mem->tmp[__FILE__][5]),
165+
hflux_srfc(args.mem->tmp[__FILE__][6][0]),
166+
tdef_sq(args.mem->tmp[__FILE__][7][0])
136167
{
137168
//assert(prandtl_num != 0);
138169
}
@@ -142,10 +173,12 @@ namespace libmpdataxx
142173
parent_t::alloc(mem, n_iters);
143174
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // rcdsn_num
144175
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // full_tht
145-
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // tdef_sq
146-
parent_t::alloc_tmp_sclr(mem, __FILE__, 1, "mix_len");
176+
parent_t::alloc_tmp_sclr(mem, __FILE__, 1, "mix_len"); // used in smg
177+
parent_t::alloc_tmp_sclr(mem, __FILE__, 1, "mix_len_h"); // used in smganiso, needs to be named; TODO: separate smg and smagni into two classes to avoid allocating too much memory
178+
parent_t::alloc_tmp_sclr(mem, __FILE__, 1, "mix_len_v"); // ditto
147179
parent_t::alloc_tmp_vctr(mem, __FILE__); // grad_tht
148180
parent_t::alloc_tmp_sclr(mem, __FILE__, 1, "", true); // hflux_srfc
181+
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // tdef_sq
149182
}
150183
};
151184
} // namespace detail

‎libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_smgani.hpp‎

Lines changed: 4 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ namespace libmpdataxx
2626
protected:
2727

2828
real_t smg_c, c_m;
29-
typename parent_t::arr_t &k_m, &tdef_sq;
29+
arrvec_t<typename parent_t::arr_t> &k_m;
3030

3131
void multiply_sgs_visc()
3232
{
@@ -36,11 +36,9 @@ namespace libmpdataxx
3636
const auto dlta_h = std::accumulate(this->dijk.begin(), this->dijk.end()-1, real_t(0.)) / (ct_params_t::n_dims-1);
3737
const auto dlta_v = this->dijk[ct_params_t::n_dims-1];
3838

39-
tdef_sq = formulae::stress::calc_tdef_sq_cmpct<ct_params_t::n_dims>(this->tau, this->ijk);
40-
4139
// Simon and Chow 2021, eqs. 9 and 10
42-
k_m[0](this->ijk) = pow(smg_c * dlta_h, 2) * tdef_sq;
43-
k_m[0](this->ijk) = pow(smg_c * dlta_v, 2) * tdef_sq;
40+
k_m[0](this->ijk) = pow(smg_c * dlta_h, 2) * formulae::stress::calc_tdef_sq_cmpct<ct_params_t::n_dims>(this->tau, this->ijk)(this->ijk); // tdef_sq could be cached, but creating a tdef_sq array without similar one in isotropic smg messes with derived classes that do cache tdef_sq also in isotropic... (e.g. boussinesq)
41+
k_m[0](this->ijk) = pow(smg_c * dlta_v, 2) * formulae::stress::calc_tdef_sq_cmpct<ct_params_t::n_dims>(this->tau, this->ijk)(this->ijk);
4442

4543
formulae::stress::multiply_tnsr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(this->tau,
4644
real_t(1.0),
@@ -67,8 +65,7 @@ namespace libmpdataxx
6765
parent_t(args, p),
6866
smg_c(p.smg_c),
6967
c_m(p.c_m),
70-
k_m(args.mem->tmp[__FILE__][0]),
71-
tdef_sq(args.mem->tmp[__FILE__][1][0])
68+
k_m(args.mem->tmp[__FILE__][0])
7269
{
7370
if (smg_c == 0) throw std::runtime_error("libmpdata++: smg_c == 0");
7471
}
@@ -79,7 +76,6 @@ namespace libmpdataxx
7976
) {
8077
parent_t::alloc(mem, n_iters);
8178
parent_t::alloc_tmp_sclr(mem, __FILE__, 2); // k_m
82-
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // tdef_sq
8379
}
8480
};
8581
} // namespace detail

‎tests/sandbox/pbl/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ add_test(pbl_smg_budget bash -c "
1616
")
1717

1818
libmpdataxx_add_test(pbl_smg_short)
19+
libmpdataxx_add_test(pbl_smgani_short)
1920
libmpdataxx_add_test(pbl_iles_short)
2021

2122
add_test(pbl_iles_short_profiles bash -c "
@@ -24,6 +25,9 @@ add_test(pbl_iles_short_profiles bash -c "
2425
add_test(pbl_smg_short_profiles bash -c "
2526
python3 ${CMAKE_CURRENT_SOURCE_DIR}/profiles.py out_pbl_smg_short
2627
")
28+
add_test(pbl_smgani_short_profiles bash -c "
29+
python3 ${CMAKE_CURRENT_SOURCE_DIR}/profiles.py out_pbl_smgani_short
30+
")
2731

2832
add_test(pbl_iles_short_prof_diff
2933
zdiff ${CMAKE_CURRENT_SOURCE_DIR}/refdata/profiles_pbl_iles_short.txt.gz profiles_pbl_iles_short.txt)
@@ -36,6 +40,9 @@ add_test(pbl_iles_short_budget bash -c "
3640
add_test(pbl_smg_short_budget bash -c "
3741
python3 ${CMAKE_CURRENT_SOURCE_DIR}/budget.py out_pbl_smg_short
3842
")
43+
add_test(pbl_smgani_short_budget bash -c "
44+
python3 ${CMAKE_CURRENT_SOURCE_DIR}/budget.py out_pbl_smgani_short
45+
")
3946

4047
add_test(pbl_iles_short_budget_diff
4148
h5diff -v2 --relative=1e-2 ${CMAKE_CURRENT_SOURCE_DIR}/refdata/budget_pbl_iles_short.h5 budget_pbl_iles_short.h5)

‎tests/sandbox/pbl/pbl.hpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,11 @@ class pbl : public libmpdataxx::output::hdf5_xdmf<libmpdataxx::solvers::boussine
2828
{
2929
parent_t::multiply_sgs_visc();
3030

31-
if (this->timestep % static_cast<int>(this->outfreq) == 0 &&
32-
static_cast<libmpdataxx::solvers::sgs_scheme_t>(ct_params_t::sgs_scheme) == libmpdataxx::solvers::smg)
33-
{
34-
tke(this->ijk) = pow2(this->k_m(this->ijk) / (this->c_m * this->mix_len(this->ijk)));
35-
}
31+
if constexpr(static_cast<libmpdataxx::solvers::sgs_scheme_t>(ct_params_t::sgs_scheme) == libmpdataxx::solvers::smg)
32+
if (this->timestep % static_cast<int>(this->outfreq) == 0)
33+
{
34+
tke(this->ijk) = pow2(this->k_m(this->ijk) / (this->c_m * this->mix_len(this->ijk)));
35+
}
3636
}
3737

3838
void vip_rhs_expl_calc()

‎tests/sandbox/pbl/pbl_test_def.hpp‎

Lines changed: 18 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
using namespace libmpdataxx;
1414

1515
struct smg_tag {};
16+
struct smgani_tag {};
1617
struct iles_tag {};
1718

1819
// set parameters specific to the sgs model
@@ -34,6 +35,12 @@ void set_sgs_specific(params_t &p, smg_tag)
3435
p.cdrag = 0.1;
3536
}
3637

38+
template <typename params_t>
39+
void set_sgs_specific(params_t &p, smgani_tag)
40+
{
41+
set_sgs_specific(p, smg_tag{});
42+
}
43+
3744
template <typename sgs_t>
3845
void test(const std::string &dirname, const int np, const int nt)
3946
{
@@ -48,7 +55,7 @@ void test(const std::string &dirname, const int np, const int nt)
4855
enum { vip_vab = solvers::impl };
4956
enum { prs_scheme = solvers::cr };
5057
enum { stress_diff = solvers::compact };
51-
enum { sgs_scheme = std::is_same<sgs_t, smg_tag>::value ? solvers::smg : solvers::iles};
58+
enum { sgs_scheme = std::is_same<sgs_t, smg_tag>::value ? solvers::smg : std::is_same<sgs_t, smgani_tag>::value ? solvers::smgani : solvers::iles};
5259
enum { impl_tht = true };
5360
struct ix { enum {
5461
u, v, w, tht,
@@ -154,17 +161,25 @@ void test(const std::string &dirname, const int np, const int nt)
154161
slv.vab_relaxed_state(1) = 0;
155162
slv.vab_relaxed_state(2) = 0;
156163

157-
if (std::is_same<sgs_t, iles_tag>::value)
164+
if constexpr (std::is_same<sgs_t, iles_tag>::value)
158165
{
159166
// iles prescribed heat flux forcing
160167
slv.sclr_array("hflux_frc") = 0.01 * 1. / p.hscale * exp(- k * p.dk / p.hscale);
161168
}
162-
else
169+
else if constexpr (std::is_same<sgs_t, smg_tag>::value)
163170
{
164171
auto dlta = (p.di + p.dj + p.dk) / 3;
165172
// smagorinsky mixing length
166173
slv.sclr_array("mix_len") = min(max(k, 1) * p.dk * 0.845, dlta);
167174
}
175+
else if constexpr (std::is_same<sgs_t, smgani_tag>::value)
176+
{
177+
// anisotropic smagorinsky mixing length
178+
auto dlta = (p.di + p.dj) / 2;
179+
slv.sclr_array("mix_len_h") = min(max(k, 1) * p.dk * 0.845, dlta);
180+
dlta = p.dk;
181+
slv.sclr_array("mix_len_v") = min(max(k, 1) * p.dk * 0.845, dlta);
182+
}
168183
}
169184
}
170185

0 commit comments

Comments
 (0)