Free-Atmosphere Damping for ABLForcing

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

@@ -5,6 +5,7 @@
 #include "amr-wind/core/SimTime.H"
 #include "amr-wind/utilities/trig_ops.H"
 #include "amr-wind/utilities/linear_interpolation.H"
+#include "amr-wind/utilities/FieldPlaneAveraging.H"
 
 namespace amr_wind::pde::icns {
 
@@ -29,12 +30,18 @@ public:
         const FieldState fstate,
         const amrex::Array4<amrex::Real>& src_term) const override;
 
+    void mean_velocity_init(const VelPlaneAveraging& /*vavg*/);
+
+    void mean_velocity_update(const VelPlaneAveraging& /*vavg*/);
+
     inline void set_target_velocities(amrex::Real ux, amrex::Real uy)
     {
         m_target_vel[0] = ux;
         m_target_vel[1] = uy;
     }
 
+    inline void set_boundary_layer_height(amrex::Real zi) { m_zi = zi; }
+
     inline void set_mean_velocities(amrex::Real ux, amrex::Real uy)
     {
         m_mean_vel[0] = ux;
@@ -71,7 +78,9 @@ public:
             outfile.open(m_force_timetable, std::ios::out | std::ios::app);
             outfile << std::setprecision(17) << nph_time << "\t"
                     << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t"
-                    << 0.0 << std::endl;
+                    << 0.0 << "\t" << m_abl_forcing[1] / m_coriolis_factor
+                    << "\t" << -m_abl_forcing[0] / m_coriolis_factor << "\t"
+                    << m_zi << std::endl;
         }
     }
 
@@ -83,9 +92,36 @@ private:
     const SimTime& m_time;
     const amrex::AmrCore& m_mesh;
 
+    amrex::Gpu::DeviceVector<amrex::Real> m_vel_ht;
+    amrex::Gpu::DeviceVector<amrex::Real> m_vel_vals;
+
+    //! Axis over which averages are computed.
+    int m_axis{2};
+
     //! Activated when water is present in domain
     bool m_use_phase_ramp{false};
 
+    //! Activate free-atmosphere damping
+    bool m_fa_damping{false};
+
+    //! Detect free-atmosphere base height (set to zi from ABLStats).
+    bool m_fa_detect_height{false};
+
+    //! Free-atmosphere base height.
+    amrex::Real m_fa_height{750.0};
+
+    //! Start time for free-atmosphere damping.
+    amrex::Real m_fa_time_start{0.0};
+
+    //! End time for free-atmosphere damping.
+    amrex::Real m_fa_time_end{1.0E9};
+
+    //! Time scale for free-atmosphere damping.
+    amrex::Real m_fa_tau{100.0};
+
+    //! Coriolis factor used for free-atmosphere damping.
+    amrex::Real m_coriolis_factor;
+
     //! Number of cells in band to prevent forcing near liquid
     int m_n_band{2};
 
@@ -97,10 +133,13 @@ private:
 
     //! 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};
 
@@ -130,6 +169,9 @@ private:
     //! Local storage of interface location
     amrex::Real m_water_level;
 
+    //! boundary layer height.
+    amrex::Real m_zi;
+
     //! VOF field, to avoid forcing on liquid above force-off height
     const Field* m_vof;
 };

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

@@ -47,16 +47,43 @@ ABLForcing::ABLForcing(const CFDSim& sim)
 
     m_write_force_timetable = pp_abl.contains("forcing_timetable_output_file");
     if (m_write_force_timetable) {
+        amrex::Print() << "Here!!!!" << std::endl;
         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()) {
+            amrex::Print() << m_force_timetable << std::endl;
             std::ofstream outfile;
             outfile.open(m_force_timetable, std::ios::out);
-            outfile << "time\tfx\tfy\tfz\n";
+            outfile << "time\tfx\tfy\tfz\tUgx\tUgy\tzi\n";
         }
     }
 
+    // If free-atmosphere damping is used, read these inputs.
+    pp_abl.query("free_atmosphere_damping", m_fa_damping);
+    if (m_fa_damping) {
+        pp_abl.query("detect_free_atmosphere_height", m_fa_detect_height);
+        if (!m_fa_detect_height) {
+            pp_abl.get("free_atmosphere_height", m_fa_height);
+        }
+        pp_abl.query("free_atmosphere_damping_start_time", m_fa_time_start);
+        pp_abl.query("free_atmosphere_damping_end_time", m_fa_time_end);
+        pp_abl.query("free_atmosphere_damping_time_scale", m_fa_tau);
+
+        amrex::ParmParse pp_coriolis("CoriolisForcing");
+        amrex::Real rot_time_period = 86164.091;
+        pp_coriolis.query("rotational_time_period", rot_time_period);
+
+        amrex::Real omega = utils::two_pi() / rot_time_period;
+
+        amrex::Real latitude = 90.0;
+        pp_coriolis.query("latitude", latitude);
+        latitude = utils::radians(latitude);
+        amrex::Real sinphi = std::sin(latitude);
+
+        m_coriolis_factor = 2.0 * omega * sinphi;
+    }
+
     for (int i = 0; i < AMREX_SPACEDIM; ++i) {
         m_mean_vel[i] = m_target_vel[i];
     }
@@ -80,6 +107,8 @@ ABLForcing::ABLForcing(const CFDSim& sim)
         // Point to something, will not be used
         m_vof = &sim.repo().get_field("velocity");
     }
+
+    mean_velocity_init(abl.abl_statistics().vel_profile_coarse());
 }
 
 ABLForcing::~ABLForcing() = default;
@@ -95,6 +124,19 @@ void ABLForcing::operator()(
     const amrex::Real dudt = m_abl_forcing[0];
     const amrex::Real dvdt = m_abl_forcing[1];
 
+    const bool fa_damping = m_fa_damping;
+    const bool fa_detect_height = m_fa_detect_height;
+    const amrex::Real fa_height = fa_detect_height ? m_zi : m_fa_height;
+    const amrex::Real fa_time_start = m_fa_time_start;
+    const amrex::Real fa_time_end = m_fa_time_end;
+    const amrex::Real fa_tau = m_fa_tau;
+    const amrex::Real fa_u = dvdt / m_coriolis_factor;
+    const amrex::Real fa_v = -dudt / m_coriolis_factor;
+
+    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 bool ph_ramp = m_use_phase_ramp;
     const int n_band = m_n_band;
     const amrex::Real wlev = m_water_level;
@@ -103,11 +145,26 @@ void ABLForcing::operator()(
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
 
+    const int idir = m_axis;
+    const amrex::Real* heights = m_vel_ht.data();
+    const amrex::Real* heights_end = m_vel_ht.end();
+    const amrex::Real* vals = m_vel_vals.data();
+
     const auto& vof = (*m_vof)(lev).const_array(mfi);
+
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+        amrex::IntVect iv(i, j, k);
+        const amrex::Real z = problo[idir] + (iv[idir] + 0.5) * dx[idir];
+
+        const amrex::Real umean =
+            amr_wind::interp::linear(heights, heights_end, vals, z, 3, 0);
+        const amrex::Real vmean =
+            amr_wind::interp::linear(heights, heights_end, vals, z, 3, 1);
+
+        // This part deals with air and water phase and only masks ABL forcing
+        // on the water phase.
         amrex::Real fac = 1.0;
         if (ph_ramp) {
-            const amrex::Real z = problo[2] + (k + 0.5) * dx[2];
             if (z - wlev < wrht0 + wrht1) {
                 if (z - wlev < wrht0) {
                     // Apply no forcing within first interval
@@ -127,11 +184,56 @@ void ABLForcing::operator()(
                 fac = 0.0;
             }
         }
-        src_term(i, j, k, 0) += fac * dudt;
-        src_term(i, j, k, 1) += fac * dvdt;
 
-        // No forcing in z-direction
+        // This part applies the free atmosphere damping, which nudges
+        // the solution in the free atmosphere toward geostrophic as
+        // computed by the wind speed controller. (i.e., the wind speed
+        // controller computes a forcing term that is really a driving
+        // pressure gradient, so there is a corresponding geostrophic
+        // wind. Nudge the free atmosphere toward that geostrophic wind.)
+        amrex::Real src_fa_damping_x = 0.0;
+        amrex::Real src_fa_damping_y = 0.0;
+        if (fa_damping &&
+            ((nph_time >= fa_time_start) && (nph_time <= fa_time_end)) &&
+            (z >= fa_height)) {
+            src_fa_damping_x = (1.0 / fa_tau) * (fa_u - umean);
+            src_fa_damping_y = (1.0 / fa_tau) * (fa_v - vmean);
+        }
+
+        // Sum up the source term
+        src_term(i, j, k, 0) += (fac * dudt) + src_fa_damping_x;
+        src_term(i, j, k, 1) += (fac * dvdt) + src_fa_damping_y;
+        src_term(i, j, k, 2) += 0.0;
     });
 }
 
+void ABLForcing::mean_velocity_init(const VelPlaneAveraging& vavg)
+{
+    m_axis = vavg.axis();
+
+    // The implementation depends the assumption that the ABL statistics class
+    // computes statistics at the cell-centeres only on level 0. If this
+    // assumption changes in future, the implementation will break... so put in
+    // a check here to catch this.
+    AMREX_ALWAYS_ASSERT(
+        m_mesh.Geom(0).Domain().length(m_axis) ==
+        static_cast<int>(vavg.line_centroids().size()));
+
+    m_vel_ht.resize(vavg.line_centroids().size());
+    m_vel_vals.resize(vavg.line_average().size());
+
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, vavg.line_centroids().begin(),
+        vavg.line_centroids().end(), m_vel_ht.begin());
+
+    mean_velocity_update(vavg);
+}
+
+void ABLForcing::mean_velocity_update(const VelPlaneAveraging& vavg)
+{
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, vavg.line_average().begin(),
+        vavg.line_average().end(), m_vel_vals.begin());
+}
+
 } // namespace amr_wind::pde::icns

amr-wind/wind_energy/ABL.cpp

@@ -173,6 +173,7 @@ void ABL::post_regrid_actions() { m_abl_anelastic->post_regrid_actions(); }
  */
 void ABL::pre_advance_work()
 {
+    const auto& zi = m_stats->zi();
     const auto& vel_pa = m_stats->vel_profile();
     m_abl_wall_func.update_umean(
         m_stats->vel_profile(), m_stats->theta_profile_fine());
@@ -201,6 +202,9 @@ void ABL::pre_advance_work()
 #endif
 
         m_abl_forcing->set_mean_velocities(vx, vy);
+        m_abl_forcing->set_boundary_layer_height(zi);
+
+        m_abl_forcing->mean_velocity_update(m_stats->vel_profile_coarse());
     }
 
     if (m_abl_mean_bous != nullptr) {

amr-wind/wind_energy/ABLStats.H

@@ -58,6 +58,8 @@ public:
     //! Compute height of capping inversion
     void compute_zi();
 
+    amrex::Real zi() const override { return m_zi; };
+
     //! Return vel plane averaging instance
     const VelPlaneAveragingFine& vel_profile() const override
     {

amr-wind/wind_energy/ABLStats.cpp

@@ -230,13 +230,29 @@ void ABLStats::post_advance_work()
 
 void ABLStats::compute_zi()
 {
-
+    // This finds the location of the capping inversion (z_i) using the
+    // method of Sullivan et al. in "Structure of the entrainment zone
+    // capping the convective atmospheric boundary layer," JAS, Vol. 55,
+    // 1998.
+    //
+    // In this method, for every x,y location in the computational domain,
+    // the maximum d(\theta)/dz is found.  You then have an x,y array
+    // of x,y local z_i.  z_i is then averaged over the x,y array to give
+    // <z_i>.
+
+    // Compute the potential temperature gradient.
     auto gradT = (this->m_sim.repo())
                      .create_scratch_field(3, m_temperature.num_grow()[0]);
     fvm::gradient(*gradT, m_temperature);
 
     // Only compute zi using coarsest level
     const int lev = 0;
+
+    // Using the AMReX ReduceToPlane function, find the maximum
+    // d(\theta)/dz over the mesh along lines in the normal direction,
+    // which typically would be z (if z is up).  We end up with this
+    // 2D "planar" array of max d(\theta)/dz over the horizontal
+    // extent of the domain.
     const int dir = m_normal_dir;
     const auto& geom = (this->m_sim.repo()).mesh().Geom(lev);
     auto const& domain_box = geom.Domain();
@@ -260,11 +276,21 @@ void ABLStats::compute_zi()
     auto& pinned_tg_fab = device_tg_fab;
 #endif
 
+    // Because of the parallel decomposition, there may be multiple
+    // of these 2D "planar" arrays local to different processes handling
+    // boxes of mesh at different heights.  So, we need to parallel
+    // reduce these by picking the max value at each x,y location to get
+    // the final resultant planar array of max temperature gradient values.
+    // This final array is collected only to the IO processor.
     amrex::ParallelReduce::Max(
         pinned_tg_fab.dataPtr(), static_cast<int>(pinned_tg_fab.size()),
         amrex::ParallelDescriptor::IOProcessorNumber(),
         amrex::ParallelDescriptor::Communicator());
 
+    // To get the planar average over the non-normal directions (typically
+    // x and y), we sum up all the max temperature gradient values.  Again,
+    // this only happens on the IO processor.  Once we have the sum, divide
+    // by the number of coarse grid points in x and y.
     if (amrex::ParallelDescriptor::IOProcessor()) {
         const auto dnval = m_dn;
         auto* p = pinned_tg_fab.dataPtr();
@@ -276,6 +302,11 @@ void ABLStats::compute_zi()
             0.0);
         m_zi /= static_cast<amrex::Real>(pinned_tg_fab.size());
     }
+
+    // It is necessary for other processors to know <z_i>, so broadcast
+    // it back out to all processors.
+    amrex::ParallelDescriptor::Bcast(
+        &m_zi, 1, amrex::ParallelDescriptor::IOProcessorNumber());
 }
 
 void ABLStats::process_output()

amr-wind/wind_energy/ABLStatsBase.H

@@ -32,6 +32,8 @@ public:
     //! Flag indicating ABL simulation mode
     virtual ABLStatsMode abl_mode() const = 0;
 
+    virtual amrex::Real zi() const = 0;
+
     //! Interpolating object for vertical velocity profile
     virtual const VelPlaneAveraging& vel_profile_coarse() const = 0;
     virtual const VelPlaneAveragingFine& vel_profile() const = 0;

amr-wind/wind_energy/ABLWallFunction.H

@@ -53,13 +53,32 @@ private:
 
     amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 
+    //! Ability to read in a table of surface temperature flux versus time.
+    std::string m_surf_temp_flux_timetable;
+    amrex::Vector<amrex::Real> m_surf_temp_flux_time;
+    amrex::Vector<amrex::Real> m_surf_temp_flux_value;
+
     //! Ability to read in a table of surface temperature versus time.
     std::string m_surf_temp_timetable;
     amrex::Vector<amrex::Real> m_surf_temp_time;
     amrex::Vector<amrex::Real> m_surf_temp_value;
 
-    bool m_tempflux{true};
-    bool m_temp_table{false};
+    //! Ability to read in a table of near-surface temperature verus time.
+    std::string m_nearsurf_temp_timetable;
+    amrex::Vector<amrex::Real> m_nearsurf_temp_time;
+    amrex::Vector<amrex::Real> m_nearsurf_temp_value;
+
+    bool m_surf_temp_flux{true};
+    bool m_surf_temp{false};
+    bool m_nearsurf_temp{false};
+
+    bool m_surf_temp_flux_use_table{false};
+
+    bool m_surf_temp_use_table{false};
+    bool m_surf_temp_use_rate{false};
+
+    bool m_nearsurf_temp_use_table{false};
+
     amrex::Real m_surf_temp_rate{0.0};
     amrex::Real m_surf_temp_rate_tstart{0.0};
     amrex::Real m_surf_temp_init{300.0};