Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Aniso smg #493

Merged
merged 2 commits into from
Aug 22, 2024
Merged

Aniso smg #493

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b493e42
add formulae for anisotropic eddy viscosity following Simon and Chow …
pdziekan Feb 14, 2024
c6eaeea
multiply_vctr_cmpct with k_ma[1] for vertical
pdziekan Feb 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 81 additions & 0 deletions libmpdata++/formulae/stress_formulae.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -492,6 +492,37 @@ namespace libmpdataxx
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
}

// multiplication of compact vector components by variable anisotropic eddy viscosity
// 3D version
// TODO: when used in tensor multiplication, it is not clear which eddy viscosity components should be used for each term; see Simon and Chow 2021 p. 17
// we opt to use the same approach as therein, i.e. horizontal for all diagonal components; this is controlled by the flag vh
template <int nd, opts_t opts, class arrvec_t, class real_t, class arr_t, class ijk_t>
inline void multiply_vctr_cmpct(const arrvec_t &av,
real_t coeff,
const arrvec_t & k_ma, // two arrays: [0] for horizontal eddy viscosity and [1] for vertical
const arr_t & rho,
const ijk_t &ijk,
const bool vh = false,
typename std::enable_if<nd == 3>::type* = 0)
{
av[0](ijk[0] + h, ijk[1], ijk[2]) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0] + 1, ijk[1], ijk[2]) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]) + G<opts, 0>(rho,ijk[0], ijk[1], ijk[2]));

av[1](ijk[0], ijk[1] + h, ijk[2]) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0], ijk[1] + 1, ijk[2]) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));

if(!vh)
av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff *
real_t(0.5) * (k_ma[1](ijk[0], ijk[1], ijk[2] + 1) + k_ma[1](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
else
av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0], ijk[1], ijk[2] + 1) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
}

// multiplication of compact tensor components by constant molecular viscosity
// 2D version
template <int nd, class arrvec_t, class real_t, class ijk_t>
Expand Down Expand Up @@ -587,6 +618,56 @@ namespace libmpdataxx
);
}

// multiplication of compact tensor components by variable anisotropic eddy viscosity
// 3D version
// TODO: it is not clear which eddy viscosity components should be used for each term; see Simon and Chow 2021 p. 17
// we opt to use the same approach as therein
template <int nd, opts_t opts, class arrvec_t, class real_t, class arr_t, class ijk_t>
inline void multiply_tnsr_cmpct(const arrvec_t &av,
const real_t coeff,
const arrvec_t &k_ma,
const arr_t &rho,
const ijk_t &ijk,
typename std::enable_if<nd == 3>::type* = 0)
{
multiply_vctr_cmpct<nd, opts>(av, coeff, k_ma, rho, ijk, true);
av[3](ijk[0] + h, ijk[1] + h, ijk[2]) *= coeff *
real_t(0.25) * ( k_ma[0](ijk[0] + 1, ijk[1] , ijk[2])
+ k_ma[0](ijk[0] , ijk[1] , ijk[2])
+ k_ma[0](ijk[0] + 1, ijk[1] + 1, ijk[2])
+ k_ma[0](ijk[0] , ijk[1] + 1, ijk[2])
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1] , ijk[2])
+ G<opts, 0>(rho, ijk[0] , ijk[1] , ijk[2])
+ G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1, ijk[2])
+ G<opts, 0>(rho, ijk[0] , ijk[1] + 1, ijk[2])
);

av[4](ijk[0] + h, ijk[1], ijk[2] + h) *= coeff *
real_t(0.25) * ( k_ma[1](ijk[0] + 1, ijk[1], ijk[2] )
+ k_ma[1](ijk[0] , ijk[1], ijk[2] )
+ k_ma[1](ijk[0] + 1, ijk[1], ijk[2] + 1)
+ k_ma[1](ijk[0] , ijk[1], ijk[2] + 1)
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] )
+ G<opts, 0>(rho, ijk[0] , ijk[1], ijk[2] )
+ G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0] , ijk[1], ijk[2] + 1)
);

av[5](ijk[0], ijk[1] + h, ijk[2] + h) *= coeff *
real_t(0.25) * ( k_ma[1](ijk[0], ijk[1] , ijk[2] + 1)
+ k_ma[1](ijk[0], ijk[1] , ijk[2] )
+ k_ma[1](ijk[0], ijk[1] + 1, ijk[2] + 1)
+ k_ma[1](ijk[0], ijk[1] + 1, ijk[2] )
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0], ijk[1] , ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0], ijk[1] , ijk[2] )
+ G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] )
);
}

// flux divergence
// 2D version
template <int nd, opts_t opts, class arrvec_t, class arr_t, class ijk_t, class dijk_t>
Expand Down
Loading