Time varying abl bodyforce (#1025)

* changes to write abl forcing to text file, add to existing reg test

* reading in body force from file

* setting up reg test

* formatting

* using start time and frequency arguments

* write to file with Geostrophic Forcing

* formatting

* reverting last change, because Geostrophic Forcing is reversible

* time-varying target velocity for geostrophic wind

* reg test for time-varying geostrophic wind

* updating to use underscores, but keeping backwards-compatible

* removing unused variable

* proliferate m_coriolis_factor (fix issue from last commit)

* input file documentation

* unit test start, better error feedback

* test to write forces file

* body forcing unit test + better precision

* geostrophic forcing unit test

* Correct source term implementations for temporal discretization

* make body force test self-contained

* add headers to time table files in reg tests, note in documentation

* unused variables

* avoiding angle interpolation issues

* Unit tests for new angle interpolation feature

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

@@ -41,20 +41,38 @@ public:
         m_mean_vel[1] = uy;
 
         const auto& current_time = m_time.current_time();
+        const auto& new_time = m_time.new_time();
+        const auto& nph_time = 0.5 * (current_time + new_time);
         const auto& dt = m_time.deltaT();
+        const auto& t_step = m_time.time_index();
 
         if (!m_vel_timetable.empty()) {
-            const amrex::Real current_spd = ::amr_wind::interp::linear(
-                m_time_table, m_speed_table, current_time);
-            const amrex::Real current_dir = ::amr_wind::interp::linear(
-                m_time_table, m_direction_table, current_time);
-
-            m_target_vel[0] = current_spd * std::cos(current_dir);
-            m_target_vel[1] = current_spd * std::sin(current_dir);
+            // Forces should be applied at n+1/2. Because ABL forcing is a
+            // delta, the difference between the target velocity (at n+1) and
+            // the current velocity (at n) puts the force term at n+1/2
+            const amrex::Real new_spd = ::amr_wind::interp::linear(
+                m_time_table, m_speed_table, new_time);
+            const amrex::Real new_dir = ::amr_wind::interp::linear_angle(
+                m_time_table, m_direction_table, new_time, 2.0 * M_PI);
+
+            m_target_vel[0] = new_spd * std::cos(new_dir);
+            m_target_vel[1] = new_spd * std::sin(new_dir);
         }
 
         m_abl_forcing[0] = (m_target_vel[0] - m_mean_vel[0]) / dt;
         m_abl_forcing[1] = (m_target_vel[1] - m_mean_vel[1]) / dt;
+
+        if (m_write_force_timetable &&
+            amrex::ParallelDescriptor::IOProcessor() &&
+            (t_step % m_force_outfreq == 0) &&
+            (current_time >= m_force_outstart)) {
+            std::ofstream outfile;
+            // Forces are recorded at n+1/2
+            outfile.open(m_force_timetable, std::ios::out | std::ios::app);
+            outfile << std::fixed << std::setprecision(15) << nph_time << "\t"
+                    << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t"
+                    << 0.0 << std::endl;
+        }
     }
 
     amrex::RealArray abl_forcing() const { return m_abl_forcing; }
@@ -77,6 +95,15 @@ private:
     //! File name for velocity forcing time table
     std::string m_vel_timetable;
 
+    //! Bool for writing forcing time table
+    bool m_write_force_timetable{false};
+    //! File name for forcing time table output
+    std::string m_force_timetable;
+    //! Output frequency for forces
+    int m_force_outfreq{1};
+    //! Output start time for force
+    amrex::Real m_force_outstart{0.0};
+
     //! Velocity forcing time table
     amrex::Vector<amrex::Real> m_time_table;
 

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

@@ -26,7 +26,9 @@ ABLForcing::ABLForcing(const CFDSim& sim)
     if (!m_vel_timetable.empty()) {
         std::ifstream ifh(m_vel_timetable, std::ios::in);
         if (!ifh.good()) {
-            amrex::Abort("Cannot find input file: " + m_vel_timetable);
+            amrex::Abort(
+                "Cannot find ABLForcing velocity_timetable file: " +
+                m_vel_timetable);
         }
         amrex::Real data_time;
         amrex::Real data_speed;
@@ -43,6 +45,18 @@ ABLForcing::ABLForcing(const CFDSim& sim)
         pp_incflo.getarr("velocity", m_target_vel);
     }
 
+    m_write_force_timetable = pp_abl.contains("forcing_timetable_output_file");
+    if (m_write_force_timetable) {
+        pp_abl.get("forcing_timetable_output_file", m_force_timetable);
+        pp_abl.query("forcing_timetable_frequency", m_force_outfreq);
+        pp_abl.query("forcing_timetable_start_time", m_force_outstart);
+        if (amrex::ParallelDescriptor::IOProcessor()) {
+            std::ofstream outfile;
+            outfile.open(m_force_timetable, std::ios::out);
+            outfile << "time\tfx\tfy\tfz\n";
+        }
+    }
+
     for (int i = 0; i < AMREX_SPACEDIM; ++i) {
         m_mean_vel[i] = m_target_vel[i];
     }

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

@@ -33,6 +33,7 @@ public:
         const amrex::Array4<amrex::Real>& src_term) const override;
 
     void read_bforce_profile(const std::string& filename);
+    void read_bforce_timetable(const std::string& filename);
 
 private:
     //! Time
@@ -44,7 +45,9 @@ private:
     amrex::Vector<amrex::Real> m_body_force{{0.0, 0.0, 0.0}};
 
     //! Body Force Type
-    std::string m_type{"constant"};
+    std::string m_type{"uniform_constant"};
+    //! Uniform time table file
+    std::string m_utt_file;
 
     //! Angular frequency used in the oscillatory forcing
     amrex::Real m_omega{0.0};
@@ -54,6 +57,12 @@ private:
     amrex::Gpu::DeviceVector<amrex::Real> m_prof_x;
     amrex::Gpu::DeviceVector<amrex::Real> m_prof_y;
     amrex::Gpu::DeviceVector<amrex::Real> m_ht;
+
+    //! Vectors for storing uniform_timetable inputs
+    amrex::Vector<amrex::Real> m_time_table;
+    amrex::Vector<amrex::Real> m_fx_table;
+    amrex::Vector<amrex::Real> m_fy_table;
+    amrex::Vector<amrex::Real> m_fz_table;
 };
 
 } // namespace amr_wind::pde::icns

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

@@ -1,6 +1,7 @@
 #include "amr-wind/equation_systems/icns/source_terms/BodyForce.H"
 #include "amr-wind/CFDSim.H"
 #include "amr-wind/utilities/trig_ops.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 #include "AMReX_ParmParse.H"
 #include "AMReX_Gpu.H"
@@ -16,18 +17,38 @@ namespace amr_wind::pde::icns {
 BodyForce::BodyForce(const CFDSim& sim) : m_time(sim.time()), m_mesh(sim.mesh())
 {
 
-    // Read the geostrophic wind speed vector (in m/s)
+    // Body Force arguments
     amrex::ParmParse pp("BodyForce");
     pp.query("type", m_type);
     m_type = amrex::toLower(m_type);
-
-    if (m_type == "height-varying") {
-        pp.get("bodyforce-file", m_bforce_file);
+    bool no_type_specified = !pp.contains("type");
+    bool file_specified = pp.contains("uniform_timetable_file");
+
+    // Prepare type of body force distribution
+    if (m_type == "height_varying" || m_type == "height-varying") {
+        // Constant in time, varies with z
+        // Using underscores is preferred, remains backwards compatible
+        pp.query("bodyforce-file", m_bforce_file);
+        if (m_bforce_file.empty()) {
+            pp.get("bodyforce_file", m_bforce_file);
+        }
         read_bforce_profile(m_bforce_file);
+    } else if (
+        m_type == "uniform_timetable" ||
+        (no_type_specified && file_specified)) {
+        // Still used if type not specified but file is
+        // Uniform in space, varies with time
+        pp.get("uniform_timetable_file", m_utt_file);
+        read_bforce_timetable(m_utt_file);
     } else {
         pp.getarr("magnitude", m_body_force);
         if (m_type == "oscillatory") {
             pp.get("angular_frequency", m_omega);
+        } else if (m_type != "uniform_constant") {
+            amrex::Abort(
+                "BodyForce type not supported. Please choose uniform_constant "
+                "(default), height_varying, oscillatory, or "
+                "uniform_timetable.\n");
         }
     }
 }
@@ -68,14 +89,36 @@ void BodyForce::read_bforce_profile(const std::string& filename)
         m_ht.begin());
 }
 
+void BodyForce::read_bforce_timetable(const std::string& filename)
+{
+    std::ifstream ifh(filename, std::ios::in);
+    if (!ifh.good()) {
+        amrex::Abort(
+            "Cannot find BodyForce uniform_timetable_file: " + filename + "\n");
+    }
+    amrex::Real data_time;
+    amrex::Real data_fx;
+    amrex::Real data_fy;
+    amrex::Real data_fz;
+    // Skip first line (header)
+    ifh.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+    while (ifh >> data_time) {
+        ifh >> data_fx >> data_fy >> data_fz;
+        m_time_table.push_back(data_time);
+        m_fx_table.push_back(data_fx);
+        m_fy_table.push_back(data_fy);
+        m_fz_table.push_back(data_fz);
+    }
+}
+
 void BodyForce::operator()(
     const int lev,
     const amrex::MFIter& /*mfi*/,
     const amrex::Box& bx,
     const FieldState /*fstate*/,
     const amrex::Array4<amrex::Real>& src_term) const
 {
-    const auto& time = m_time.current_time();
+    const auto& nph_time = 0.5 * (m_time.current_time() + m_time.new_time());
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
     const int lp1 = lev + 1;
@@ -85,7 +128,7 @@ void BodyForce::operator()(
     const amrex::Real* force_x = m_prof_x.data();
     const amrex::Real* force_y = m_prof_y.data();
 
-    if (m_type == "height-varying") {
+    if (m_type == "height_varying" || m_type == "height-varying") {
 
         amrex::ParallelFor(
             bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -113,8 +156,18 @@ void BodyForce::operator()(
         amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> forcing{
             {m_body_force[0], m_body_force[1], m_body_force[2]}};
 
+        if (!m_utt_file.empty()) {
+            // Populate forcing from file if supplied
+            forcing[0] =
+                amr_wind::interp::linear(m_time_table, m_fx_table, nph_time);
+            forcing[1] =
+                amr_wind::interp::linear(m_time_table, m_fy_table, nph_time);
+            forcing[2] =
+                amr_wind::interp::linear(m_time_table, m_fz_table, nph_time);
+        }
+
         amrex::Real coeff =
-            (m_type == "oscillatory") ? std::cos(m_omega * time) : 1.0;
+            (m_type == "oscillatory") ? std::cos(m_omega * nph_time) : 1.0;
         amrex::ParallelFor(
             bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                 src_term(i, j, k, 0) += coeff * forcing[0];

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

@@ -26,8 +26,24 @@ public:
         const amrex::Array4<amrex::Real>& src_term) const override;
 
 private:
+    const SimTime& m_time;
     const amrex::AmrCore& m_mesh;
 
+    //! File name for target velocity time table
+    std::string m_vel_timetable;
+
+    //! Velocity forcing time table
+    amrex::Vector<amrex::Real> m_time_table;
+
+    //! Velocity forcing speed table
+    amrex::Vector<amrex::Real> m_speed_table;
+
+    //! Velocity forcing direction table
+    amrex::Vector<amrex::Real> m_direction_table;
+
+    //! Coriolis factor
+    amrex::Real m_coriolis_factor;
+
     //! Activated when water is present in domain
     bool m_use_phase_ramp{false};
 

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

@@ -4,6 +4,7 @@
 #include "amr-wind/core/vs/vstraits.H"
 #include "amr-wind/physics/multiphase/MultiPhase.H"
 #include "amr-wind/equation_systems/vof/volume_fractions.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 #include "AMReX_ParmParse.H"
 #include "AMReX_Gpu.H"
@@ -22,10 +23,9 @@ namespace amr_wind::pde::icns {
  *    GeostrophicForcing namespace
  *
  */
-GeostrophicForcing::GeostrophicForcing(const CFDSim& sim) : m_mesh(sim.mesh())
+GeostrophicForcing::GeostrophicForcing(const CFDSim& sim)
+    : m_time(sim.time()), m_mesh(sim.mesh())
 {
-    amrex::Real coriolis_factor = 0.0;
-
     // Read the rotational time period (in seconds)
     amrex::ParmParse ppc("CoriolisForcing");
     // Read the rotational time period (in seconds) -- This is 23hrs and 56
@@ -39,18 +39,40 @@ GeostrophicForcing::GeostrophicForcing(const CFDSim& sim) : m_mesh(sim.mesh())
     latitude = utils::radians(latitude);
     amrex::Real sinphi = std::sin(latitude);
 
-    coriolis_factor = 2.0 * omega * sinphi;
+    m_coriolis_factor = 2.0 * omega * sinphi;
     ppc.query("is_horizontal", m_is_horizontal);
     amrex::Print() << "Geostrophic forcing: Coriolis factor = "
-                   << coriolis_factor << std::endl;
+                   << m_coriolis_factor << std::endl;
 
     // Read the geostrophic wind speed vector (in m/s)
     amrex::ParmParse ppg("GeostrophicForcing");
-    ppg.getarr("geostrophic_wind", m_target_vel);
+    ppg.query("geostrophic_wind_timetable", m_vel_timetable);
+    if (!m_vel_timetable.empty()) {
+        std::ifstream ifh(m_vel_timetable, std::ios::in);
+        if (!ifh.good()) {
+            amrex::Abort(
+                "Cannot find GeostrophicForcing geostrophic_wind_timetable "
+                "file: " +
+                m_vel_timetable);
+        }
+        amrex::Real data_time;
+        amrex::Real data_speed;
+        amrex::Real data_deg;
+        ifh.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+        while (ifh >> data_time) {
+            ifh >> data_speed >> data_deg;
+            amrex::Real data_rad = utils::radians(data_deg);
+            m_time_table.push_back(data_time);
+            m_speed_table.push_back(data_speed);
+            m_direction_table.push_back(data_rad);
+        }
+    } else {
+        ppg.getarr("geostrophic_wind", m_target_vel);
+    }
 
     m_g_forcing = {
-        -coriolis_factor * m_target_vel[1], coriolis_factor * m_target_vel[0],
-        0.0};
+        -m_coriolis_factor * m_target_vel[1],
+        m_coriolis_factor * m_target_vel[0], 0.0};
 
     // Set up relaxation toward 0 forcing near the air-water interface
     if (sim.repo().field_exists("vof")) {
@@ -83,6 +105,8 @@ void GeostrophicForcing::operator()(
     const amrex::Array4<amrex::Real>& src_term) const
 {
     amrex::Real hfac = (m_is_horizontal) ? 0. : 1.;
+    // Forces applied at n+1/2
+    const auto& nph_time = 0.5 * (m_time.current_time() + m_time.new_time());
 
     const bool ph_ramp = m_use_phase_ramp;
     const int n_band = m_n_band;
@@ -94,6 +118,21 @@ void GeostrophicForcing::operator()(
     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> forcing{
         {m_g_forcing[0], m_g_forcing[1], m_g_forcing[2]}};
 
+    // Calculate forcing values if target velocity is a function of time
+    if (!m_vel_timetable.empty()) {
+        const amrex::Real nph_spd =
+            amr_wind::interp::linear(m_time_table, m_speed_table, nph_time);
+        const amrex::Real nph_dir = amr_wind::interp::linear_angle(
+            m_time_table, m_direction_table, nph_time, 2.0 * M_PI);
+
+        const amrex::Real target_u = nph_spd * std::cos(nph_dir);
+        const amrex::Real target_v = nph_spd * std::sin(nph_dir);
+
+        forcing[0] = -m_coriolis_factor * target_v;
+        forcing[1] = m_coriolis_factor * target_u;
+        forcing[2] = 0.0;
+    }
+
     const auto& vof = (*m_vof)(lev).const_array(mfi);
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::Real wfac = 1.0;

amr-wind/equation_systems/temperature/source_terms/BodyForce.H

@@ -33,7 +33,7 @@ private:
     const SimTime& m_time;
     const amrex::AmrCore& m_mesh;
 
-    std::string m_type{"height-varying"};
+    std::string m_type{"height_varying"};
     std::string m_bforce_file;
     size_t m_bforce_profile_nhts;
 

amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp

@@ -23,8 +23,12 @@ BodyForce::BodyForce(const CFDSim& sim) : m_time(sim.time()), m_mesh(sim.mesh())
     pp.query("type", m_type);
     m_type = amrex::toLower(m_type);
 
-    if (m_type == "height-varying") {
-        pp.get("bodyforce-file", m_bforce_file);
+    // Updated to underscores (more common) but still backwards-compatible
+    if (m_type == "height_varying" || m_type == "height-varying") {
+        pp.query("bodyforce-file", m_bforce_file);
+        if (m_bforce_file.empty()) {
+            pp.get("bodyforce_file", m_bforce_file);
+        }
         read_bforce_profile(m_bforce_file);
     }
 }
@@ -74,7 +78,7 @@ void BodyForce::operator()(
     const amrex::Real* force_ht = m_ht.data();
     const amrex::Real* force_theta = m_prof_theta.data();
 
-    if (m_type == "height-varying") {
+    if (m_type == "height_varying" || m_type == "height-varying") {
 
         amrex::ParallelFor(
             bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {