Skip to content

Commit

Permalink
Refactor ConstantAlpha to allow running with turbulence model (erf-mo…
Browse files Browse the repository at this point in the history 
…del#1391)

* Enable ConstantAlpha with LES/PBL

* Throw error if Smagorinsky selected w/o setting Cs

* Fix bug in incompressible path

* Cleanup: variable names for clarity and consistency, update comments

* Handle ConstantAlpha case when calling ComputeStress*

- ComputeStress*() now have cell_data as an argument
- cell_data is set if l_use_constAlpha, otherwise a nullptr is passed
- mu_eff is now converted to alpha prior to being passed to
  ComputeStress*()

* Apply the same changes for ConstantAlpha as ERF_slow_rhs_pre

* No need to scale by rho here in the ConstantAlpha solver path

For generality, rho scaling is applied when tau is calculated. This
enables ConstantAlpha to be used with a turbulence model.

* Make interpolation of rhoAlpha consistent with mu_turb

* Calculate rhoalpha for ghost cells in xyz

* New gold data after fixing indexing and interpolation

* Average rho to faces, following mu_turb; update gold data

* Removed fixed_mri_dt_ratio

Otherwise ref solution crashes on step 938 (t=234.25). Orig fixed dt ratio was
4; allowing a variable dt ratio results in a constant dt ratio of 6.

---------

Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
@ewquon @AMLattanzi
ewquon and AMLattanzi authored Jan 22, 2024
1 parent 88ea1f0 commit 9ea9813
Show file tree
Hide file tree
Showing 37 changed files with 1,022 additions and 904 deletions.
16 changes: 10 additions & 6 deletions Exec/RegTests/DensityCurrent/inputs_refsoln
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,11 @@ zlo.type = "SlipWall"
zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.fixed_dt = 0.25 # fixed time step [s] -- Straka et al 1993
erf.fixed_fast_dt = 0.0625 # fixed time step [s] -- Straka et al 1993
# Straka et al 1993:
#erf.fixed_dt = 1.5625e-2 # [s]
#erf.no_substepping = 1

erf.fixed_dt = 0.25 # large time step [s] -- results in dt ratio = 6

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
Expand All @@ -31,11 +34,12 @@ amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = 57600 # number of timesteps between checkpoints
erf.check_int = -1 # number of timesteps between checkpoints

# PLOTFILES
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 3840 # number of timesteps between plotfiles
erf.plot_int_1 = 19200 # number of timesteps between plotfiles
erf.plot_int_1 = 1200 # number of timesteps between plotfiles
erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens

# SOLVER CHOICE
Expand All @@ -47,9 +51,9 @@ erf.les_type = "None"
#
# diffusion coefficient from Straka, K = 75 m^2/s
#
erf.molec_diff_type = "ConstantAlpha"
erf.molec_diff_type = "ConstantAlpha" # where alpha == "K" in Straka et al 1993
erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
erf.dynamicViscosity = 75.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
erf.dynamicViscosity = 75.0 # [kg/(m-s)] ==> alpha = 75.0 m^2/s
erf.alpha_T = 75.0 # [m^2/s]

erf.c_p = 1004.0
Expand Down
7 changes: 3 additions & 4 deletions Exec/RegTests/DensityCurrent/runscript_refsoln
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,9 @@
#!/bin/bash
#SBATCH --account=erf
#SBATCH --time=1:30:00
#SBATCH --time=0:30:00
#SBATCH --job-name=DensityCurrent_ref_soln
#SBATCH --ntasks=256
#SBATCH --mail-user [email protected]
#SBATCH --mail-type BEGIN,END,FAIL
#SBATCH --ntasks=288
# Timing with 288 CPU cores on NREL Eagle HPC: 427.1469369 s

# load environment
. ~/.bash_profile
Expand Down
9 changes: 0 additions & 9 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,15 +210,6 @@ struct SolverChoice {
turbChoice[lev].init_params(lev,max_level);
}

// If running LES/PBL then molecular diffusion must be "Constant" or "None"
for (int lev = 0; lev <= max_level; lev++) {
if (turbChoice[lev].les_type != LESType::None) {
if ( diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha ) {
amrex::Error("We don't allow LES with MolecDiffType::ConstantAlpha");
}
}
}

// Which type of refinement
static std::string coupling_type_string = "OneWay";
pp.query("coupling_type",coupling_type_string);
Expand Down
7 changes: 7 additions & 0 deletions Source/DataStructs/TurbStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,13 @@ struct TurbChoice {

pp.query("theta_ref", theta_ref);

// Validate inputs
if (les_type == LESType::Smagorinsky) {
if (Cs == 0) {
amrex::Error("Need to specify Cs for Smagorsinky LES");
}
}

// Here we cover the case where either one is > 1 and the other is 0, or they both are the same and > 1
} else {

Expand Down
174 changes: 128 additions & 46 deletions Source/Diffusion/ComputeStress_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ using namespace amrex;
* @param[in] tbxxz nodal xz box for tau_13
* @param[in] tbxyz nodal yz box for tau_23
* @param[in] mu_eff constant molecular viscosity
* @param[in] cell_data to access rho if ConstantAlpha
* @param[in,out] tau11 11 strain -> stress
* @param[in,out] tau22 22 strain -> stress
* @param[in,out] tau33 33 strain -> stress
Expand All @@ -20,32 +21,66 @@ using namespace amrex;
*/
void
ComputeStressConsVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff,
const Array4<const Real>& cell_data,
Array4<Real>& tau11, Array4<Real>& tau22, Array4<Real>& tau33,
Array4<Real>& tau12, Array4<Real>& tau13, Array4<Real>& tau23,
const Array4<const Real>& er_arr)
{
Real OneThird = (1./3.);

// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
if (cell_data)
// constant alpha (stored in mu_eff)
{
tau11(i,j,k) = -mu_eff * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -mu_eff * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -mu_eff * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});
// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real rhoAlpha = cell_data(i, j, k, Rho_comp) * mu_eff;
tau11(i,j,k) = -rhoAlpha * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -rhoAlpha * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -rhoAlpha * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});

// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) *= -mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) *= -mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) *= -mu_eff;
});
// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i-1, j , k, Rho_comp) + cell_data(i, j , k, Rho_comp)
+ cell_data(i-1, j-1, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
tau12(i,j,k) *= -rho_bar * mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i-1, j, k , Rho_comp) + cell_data(i, j, k , Rho_comp)
+ cell_data(i-1, j, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
tau13(i,j,k) *= -rho_bar * mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i, j-1, k , Rho_comp) + cell_data(i, j, k , Rho_comp)
+ cell_data(i, j-1, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
tau23(i,j,k) *= -rho_bar * mu_eff;
});
}
else
// constant mu_eff
{
// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
tau11(i,j,k) = -mu_eff * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -mu_eff * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -mu_eff * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});

// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) *= -mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau13(i,j,k) *= -mu_eff;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau23(i,j,k) *= -mu_eff;
});
}
}

/**
Expand All @@ -57,6 +92,7 @@ ComputeStressConsVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff,
* @param[in] tbxyz nodal yz box for tau_23
* @param[in] mu_eff constant molecular viscosity
* @param[in] mu_turb variable turbulent viscosity
* @param[in] cell_data to access rho if ConstantAlpha
* @param[in,out] tau11 11 strain -> stress
* @param[in,out] tau22 22 strain -> stress
* @param[in,out] tau33 33 strain -> stress
Expand All @@ -68,40 +104,86 @@ ComputeStressConsVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff,
void
ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff,
const Array4<const Real>& mu_turb,
const Array4<const Real>& cell_data,
Array4<Real>& tau11, Array4<Real>& tau22, Array4<Real>& tau33,
Array4<Real>& tau12, Array4<Real>& tau13, Array4<Real>& tau23,
const Array4<const Real>& er_arr)
{
Real OneThird = (1./3.);

// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_11 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_h);
Real mu_22 = mu_11;
Real mu_33 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v);
tau11(i,j,k) = -mu_11 * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -mu_22 * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -mu_33 * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});
if (cell_data)
// constant alpha (stored in mu_eff)
{
// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rhoAlpha = cell_data(i, j, k, Rho_comp) * mu_eff;
Real mu_11 = rhoAlpha + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_h);
Real mu_22 = mu_11;
Real mu_33 = rhoAlpha + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v);
tau11(i,j,k) = -mu_11 * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -mu_22 * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -mu_33 * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});

// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i-1, j , k, Rho_comp) + cell_data(i, j , k, Rho_comp)
+ cell_data(i-1, j-1, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) );
Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
+ mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
Real mu_12 = rho_bar*mu_eff + 2.0*mu_bar;
tau12(i,j,k) *= -mu_12;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i-1, j, k , Rho_comp) + cell_data(i, j, k , Rho_comp)
+ cell_data(i-1, j, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_13 = rho_bar*mu_eff + 2.0*mu_bar;
tau13(i,j,k) *= -mu_13;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real rho_bar = 0.25*( cell_data(i, j-1, k , Rho_comp) + cell_data(i, j, k , Rho_comp)
+ cell_data(i, j-1, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_23 = rho_bar*mu_eff + 2.0*mu_bar;
tau23(i,j,k) *= -mu_23;
});
}
else
// constant mu_eff
{
// Cell centered strains
ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_11 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_h);
Real mu_22 = mu_11;
Real mu_33 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v);
tau11(i,j,k) = -mu_11 * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) );
tau22(i,j,k) = -mu_22 * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) );
tau33(i,j,k) = -mu_33 * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) );
});

// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
+ mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
Real mu_12 = mu_eff + 2.0*mu_bar;
tau12(i,j,k) *= -mu_12;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_13 = mu_eff + 2.0*mu_bar;
tau13(i,j,k) *= -mu_13;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_23 = mu_eff + 2.0*mu_bar;
tau23(i,j,k) *= -mu_23;
});
// Off-diagonal strains
ParallelFor(tbxxy,tbxxz,tbxyz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i-1, j , k, EddyDiff::Mom_h) + mu_turb(i, j , k, EddyDiff::Mom_h)
+ mu_turb(i-1, j-1, k, EddyDiff::Mom_h) + mu_turb(i, j-1, k, EddyDiff::Mom_h) );
Real mu_12 = mu_eff + 2.0*mu_bar;
tau12(i,j,k) *= -mu_12;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_13 = mu_eff + 2.0*mu_bar;
tau13(i,j,k) *= -mu_13;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v)
+ mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) );
Real mu_23 = mu_eff + 2.0*mu_bar;
tau23(i,j,k) *= -mu_23;
});
}
}
Loading

0 comments on commit 9ea9813

Please sign in to comment.