Skip to content

Commit

Permalink
Add option to input a constant geostrophic wind profile, update docs (e…
Browse files Browse the repository at this point in the history
…rf-model#1720)

* Cleanup comments

* Add additional abl_geo_wind_table param and associated function stub

* Update function signature, catch undefined terrain use case

* Allow "terrain" for grid stretching

* Implement init_geo_wind_profile

* Provide a bit of verbose output

* Describe large-scale forcing options in docs

* Make compiler happy

* Fix name conflict.

* missed a name conflict.

---------

Co-authored-by: Aaron M. Lattanzi <[email protected]>
Co-authored-by: AMLattanzi <[email protected]>
  • Loading branch information
3 people authored Aug 6, 2024
1 parent 1330413 commit 6b1a193
Show file tree
Hide file tree
Showing 8 changed files with 145 additions and 25 deletions.
31 changes: 31 additions & 0 deletions Docs/sphinx_doc/Inputs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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, | | |
Expand Down Expand Up @@ -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
==============
Expand Down
23 changes: 16 additions & 7 deletions Source/DataStructs/DataStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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<amrex::Real> 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;
}
}
}

Expand Down Expand Up @@ -544,6 +551,8 @@ struct SolverChoice {
ABLDriverType abl_driver_type;
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_pressure_grad;
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_geo_forcing;
std::string abl_geo_wind_table;
bool have_geo_wind_profile {false};

int ave_plane {2};
// Microphysics params
Expand Down
9 changes: 9 additions & 0 deletions Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -982,6 +982,15 @@ private:
amrex::Vector< amrex::Vector<amrex::Real> > h_v_geos;
amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real> > 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<amrex::Real>& u_geos,
amrex::Gpu::DeviceVector<amrex::Real>& u_geos_d,
amrex::Vector<amrex::Real>& v_geos,
amrex::Gpu::DeviceVector<amrex::Real>& v_geos_d,
const amrex::Geometry& geom,
const amrex::Vector<amrex::Real>& zlev_stag);

// This is a vector over levels of vectors across quantities of Vectors
amrex::Vector<amrex::Vector< amrex::Vector<amrex::Real> > > h_rayleigh_ptrs;

Expand Down
21 changes: 16 additions & 5 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real>(0));
d_u_geos.resize(max_level+1, Gpu::DeviceVector<Real>(0));
Expand All @@ -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);
}
}
}

Expand Down
58 changes: 58 additions & 0 deletions Source/Initialization/ERF_init1d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
#include <prob_common.H>
#include <Utils/ParFunctions.H>

#include <Interpolation_1D.H>

using namespace amrex;

/**
Expand Down Expand Up @@ -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<Real>& u_geos,
Gpu::DeviceVector<Real>& u_geos_d,
Vector<Real>& v_geos,
Gpu::DeviceVector<Real>& v_geos_d,
const Geometry& geom,
const Vector<Real>& 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<Real> 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();
}
2 changes: 1 addition & 1 deletion Source/SourceTerms/ERF_make_mom_sources.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions Source/TimeIntegration/TI_slow_rhs_fun.H
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
18 changes: 10 additions & 8 deletions Source/prob_common.H
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
}
}
Expand All @@ -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
Expand All @@ -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;
}
}
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down

0 comments on commit 6b1a193

Please sign in to comment.