Refactor ConstantAlpha to allow running with turbulence model (#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]>
erf-model · Jan 22, 2024 · 9ea9813 · 9ea9813
1 parent 88ea1f0
commit 9ea9813
Show file tree

Hide file tree

Showing 37 changed files with 1,022 additions and 904 deletions.
diff --git a/Exec/RegTests/DensityCurrent/inputs_refsoln b/Exec/RegTests/DensityCurrent/inputs_refsoln
@@ -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
@@ -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
@@ -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

diff --git a/Exec/RegTests/DensityCurrent/runscript_refsoln b/Exec/RegTests/DensityCurrent/runscript_refsoln
@@ -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

diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H
@@ -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);

diff --git a/Source/DataStructs/TurbStruct.H b/Source/DataStructs/TurbStruct.H
@@ -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 {
 

diff --git a/Source/Diffusion/ComputeStress_N.cpp b/Source/Diffusion/ComputeStress_N.cpp
@@ -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
@@ -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;
+ });
+ }
 }
 
 /**
@@ -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
@@ -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;
+ });
+ }
 }