Correcting Coriolis and geostrophic wind forcing to allow for any lat…

…itude (#897)

* changing Coriolis and Geostrophic forces to admit a two-component formulation

* making changes to make sure that it is 2\pi instead of 2* 2\pi

* changed naming convention from is_2d to is_horizontal + added fac to geostrophic wind to allow for a vertical component

* small fix

amr-wind/equation_systems/icns/source_terms/CoriolisForcing.H

@@ -44,8 +44,13 @@ private:
     amrex::Real m_sinphi{0.0};
     amrex::Real m_cosphi{0.0};
 
-    //! `2.0 * \Omega`
+    //! \Omega --> Earth’s rotation angular speed
+    amrex::Real m_omega{7.3e-5};
+
+    //! `2.0 * \Omega
     amrex::Real m_coriolis_factor{0.0};
+
+    bool m_is_horizontal{false};
 };
 
 } // namespace amr_wind::pde::icns

amr-wind/equation_systems/icns/source_terms/CoriolisForcing.cpp

@@ -36,16 +36,21 @@ CoriolisForcing::CoriolisForcing(const CFDSim& sim)
     m_sinphi = std::sin(m_latitude);
     m_cosphi = std::cos(m_latitude);
 
-    // Read the rotational time period (in seconds)
-    amrex::Real rot_time_period = 86400.0;
+    // Read the rotational time period (in seconds) -- This is 23hrs and 56
+    // minutes and 4.091 seconds
+    amrex::Real rot_time_period = 86164.091;
     pp.query("rotational_time_period", rot_time_period);
-    m_coriolis_factor = 2.0 * utils::two_pi() / rot_time_period;
+
+    m_omega = utils::two_pi() / rot_time_period;
+    m_coriolis_factor = 2.0 * m_omega;
 
     pp.queryarr("east_vector", m_east, 0, AMREX_SPACEDIM);
     pp.queryarr("north_vector", m_north, 0, AMREX_SPACEDIM);
     utils::vec_normalize(m_east.data());
     utils::vec_normalize(m_north.data());
     utils::cross_prod(m_east.data(), m_north.data(), m_up.data());
+
+    pp.query("is_horizontal", m_is_horizontal);
 }
 
 CoriolisForcing::~CoriolisForcing() = default;
@@ -70,6 +75,8 @@ void CoriolisForcing::operator()(
     const auto& vel =
         m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
 
+    amrex::Real fac = (m_is_horizontal) ? 0. : 1.;
+
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const amrex::Real ue = east[0] * vel(i, j, k, 0) +
                                east[1] * vel(i, j, k, 1) +
@@ -81,9 +88,9 @@ void CoriolisForcing::operator()(
                                up[1] * vel(i, j, k, 1) +
                                up[2] * vel(i, j, k, 2);
 
-        const amrex::Real ae = +corfac * (un * sinphi - uu * cosphi);
+        const amrex::Real ae = +corfac * (un * sinphi - fac * uu * cosphi);
         const amrex::Real an = -corfac * ue * sinphi;
-        const amrex::Real au = +corfac * ue * cosphi;
+        const amrex::Real au = +fac * corfac * ue * cosphi;
 
         const amrex::Real ax = ae * east[0] + an * north[0] + au * up[0];
         const amrex::Real ay = ae * east[1] + an * north[1] + au * up[1];

amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H

@@ -31,6 +31,8 @@ private:
 
     //! Forcing source term (pressure gradient)
     amrex::Vector<amrex::Real> m_g_forcing{{0.0, 0.0, 0.0}};
+
+    bool m_is_horizontal{false};
 };
 
 } // namespace amr_wind::pde::icns

amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.cpp

@@ -22,28 +22,29 @@ namespace amr_wind::pde::icns {
  */
 GeostrophicForcing::GeostrophicForcing(const CFDSim& /*unused*/)
 {
-    amrex::Real coriolis_factor;
-    {
-        // Read the rotational time period (in seconds)
-        amrex::ParmParse pp("CoriolisForcing");
-        amrex::Real rot_time_period = 86400.0;
-        pp.query("rotational_time_period", rot_time_period);
-        coriolis_factor = 2.0 * utils::two_pi() / rot_time_period;
-        amrex::Print() << "Geostrophic forcing: Coriolis factor = "
-                       << coriolis_factor << std::endl;
+    amrex::Real coriolis_factor = 0.0;
 
-        amrex::Real latitude = 90.0;
-        pp.query("latitude", latitude);
-        AMREX_ALWAYS_ASSERT(
-            std::abs(latitude - 90.0) <
-            static_cast<amrex::Real>(vs::DTraits<float>::eps()));
-    }
+    // Read the rotational time period (in seconds)
+    amrex::ParmParse ppc("CoriolisForcing");
+    // Read the rotational time period (in seconds) -- This is 23hrs and 56
+    // minutes and 4.091 seconds
+    amrex::Real rot_time_period = 86164.091;
+    ppc.query("rotational_time_period", rot_time_period);
 
-    {
-        // Read the geostrophic wind speed vector (in m/s)
-        amrex::ParmParse pp("GeostrophicForcing");
-        pp.getarr("geostrophic_wind", m_target_vel);
-    }
+    amrex::Real omega = utils::two_pi() / rot_time_period;
+    amrex::Real latitude = 90;
+    ppc.get("latitude", latitude);
+    latitude = utils::radians(latitude);
+    amrex::Real sinphi = std::sin(latitude);
+
+    coriolis_factor = 2.0 * omega * sinphi;
+    ppc.query("is_horizontal", m_is_horizontal);
+    amrex::Print() << "Geostrophic forcing: Coriolis factor = "
+                   << coriolis_factor << std::endl;
+
+    // Read the geostrophic wind speed vector (in m/s)
+    amrex::ParmParse ppg("GeostrophicForcing");
+    ppg.getarr("geostrophic_wind", m_target_vel);
 
     m_g_forcing = {
         -coriolis_factor * m_target_vel[1], coriolis_factor * m_target_vel[0],
@@ -59,12 +60,13 @@ void GeostrophicForcing::operator()(
     const FieldState /*fstate*/,
     const amrex::Array4<amrex::Real>& src_term) const
 {
+    amrex::Real fac = (m_is_horizontal) ? 0. : 1.;
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> forcing{
         {m_g_forcing[0], m_g_forcing[1], m_g_forcing[2]}};
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         src_term(i, j, k, 0) += forcing[0];
         src_term(i, j, k, 1) += forcing[1];
-        // No forcing in z-direction
+        src_term(i, j, k, 1) += fac * forcing[2];
     });
 }
 

unit_tests/wind_energy/test_abl_src.cpp

@@ -178,7 +178,7 @@ TEST_F(ABLMeshTest, geostrophic_forcing)
     populate_parameters();
 
     amrex::ParmParse pp("CoriolisForcing");
-    pp.add("latitude", 90.0);
+    pp.add("latitude", 54.0);
 
     initialize_mesh();
 
@@ -199,7 +199,8 @@ TEST_F(ABLMeshTest, geostrophic_forcing)
         geostrophic_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr);
     });
 
-    constexpr amrex::Real corfac = 2.0 * amr_wind::utils::two_pi() / 86400.0;
+    constexpr amrex::Real corfac =
+        2.0 * amr_wind::utils::two_pi() / 86164.091 * 0.80901699437;
     const amrex::Array<amrex::Real, AMREX_SPACEDIM> golds{
         {-corfac * 6.0, corfac * 10.0, 0.0}};
     for (int i = 0; i < AMREX_SPACEDIM; ++i) {
@@ -347,7 +348,7 @@ TEST_F(ABLMeshTest, hurricane_forcing)
 TEST_F(ABLMeshTest, coriolis_const_vel)
 {
     constexpr amrex::Real tol = 1.0e-12;
-    constexpr amrex::Real corfac = 2.0 * amr_wind::utils::two_pi() / 86400.0;
+    constexpr amrex::Real corfac = 2.0 * amr_wind::utils::two_pi() / 86164.091;
     // Latitude is set to 45 degrees in the input file so sinphi = cosphi
     const amrex::Real latfac = std::sin(amr_wind::utils::radians(45.0));
     // Initialize a random value for the velocity component
@@ -429,7 +430,7 @@ TEST_F(ABLMeshTest, coriolis_height_variation)
     // ABL unit test mesh has 64 cells in z
     constexpr int kdim = 63;
     constexpr amrex::Real tol = 1.0e-12;
-    constexpr amrex::Real corfac = 2.0 * amr_wind::utils::two_pi() / 86400.0;
+    constexpr amrex::Real corfac = 2.0 * amr_wind::utils::two_pi() / 86164.091;
     // Latitude is set to 45 degrees in the input file so sinphi = cosphi
     const amrex::Real latfac = std::sin(amr_wind::utils::radians(45.0));