From e47a0f3418ddce1d7da0a2bd57390c64673e3ffe Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 7 Feb 2024 10:03:55 -0700 Subject: [PATCH] Properly initialize reference values for Rayleigh damping if using Rayleigh + input_sounding + terrain (#1429) * Remove duplicative print statement * Clarify comment * Properly initialize rayleigh from sounding when there's terrain * Remove out-of-date notification --- Source/ERF.cpp | 8 +++----- Source/Initialization/ERF_init1d.cpp | 28 ++++++++++++++++++++-------- Source/Prob/init_rayleigh_damping.H | 18 +++--------------- Source/Utils/Make.package | 2 ++ Source/Utils/ParFunctions.H | 26 ++++++++++++++++++++++++++ 5 files changed, 54 insertions(+), 28 deletions(-) create mode 100644 Source/Utils/ParFunctions.H diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 307b03f22..fa5204cef 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -506,10 +506,6 @@ ERF::InitData () if (solverChoice.use_terrain) { if (init_type == "ideal") { amrex::Abort("We do not currently support init_type = ideal with terrain"); - } else if (init_type == "input_sounding") { - amrex::Print() << "Note: use_terrain==true with input_sounding" - << " -- the expected use case is to enable grid stretching" - << std::endl; } } @@ -665,7 +661,9 @@ ERF::InitData () initRayleigh(); if (init_type == "input_sounding") { - // overwrite Ubar, Tbar, and thetabar with input profiles + // Overwrite ubar, vbar, and thetabar with input profiles; wbar is + // assumed to be 0. Note: the tau coefficient set by + // prob->erf_init_rayleigh() is still used bool restarting = (!restart_chkfile.empty()); setRayleighRefFromSounding(restarting); } diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 11527b531..8b67f1d18 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -5,6 +5,7 @@ #include #include #include +#include using namespace amrex; @@ -72,18 +73,29 @@ ERF::setRayleighRefFromSounding (bool restarting) for (int lev = 0; lev <= finest_level; lev++) { const int khi = geom[lev].Domain().bigEnd()[2]; - const auto *const prob_lo = geom[lev].ProbLo(); - const auto *const dx = geom[lev].CellSize(); + Vector zcc(khi+1); + + if (z_phys_cc[lev]) { + // use_terrain=1 + // calculate the damping strength based on the max height at each k + reduce_to_max_per_level(zcc, z_phys_cc[lev]); + } else { + const auto *const prob_lo = geom[lev].ProbLo(); + const auto *const dx = geom[lev].CellSize(); + for (int k = 0; k <= khi; k++) + { + zcc[k] = prob_lo[2] + (k+0.5) * dx[2]; + } + } for (int k = 0; k <= khi; k++) { - const Real z = prob_lo[2] + (k + 0.5) * dx[2]; - h_rayleigh_ptrs[lev][Rayleigh::ubar][k] = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size); - h_rayleigh_ptrs[lev][Rayleigh::vbar][k] = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size); - h_rayleigh_ptrs[lev][Rayleigh::wbar][k] = Real(0.0); - h_rayleigh_ptrs[lev][Rayleigh::thetabar][k] = interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size); + h_rayleigh_ptrs[lev][Rayleigh::ubar][k] = interpolate_1d(z_inp_sound, U_inp_sound, zcc[k], inp_sound_size); + h_rayleigh_ptrs[lev][Rayleigh::vbar][k] = interpolate_1d(z_inp_sound, V_inp_sound, zcc[k], inp_sound_size); + h_rayleigh_ptrs[lev][Rayleigh::wbar][k] = Real(0.0); + h_rayleigh_ptrs[lev][Rayleigh::thetabar][k] = interpolate_1d(z_inp_sound, theta_inp_sound, zcc[k], inp_sound_size); if (h_rayleigh_ptrs[lev][Rayleigh::tau][k] > 0) { - amrex::Print() << z << ":" << " tau=" << h_rayleigh_ptrs[lev][Rayleigh::tau][k]; + amrex::Print() << zcc[k] << ":" << " tau=" << h_rayleigh_ptrs[lev][Rayleigh::tau][k]; if (solverChoice.rayleigh_damp_U) amrex::Print() << " ubar = " << h_rayleigh_ptrs[lev][Rayleigh::ubar][k]; if (solverChoice.rayleigh_damp_V) amrex::Print() << " vbar = " << h_rayleigh_ptrs[lev][Rayleigh::vbar][k]; if (solverChoice.rayleigh_damp_W) amrex::Print() << " wbar = " << h_rayleigh_ptrs[lev][Rayleigh::wbar][k]; diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index 110e65deb..5c0d824f4 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -1,3 +1,5 @@ +#include + /** * Initialize a Rayleigh damping layer with the same structure as in WRF, based * on Durran and Klemp 1983 @@ -13,20 +15,7 @@ erf_init_rayleigh (amrex::Vector >& rayleigh_ptrs, if (z_phys_cc) { // use_terrain=1 // calculate the damping strength based on the max height at each k - amrex::MultiArray4 const& ma_z_phys = z_phys_cc->const_arrays(); - for (int k = 0; k <= khi; k++) { - zcc[k] = amrex::ParReduce(amrex::TypeList{}, - amrex::TypeList{}, - *z_phys_cc, - amrex::IntVect(0), - [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept - -> amrex::GpuTuple - { - return { ma_z_phys[box_no](i,j,k) }; - }); - } - amrex::ParallelDescriptor::ReduceRealMax(zcc.data(), zcc.size()); - + reduce_to_max_per_level(zcc, z_phys_cc); } else { const amrex::Real* prob_lo = geom.ProbLo(); const auto dx = geom.CellSize(); @@ -52,7 +41,6 @@ erf_init_rayleigh (amrex::Vector >& rayleigh_ptrs, rayleigh_ptrs[Rayleigh::vbar][k] = parms.V_0; rayleigh_ptrs[Rayleigh::wbar][k] = parms.W_0; rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; - amrex::Print() << " " << zcc[k] << " " << sinefac*sinefac << std::endl; } else { diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 90df9ccb6..f6a9017f9 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -13,6 +13,8 @@ CEXE_headers += Interpolation.H CEXE_headers += Interpolation_1D.H CEXE_sources += TerrainMetrics.cpp +CEXE_headers += ParFunctions.H + CEXE_headers += Sat_methods.H CEXE_headers += Water_vapor_saturation.H CEXE_headers += DirectionSelector.H diff --git a/Source/Utils/ParFunctions.H b/Source/Utils/ParFunctions.H new file mode 100644 index 000000000..d88d0edf0 --- /dev/null +++ b/Source/Utils/ParFunctions.H @@ -0,0 +1,26 @@ +#ifndef ParFunctions_H +#define ParFunctions_H + +/** + * Reduce a multifab to a vector of values at height levels + */ +AMREX_FORCE_INLINE +void reduce_to_max_per_level (amrex::Vector& v, + std::unique_ptr& mf) +{ + amrex::MultiArray4 const& ma = mf->const_arrays(); + for (int k = 0; k < v.size(); k++) { + v[k] = amrex::ParReduce(amrex::TypeList{}, + amrex::TypeList{}, + *mf, + amrex::IntVect(0), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept + -> amrex::GpuTuple + { + return { ma[box_no](i,j,k) }; + }); + } + amrex::ParallelDescriptor::ReduceRealMax(v.data(), v.size()); +} + +#endif /* ParFunctions_H */