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

QKE src fixes #1252

Merged
merged 4 commits into from
Sep 19, 2023
+15 −15
Merged

QKE src fixes #1252

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
50551e7
This runs multiple steps without a negative QKE and the field looks r…
AMLattanzi Sep 19, 2023
534c026
Fix debug modification.
AMLattanzi Sep 19, 2023
f50fd2a
fix sqrt in brunt vaisala freq
baperry2 Sep 19, 2023
89feda7
Merge pull request #1 from baperry2/sqrt-brunt-vaisala
AMLattanzi Sep 19, 2023
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
18 changes: 9 additions & 9 deletions Exec/DevTests/MiguelDev/inputs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ zhi.type = "SlipWall"

# TIME STEP CONTROL
erf.no_substepping = 0
erf.fixed_dt = 1 # fixed time step
erf.fixed_dt = 0.5 # fixed time step

# DIAGNOSTICS & VERBOSITY
erf.sum_interval = 1 # timesteps between computing mass
Expand All @@ -40,8 +40,8 @@ erf.check_int =-100 # number of timesteps between checkpoints
# PLOTFILES
erf.plotfile_type = netcdf
erf.plot_file_1 = plt # prefix of plotfile name
erf.plot_int_1 = 100 # number of timesteps between plotfiles
erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure temp theta pres_hse
erf.plot_int_1 = 100 # number of timesteps between plotfiles
erf.plot_vars_1 = density rhotheta rhoQKE x_velocity y_velocity z_velocity pressure temp theta QKE pres_hse

# SOLVER CHOICE
erf.alpha_T = 0.0
Expand All @@ -55,7 +55,7 @@ erf.dynamicViscosity = 0.0 # 5.0
erf.Cs = 0.16

erf.dycore_horiz_adv_type = "Upwind_3rd"
erf.dycore_vert_adv_type = "Upwind_5th"
erf.dycore_vert_adv_type = "Upwind_5th"

erf.use_coriolis = true
erf.abl_driver_type = "GeostrophicWind"
Expand All @@ -66,11 +66,11 @@ erf.use_rayleigh_damping = true #true

# PROBLEM PARAMETERS (optional)
prob.rho_0 = 1.0
prob.T_0 = 300.0
prob.Xc_0 = 100000
prob.Yc_0 = 100000
prob.VMAX = 15
prob.RMAX = 20000
prob.T_0 = 300.0
prob.Xc_0 = 100000
prob.Yc_0 = 100000
prob.VMAX = 15
prob.RMAX = 20000
prob.RZERO = 80000
prob.ZZERO = 3000
prob.dampcoef = 0.2
Expand Down
8 changes: 4 additions & 4 deletions Source/Diffusion/ComputeQKESourceTerm.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,16 +48,16 @@ ComputeQKESourceTerms (int i, int j, int k,
dvdz = 0.25*(vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1))*dz_inv;
}

// Bouyancy
source_term += 2*CONST_GRAV/theta_mean*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;

// Production
source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);

// Bouyancy
source_term -= 2*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;

// Dissipation
amrex::Real qke = cell_prim(i,j,k,PrimQKE_comp);
if (std::abs(qke) > 0.0) {
source_term += 2.0 * cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
source_term -= 2.0 * cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
(pbl_B1_l * K_turb(i,j,k,EddyDiff::PBL_lengthscale));
}

Expand Down
4 changes: 2 additions & 2 deletions Source/Diffusion/PBLModels.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
// Third Length Scale
amrex::Real l_B;
if (dthetadz > 0) {
amrex::Real N_brunt_vaisala = CONST_GRAV/theta0 * std::sqrt(dthetadz);
amrex::Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
if (zeta < 0) {
amrex::Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T;
qc = std::pow(qc,1.0/3.0);
Expand All @@ -172,7 +172,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
// Compute non-dimensional parameters
amrex::Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
amrex::Real GM = l2_over_q2 * (dudz*dudz + dvdz*dvdz);
amrex::Real GH = -l2_over_q2 / theta0 * dthetadz;
amrex::Real GH = -l2_over_q2 * (CONST_GRAV/theta0) * dthetadz;
amrex::Real E1 = 1.0 + 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH;
amrex::Real E2 = -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH;
amrex::Real E3 = 6.0*A2*A1*GM;
Expand Down