Skip to content

Commit

Permalink
Merge pull request #493 from pdziekan/aniso_smg
Browse files Browse the repository at this point in the history 
Aniso smg
  • Loading branch information
@pdziekan
pdziekan authored Aug 22, 2024
2 parents d53ff5e + c6eaeea commit df9e0b2
Showing 1 changed file with 81 additions and 0 deletions.
81 changes: 81 additions & 0 deletions libmpdata++/formulae/stress_formulae.hpp
View file
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

0 comments on commit df9e0b2

Please sign in to comment.