Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor ConstantAlpha to allow running with turbulence model #1391

Merged
merged 16 commits into from
Jan 22, 2024
Merged

Refactor ConstantAlpha to allow running with turbulence model #1391

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
e1bdd14
Enable ConstantAlpha with LES/PBL
ewquon Jan 19, 2024
3db98da
Throw error if Smagorinsky selected w/o setting Cs
ewquon Jan 19, 2024
c012775
Fix bug in incompressible path
ewquon Jan 19, 2024
6fa9393
Cleanup: variable names for clarity and consistency, update comments
ewquon Jan 19, 2024
6200a0e
Handle ConstantAlpha case when calling ComputeStress*
ewquon Jan 19, 2024
231c9c0
Apply the same changes for ConstantAlpha as ERF_slow_rhs_pre
ewquon Jan 19, 2024
00b3352
No need to scale by rho here in the ConstantAlpha solver path
ewquon Jan 19, 2024
243f252
Merge branch 'development' into constalpha_with_turb
AMLattanzi Jan 19, 2024
c308955
Merge branch 'development' into constalpha_with_turb
AMLattanzi Jan 20, 2024
71f1c31
Update golds after ConstantAlpha refactor
ewquon Jan 20, 2024
1834373
Make interpolation of rhoAlpha consistent with mu_turb
ewquon Jan 21, 2024
23d55a5
Calculate rhoalpha for ghost cells in xyz
ewquon Jan 21, 2024
3792202
Revert "Update golds after ConstantAlpha refactor"
ewquon Jan 22, 2024
c290c0a
New gold data after fixing indexing and interpolation
ewquon Jan 22, 2024
1ebf0bd
Average rho to faces, following mu_turb; update gold data
ewquon Jan 22, 2024
0d4d6b7
Removed fixed_mri_dt_ratio
ewquon Jan 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 10 additions & 6 deletions Exec/RegTests/DensityCurrent/inputs_refsoln
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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;
});
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the off-diagonal stresses, shouldn't the cell_data be averaged to the location of the strain (this would hold for other places in the code as well)?
image

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, you're right. I have a fix coming.

}
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
Loading