From 78d41f70755298394b572b013ea0752b4bfcaeeb Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Thu, 19 Oct 2023 15:49:59 -0700 Subject: [PATCH 01/15] Fix negative K_turb below surface with PBL. --- Source/Diffusion/PBLModels.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index b6e7d00b3..1e7ac7b6a 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -57,15 +57,15 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, for ( amrex::MFIter mfi(eddyViscosity,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { const amrex::Box &bx = mfi.growntilebox(1); - const amrex::Array4 &cell_data = cons_in.array(mfi); - const amrex::Array4 &cell_data_old = cons_old.array(mfi); - const amrex::Array4 &K_turb = eddyViscosity.array(mfi); + const amrex::Array4 &cell_data = cons_in.array(mfi); + const amrex::Array4 &cell_data_old = cons_old.array(mfi); + const amrex::Array4 &K_turb = eddyViscosity.array(mfi); const amrex::Array4 &uvel = xvel.array(mfi); const amrex::Array4 &vvel = yvel.array(mfi); // Compute some quantities that are constant in each column // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals - // NOTE: Here we requite that sbx covers the entire vertical domain + // Box includes one ghost cell in each direction const amrex::Box &dbx = geom.Domain(); amrex::Box sbx(bx.smallEnd(), bx.bigEnd()); sbx.grow(2,-1); @@ -163,7 +163,8 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, // First Length Scale AMREX_ASSERT(l_obukhov != 0); - const amrex::Real zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2); + int lk = amrex::max(k,0); + const amrex::Real zval = gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2); const amrex::Real zeta = zval/l_obukhov; amrex::Real l_S; if (zeta >= 1.0) { From a41babd78cce664d9406d74782e7407c1fb315ac Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Mon, 16 Oct 2023 20:54:44 -0600 Subject: [PATCH 02/15] Add optional most.surf_heating_rate param --- Source/BoundaryConditions/ABLMost.H | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 77d8780a2..92260cd48 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -54,8 +54,15 @@ public: auto erf_st = pp.query("most.surf_temp", surf_temp); if (erf_st) { theta_type = ThetaCalcType::SURFACE_TEMPERATURE; + pp.query("most.surf_heating_rate", surf_heating_rate); + if (pp.query("most.surf_temp_flux", surf_temp_flux)) { + amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate"); + } } else { pp.query("most.surf_temp_flux", surf_temp_flux); + if (pp.query("most.surf_heating_rate", surf_heating_rate)) { + amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate"); + } if (std::abs(surf_temp_flux) > std::numeric_limits::epsilon()) { theta_type = ThetaCalcType::HEAT_FLUX; } else { @@ -192,6 +199,7 @@ public: private: amrex::Real z0_const; amrex::Real surf_temp; + amrex::Real surf_heating_rate{0}; amrex::Real surf_temp_flux{0}; amrex::Real cnk_a{0.0185}; amrex::Real depth{30.0}; From 08b4e3e1f99c23ce625d504b02d0e2abb50635c1 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Mon, 16 Oct 2023 23:11:22 -0600 Subject: [PATCH 03/15] Update t_surf multifab t_surf = surf_temp (initial input) + surf_heating_rate*cur_time Note: Need to call update_fluxes (for MOST) w/ current sim time so that it can be passed to update_surf_temp --- Source/BoundaryConditions/ABLMost.H | 23 +++++++++++++++++++++-- Source/BoundaryConditions/ABLMost.cpp | 2 ++ Source/ERF.cpp | 2 +- Source/TimeIntegration/ERF_Advance.cpp | 2 +- 4 files changed, 25 insertions(+), 4 deletions(-) diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 92260cd48..ac42e363b 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -54,7 +54,8 @@ public: auto erf_st = pp.query("most.surf_temp", surf_temp); if (erf_st) { theta_type = ThetaCalcType::SURFACE_TEMPERATURE; - pp.query("most.surf_heating_rate", surf_heating_rate); + pp.query("most.surf_heating_rate", surf_heating_rate); // [K/h] + surf_heating_rate = surf_heating_rate / 3600.0; // [K/s] if (pp.query("most.surf_temp_flux", surf_temp_flux)) { amrex::Abort("Can only specify one of surf_temp_flux or surf_heating_rate"); } @@ -136,6 +137,22 @@ public: }// lev } + void + update_surf_temp (amrex::Real cur_time) + { + if (surf_heating_rate != 0) { + int nlevs = m_geom.size(); + for (int lev = 0; lev < nlevs; lev++) + { + t_surf[lev]->setVal(surf_temp + surf_heating_rate*cur_time); + amrex::Print() << "Surface temp at t=" << cur_time + << ": " + << surf_temp + surf_heating_rate*cur_time + << std::endl; + } + } + } + void impose_most_bcs (const int lev, const amrex::Vector& mfs, @@ -149,7 +166,9 @@ public: const FluxCalc& flux_comp); void - update_fluxes (int lev, int max_iters = 25); + update_fluxes (int lev, + amrex::Real cur_time, + int max_iters = 25); template void diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index bc270cb8c..7ce2a6b3d 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -11,6 +11,7 @@ using namespace amrex; */ void ABLMost::update_fluxes (int lev, + Real cur_time, int max_iters) { // Compute plane averages for all vars (regardless of flux type) @@ -30,6 +31,7 @@ ABLMost::update_fluxes (int lev, compute_fluxes(lev, max_iters, most_flux); } } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) { + update_surf_temp(cur_time); if (rough_type == RoughCalcType::CONSTANT) { surface_temp most_flux(m_ma.get_zref(), surf_temp_flux); compute_fluxes(lev, max_iters, most_flux); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2e54eded7..e3a875477 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -540,7 +540,7 @@ ERF::InitData () MultiFab::Copy( *Theta_prim[lev], S, Cons::RhoTheta, 0, 1, ng); MultiFab::Divide(*Theta_prim[lev], S, Cons::Rho , 0, 1, ng); m_most->update_mac_ptrs(lev, vars_new, Theta_prim); - m_most->update_fluxes(lev); + m_most->update_fluxes(lev,t_new[lev]); } } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 9f302b23c..207440853 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -42,7 +42,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ // NOTE: std::swap above causes the field ptrs to be out of date. // Reassign the field ptrs for MAC avg computation. m_most->update_mac_ptrs(lev, vars_old, Theta_prim); - m_most->update_fluxes(lev); + m_most->update_fluxes(lev, time); } } From 7a5cbbe07f61e96daa690d5f1fa57baaef03c379 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Mon, 16 Oct 2023 23:16:56 -0600 Subject: [PATCH 04/15] Add option to specify initial uniform QKE for PBL cases --- Source/ERF.cpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index e3a875477..b3e40f8e9 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -459,6 +459,24 @@ ERF::InitData () } // lev } + Real QKE_0; + if (pp.query("QKE_0", QKE_0)) { + for (int lev = 0; lev <= finest_level; lev++) { + auto& lev_new = vars_new[lev]; + for (MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box &bx = mfi.tilebox(); + const auto &cons_arr = lev_new[Vars::cons].array(mfi); + // We want to set the lateral BC values, too + Box gbx = bx; // Copy constructor + gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + cons_arr(i,j,k,RhoQKE_comp) = cons_arr(i,j,k,Rho_comp) * QKE_0; + }); + } // mfi + } + } + + if (coupling_type == "TwoWay") { AverageDown(); } From a4fa245a834945da1e598e4adedca87e0873836e Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 09:53:20 -0600 Subject: [PATCH 05/15] Handle diabatic case w/ specified temperature --- Source/BoundaryConditions/MOSTStress.H | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 7c021d7d9..c3d5089a8 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -649,7 +649,8 @@ struct moeng_flux amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real num1 = (theta-theta_mean)*wsp_mean; amrex::Real num2 = (theta_mean-theta_surf)*wsp; - amrex::Real moflux = tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean); + amrex::Real moflux = (std::abs(tstar) > EPS) ? + tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0; amrex::Real deltaz = dz * (zlo - k); dest_arr(i,j,k,icomp+n) = rho*(theta - moflux*rho/eta*deltaz); @@ -782,6 +783,7 @@ struct moeng_flux private: int zlo; amrex::Real dz; + static constexpr amrex::Real EPS = 1e-16; }; From eca8e362d6a2d9fe853836ba1b91f728c3c89f9e Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 10:23:44 -0600 Subject: [PATCH 06/15] Add GABLS1 PBL column inputs --- Exec/ABL/inputs_gabls1_mynn25_smag | 74 ++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 Exec/ABL/inputs_gabls1_mynn25_smag diff --git a/Exec/ABL/inputs_gabls1_mynn25_smag b/Exec/ABL/inputs_gabls1_mynn25_smag new file mode 100644 index 000000000..f19ef8c7e --- /dev/null +++ b/Exec/ABL/inputs_gabls1_mynn25_smag @@ -0,0 +1,74 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY (Cuxart et al. 2006) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 + +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.z0 = 0.1 # from Cuxart et al. 2006 +erf.most.zref = 3.125 # == dz/2 +erf.most.surf_temp_flux = 0.0 +#erf.most.surf_temp = 265.0 # initial value, should match input_sounding +#erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 + +zhi.type = "SlipWall" +zhi.theta_grad = 0.01 # [K/m] to match the input sounding + +# INITIALIZATION (Cuxart et al. 2006) +erf.init_type = "input_sounding" +erf.init_sounding_ideal = 1 + +# TIME STEP CONTROL +erf.fixed_dt = 10.0 # fixed time step (Cuxart et al. 2006) +erf.no_substepping = 1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 100 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 1 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta + + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = true + +# Coriolis (Cuxart et al. 2006) +# f = 1.39e-4 s^-1 +erf.use_coriolis = true +erf.latitude = 73.0 +erf.rotational_time_period = 86455.2516813368 + +# Geostrophic wind (Cuxart et al. 2006) +erf.abl_driver_type = "GeostrophicWind" +erf.abl_geo_wind = 8.0 0.0 0.0 + +# Turbulence closure +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 +erf.rho0_trans = 1.3223 # from Cuxart et al. 2006 +erf.theta_ref = 263.5 # from Cuxart et al. 2006 + +erf.pbl_type = "MYNN2.5" +#erf.QKE_0 = 0.8 # == 2*KE_0; note: Cuxart et al. 2006 has an initial TKE profile From 5d827d97bed3a56099ce3622afd320246633a31d Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 12:58:55 -0600 Subject: [PATCH 07/15] Update comments, add code for testing --- Source/Initialization/InputSoundingData.H | 30 +++++++++++++++++++---- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/Source/Initialization/InputSoundingData.H b/Source/Initialization/InputSoundingData.H index 046567b4a..3e292409f 100644 --- a/Source/Initialization/InputSoundingData.H +++ b/Source/Initialization/InputSoundingData.H @@ -129,10 +129,9 @@ public: * density and pressure, then back down to get the dry density and * pressure. * - * This deviates from that code slightly because we calculate the - * correct moist theta, equal to virtual potential temperature, and - * also remove the factor (1+qv) for consistency with the the surface - * density. + * This deviates from that code slightly because we calculate a moist + * theta that is equal to virtual potential temperature, and also + * remove the factor of (1+qv) for consistency with the surface density. */ const int maxiter = 10; const int Ninp = size(); @@ -158,6 +157,18 @@ public: // integrate from surface to domain top amrex::Real qvf, dz; +#if 0 + // In this absence of moisture, this moist profile will match the + // following dry profile + amrex::Print() << "z p_m rho_m theta U V" << std::endl; + amrex::Print() << z_inp_sound[0] + << " " << pm_integ[0] + << " " << rhom_integ[0] + << " " << theta_inp_sound[0] + << " " << U_inp_sound[0] + << " " << V_inp_sound[0] + << std::endl; +#endif for (int k=1; k < size(); ++k) { qvf = 1 + (R_v/R_d-1) * qv_inp_sound[k]; @@ -171,9 +182,18 @@ public: rhom_integ[k] = 1 / ( R_d/p_0 * theta_inp_sound[k]*qvf * std::pow(pm_integ[k]/p_0, -iGamma)); } +#if 0 + amrex::Print() << z_inp_sound[k] + << " " << pm_integ[k] + << " " << rhom_integ[k] + << " " << theta_inp_sound[k] + << " " << U_inp_sound[k] + << " " << V_inp_sound[k] + << std::endl; +#endif } - // now, integrate from the top of the sounding (where it's dry) back + // now, integrate from the top of the sounding (where it's assumed dry) back // down to get the dry air column properties pd_integ[Ninp-1] = pm_integ[Ninp-1]; rhod_integ[Ninp-1] = 1 / ( From 9efdb2c5d416c0833f6af2cfb7d018352faf9aee Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 13:04:03 -0600 Subject: [PATCH 08/15] Add input_sounding input file --- Exec/ABL/input_sounding_GABLS1 | 4 ++++ Exec/ABL/inputs_gabls1_mynn25_smag | 1 + 2 files changed, 5 insertions(+) create mode 100644 Exec/ABL/input_sounding_GABLS1 diff --git a/Exec/ABL/input_sounding_GABLS1 b/Exec/ABL/input_sounding_GABLS1 new file mode 100644 index 000000000..fc02ca975 --- /dev/null +++ b/Exec/ABL/input_sounding_GABLS1 @@ -0,0 +1,4 @@ +1013.2 265.0 0.0 + 0.0 265.0 0.0 8.0 0.0 + 100.0 265.0 0.0 8.0 0.0 + 400.0 268.0 0.0 8.0 0.0 diff --git a/Exec/ABL/inputs_gabls1_mynn25_smag b/Exec/ABL/inputs_gabls1_mynn25_smag index f19ef8c7e..ea1cad0d4 100644 --- a/Exec/ABL/inputs_gabls1_mynn25_smag +++ b/Exec/ABL/inputs_gabls1_mynn25_smag @@ -25,6 +25,7 @@ zhi.theta_grad = 0.01 # [K/m] to match the input sounding # INITIALIZATION (Cuxart et al. 2006) erf.init_type = "input_sounding" erf.init_sounding_ideal = 1 +erf.input_sounding_file = "input_sounding_GABLS1" # TIME STEP CONTROL erf.fixed_dt = 10.0 # fixed time step (Cuxart et al. 2006) From ffeae63039cc1a11287214048de29e3004990b64 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 15:47:20 -0600 Subject: [PATCH 09/15] Manually calculate geostrophic forcing to neglect z component --- Exec/ABL/inputs_gabls1_mynn25_smag | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Exec/ABL/inputs_gabls1_mynn25_smag b/Exec/ABL/inputs_gabls1_mynn25_smag index ea1cad0d4..5012adedb 100644 --- a/Exec/ABL/inputs_gabls1_mynn25_smag +++ b/Exec/ABL/inputs_gabls1_mynn25_smag @@ -46,7 +46,7 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure temp theta # SOLVER CHOICE @@ -54,15 +54,15 @@ erf.alpha_T = 0.0 erf.alpha_C = 1.0 erf.use_gravity = true -# Coriolis (Cuxart et al. 2006) +# Geostrophic wind (Cuxart et al. 2006) # f = 1.39e-4 s^-1 erf.use_coriolis = true -erf.latitude = 73.0 -erf.rotational_time_period = 86455.2516813368 - -# Geostrophic wind (Cuxart et al. 2006) -erf.abl_driver_type = "GeostrophicWind" -erf.abl_geo_wind = 8.0 0.0 0.0 +#erf.latitude = 73.0 +#erf.rotational_time_period = 86455.2516813368 +#erf.abl_driver_type = "GeostrophicWind" +#erf.abl_geo_wind = 8.0 0.0 0.0 +# manually calculate: +erf.abl_geo_forcing = 0.0 0.001112 0.0 # neglect vertical component # Turbulence closure erf.molec_diff_type = "None" From 24c5a634f7fe0c8085d07be90f0669f01c06ab19 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 15:49:04 -0600 Subject: [PATCH 10/15] Print out specified abl_geo_forcing --- Source/DataStructs/DataStruct.H | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index b52ea0556..24be278c1 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -112,8 +112,14 @@ struct SolverChoice { for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i]; amrex::Vector abl_geo_forcing_in = {0.0, 0.0, 0.0}; - pp.queryarr("abl_geo_forcing",abl_geo_forcing_in); - for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_geo_forcing[i] = abl_geo_forcing_in[i]; + if(pp.queryarr("abl_geo_forcing",abl_geo_forcing_in)) { + amrex::Print() << "Specified abl_geo_forcing: ("; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + abl_geo_forcing[i] = abl_geo_forcing_in[i]; + amrex::Print() << abl_geo_forcing[i] << " "; + } + amrex::Print() << ")" << std::endl; + } if (use_coriolis) { From 739cdd3ec48ebdf1a3d680c4e038a629d31c4bb8 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Tue, 17 Oct 2023 23:36:14 -0600 Subject: [PATCH 11/15] Move GABLS1 case to DevTests; init QKE following Cuxart et al 2006 Min value of QKE taken from WRF/phys/module_bl_mynn.F Note: Case description provides TKE, we need to set rho * QKE = rho_hse * (2*TKE) --- Exec/CMakeLists.txt | 1 + Exec/RegTests/GABLS1/CMakeLists.txt | 12 ++ Exec/RegTests/GABLS1/GNUmakefile | 37 ++++ Exec/RegTests/GABLS1/Make.package | 2 + .../GABLS1}/input_sounding_GABLS1 | 0 .../GABLS1}/inputs_gabls1_mynn25_smag | 0 Exec/RegTests/GABLS1/prob.H | 86 +++++++++ Exec/RegTests/GABLS1/prob.cpp | 180 ++++++++++++++++++ 8 files changed, 318 insertions(+) create mode 100644 Exec/RegTests/GABLS1/CMakeLists.txt create mode 100644 Exec/RegTests/GABLS1/GNUmakefile create mode 100644 Exec/RegTests/GABLS1/Make.package rename Exec/{ABL => RegTests/GABLS1}/input_sounding_GABLS1 (100%) rename Exec/{ABL => RegTests/GABLS1}/inputs_gabls1_mynn25_smag (100%) create mode 100644 Exec/RegTests/GABLS1/prob.H create mode 100644 Exec/RegTests/GABLS1/prob.cpp diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 38fe84613..6474d52a0 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -14,6 +14,7 @@ else () add_subdirectory(RegTests/EkmanSpiral_custom) add_subdirectory(RegTests/EkmanSpiral_ideal) add_subdirectory(RegTests/EkmanSpiral_input_sounding) + add_subdirectory(RegTests/GABLS1) add_subdirectory(RegTests/IsentropicVortex) add_subdirectory(RegTests/PoiseuilleFlow) add_subdirectory(RegTests/ScalarAdvDiff) diff --git a/Exec/RegTests/GABLS1/CMakeLists.txt b/Exec/RegTests/GABLS1/CMakeLists.txt new file mode 100644 index 000000000..6d0ec2792 --- /dev/null +++ b/Exec/RegTests/GABLS1/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name erf_gabls1) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) diff --git a/Exec/RegTests/GABLS1/GNUmakefile b/Exec/RegTests/GABLS1/GNUmakefile new file mode 100644 index 000000000..dd835a830 --- /dev/null +++ b/Exec/RegTests/GABLS1/GNUmakefile @@ -0,0 +1,37 @@ +# AMReX +COMP = gnu +PRECISION = DOUBLE + +# Profiling +PROFILE = FALSE +TINY_PROFILE = TRUE +COMM_PROFILE = FALSE +TRACE_PROFILE = FALSE +MEM_PROFILE = FALSE +USE_GPROF = FALSE + +# Performance +USE_MPI = TRUE +USE_OMP = FALSE + +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +# Physics +USE_MOISTURE = FALSE + +# Debugging +DEBUG = FALSE + +TEST = TRUE +USE_ASSERTION = TRUE + +#USE_POISSON_SOLVE = TRUE + +# GNU Make +Bpack := ./Make.package +Blocs := . +ERF_HOME := ../.. +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/RegTests/GABLS1 +include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/RegTests/GABLS1/Make.package b/Exec/RegTests/GABLS1/Make.package new file mode 100644 index 000000000..aeacb72f0 --- /dev/null +++ b/Exec/RegTests/GABLS1/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += prob.H +CEXE_sources += prob.cpp diff --git a/Exec/ABL/input_sounding_GABLS1 b/Exec/RegTests/GABLS1/input_sounding_GABLS1 similarity index 100% rename from Exec/ABL/input_sounding_GABLS1 rename to Exec/RegTests/GABLS1/input_sounding_GABLS1 diff --git a/Exec/ABL/inputs_gabls1_mynn25_smag b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag similarity index 100% rename from Exec/ABL/inputs_gabls1_mynn25_smag rename to Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag diff --git a/Exec/RegTests/GABLS1/prob.H b/Exec/RegTests/GABLS1/prob.H new file mode 100644 index 000000000..e4a8ee5be --- /dev/null +++ b/Exec/RegTests/GABLS1/prob.H @@ -0,0 +1,86 @@ +#ifndef _PROB_H_ +#define _PROB_H_ + +#include + +#include "AMReX_REAL.H" + +#include "prob_common.H" + +struct ProbParm : ProbParmDefaults { + amrex::Real rho_0 = 0.0; + amrex::Real T_0 = 0.0; + amrex::Real A_0 = 1.0; + amrex::Real QKE_0 = 1e-12; // min value + + amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; + amrex::Real W_0 = 0.0; + + // random initial perturbations (legacy code) + amrex::Real U_0_Pert_Mag = 0.0; + amrex::Real V_0_Pert_Mag = 0.0; + amrex::Real W_0_Pert_Mag = 0.0; + amrex::Real T_0_Pert_Mag = 0.0; // perturbation to rho*Theta + + // divergence-free initial perturbations + amrex::Real pert_deltaU = 0.0; + amrex::Real pert_deltaV = 0.0; + amrex::Real pert_periods_U = 5.0; + amrex::Real pert_periods_V = 5.0; + amrex::Real pert_ref_height = 100.0; + + // rayleigh damping + amrex::Real dampcoef = 0.2; // inverse time scale [1/s] + amrex::Real zdamp = 500.0; // damping depth [m] from model top + + // helper vars + amrex::Real aval; + amrex::Real bval; + amrex::Real ufac; + amrex::Real vfac; +}; // namespace ProbParm + +class Problem : public ProblemBase +{ +public: + Problem(const amrex::Real* problo, const amrex::Real* probhi); + +#include "Prob/init_constant_density_hse.H" +#include "Prob/init_rayleigh_damping.H" + + void init_custom_pert ( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& x_vel, + amrex::Array4 const& y_vel, + amrex::Array4 const& z_vel, + amrex::Array4 const& r_hse, + amrex::Array4 const& p_hse, + amrex::Array4 const& z_nd, + amrex::Array4 const& z_cc, +#if defined(ERF_USE_MOISTURE) + amrex::Array4 const& qv, + amrex::Array4 const& qc, + amrex::Array4 const& qi, +#elif defined(ERF_USE_WARM_NO_PRECIP) + amrex::Array4 const& qv, + amrex::Array4 const& qc, +#endif + amrex::GeometryData const& geomdata, + amrex::Array4 const& mf_m, + amrex::Array4 const& mf_u, + amrex::Array4 const& mf_v, + const SolverChoice& sc) override; + +protected: + std::string name() override { return "ABL"; } + +private: + ProbParm parms; +}; + +#endif diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp new file mode 100644 index 000000000..25d70f4f0 --- /dev/null +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -0,0 +1,180 @@ +#include "prob.H" +#include "AMReX_Random.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit(const amrex_real* problo, const amrex_real* probhi) +{ + return std::make_unique(problo, probhi); +} + +Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi) +{ + // Parse params + ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("T_0", parms.T_0); + pp.query("A_0", parms.A_0); + pp.query("QKE_0", parms.QKE_0); + + pp.query("U_0", parms.U_0); + pp.query("V_0", parms.V_0); + pp.query("W_0", parms.W_0); + pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag); + pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag); + pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag); + pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag); + + pp.query("pert_deltaU", parms.pert_deltaU); + pp.query("pert_deltaV", parms.pert_deltaV); + pp.query("pert_periods_U", parms.pert_periods_U); + pp.query("pert_periods_V", parms.pert_periods_V); + pp.query("pert_ref_height", parms.pert_ref_height); + parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]); + parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]); + parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height; + parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height; + + pp.query("dampcoef", parms.dampcoef); + pp.query("zdamp", parms.zdamp); + + init_base_parms(parms.rho_0, parms.T_0); +} + +void +Problem::init_custom_pert( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& state, + amrex::Array4 const& x_vel, + amrex::Array4 const& y_vel, + amrex::Array4 const& z_vel, + amrex::Array4 const& r_hse, + amrex::Array4 const& /*p_hse*/, + amrex::Array4 const& /*z_nd*/, + amrex::Array4 const& /*z_cc*/, +#if defined(ERF_USE_MOISTURE) + amrex::Array4 const& /*qv*/, + amrex::Array4 const& /*qc*/, + amrex::Array4 const& /*qi*/, +#elif defined(ERF_USE_WARM_NO_PRECIP) + amrex::Array4 const& /*qv*/, + amrex::Array4 const& /*qc*/, +#endif + amrex::GeometryData const& geomdata, + amrex::Array4 const& /*mf_m*/, + amrex::Array4 const& /*mf_u*/, + amrex::Array4 const& /*mf_v*/, + const SolverChoice& /*sc*/) +{ + ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]); + const Real yc = 0.5 * (prob_lo[1] + prob_hi[1]); + const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]); + + const Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + + // Add temperature perturbations + if ((z <= parms.pert_ref_height) && (parms.T_0_Pert_Mag != 0.0)) { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + state(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; + } + + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain + state(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r*r); + + // Set an initial value for QKE = 2*e + if (z < 250) { + state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3); + } else { + state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0; + } + +#if defined(ERF_USE_MOISTURE) + state(i, j, k, RhoQt_comp) = 0.0; + state(i, j, k, RhoQp_comp) = 0.0; +#elif defined(ERF_USE_WARM_NO_PRECIP) + state(i, j, k, RhoQv_comp) = 0.0; + state(i, j, k, RhoQc_comp) = 0.0; +#endif + }); + + // Set the x-velocity + ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the x-velocity + x_vel(i, j, k) = parms.U_0; + if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag; + x_vel(i, j, k) += x_vel_prime; + } + if (parms.pert_deltaU != 0.0) + { + const amrex::Real yl = y - prob_lo[1]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + x_vel(i, j, k) += parms.ufac * damp * z * std::cos(parms.aval * yl); + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the y-velocity + y_vel(i, j, k) = parms.V_0; + if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag; + y_vel(i, j, k) += y_vel_prime; + } + if (parms.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl); + } + }); + + // Set the z-velocity + ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const int dom_lo_z = geomdata.Domain().smallEnd()[2]; + const int dom_hi_z = geomdata.Domain().bigEnd()[2]; + + // Set the z-velocity + if (k == dom_lo_z || k == dom_hi_z+1) + { + z_vel(i, j, k) = 0.0; + } + else if (parms.W_0_Pert_Mag != 0.0) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag; + z_vel(i, j, k) = parms.W_0 + z_vel_prime; + } + }); +} + From 19a89d8ae52ceaa8b74e0cff93640aeb0b7922bf Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 18 Oct 2023 14:05:19 -0600 Subject: [PATCH 12/15] Add units to Rayleigh damping screen info --- Source/Prob/init_rayleigh_damping.H | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index f0c550d51..b10d4271c 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -11,8 +11,9 @@ erf_init_rayleigh (amrex::Vector& tau, amrex::Geometry const& geom) override { const auto ztop = geom.ProbHi()[2]; - amrex::Print() << "Rayleigh damping (coef="< Date: Sun, 22 Oct 2023 13:18:12 -0600 Subject: [PATCH 13/15] Clip QKE --- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index 62f7e2020..32579dbfb 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -441,7 +441,8 @@ void erf_slow_rhs_post (int level, // NOTE: we don't include additional source terms when terrain is moving Real temp_val = detJ_arr(i,j,k) * old_cons(i,j,k,n) + dt * detJ_arr(i,j,k) * cell_rhs(i,j,k,n); cur_cons(i,j,k,n) = temp_val / detJ_new_arr(i,j,k); -#ifndef AMREX_USE_GPU + cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); +#if 0 if (cur_cons(i,j,k,n) < Real(0.)) { amrex::AllPrint() << "MAKING NEGATIVE QKE " << IntVect(i,j,k) << " NEW / OLD " << cur_cons(i,j,k,n) << " " << old_cons(i,j,k,n) << std::endl; @@ -482,7 +483,8 @@ void erf_slow_rhs_post (int level, const int n = start_comp + nn; cell_rhs(i,j,k,n) += src_arr(i,j,k,n); cur_cons(i,j,k,n) = old_cons(i,j,k,n) + dt * cell_rhs(i,j,k,n); -#ifndef AMREX_USE_GPU + cur_cons(i,j,k,n) = amrex::max(cur_cons(i,j,k,n), 1e-12); +#if 0 if (cur_cons(i,j,k,n) < Real(0.)) { amrex::AllPrint() << "MAKING NEGATIVE QKE " << IntVect(i,j,k) << " NEW / OLD " << cur_cons(i,j,k,n) << " " << old_cons(i,j,k,n) << std::endl; From a60b65ec47b19ed8e6ce84f4d557a02998d5ff60 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Sun, 22 Oct 2023 13:20:57 -0600 Subject: [PATCH 14/15] Add output options for diffusivities --- Docs/sphinx_doc/Plotfiles.rst | 16 ++++++++++++++++ Source/ERF.H | 7 ++++++- Source/IO/Plotfile.cpp | 17 +++++++++++++++++ 3 files changed, 39 insertions(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst index 4626358b2..59d18c183 100644 --- a/Docs/sphinx_doc/Plotfiles.rst +++ b/Docs/sphinx_doc/Plotfiles.rst @@ -209,6 +209,22 @@ Output Options | | | | | | +-----------------------------+------------------+ +| **Kmv** | Vertical | +| | Eddy Diffusivity | +| | of Momentum | ++-----------------------------+------------------+ +| **Kmh** | Horizontal | +| | Eddy Diffusivity | +| | of Momentum | ++-----------------------------+------------------+ +| **Khv** | Vertical | +| | Eddy Diffusivity | +| | of Heat | ++-----------------------------+------------------+ +| **Khh** | Horizontal | +| | Eddy Diffusivity | +| | of Heat | ++-----------------------------+------------------+ | **qt** | Nonprecipitating | | | water (qv + qc + | | | qi) | diff --git a/Source/ERF.H b/Source/ERF.H index 02beaddb7..efcfd4694 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -651,7 +651,12 @@ private: // Note that the order of variable names here must match the order in Derive.cpp const amrex::Vector derived_names {"pressure", "soundspeed", "temp", "theta", "KE", "QKE", "scalar", "pres_hse", "dens_hse", "pert_pres", "pert_dens", - "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac" + "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", + "z_phys", "detJ" , "mapfac", + // eddy viscosity + "Kmv","Kmh", + // eddy diffusivity of heat + "Khv","Khh" #if defined(ERF_USE_MOISTURE) ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup" #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 144599879..5763a585a 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -544,6 +544,23 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } + if (containerHasElement(plot_var_names, "Kmv")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_v,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Kmh")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Mom_h,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Khv")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_v,mf_comp,1,0); + mf_comp ++; + } + if (containerHasElement(plot_var_names, "Khh")) { + MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_h,mf_comp,1,0); + mf_comp ++; + } + #if defined(ERF_USE_MOISTURE) calculate_derived("qt", derived::erf_derQt); calculate_derived("qp", derived::erf_derQp); From 14da0ba0b8823e679147345b63bd393f94166486 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Sun, 22 Oct 2023 13:32:33 -0600 Subject: [PATCH 15/15] Update inputs for GABLS1 case --- .../RegTests/GABLS1/inputs_gabls1_mynn25_smag | 42 ++++++++++--------- Exec/RegTests/GABLS1/prob.cpp | 1 + Source/Diffusion/PBLModels.cpp | 2 +- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag index 5012adedb..f31d16159 100644 --- a/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag +++ b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag @@ -15,9 +15,8 @@ geometry.is_periodic = 1 1 0 zlo.type = "Most" erf.most.z0 = 0.1 # from Cuxart et al. 2006 erf.most.zref = 3.125 # == dz/2 -erf.most.surf_temp_flux = 0.0 -#erf.most.surf_temp = 265.0 # initial value, should match input_sounding -#erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 +erf.most.surf_temp = 265.0 # initial value, should match input_sounding +erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 zhi.type = "SlipWall" zhi.theta_grad = 0.01 # [K/m] to match the input sounding @@ -28,8 +27,10 @@ erf.init_sounding_ideal = 1 erf.input_sounding_file = "input_sounding_GABLS1" # TIME STEP CONTROL -erf.fixed_dt = 10.0 # fixed time step (Cuxart et al. 2006) -erf.no_substepping = 1 +#erf.fixed_dt = 10.0 # fixed time step (Cuxart et al. 2006) +#erf.no_substepping = 1 +erf.fixed_dt = 1.0 # largest stable low Mach dt +erf.fixed_mri_dt_ratio = 6 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -41,35 +42,36 @@ amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file -erf.check_int = 100 # number of timesteps between checkpoints +erf.check_int = 600 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure temp theta +erf.plot_int_1 = 600 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv # SOLVER CHOICE -erf.alpha_T = 0.0 -erf.alpha_C = 1.0 +erf.dycore_vert_adv_type = "Upwind_3rd" + +erf.molec_diff_type = "None" +#erf.alpha_T = 0.0 +#erf.alpha_C = 1.0 + erf.use_gravity = true -# Geostrophic wind (Cuxart et al. 2006) -# f = 1.39e-4 s^-1 +# Coriolis parameter f = 1.39e-4 s^-1 (Cuxart et al. 2006) erf.use_coriolis = true -#erf.latitude = 73.0 -#erf.rotational_time_period = 86455.2516813368 -#erf.abl_driver_type = "GeostrophicWind" -#erf.abl_geo_wind = 8.0 0.0 0.0 -# manually calculate: -erf.abl_geo_forcing = 0.0 0.001112 0.0 # neglect vertical component +erf.latitude = 73.0 +erf.rotational_time_period = 86455.2516813368 + +# Geostrophic wind (Cuxart et al. 2006) +erf.abl_driver_type = "GeostrophicWind" +erf.abl_geo_wind = 8.0 0.0 0.0 # Turbulence closure -erf.molec_diff_type = "None" erf.les_type = "Smagorinsky" erf.Cs = 0.1 erf.rho0_trans = 1.3223 # from Cuxart et al. 2006 erf.theta_ref = 263.5 # from Cuxart et al. 2006 erf.pbl_type = "MYNN2.5" -#erf.QKE_0 = 0.8 # == 2*KE_0; note: Cuxart et al. 2006 has an initial TKE profile diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp index 25d70f4f0..c9bffb30f 100644 --- a/Exec/RegTests/GABLS1/prob.cpp +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -97,6 +97,7 @@ Problem::init_custom_pert( // Set an initial value for QKE = 2*e if (z < 250) { + // Following Cuxart et al. 2006 state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3); } else { state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0; diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 1e7ac7b6a..c890a59d4 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -227,7 +227,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel, amrex::Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4); amrex::Real SH = (R1*E3 - R2*E1)/(E2*E3 - E1*E4); - amrex::Real SQ = 3.0 * SM; + amrex::Real SQ = 3.0 * SM; // Nakanishi & Niino 2009 // Finally, compute the eddy viscosity/diffusivities const amrex::Real rho = cell_data(i,j,k,Rho_comp);