From 50551e7d420a995d7ea9fbd939c0ec5fc6a0ff98 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 19 Sep 2023 12:39:14 -0700 Subject: [PATCH 1/3] This runs multiple steps without a negative QKE and the field looks reasonable. --- Exec/DevTests/MiguelDev/inputs | 18 +++++++++--------- Source/Diffusion/ComputeQKESourceTerm.H | 8 ++++---- Source/Diffusion/PBLModels.cpp | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/Exec/DevTests/MiguelDev/inputs b/Exec/DevTests/MiguelDev/inputs index da357a28c..69a5771df 100644 --- a/Exec/DevTests/MiguelDev/inputs +++ b/Exec/DevTests/MiguelDev/inputs @@ -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 @@ -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 @@ -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" @@ -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 diff --git a/Source/Diffusion/ComputeQKESourceTerm.H b/Source/Diffusion/ComputeQKESourceTerm.H index f592f7921..1812f0d46 100644 --- a/Source/Diffusion/ComputeQKESourceTerm.H +++ b/Source/Diffusion/ComputeQKESourceTerm.H @@ -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,3.0) / (pbl_B1_l * K_turb(i,j,k,EddyDiff::PBL_lengthscale)); } diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index a29b26fc2..a3f96e7cc 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -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; From 534c0268d63ff06e8c16c594dfbf264b66518e3d Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 19 Sep 2023 12:58:16 -0700 Subject: [PATCH 2/3] Fix debug modification. --- Source/Diffusion/ComputeQKESourceTerm.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Diffusion/ComputeQKESourceTerm.H b/Source/Diffusion/ComputeQKESourceTerm.H index 1812f0d46..965c9cce4 100644 --- a/Source/Diffusion/ComputeQKESourceTerm.H +++ b/Source/Diffusion/ComputeQKESourceTerm.H @@ -57,7 +57,7 @@ ComputeQKESourceTerms (int i, int j, int k, // 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,3.0) / + 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)); } From f50fd2ac06789f863b92f351d4940137d88a62d6 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 19 Sep 2023 14:25:14 -0600 Subject: [PATCH 3/3] fix sqrt in brunt vaisala freq --- Source/Diffusion/PBLModels.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index a3f96e7cc..66cfa072f 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -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);