/
mpdata_rhs_vip_prs_sgs_smg.hpp
89 lines (74 loc) · 2.83 KB
/
mpdata_rhs_vip_prs_sgs_smg.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
/**
* @file
* @copyright University of Warsaw
* @section LICENSE
* GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
*/
#pragma once
#include <libmpdata++/solvers/detail/mpdata_rhs_vip_prs_sgs_common.hpp>
namespace libmpdataxx
{
namespace solvers
{
namespace detail
{
template <class ct_params_t, int minhalo>
class mpdata_rhs_vip_prs_sgs_smg : public detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>
{
using parent_t = detail::mpdata_rhs_vip_prs_sgs_common<ct_params_t, minhalo>;
public:
using real_t = typename ct_params_t::real_t;
protected:
real_t smg_c, c_m;
typename parent_t::arr_t &k_m;
void multiply_sgs_visc()
{
const auto dlta = std::accumulate(this->dijk.begin(), this->dijk.end(), real_t(0.)) / 3;
if (static_cast<stress_diff_t>(ct_params_t::stress_diff) == compact)
{
k_m(this->ijk) = pow(smg_c * dlta, 2) * formulae::stress::calc_tdef_sq_cmpct<ct_params_t::n_dims>(this->tau, this->ijk);
formulae::stress::multiply_tnsr_cmpct<ct_params_t::n_dims, ct_params_t::opts>(this->tau,
real_t(1.0),
this->k_m,
*this->mem->G,
this->ijkm_sep);
this->xchng_sgs_tnsr_diag(this->tau, this->vips()[ct_params_t::n_dims - 1], this->vip_div, this->ijk);
this->xchng_sgs_tnsr_offdiag(this->tau, this->tau_srfc, this->ijk, this->ijkm);
}
else
{
k_m(this->ijk) = pow(smg_c * dlta, 2) * formulae::stress::calc_tdef_sq<ct_params_t::n_dims>(this->tau, this->ijk);
for (auto& t : this->tau)
{
t(this->ijk) *= k_m(this->ijk);
}
}
}
public:
struct rt_params_t : parent_t::rt_params_t
{
real_t smg_c, c_m;
};
// ctor
mpdata_rhs_vip_prs_sgs_smg(
typename parent_t::ctor_args_t args,
const rt_params_t &p
) :
parent_t(args, p),
smg_c(p.smg_c),
c_m(p.c_m),
k_m(args.mem->tmp[__FILE__][0][0])
{
if (smg_c == 0) throw std::runtime_error("libmpdata++: smg_c == 0");
}
static void alloc(
typename parent_t::mem_t *mem,
const int &n_iters
) {
parent_t::alloc(mem, n_iters);
parent_t::alloc_tmp_sclr(mem, __FILE__, 1); // k_m
}
};
} // namespace detail
} // namespace solvers
} // namespace libmpdataxx