From b493e423f868169ead6c1cf03d8ebc1e5e13a6f8 Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Wed, 14 Feb 2024 16:49:34 +0100 Subject: [PATCH 1/2] add formulae for anisotropic eddy viscosity following Simon and Chow 2021 --- libmpdata++/formulae/stress_formulae.hpp | 75 ++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/libmpdata++/formulae/stress_formulae.hpp b/libmpdata++/formulae/stress_formulae.hpp index b009dccf..249f7cbf 100644 --- a/libmpdata++/formulae/stress_formulae.hpp +++ b/libmpdata++/formulae/stress_formulae.hpp @@ -492,6 +492,31 @@ namespace libmpdataxx real_t(0.5) * (G(rho, ijk[0], ijk[1], ijk[2] + 1) + G(rho, ijk[0], ijk[1], ijk[2])); } + // multiplication of compact vector 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, i.e. horizontal for all diagonal components + template + 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, + typename std::enable_if::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(rho, ijk[0] + 1, ijk[1], ijk[2]) + G(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(rho, ijk[0], ijk[1] + 1, ijk[2]) + G(rho, ijk[0], ijk[1], ijk[2])); + + 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(rho, ijk[0], ijk[1], ijk[2] + 1) + G(rho, ijk[0], ijk[1], ijk[2])); + } + // multiplication of compact tensor components by constant molecular viscosity // 2D version template @@ -587,6 +612,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 + 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::type* = 0) + { + multiply_vctr_cmpct(av, coeff, k_ma, rho, ijk); + 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(rho, ijk[0] + 1, ijk[1] , ijk[2]) + + G(rho, ijk[0] , ijk[1] , ijk[2]) + + G(rho, ijk[0] + 1, ijk[1] + 1, ijk[2]) + + G(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(rho, ijk[0] + 1, ijk[1], ijk[2] ) + + G(rho, ijk[0] , ijk[1], ijk[2] ) + + G(rho, ijk[0] + 1, ijk[1], ijk[2] + 1) + + G(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(rho, ijk[0], ijk[1] , ijk[2] + 1) + + G(rho, ijk[0], ijk[1] , ijk[2] ) + + G(rho, ijk[0], ijk[1] + 1, ijk[2] + 1) + + G(rho, ijk[0], ijk[1] + 1, ijk[2] ) + ); + } + // flux divergence // 2D version template From c6eaeeac3b1f4d5219806b226ff77fbab0a5e3de Mon Sep 17 00:00:00 2001 From: Piotr Dziekan Date: Thu, 15 Feb 2024 15:48:04 +0100 Subject: [PATCH 2/2] multiply_vctr_cmpct with k_ma[1] for vertical --- libmpdata++/formulae/stress_formulae.hpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/libmpdata++/formulae/stress_formulae.hpp b/libmpdata++/formulae/stress_formulae.hpp index 249f7cbf..7724a73c 100644 --- a/libmpdata++/formulae/stress_formulae.hpp +++ b/libmpdata++/formulae/stress_formulae.hpp @@ -494,14 +494,15 @@ namespace libmpdataxx // multiplication of compact vector 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, i.e. horizontal for all diagonal components + // 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 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::type* = 0) { av[0](ijk[0] + h, ijk[1], ijk[2]) *= coeff * @@ -512,7 +513,12 @@ namespace libmpdataxx 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(rho, ijk[0], ijk[1] + 1, ijk[2]) + G(rho, ijk[0], ijk[1], ijk[2])); - av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff * + 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(rho, ijk[0], ijk[1], ijk[2] + 1) + G(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(rho, ijk[0], ijk[1], ijk[2] + 1) + G(rho, ijk[0], ijk[1], ijk[2])); } @@ -624,7 +630,7 @@ namespace libmpdataxx const ijk_t &ijk, typename std::enable_if::type* = 0) { - multiply_vctr_cmpct(av, coeff, k_ma, rho, ijk); + multiply_vctr_cmpct(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])