From 669992efe1c8314fdee3f1d092c18790dcf00ed8 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 24 Jan 2024 20:59:48 -0700 Subject: [PATCH 1/5] Pass terrain heights to erf_init_rayleigh --- Exec/DevTests/MovingTerrain/prob.H | 3 ++- Exec/DevTests/MovingTerrain/prob.cpp | 3 ++- Exec/Radiation/prob.H | 3 ++- Exec/Radiation/prob.cpp | 3 ++- Exec/RegTests/ScalarAdvDiff/prob.H | 3 ++- Exec/RegTests/ScalarAdvDiff/prob.cpp | 3 ++- Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence | 3 ++- Exec/SquallLine_2D/prob.H | 3 ++- Exec/SquallLine_2D/prob.cpp | 3 ++- Exec/SuperCell/prob.H | 3 ++- Exec/SuperCell/prob.cpp | 3 ++- Source/Initialization/ERF_init1d.cpp | 3 ++- Source/Prob/init_rayleigh_damping.H | 3 ++- Source/prob_common.H | 14 ++++++++------ 14 files changed, 34 insertions(+), 19 deletions(-) diff --git a/Exec/DevTests/MovingTerrain/prob.H b/Exec/DevTests/MovingTerrain/prob.H index fea7bd7f9..c56b3bb7b 100644 --- a/Exec/DevTests/MovingTerrain/prob.H +++ b/Exec/DevTests/MovingTerrain/prob.H @@ -49,7 +49,8 @@ public: amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override; + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override; protected: std::string name() override { return "MovingTerrain"; } diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index 9a25dc3fc..dc1859077 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -137,7 +137,8 @@ Problem::erf_init_rayleigh( amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { // amrex::Error("Should never get here for MovingTerrain problem"); const int khi = geom.Domain().bigEnd()[2]; diff --git a/Exec/Radiation/prob.H b/Exec/Radiation/prob.H index 04025026f..5c20a28c2 100644 --- a/Exec/Radiation/prob.H +++ b/Exec/Radiation/prob.H @@ -52,7 +52,8 @@ public: amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override; + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override; protected: std::string name() override { return "Supercell"; } diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp index f1fd2c10d..912d82b54 100644 --- a/Exec/Radiation/prob.cpp +++ b/Exec/Radiation/prob.cpp @@ -219,7 +219,8 @@ Problem::erf_init_rayleigh( amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; diff --git a/Exec/RegTests/ScalarAdvDiff/prob.H b/Exec/RegTests/ScalarAdvDiff/prob.H index 00cbd5f83..4208fc8d5 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.H +++ b/Exec/RegTests/ScalarAdvDiff/prob.H @@ -58,7 +58,8 @@ public: amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override; + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override; protected: std::string name() override { return "Scalar Advection/Diffusion"; } diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index f161da47a..221c3b2ce 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -41,7 +41,8 @@ Problem::erf_init_rayleigh( amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) { const int khi = geom.Domain().bigEnd()[2]; diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence b/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence index 18ee86f46..460f96575 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence @@ -16,7 +16,8 @@ erf_init_rayleigh(amrex::Vector& tau, amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) { const int khi = geom.Domain().bigEnd()[2]; diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 4035fcfea..9f2fd7122 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -67,7 +67,8 @@ public: amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override; + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override; amrex::Real compute_theta (amrex::Real z); diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 469a24477..8c4ccdf81 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -407,7 +407,8 @@ Problem::erf_init_rayleigh (amrex::Vector& tau, amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; diff --git a/Exec/SuperCell/prob.H b/Exec/SuperCell/prob.H index 04025026f..5c20a28c2 100644 --- a/Exec/SuperCell/prob.H +++ b/Exec/SuperCell/prob.H @@ -52,7 +52,8 @@ public: amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override; + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override; protected: std::string name() override { return "Supercell"; } diff --git a/Exec/SuperCell/prob.cpp b/Exec/SuperCell/prob.cpp index ccd1d6242..128fa7471 100644 --- a/Exec/SuperCell/prob.cpp +++ b/Exec/SuperCell/prob.cpp @@ -219,7 +219,8 @@ Problem::erf_init_rayleigh( amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) + amrex::Geometry const& geom, + std::unique_ptr& /*z_phys_cc*/) { const int khi = geom.Domain().bigEnd()[2]; diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index f3bdd2cec..b2543e7b7 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -44,7 +44,8 @@ ERF::initRayleigh () d_rayleigh_thetabar[lev].resize(zlen_rayleigh, 0.0_rt); prob->erf_init_rayleigh(h_rayleigh_tau[lev], h_rayleigh_ubar[lev], h_rayleigh_vbar[lev], - h_rayleigh_wbar[lev], h_rayleigh_thetabar[lev], geom[lev]); + h_rayleigh_wbar[lev], h_rayleigh_thetabar[lev], geom[lev], + z_phys_cc[lev]); // Copy from host version to device version amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_rayleigh_tau[lev].begin(), h_rayleigh_tau[lev].end(), diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index b10d4271c..a850e5d84 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -8,7 +8,8 @@ erf_init_rayleigh (amrex::Vector& tau, amrex::Vector& vbar, amrex::Vector& wbar, amrex::Vector& thetabar, - amrex::Geometry const& geom) override + amrex::Geometry const& geom, + std::unique_ptr& z_phys_cc) override { const auto ztop = geom.ProbHi()[2]; amrex::Print() << "Rayleigh damping (coef="<& /*tau*/, @@ -163,7 +164,8 @@ public: amrex::Vector& /*vbar*/, amrex::Vector& /*wbar*/, amrex::Vector& /*thetabar*/, - amrex::Geometry const& /*geom*/) + amrex::Geometry const& /*geom*/, + std::unique_ptr& /*z_phys_cc*/) { amrex::Error("Should never call erf_init_rayleigh for "+name()+" problem"); } From 333d5fdd85c9c4cfcfa898c516d8968cbbad6ce7 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 24 Jan 2024 22:52:17 -0700 Subject: [PATCH 2/5] Add terrain support to default init_rayleigh_damping function --- Source/Prob/init_rayleigh_damping.H | 37 +++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index a850e5d84..ba5298d0d 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -11,19 +11,41 @@ erf_init_rayleigh (amrex::Vector& tau, amrex::Geometry const& geom, std::unique_ptr& z_phys_cc) override { - const auto ztop = geom.ProbHi()[2]; + const int khi = geom.Domain().bigEnd()[2]; + amrex::Vector zcc(khi); + + 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) }; + }); + } + } else { + const amrex::Real* prob_lo = geom.ProbLo(); + const auto dx = geom.CellSize(); + for (int k = 0; k <= khi; k++) + { + zcc[k] = prob_lo[2] + (k+0.5) * dx[2]; + } + } + + const auto ztop = zcc[khi-1]; amrex::Print() << "Rayleigh damping (coef="<= 0) { const amrex::Real sinefac = std::sin(PIoTwo*zfrac); @@ -32,6 +54,7 @@ erf_init_rayleigh (amrex::Vector& tau, vbar[k] = parms.V_0; wbar[k] = parms.W_0; thetabar[k] = parms.T_0; + amrex::Print() << " " << zcc[k] << " " << sinefac*sinefac << std::endl; } else { From 6e38097b31efcbdd3720c42081d26067e5e9d09a Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Thu, 25 Jan 2024 00:05:34 -0700 Subject: [PATCH 3/5] Fix indexing, parallel behavior Note: The Rayleigh damping implementation will differ slightly from before in the case without terrain. Previously, ztop was set to the top of the domain; now, ztop is set to max(z_phys_cc[:,:,khi]) -- where ztop is the value at which the sine-squared damping factor is one. Therefore, the top layer of cells will now have a factor of exactly one instead of nearly one. --- Source/Prob/init_rayleigh_damping.H | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Source/Prob/init_rayleigh_damping.H b/Source/Prob/init_rayleigh_damping.H index ba5298d0d..d1fd58afe 100644 --- a/Source/Prob/init_rayleigh_damping.H +++ b/Source/Prob/init_rayleigh_damping.H @@ -12,12 +12,12 @@ erf_init_rayleigh (amrex::Vector& tau, std::unique_ptr& z_phys_cc) override { const int khi = geom.Domain().bigEnd()[2]; - amrex::Vector zcc(khi); + amrex::Vector zcc(khi+1); 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(); + 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{}, @@ -29,6 +29,8 @@ erf_init_rayleigh (amrex::Vector& tau, return { ma_z_phys[box_no](i,j,k) }; }); } + amrex::ParallelDescriptor::ReduceRealMax(zcc.data(), zcc.size()); + } else { const amrex::Real* prob_lo = geom.ProbLo(); const auto dx = geom.CellSize(); @@ -38,7 +40,7 @@ erf_init_rayleigh (amrex::Vector& tau, } } - const auto ztop = zcc[khi-1]; + const auto ztop = zcc[khi]; amrex::Print() << "Rayleigh damping (coef="< Date: Thu, 25 Jan 2024 00:24:44 -0700 Subject: [PATCH 4/5] Enable Rayleigh damping for WoA --- Exec/RegTests/WitchOfAgnesi/prob.H | 5 +++++ Exec/RegTests/WitchOfAgnesi/prob.cpp | 4 ++++ 2 files changed, 9 insertions(+) diff --git a/Exec/RegTests/WitchOfAgnesi/prob.H b/Exec/RegTests/WitchOfAgnesi/prob.H index 922c5f19f..e7f9de6bd 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.H +++ b/Exec/RegTests/WitchOfAgnesi/prob.H @@ -9,6 +9,10 @@ struct ProbParm : ProbParmDefaults { amrex::Real U_0 = 0.0; + amrex::Real V_0 = 0.0; + amrex::Real W_0 = 0.0; + amrex::Real dampcoef = 0.003; + amrex::Real zdamp = 1.0; }; // namespace ProbParm class Problem : public ProblemBase @@ -17,6 +21,7 @@ public: Problem(); #include "Prob/init_density_hse_dry_terrain.H" +#include "Prob/init_rayleigh_damping.H" void init_custom_pert ( const amrex::Box& bx, diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index cf14e4aae..a86a7218b 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -16,6 +16,10 @@ Problem::Problem() // Parse params amrex::ParmParse pp("prob"); pp.query("U_0", parms.U_0); + pp.query("V_0", parms.V_0); + pp.query("W_0", parms.W_0); + pp.query("dampcoef", parms.dampcoef); + pp.query("zdamp", parms.zdamp); init_base_parms(parms.rho_0, parms.T_0); } From 193dbd9c66967e19387d0cc6c40bacb00c9811f6 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Thu, 25 Jan 2024 00:26:14 -0700 Subject: [PATCH 5/5] Update WoA input to demo Rayleigh damping with terrain --- Exec/RegTests/WitchOfAgnesi/inputs | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Exec/RegTests/WitchOfAgnesi/inputs b/Exec/RegTests/WitchOfAgnesi/inputs index fb45433da..1a4344c21 100644 --- a/Exec/RegTests/WitchOfAgnesi/inputs +++ b/Exec/RegTests/WitchOfAgnesi/inputs @@ -41,13 +41,12 @@ erf.check_int = -57600 # 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_int_1 = 10 # number of timesteps between plotfiles erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y # SOLVER CHOICE erf.use_gravity = true erf.use_coriolis = false -erf.use_rayleigh_damping = false erf.les_type = "None" # TERRRAIN GRID TYPE @@ -62,5 +61,9 @@ erf.alpha_T = 0.0 # [m^2/s] # PROBLEM PARAMETERS (optional) prob.T_0 = 300.0 -prob.U_0 = 0.0 +prob.U_0 = 1.0 prob.rho_0 = 1.16 + +erf.use_rayleigh_damping = true +prob.dampcoef = 2000. # ==> time scale = 0.0005 s (for testing) +prob.zdamp = 1.0