Initialize anelastic with isentropic hydrostatic equilibrium (Exawind#1484)

marchdf · web-flow · commit d2b200d068c2 · 2025-02-13T16:13:28.000-07:00
diff --git a/amr-wind/CMakeLists.txt b/amr-wind/CMakeLists.txt
@@ -33,6 +33,7 @@ add_subdirectory(overset)
 add_subdirectory(immersed_boundary)
 add_subdirectory(mesh_mapping_models)
 add_subdirectory(ocean_waves)
+add_subdirectory(eos_models)
 
 include(AMReXBuildInfo)
 generate_buildinfo(${amr_wind_lib_name} ${CMAKE_SOURCE_DIR})
diff --git a/amr-wind/eos_models/CMakeLists.txt b/amr-wind/eos_models/CMakeLists.txt
@@ -0,0 +1 @@
+
diff --git a/amr-wind/eos_models/EOSModel.H b/amr-wind/eos_models/EOSModel.H
@@ -0,0 +1,87 @@
+#ifndef EOSMODEL_H
+#define EOSMODEL_H
+
+#include "AMReX_REAL.H"
+#include "amr-wind/utilities/constants.H"
+
+namespace amr_wind::eos {
+
+/**
+ *  \defgroup eos Equation of state models
+ *
+ *  AMR-Wind representation of equation of state models. These
+ *  expression are adapted from those found in ERF
+ *  https://github.com/erf-model/ERF.
+ *
+ */
+
+/** Gamma law equation of state
+ *  \ingroup eos
+ */
+struct GammaLaw
+{
+    using eos_type = GammaLaw;
+    static std::string identifier() { return "GammaLaw"; }
+
+    /**
+     * Return pressure given density and potential temperature
+     *
+     * \param [in] rho    density
+     * \param [in] theta  potential temperature
+     * \param [in] qv     water vapor
+     * \return            pressure
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real p_rth(
+        const amrex::Real rho,
+        const amrex::Real theta,
+        const amrex::Real qv = 0.0) const
+    {
+        const auto rhotheta = rho * theta;
+        return m_p0 * std::pow(
+                          m_air_gas_constant * rhotheta *
+                              (1.0 + (m_water_vapor_gas_constant /
+                                      m_air_gas_constant) *
+                                         qv) *
+                              m_ip0,
+                          m_gamma);
+    }
+
+    /**
+     * Return dP/drho at constant theta
+     *
+     * \param [in] rho   density
+     * \param [in] theta potential temperature
+     * \param [in] qv    water vapor
+     * \return           pressure
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real dp_constanttheta(
+        const amrex::Real rho,
+        const amrex::Real theta,
+        const amrex::Real qv = 0.0) const
+    {
+        return m_gamma * m_p0 *
+               std::pow(
+                   (m_air_gas_constant * theta *
+                    (1.0 +
+                     m_water_vapor_gas_constant / m_air_gas_constant * qv) *
+                    m_ip0),
+                   m_gamma) *
+               std::pow(rho, m_gamma - 1.0);
+    }
+
+    template <class... Args>
+    AMREX_GPU_HOST_DEVICE explicit GammaLaw(const amrex::Real p0 = 1e-5)
+        : m_p0{p0}, m_ip0{1.0 / p0}
+    {}
+
+    const amrex::Real m_gamma{constants::HEAT_CAPACITY_RATIO};
+    const amrex::Real m_p0{1e5};
+    const amrex::Real m_ip0{1e-5};
+    const amrex::Real m_air_gas_constant{
+        constants::UNIVERSAL_GAS_CONSTANT / constants::MOLAR_MASS_AIR};
+    const amrex::Real m_water_vapor_gas_constant{
+        constants::UNIVERSAL_GAS_CONSTANT / constants::MOLAR_MASS_WATER_VAPOR};
+};
+} // namespace amr_wind::eos
+
+#endif /* EOSMODEL_H */
diff --git a/amr-wind/utilities/constants.H b/amr-wind/utilities/constants.H
@@ -46,5 +46,14 @@ static constexpr amrex::Real AVOGADRO_CONSTANT = 6.02214076 * 1e23;
 static constexpr amrex::Real UNIVERSAL_GAS_CONSTANT =
     AVOGADRO_CONSTANT * BOLTZMANN_CONSTANT;
 
+//! Heat capacity ratio
+static constexpr amrex::Real HEAT_CAPACITY_RATIO = 1.4;
+
+//! Molar mass of air (kg/mol)
+static constexpr amrex::Real MOLAR_MASS_AIR = 0.02896492;
+
+//! Molar mass of water vapor (kg/mol)
+static constexpr amrex::Real MOLAR_MASS_WATER_VAPOR = 0.01801528;
+
 } // namespace amr_wind::constants
 #endif
diff --git a/amr-wind/wind_energy/ABLAnelastic.H b/amr-wind/wind_energy/ABLAnelastic.H
@@ -1,8 +1,10 @@
 #ifndef ABLANELASTIC_H
 #define ABLANELASTIC_H
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/utilities/constants.H"
 #include "amr-wind/utilities/MultiLevelVector.H"
+#include "amr-wind/eos_models/EOSModel.H"
 #include "AMReX_ParmParse.H"
 
 namespace amr_wind {
@@ -26,20 +28,24 @@ public:
     bool is_anelastic() const { return m_is_anelastic; }
 
 private:
+    //! Function to calculate the hydrostatic density and pressure
+    void initialize_isentropic_hse();
+
     const CFDSim& m_sim;
 
     bool m_is_anelastic{false};
 
     amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 
-    amrex::Real m_rho0_const{1.0};
+    amrex::Real m_reference_density_constant{1.0};
 
-    amrex::Real m_bottom_reference_pressure{1.01325e5};
+    amrex::Real m_bottom_reference_pressure{1.0e5};
 
     const int m_axis{2};
 
     MultiLevelVector m_density{FieldLoc::CELL};
     MultiLevelVector m_pressure{FieldLoc::NODE};
+    MultiLevelVector m_theta{FieldLoc::CELL};
 };
 } // namespace amr_wind
 #endif
diff --git a/amr-wind/wind_energy/ABLAnelastic.cpp b/amr-wind/wind_energy/ABLAnelastic.cpp
@@ -14,7 +14,7 @@ ABLAnelastic::ABLAnelastic(CFDSim& sim) : m_sim(sim)
     {
         amrex::ParmParse pp("incflo");
         pp.queryarr("gravity", m_gravity);
-        pp.query("density", m_rho0_const);
+        pp.query("density", m_reference_density_constant);
         pp.query("godunov_type", godunov_type);
         pp.query("icns_conserv", conserv);
     }
@@ -35,8 +35,8 @@ ABLAnelastic::ABLAnelastic(CFDSim& sim) : m_sim(sim)
             "reference_density", 1, density.num_grow()[0], 1);
         ref_density.set_default_fillpatch_bc(m_sim.time());
         m_sim.repo().declare_nd_field("reference_pressure", 1, 0, 1);
-        auto& ref_theta = m_sim.repo().declare_field(
-            "reference_temperature", 1, density.num_grow()[0], 1);
+        auto& ref_theta =
+            m_sim.repo().declare_field("reference_temperature", 1, 1, 1);
         ref_theta.set_default_fillpatch_bc(m_sim.time());
     }
 }
@@ -62,45 +62,108 @@ void ABLAnelastic::initialize_data()
 
     m_density.resize(m_axis, m_sim.mesh().Geom());
     m_pressure.resize(m_axis, m_sim.mesh().Geom());
+    m_theta.resize(m_axis, m_sim.mesh().Geom());
 
     AMREX_ALWAYS_ASSERT(m_sim.repo().num_active_levels() == m_density.size());
     AMREX_ALWAYS_ASSERT(m_sim.repo().num_active_levels() == m_pressure.size());
+    AMREX_ALWAYS_ASSERT(m_sim.repo().num_active_levels() == m_theta.size());
+
+    initialize_isentropic_hse();
+
+    auto& density0_field = m_sim.repo().get_field("reference_density");
+    auto& p0 = m_sim.repo().get_field("reference_pressure");
+    auto& temp0 = m_sim.repo().get_field("reference_temperature");
+    m_density.copy_to_field(density0_field);
+    m_pressure.copy_to_field(p0);
+    m_theta.copy_to_field(temp0);
+    density0_field.fillpatch(m_sim.time().current_time());
+    temp0.fillpatch(m_sim.time().current_time());
+}
+
+void ABLAnelastic::initialize_isentropic_hse()
+{
+    const int max_iterations = 10;
+    const auto ref_theta = m_sim.transport_model().reference_temperature();
+    const auto eos = amr_wind::eos::GammaLaw(m_bottom_reference_pressure);
 
     for (int lev = 0; lev < m_density.size(); lev++) {
         auto& dens = m_density.host_data(lev);
         auto& pres = m_pressure.host_data(lev);
-        const auto& dx = m_sim.mesh().Geom(lev).CellSizeArray();
-        dens.assign(dens.size(), m_rho0_const);
-        pres[0] = m_bottom_reference_pressure;
-        for (int k = 0; k < pres.size() - 1; k++) {
-            pres[k + 1] = pres[k] - dens[k] * m_gravity[m_axis] * dx[m_axis];
+        auto& theta = m_theta.host_data(lev);
+        theta.assign(theta.size(), ref_theta);
+
+        const auto& dx = m_sim.mesh().Geom(lev).CellSizeArray()[m_axis];
+        const amrex::Real half_dx = 0.5 * dx;
+
+        // Initial guess
+        dens[0] = m_reference_density_constant;
+        pres[0] =
+            m_bottom_reference_pressure + half_dx * dens[0] * m_gravity[m_axis];
+
+        // We do a Newton iteration to satisfy the EOS & HSE (with constant
+        // theta) at the surface
+        bool converged_hse = false;
+        amrex::Real p_hse = 0.0;
+        amrex::Real p_eos = 0.0;
+
+        for (int iter = 0; (iter < max_iterations) && (!converged_hse);
+             iter++) {
+            p_hse = m_bottom_reference_pressure +
+                    half_dx * dens[0] * m_gravity[m_axis];
+            p_eos = eos.p_rth(dens[0], ref_theta);
+
+            const amrex::Real p_diff = p_hse - p_eos;
+            const amrex::Real dpdr = eos.dp_constanttheta(dens[0], ref_theta);
+            const amrex::Real ddens =
+                p_diff / (dpdr - half_dx * m_gravity[m_axis]);
+
+            dens[0] = dens[0] + ddens;
+            pres[0] = eos.p_rth(dens[0], ref_theta);
+
+            if (std::abs(ddens) < constants::LOOSE_TOL) {
+                converged_hse = true;
+                break;
+            }
+        }
+        if (!converged_hse) {
+            amrex::Abort("Initializing with hydrostatic equilibrium failed");
         }
-    }
-    m_density.copy_host_to_device();
-    m_pressure.copy_host_to_device();
 
-    auto& rho0 = m_sim.repo().get_field("reference_density");
-    auto& p0 = m_sim.repo().get_field("reference_pressure");
-    m_density.copy_to_field(rho0);
-    m_pressure.copy_to_field(p0);
+        // To get values at k > 0 we do a Newton iteration to satisfy the EOS
+        // (with constant theta)
+        for (int k = 1; k < dens.size(); k++) {
+            converged_hse = false;
 
-    rho0.fillpatch(m_sim.time().current_time());
+            dens[k] = dens[k - 1];
+            for (int iter = 0; (iter < max_iterations) && (!converged_hse);
+                 iter++) {
+                const amrex::Real dens_avg = 0.5 * (dens[k - 1] + dens[k]);
+                p_hse = pres[k - 1] + dx * dens_avg * m_gravity[m_axis];
+                p_eos = eos.p_rth(dens[k], ref_theta);
 
-    auto& temp0 = m_sim.repo().get_field("reference_temperature");
-    const amrex::Real air_molar_mass = 0.02896492; // kg/mol
-    const amrex::Real Rair = constants::UNIVERSAL_GAS_CONSTANT / air_molar_mass;
-    for (int lev = 0; lev < m_sim.repo().num_active_levels(); ++lev) {
-        const auto& rho0_arrs = rho0(lev).const_arrays();
-        const auto& p0_arrs = p0(lev).const_arrays();
-        const auto& temp0_arrs = temp0(lev).arrays();
-        amrex::ParallelFor(
-            temp0(lev), amrex::IntVect(0),
-            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
-                temp0_arrs[nbx](i, j, k) =
-                    p0_arrs[nbx](i, j, k) / (Rair * rho0_arrs[nbx](i, j, k));
-            });
+                const amrex::Real p_diff = p_hse - p_eos;
+                const amrex::Real dpdr =
+                    eos.dp_constanttheta(dens[k], ref_theta);
+                const amrex::Real ddens =
+                    p_diff / (dpdr - dx * m_gravity[m_axis]);
+
+                dens[k] = dens[k] + ddens;
+                pres[k] = eos.p_rth(dens[k], ref_theta);
+
+                if (std::abs(ddens) < constants::LOOSE_TOL * dens[k - 1]) {
+                    converged_hse = true;
+                    break;
+                }
+            }
+            if (!converged_hse) {
+                amrex::Abort(
+                    "Initializing with hydrostatic equilibrium failed");
+            }
+        }
     }
-    amrex::Gpu::synchronize();
-    temp0.fillpatch(m_sim.time().current_time());
+    m_density.copy_host_to_device();
+    m_pressure.copy_host_to_device();
+    m_theta.copy_host_to_device();
 }
+
 } // namespace amr_wind
diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt
@@ -16,6 +16,7 @@ add_subdirectory(multiphase)
 add_subdirectory(ocean_waves)
 add_subdirectory(projection)
 add_subdirectory(boundary_conditions)
+add_subdirectory(eos_models)
 
 if(AMR_WIND_ENABLE_MASA)
   add_subdirectory(mms)
diff --git a/unit_tests/eos_models/CMakeLists.txt b/unit_tests/eos_models/CMakeLists.txt
@@ -0,0 +1,5 @@
+target_sources(${amr_wind_unit_test_exe_name}
+  PRIVATE
+
+  test_eos.cpp
+  )
diff --git a/unit_tests/eos_models/test_eos.cpp b/unit_tests/eos_models/test_eos.cpp
@@ -0,0 +1,27 @@
+#include <AMReX_Gpu.H>
+#include "aw_test_utils/AmrexTest.H"
+#include "amr-wind/eos_models/EOSModel.H"
+
+namespace amr_wind_tests {
+
+void test_eos_impl()
+{
+    const auto eos = amr_wind::eos::GammaLaw(1.01325e5);
+    amrex::Gpu::DeviceScalar<amrex::Real> val(0.0);
+    auto* d_val = val.dataPtr();
+
+    amrex::ParallelFor(1, [=] AMREX_GPU_DEVICE(int /*unused*/) {
+        d_val[0] = eos.p_rth(1.225, 300.0, 0.5);
+    });
+    EXPECT_NEAR(
+        val.dataValue(), 244859.65251771925, amr_wind::constants::LOOSE_TOL);
+
+    amrex::ParallelFor(1, [=] AMREX_GPU_DEVICE(int /*unused*/) {
+        d_val[0] = eos.dp_constanttheta(1.225, 300.0, 0.5);
+    });
+    EXPECT_NEAR(
+        val.dataValue(), 279839.60287739342, amr_wind::constants::LOOSE_TOL);
+}
+
+TEST(EOSModelTest, test_eos) { test_eos_impl(); }
+} // namespace amr_wind_tests
