diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 01c4e14ec..03de35889 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -919,6 +919,15 @@ List of Parameters | | abl.driver_type = | | | | | GeostrophicWind) | | | +----------------------------------+-------------------+-------------------+-------------+ +| **erf.abl_geo_wind_table** | Path to text file | String | None | +| | containing a | | | +| | geostrophic wind | | | +| | profile | | | +| | (with z, Ug, and | | | +| | Vg whtiespace | | | +| | delimited | | | +| | columns) | | | ++----------------------------------+-------------------+-------------------+-------------+ | **erf.use_gravity** | Include gravity | true / false | false | | | in momentum | | | | | update? If true, | | | @@ -955,6 +964,28 @@ function(s). | | temperature source| | | | | term | | | +--------------------------------------------+-------------------+-------------------+-------------+ +| **erf.add_custom_moisture_forcing** | Apply the | true or false | false | +| | user-defined | | | +| | qv source | | | +| | term | | | ++--------------------------------------------+-------------------+-------------------+-------------+ +| **erf.add_custom_w_subsidence** | Apply the | true or false | false | +| | user-defined | | | +| | vertical velocity | | | +| | profile for use in| | | +| | calculating | | | +| | subsidence source | | | +| | terms | | | ++--------------------------------------------+-------------------+-------------------+-------------+ +| **erf.add_custom_geostrophic_profile** | Apply the | true or false | false | +| | user-defined | | | +| | geostrophic wind | | | +| | profile | | | ++--------------------------------------------+-------------------+-------------------+-------------+ + +Note that ``erf.add_custom_geostrophic_profile`` cannot be used in combination +with an ``erf.abl_geo_wind_table``. + Initialization ============== diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index a12eeea65..5ac155f58 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -279,6 +279,10 @@ struct SolverChoice { pp.query("add_custom_geostrophic_profile", custom_geostrophic_profile); pp.query("custom_forcing_uses_primitive_vars", custom_forcing_prim_vars); + have_geo_wind_profile = (!abl_geo_wind_table.empty() || custom_geostrophic_profile); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(!abl_geo_wind_table.empty() && custom_geostrophic_profile), + "Should not have both abl_geo_wind_table and custom_geostrophic_profile set."); + pp.query("Ave_Plane", ave_plane); pp.query("mp_clouds", do_cloud); @@ -441,18 +445,21 @@ struct SolverChoice { amrex::Print() << "Coriolis frequency, f = " << coriolis_factor * sinphi << " 1/s" << std::endl; - if (abl_driver_type == ABLDriverType::GeostrophicWind) - { + if (abl_driver_type == ABLDriverType::GeostrophicWind) { // Read in the geostrophic wind -- we only use this to construct // the forcing term so no need to keep it amrex::Vector abl_geo_wind(3); pp.queryarr("abl_geo_wind",abl_geo_wind); - abl_geo_forcing = { - -coriolis_factor * (abl_geo_wind[1]*sinphi - abl_geo_wind[2]*cosphi), - coriolis_factor * abl_geo_wind[0]*sinphi, - -coriolis_factor * abl_geo_wind[0]*cosphi - }; + if(!pp.query("abl_geo_wind_table",abl_geo_wind_table)) { + abl_geo_forcing = { + -coriolis_factor * (abl_geo_wind[1]*sinphi - abl_geo_wind[2]*cosphi), + coriolis_factor * abl_geo_wind[0]*sinphi, + -coriolis_factor * abl_geo_wind[0]*cosphi + }; + } else { + amrex::Print() << "NOTE: abl_geo_wind_table provided, ignoring input abl_geo_wind" << std::endl; + } } } @@ -544,6 +551,8 @@ struct SolverChoice { ABLDriverType abl_driver_type; amrex::GpuArray abl_pressure_grad; amrex::GpuArray abl_geo_forcing; + std::string abl_geo_wind_table; + bool have_geo_wind_profile {false}; int ave_plane {2}; // Microphysics params diff --git a/Source/ERF.H b/Source/ERF.H index 01b6d2eeb..95230c662 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -982,6 +982,15 @@ private: amrex::Vector< amrex::Vector > h_v_geos; amrex::Vector > d_v_geos; + // Function to read and populate above host vectors (if input file exists) + void init_geo_wind_profile (const std::string input_file, + amrex::Vector& u_geos, + amrex::Gpu::DeviceVector& u_geos_d, + amrex::Vector& v_geos, + amrex::Gpu::DeviceVector& v_geos_d, + const amrex::Geometry& geom, + const amrex::Vector& zlev_stag); + // This is a vector over levels of vectors across quantities of Vectors amrex::Vector > > h_rayleigh_ptrs; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 0d2193f07..71c0cf43d 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -750,7 +750,7 @@ ERF::InitData () } } - if (solverChoice.custom_geostrophic_profile) + if (solverChoice.have_geo_wind_profile) { h_u_geos.resize(max_level+1, Vector(0)); d_u_geos.resize(max_level+1, Gpu::DeviceVector(0)); @@ -762,10 +762,21 @@ ERF::InitData () d_u_geos[lev].resize(domlen, 0.0_rt); h_v_geos[lev].resize(domlen, 0.0_rt); d_v_geos[lev].resize(domlen, 0.0_rt); - prob->update_geostrophic_profile(t_new[0], - h_u_geos[lev], d_u_geos[lev], - h_v_geos[lev], d_v_geos[lev], - geom[lev], z_phys_cc[lev]); + if (solverChoice.custom_geostrophic_profile) { + prob->update_geostrophic_profile(t_new[0], + h_u_geos[lev], d_u_geos[lev], + h_v_geos[lev], d_v_geos[lev], + geom[lev], z_phys_cc[lev]); + } else { + if (solverChoice.use_terrain > 0) { + amrex::Print() << "Note: 1-D geostrophic wind profile input is only defined for grid stretching, not real terrain" << std::endl; + } + init_geo_wind_profile(solverChoice.abl_geo_wind_table, + h_u_geos[lev], d_u_geos[lev], + h_v_geos[lev], d_v_geos[lev], + geom[lev], + zlevels_stag); + } } } diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index f86cfb99d..c87ce0bb0 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -7,6 +7,8 @@ #include #include +#include + using namespace amrex; /** @@ -403,3 +405,59 @@ ERF::erf_enforce_hse (int lev, dens.FillBoundary(geom[lev].periodicity()); pres.FillBoundary(geom[lev].periodicity()); } + +void ERF::init_geo_wind_profile(const std::string input_file, + Vector& u_geos, + Gpu::DeviceVector& u_geos_d, + Vector& v_geos, + Gpu::DeviceVector& v_geos_d, + const Geometry& geom, + const Vector& zlev_stag) +{ + const int klo = 0; + const int khi = geom.Domain().bigEnd()[AMREX_SPACEDIM-1]; + const amrex::Real dz = geom.CellSize()[AMREX_SPACEDIM-1]; + + const bool grid_stretch = (zlev_stag.size() > 0); + const Real zbot = (grid_stretch) ? zlev_stag[klo] : geom.ProbLo(AMREX_SPACEDIM-1); + const Real ztop = (grid_stretch) ? zlev_stag[khi+1] : geom.ProbHi(AMREX_SPACEDIM-1); + + amrex::Print() << "Reading geostrophic wind profile from " << input_file << std::endl; + std::ifstream profile_reader(input_file); + if(!profile_reader.is_open()) { + amrex::Error("Error opening the abl_geo_wind_table\n"); + } + + // First, read the input data into temp vectors + std::string line; + Vector z_inp, Ug_inp, Vg_inp; + Real z, Ug, Vg; + amrex::Print() << "z Ug Vg" << std::endl; + while(std::getline(profile_reader, line)) { + std::istringstream iss(line); + iss >> z >> Ug >> Vg; + amrex::Print() << z << " " << Ug << " " << Vg << std::endl; + z_inp.push_back(z); + Ug_inp.push_back(Ug); + Vg_inp.push_back(Vg); + if (z >= ztop) break; + } + + const int Ninp = z_inp.size(); + AMREX_ALWAYS_ASSERT(z_inp[0] <= zbot); + AMREX_ALWAYS_ASSERT(z_inp[Ninp-1] >= ztop); + + // Now, interpolate vectors to the cell centers + for (int k = 0; k <= khi; k++) { + z = (grid_stretch) ? 0.5 * (zlev_stag[k] + zlev_stag[k+1]) + : zbot + (k + 0.5) * dz; + u_geos[k] = interpolate_1d(z_inp.dataPtr(), Ug_inp.dataPtr(), z, Ninp); + v_geos[k] = interpolate_1d(z_inp.dataPtr(), Vg_inp.dataPtr(), z, Ninp); + } + + // Copy from host version to device version + Gpu::copy(Gpu::hostToDevice, u_geos.begin(), u_geos.end(), u_geos_d.begin()); + Gpu::copy(Gpu::hostToDevice, v_geos.begin(), v_geos.end(), v_geos_d.begin()); + + profile_reader.close(); +} diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index dc65a054c..563f94dc1 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -93,7 +93,7 @@ void make_mom_sources (int /*level*/, // Flag for Geostrophic forcing // ***************************************************************************** auto abl_geo_forcing = solverChoice.abl_geo_forcing; - auto geo_wind_profile = solverChoice.custom_geostrophic_profile; + auto geo_wind_profile = solverChoice.have_geo_wind_profile; // ***************************************************************************** // Data for Rayleigh damping diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index dc42e9082..393b9a3ad 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -33,8 +33,8 @@ } } - Real* dptr_u_geos = solverChoice.custom_geostrophic_profile ? d_u_geos[level].data(): nullptr; - Real* dptr_v_geos = solverChoice.custom_geostrophic_profile ? d_v_geos[level].data(): nullptr; + Real* dptr_u_geos = solverChoice.have_geo_wind_profile ? d_u_geos[level].data(): nullptr; + Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr; // Construct the source terms for the cell-centered (conserved) variables make_sources(level, nrk, slow_dt, S_data, S_prim, cc_src, @@ -404,8 +404,8 @@ } } - Real* dptr_u_geos = solverChoice.custom_geostrophic_profile ? d_u_geos[level].data(): nullptr; - Real* dptr_v_geos = solverChoice.custom_geostrophic_profile ? d_v_geos[level].data(): nullptr; + Real* dptr_u_geos = solverChoice.have_geo_wind_profile ? d_u_geos[level].data(): nullptr; + Real* dptr_v_geos = solverChoice.have_geo_wind_profile ? d_v_geos[level].data(): nullptr; make_sources(level, nrk, slow_dt, S_data, S_prim, cc_src, #if defined(ERF_USE_RRTMGP) diff --git a/Source/prob_common.H b/Source/prob_common.H index 68e49196e..8b65314d9 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -142,7 +142,7 @@ public: } /** - * Function to update user-specified temperature source terms that can + * Function to update user-specified moisture source terms that can * vary with time and height. * * @param[in] time current time @@ -169,7 +169,7 @@ public: for (int k = 0; k <= khi; k++) { // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2]; - // set RHS term of RhoTheta equation based on time, z_cc here + // set RHS term of RhoQ1 equation based on time, z_cc here qsrc[k] = 0.0; } } @@ -179,8 +179,11 @@ public: } /** - * Function to update user-specified temperature source terms that can - * vary with time and height. + * Function to update the vertical velocity profile, used to add subsidence + * source terms for x-mom, y-mom, rho*theta, rho*Q1, and rho*Q2. + * + * TODO: Currently, this is only called by InitData, so there is no time + * dependence. * * @param[in] time current time * @param[out] wbar w vel forcing profile @@ -206,7 +209,7 @@ public: for (int k = 0; k <= khi; k++) { // const amrex::Real z_cc = prob_lo[2] + (k+0.5)* dx[2]; - // set RHS term of RhoTheta equation based on time, z_cc here + // set vertical velocity profile based on time, z_cc here wbar[k] = 0.0; } } @@ -216,8 +219,7 @@ public: } /** - * Function to update user-specified temperature source terms that can - * vary with time and height. + * Function to update user-specified geostrophic wind profile. * * @param[in] time current time * @param[out] u_geos geostrophic wind profile @@ -257,7 +259,7 @@ public: amrex::Gpu::copy(amrex::Gpu::hostToDevice, v_geos.begin(), v_geos.end(), d_v_geos.begin()); } - /** + /** * Function to perform custom initialization of terrain * * Note: Terrain functionality can also be used to provide grid stretching.