Skip to content

Commit

Permalink
Implement terrain-aware input-sounding initialization
Browse files Browse the repository at this point in the history 
Notes:
- Nominal z levels (e.g., from input terrain_z_levels) that are used for
  terrain smoothing should range from 0 (on the surface) to ztop
- The same levels are used to calculate the base state by integrating
  the hydrostatic equation through a column of air
- In the case of terrain or grid stretching, the base state profiles
  may have non-uniform vertical grid spacing
- These 1-D profiles are then interpolated to the 3-D terrain grid; the
  state and hse fab arrays are interpolated at z_phys_cc whereas the
  velocity fab arrays are interpolated at face centers calculated from
  z_phys_nd
  • Loading branch information
@ewquon
ewquon committed Jan 9, 2024
1 parent 343ebff commit 302d57e
Show file tree
Hide file tree
Showing 4 changed files with 62 additions and 26 deletions.
2 changes: 1 addition & 1 deletion Source/Initialization/ERF_init1d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ ERF::setRayleighRefFromSounding (bool restarting)
// so we need to read it here
// TODO: should we store this information in the checkpoint file instead?
if (restarting) {
input_sounding_data.read_from_file(input_sounding_file, geom[0]);
input_sounding_data.read_from_file(input_sounding_file, geom[0], zlevels_stag);
}

const Real* z_inp_sound = input_sounding_data.z_inp_sound.dataPtr();
Expand Down
1 change: 1 addition & 0 deletions Source/Initialization/ERF_init_custom.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@

/**
* \file ERF_init_custom.cpp
*/
Expand Down
69 changes: 49 additions & 20 deletions Source/Initialization/ERF_init_from_input_sounding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ void
init_bx_scalars_from_input_sounding (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &state,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_cc_arr,
const bool& l_moist,
InputSoundingData const &inputSoundingData);
void
Expand All @@ -23,6 +24,7 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &p_hse_arr,
amrex::Array4<amrex::Real> const &pi_hse_arr,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_cc_arr,
const amrex::Real& l_gravity,
const amrex::Real& l_rdOcp,
const bool& l_moist,
Expand All @@ -34,6 +36,7 @@ init_bx_velocities_from_input_sounding (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &y_vel,
amrex::Array4<amrex::Real> const &z_vel,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_nd_arr,
InputSoundingData const &inputSoundingData);

/**
Expand All @@ -50,8 +53,12 @@ ERF::init_from_input_sounding (int lev)
if (input_sounding_file.empty())
amrex::Error("input_sounding file name must be provided via input");

input_sounding_data.read_from_file(input_sounding_file, geom[lev]);
// this will interpolate the input profiles to the nominal height levels
// (ranging from 0 to the domain top)
input_sounding_data.read_from_file(input_sounding_file, geom[lev], zlevels_stag);

// this will calculate the hydrostatically balanced density and pressure
// profiles following WRF ideal.exe
if (init_sounding_ideal) input_sounding_data.calc_rho_p();
}

Expand Down Expand Up @@ -79,25 +86,32 @@ ERF::init_from_input_sounding (int lev)
Array4<Real> p_hse_arr = p_hse.array(mfi);
Array4<Real> pi_hse_arr = pi_hse.array(mfi);

Array4<Real const> z_cc_arr = (solverChoice.use_terrain) ? z_phys_cc[lev]->const_array(mfi) : Array4<Real const>{};
Array4<Real const> z_nd_arr = (solverChoice.use_terrain) ? z_phys_nd[lev]->const_array(mfi) : Array4<Real const>{};

if (init_sounding_ideal)
{
// HSE will be calculated here, following WRF ideal.exe
// HSE will be initialized here, interpolated from values previously
// calculated by calc_rho_p()
init_bx_scalars_from_input_sounding_hse(
bx, cons_arr,
r_hse_arr, p_hse_arr, pi_hse_arr,
geom[lev].data(), l_gravity, l_rdOcp, l_moist, input_sounding_data);
geom[lev].data(), z_cc_arr,
l_gravity, l_rdOcp, l_moist, input_sounding_data);
}
else
{
// HSE will be calculated later with call to initHSE
init_bx_scalars_from_input_sounding(
bx, cons_arr,
geom[lev].data(), l_moist, input_sounding_data);
geom[lev].data(), z_cc_arr,
l_moist, input_sounding_data);
}

init_bx_velocities_from_input_sounding(
bx, xvel_arr, yvel_arr, zvel_arr,
geom[lev].data(), input_sounding_data);
geom[lev].data(), z_nd_arr,
input_sounding_data);

} //mfi
}
Expand All @@ -115,6 +129,7 @@ void
init_bx_scalars_from_input_sounding (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &state,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_cc_arr,
const bool& l_moist,
InputSoundingData const &inputSoundingData)
{
Expand All @@ -123,15 +138,17 @@ init_bx_scalars_from_input_sounding (const amrex::Box &bx,
const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d.dataPtr();
const int inp_sound_size = inputSoundingData.size();

// Geometry
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();

// 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 {
// Geometry
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();
const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
const amrex::Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
: prob_lo[2] + (k + 0.5) * dx[2];

amrex::Real rho_0 = 1.0;

Expand Down Expand Up @@ -172,6 +189,7 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &p_hse_arr,
amrex::Array4<amrex::Real> const &pi_hse_arr,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_cc_arr,
const amrex::Real& l_gravity,
const amrex::Real& l_rdOcp,
const bool& l_moist,
Expand All @@ -183,16 +201,18 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx,
const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d.dataPtr();
const int inp_sound_size = inputSoundingData.size();

// Geometry
int ktop = bx.bigEnd(2);
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();

// 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 {
// Geometry
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();
const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
int ktop = bx.bigEnd(2);
const amrex::Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
: prob_lo[2] + (k + 0.5) * dx[2];

Real rho_k, rhoTh_k;

Expand Down Expand Up @@ -253,13 +273,18 @@ init_bx_velocities_from_input_sounding (const amrex::Box &bx,
amrex::Array4<amrex::Real> const &y_vel,
amrex::Array4<amrex::Real> const &z_vel,
amrex::GeometryData const &geomdata,
amrex::Array4<const amrex::Real> const &z_nd_arr,
InputSoundingData const &inputSoundingData)
{
const Real* z_inp_sound = inputSoundingData.z_inp_sound_d.dataPtr();
const Real* U_inp_sound = inputSoundingData.U_inp_sound_d.dataPtr();
const Real* V_inp_sound = inputSoundingData.V_inp_sound_d.dataPtr();
const int inp_sound_size = inputSoundingData.size();

// Geometry
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();

// 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
Expand All @@ -275,18 +300,22 @@ init_bx_velocities_from_input_sounding (const amrex::Box &bx,
amrex::ParallelFor(xbx, ybx, zbx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
// Note that this is called on a box of x-faces
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();
const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
const amrex::Real z = (z_nd_arr) ? 0.25*( z_nd_arr(i,j ,k )
+ z_nd_arr(i,j+1,k )
+ z_nd_arr(i,j ,k+1)
+ z_nd_arr(i,j+1,k+1))
: prob_lo[2] + (k + 0.5) * dx[2];

// Set the x-velocity
x_vel(i, j, k) = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
// Note that this is called on a box of y-faces
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();
const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
const amrex::Real z = (z_nd_arr) ? 0.25*( z_nd_arr(i ,j,k )
+ z_nd_arr(i+1,j,k )
+ z_nd_arr(i ,j,k+1)
+ z_nd_arr(i+1,j,k+1))
: prob_lo[2] + (k + 0.5) * dx[2];

// Set the y-velocity
y_vel(i, j, k) = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size);
Expand Down
16 changes: 11 additions & 5 deletions Source/Initialization/InputSoundingData.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,17 @@ public:
InputSoundingData () {}

void read_from_file (const std::string input_sounding_file,
amrex::Geometry const &geom)
const amrex::Geometry &geom,
const amrex::Vector<amrex::Real>& zlevels_stag)
{
const amrex::Real zbot = geom.ProbLo(AMREX_SPACEDIM-1);
const amrex::Real ztop = geom.ProbHi(AMREX_SPACEDIM-1);
const amrex::Real dz = geom.CellSize(AMREX_SPACEDIM-1);
const int klo = 0;
const int khi = geom.Domain().bigEnd()[AMREX_SPACEDIM-1];
const int Nz = geom.Domain().size()[AMREX_SPACEDIM-1];
const amrex::Real dz = geom.CellSize()[AMREX_SPACEDIM-1];

const bool use_terrain = (zlevels_stag.size() > 0);
const amrex::Real zbot = (use_terrain) ? zlevels_stag[klo] : geom.ProbLo(AMREX_SPACEDIM-1);
const amrex::Real ztop = (use_terrain) ? zlevels_stag[khi+1] : geom.ProbHi(AMREX_SPACEDIM-1);

z_inp_sound.resize(Nz+2);
theta_inp_sound.resize(Nz+2);
Expand Down Expand Up @@ -95,7 +100,8 @@ public:
U_inp_sound[0] = U_inp_sound_tmp[0];
V_inp_sound[0] = V_inp_sound_tmp[0];
for (int k=0; k < Nz; ++k) {
z_inp_sound[k+1] = zbot + (k + 0.5) * dz;
z_inp_sound[k+1] = (use_terrain) ? 0.5 * (zlevels_stag[k] + zlevels_stag[k+1])
: zbot + (k + 0.5) * dz;
theta_inp_sound[k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), theta_inp_sound_tmp.dataPtr(), z_inp_sound[k+1], Ninp);
qv_inp_sound[k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), qv_inp_sound_tmp.dataPtr(), z_inp_sound[k+1], Ninp);
U_inp_sound[k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), U_inp_sound_tmp.dataPtr(), z_inp_sound[k+1], Ninp);
Expand Down

0 comments on commit 302d57e

Please sign in to comment.