Skip to content

Commit

Permalink
Templated strainrate to reduce code duplication (Exawind#25)
Browse files Browse the repository at this point in the history 
* templated strainrate to reduce code duplication, thanks Shreyas for the suggestions!

* removed a const
  • Loading branch information
@michaeljbrazell
michaeljbrazell authored Mar 25, 2020
1 parent b3f79e5 commit 199261b
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 162 deletions.
238 changes: 83 additions & 155 deletions src/derive/derive_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,174 +4,102 @@
#include <AMReX_FArrayBox.H>
#include <cmath>

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
struct StencilInterior
{
using namespace amrex;

const Real ux = 0.5 * (vel(i+1,j,k,0) - vel(i-1,j,k,0)) * idx;
const Real vx = 0.5 * (vel(i+1,j,k,1) - vel(i-1,j,k,1)) * idx;
const Real wx = 0.5 * (vel(i+1,j,k,2) - vel(i-1,j,k,2)) * idx;

const Real uy = 0.5 * (vel(i,j+1,k,0) - vel(i,j-1,k,0)) * idy;
const Real vy = 0.5 * (vel(i,j+1,k,1) - vel(i,j-1,k,1)) * idy;
const Real wy = 0.5 * (vel(i,j+1,k,2) - vel(i,j-1,k,2)) * idy;

const Real uz = 0.5 * (vel(i,j,k+1,0) - vel(i,j,k-1,0)) * idz;
const Real vz = 0.5 * (vel(i,j,k+1,1) - vel(i,j,k-1,1)) * idz;
const Real wz = 0.5 * (vel(i,j,k+1,2) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// i = 0 first cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_ilo (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5}
};
};

struct StencilILO
{
using namespace amrex;

// average velocity to face i+1/2 and take derivative about cell center between interior face and boundary face i-1
const Real ux = (0.5 * (vel(i+1,j,k,0) + vel(i,j,k,0)) - vel(i-1,j,k,0)) * idx;
const Real vx = (0.5 * (vel(i+1,j,k,1) + vel(i,j,k,1)) - vel(i-1,j,k,1)) * idx;
const Real wx = (0.5 * (vel(i+1,j,k,2) + vel(i,j,k,2)) - vel(i-1,j,k,2)) * idx;

const Real uy = 0.5 * (vel(i,j+1,k,0) - vel(i,j-1,k,0)) * idy;
const Real vy = 0.5 * (vel(i,j+1,k,1) - vel(i,j-1,k,1)) * idy;
const Real wy = 0.5 * (vel(i,j+1,k,2) - vel(i,j-1,k,2)) * idy;

const Real uz = 0.5 * (vel(i,j,k+1,0) - vel(i,j,k-1,0)) * idz;
const Real vz = 0.5 * (vel(i,j,k+1,1) - vel(i,j,k-1,1)) * idz;
const Real wz = 0.5 * (vel(i,j,k+1,2) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// j = 0 first cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_jlo (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.5, -1.0},
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5}
};
};

struct StencilJLO
{
using namespace amrex;

const Real ux = 0.5 * (vel(i+1,j,k,0) - vel(i-1,j,k,0)) * idx;
const Real vx = 0.5 * (vel(i+1,j,k,1) - vel(i-1,j,k,1)) * idx;
const Real wx = 0.5 * (vel(i+1,j,k,2) - vel(i-1,j,k,2)) * idx;

// average velocity to face k+1/2 and take derivative about cell center between interior face and boundary face k-1
const Real uy = (0.5 * (vel(i,j+1,k,0) + vel(i,j,k,0)) - vel(i,j-1,k,0)) * idy;
const Real vy = (0.5 * (vel(i,j+1,k,1) + vel(i,j,k,1)) - vel(i,j-1,k,1)) * idy;
const Real wy = (0.5 * (vel(i,j+1,k,2) + vel(i,j,k,2)) - vel(i,j-1,k,2)) * idy;

const Real uz = 0.5 * (vel(i,j,k+1,0) - vel(i,j,k-1,0)) * idz;
const Real vz = 0.5 * (vel(i,j,k+1,1) - vel(i,j,k-1,1)) * idz;
const Real wz = 0.5 * (vel(i,j,k+1,2) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// k = 0 first cell
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_klo (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.0, -0.5},
{0.5, 0.5, -1.0},
{0.5, 0.0, -0.5}
};
};

struct StencilKLO
{
using namespace amrex;

const Real ux = 0.5 * (vel(i+1,j,k,0) - vel(i-1,j,k,0)) * idx;
const Real vx = 0.5 * (vel(i+1,j,k,1) - vel(i-1,j,k,1)) * idx;
const Real wx = 0.5 * (vel(i+1,j,k,2) - vel(i-1,j,k,2)) * idx;

const Real uy = 0.5 * (vel(i,j+1,k,0) - vel(i,j-1,k,0)) * idy;
const Real vy = 0.5 * (vel(i,j+1,k,1) - vel(i,j-1,k,1)) * idy;
const Real wy = 0.5 * (vel(i,j+1,k,2) - vel(i,j-1,k,2)) * idy;

// average velocity to face k+1/2 and take derivative about cell center between interior face and boundary face k-1
const Real uz = (0.5 * (vel(i,j,k+1,0) + vel(i,j,k,0)) - vel(i,j,k-1,0)) * idz;
const Real vz = (0.5 * (vel(i,j,k+1,1) + vel(i,j,k,1)) - vel(i,j,k-1,1)) * idz;
const Real wz = (0.5 * (vel(i,j,k+1,2) + vel(i,j,k,2)) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// i = ncell-1
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_ihi (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5},
{0.5, 0.5, -1.0}
};
};

struct StencilIHI
{
using namespace amrex;

// average velocity to face i-1/2 and take derivative about cell center between interior face and boundary face i+1
const Real ux = (vel(i+1,j,k,0) - 0.5 * (vel(i,j,k,0) + vel(i-1,j,k,0))) * idx;
const Real vx = (vel(i+1,j,k,1) - 0.5 * (vel(i,j,k,1) + vel(i-1,j,k,1))) * idx;
const Real wx = (vel(i+1,j,k,2) - 0.5 * (vel(i,j,k,2) + vel(i-1,j,k,2))) * idx;

const Real uy = 0.5 * (vel(i,j+1,k,0) - vel(i,j-1,k,0)) * idy;
const Real vy = 0.5 * (vel(i,j+1,k,1) - vel(i,j-1,k,1)) * idy;
const Real wy = 0.5 * (vel(i,j+1,k,2) - vel(i,j-1,k,2)) * idy;

const Real uz = 0.5 * (vel(i,j,k+1,0) - vel(i,j,k-1,0)) * idz;
const Real vz = 0.5 * (vel(i,j,k+1,1) - vel(i,j,k-1,1)) * idz;
const Real wz = 0.5 * (vel(i,j,k+1,2) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// j = ncell-1
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_jhi (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
static constexpr amrex::Real coeffs[3][3] = {
{1.0, -0.5, -0.5},
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5}
};
};

struct StencilJHI
{
using namespace amrex;

const Real ux = 0.5 * (vel(i+1,j,k,0) - vel(i-1,j,k,0)) * idx;
const Real vx = 0.5 * (vel(i+1,j,k,1) - vel(i-1,j,k,1)) * idx;
const Real wx = 0.5 * (vel(i+1,j,k,2) - vel(i-1,j,k,2)) * idx;

// average velocity to face j-1/2 and take derivative about cell center between interior face and boundary face j+1
const Real uy = (vel(i,j+1,k,0) - 0.5 * (vel(i,j,k,0) + vel(i,j-1,k,0))) * idy;
const Real vy = (vel(i,j+1,k,1) - 0.5 * (vel(i,j,k,1) + vel(i,j-1,k,1))) * idy;
const Real wy = (vel(i,j+1,k,2) - 0.5 * (vel(i,j,k,2) + vel(i,j-1,k,2))) * idy;

const Real uz = 0.5 * (vel(i,j,k+1,0) - vel(i,j,k-1,0)) * idz;
const Real vz = 0.5 * (vel(i,j,k+1,1) - vel(i,j,k-1,1)) * idz;
const Real wz = 0.5 * (vel(i,j,k+1,2) - vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
}

// k = ncell-1
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.0, -0.5},
{1.0, -0.5, -0.5},
{0.5, 0.0, -0.5}
};
};

struct StencilKHI
{
static constexpr amrex::Real coeffs[3][3] = {
{0.5, 0.0, -0.5},
{0.5, 0.0, -0.5},
{1.0, -0.5, -0.5}
};
};

template<typename Stencil>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real incflo_strainrate_khi (int i, int j, int k,
amrex::Real incflo_strainrate (int i, int j, int k,
amrex::Real idx, amrex::Real idy, amrex::Real idz,
amrex::Array4<amrex::Real const> const& vel) noexcept
{
using namespace amrex;

Real cp1, c, cm1;

cp1 = Stencil::coeffs[0][0];
c = Stencil::coeffs[0][1];
cm1 = Stencil::coeffs[0][2];

const Real ux = (cp1*vel(i+1,j,k,0) + c*vel(i,j,k,0) + cm1*vel(i-1,j,k,0)) * idx;
const Real vx = (cp1*vel(i+1,j,k,1) + c*vel(i,j,k,1) + cm1*vel(i-1,j,k,1)) * idx;
const Real wx = (cp1*vel(i+1,j,k,2) + c*vel(i,j,k,2) + cm1*vel(i-1,j,k,2)) * idx;

const Real ux = 0.5 * (vel(i+1,j,k,0) - vel(i-1,j,k,0)) * idx;
const Real vx = 0.5 * (vel(i+1,j,k,1) - vel(i-1,j,k,1)) * idx;
const Real wx = 0.5 * (vel(i+1,j,k,2) - vel(i-1,j,k,2)) * idx;

const Real uy = 0.5 * (vel(i,j+1,k,0) - vel(i,j-1,k,0)) * idy;
const Real vy = 0.5 * (vel(i,j+1,k,1) - vel(i,j-1,k,1)) * idy;
const Real wy = 0.5 * (vel(i,j+1,k,2) - vel(i,j-1,k,2)) * idy;
cp1 = Stencil::coeffs[1][0];
c = Stencil::coeffs[1][1];
cm1 = Stencil::coeffs[1][2];
const Real uy = (cp1*vel(i,j+1,k,0) + c*vel(i,j,k,0) + cm1*vel(i,j-1,k,0)) * idy;
const Real vy = (cp1*vel(i,j+1,k,1) + c*vel(i,j,k,1) + cm1*vel(i,j-1,k,1)) * idy;
const Real wy = (cp1*vel(i,j+1,k,2) + c*vel(i,j,k,2) + cm1*vel(i,j-1,k,2)) * idy;

// average velocity to face k-1/2 and take derivative about cell center between interior face and boundary face k+1
const Real uz = (vel(i,j,k+1,0) - 0.5 * (vel(i,j,k,0) + vel(i,j,k-1,0))) * idz;
const Real vz = (vel(i,j,k+1,1) - 0.5 * (vel(i,j,k,1) + vel(i,j,k-1,1))) * idz;
const Real wz = (vel(i,j,k+1,2) - 0.5 * (vel(i,j,k,2) + vel(i,j,k-1,2))) * idz;
cp1 = Stencil::coeffs[2][0];
c = Stencil::coeffs[2][1];
cm1 = Stencil::coeffs[2][2];

const Real uz = (cp1*vel(i,j,k+1,0) + c*vel(i,j,k,0) + cm1*vel(i,j,k-1,0)) * idz;
const Real vz = (cp1*vel(i,j,k+1,1) + c*vel(i,j,k,1) + cm1*vel(i,j,k-1,1)) * idz;
const Real wz = (cp1*vel(i,j,k+1,2) + c*vel(i,j,k,2) + cm1*vel(i,j,k-1,2)) * idz;

return std::sqrt(2.0 * ux*ux + 2.0 * vy*vy + 2.0 * wz*wz
+ (uy+vx)*(uy+vx) + (vz+wy)*(vz+wy) + (wx+uz)*(wx+uz));
Expand Down
14 changes: 7 additions & 7 deletions src/rheology/incflo_rheology.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// fixme strainrate does not account for one sided stencils near walls
Real sr = incflo_strainrate(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilInterior>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);

Expand All @@ -157,7 +157,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxlo,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_ilo(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilILO>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand All @@ -176,7 +176,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxhi,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_ihi(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilIHI>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand All @@ -198,7 +198,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxlo,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_jlo(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilJLO>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand All @@ -217,7 +217,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxhi,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_jhi(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilJHI>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand All @@ -239,7 +239,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxlo,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_klo(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilKLO>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand All @@ -258,7 +258,7 @@ void incflo::compute_viscosity (Vector<MultiFab*> const& vel_eta,
amrex::ParallelFor(bxhi,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real sr = incflo_strainrate_khi(i,j,k,idx,idy,idz,vel_arr);
Real sr = incflo_strainrate<StencilKHI>(i,j,k,idx,idy,idz,vel_arr);
Real den = rho_arr(i,j,k);
eta_arr(i,j,k) = non_newtonian_viscosity(sr,den,ds);
});
Expand Down

0 comments on commit 199261b

Please sign in to comment.