Changes for coupling with ERF using MultiBlock

CMakeLists.txt

@@ -3,7 +3,7 @@
 cmake_minimum_required (VERSION 3.20 FATAL_ERROR)
 project(AMR-Wind CXX C)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON CACHE INTERNAL "")
-list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
+list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake")
 include(CMakePackageConfigHelpers)
 include(GNUInstallDirs)
 include(amr-wind-utils)
@@ -129,6 +129,10 @@ endif()
 
 include(set_compile_flags)
 
+if(ERF_AW_MULTIBLOCK)
+    target_compile_definitions(${amr_wind_lib_name} PUBLIC ERF_AW_MULTIBLOCK)
+endif()
+
 if (AMR_WIND_ENABLE_W2A)
   # Turn on the use of the library via bool, it has already been included above
   target_compile_definitions(${amr_wind_lib_name} PUBLIC AMR_WIND_USE_W2A)
@@ -185,7 +189,7 @@ if(AMR_WIND_ENABLE_ASCENT)
   target_compile_definitions(${amr_wind_lib_name} PRIVATE AMR_WIND_USE_ASCENT)
 endif()
 
-# Link with HELICS module 
+# Link with HELICS module
 if(AMR_WIND_ENABLE_HELICS)
   set(CMAKE_PREFIX_PATH ${HELICS_DIR} ${CMAKE_PREFIX_PATH})
   find_package(HELICS 3 REQUIRED)

amr-wind/CFDSim.H

@@ -9,6 +9,11 @@
 #include "amr-wind/core/MeshMap.H"
 #include "amr-wind/helics.H"
 
+#ifdef ERF_AW_MULTIBLOCK
+#include "amr-wind/wind_energy/ABLReadERFFcn.H"
+class MultiBlockContainer;
+#endif
+
 /** AMR-Wind
  *
  *  All C++ code in AMR-Wind is organized within the amr_wind namespace.
@@ -58,6 +63,14 @@ public:
     SimTime& time() { return m_time; }
     const SimTime& time() const { return m_time; }
 
+#ifdef ERF_AW_MULTIBLOCK
+    void set_mbc(MultiBlockContainer* mbc) { m_mbc = mbc; }
+    MultiBlockContainer** mbc() { return &m_mbc; }
+
+    void set_read_erf(ReadERFFcn fcn) { m_read_erf = fcn; }
+    ReadERFFcn* get_read_erf() { return &m_read_erf; }
+#endif
+
     //! Return the field repository
     FieldRepo& repo() { return m_repo; }
     FieldRepo& repo() const { return m_repo; }
@@ -145,6 +158,11 @@ private:
     std::unique_ptr<HelicsStorage> m_helics;
 
     bool m_mesh_mapping{false};
+
+#ifdef ERF_AW_MULTIBLOCK
+    MultiBlockContainer* m_mbc = nullptr;
+    ReadERFFcn m_read_erf {nullptr};
+#endif
 };
 
 } // namespace amr_wind

amr-wind/CMakeLists.txt

@@ -38,7 +38,7 @@ include(AMReXBuildInfo)
 generate_buildinfo(${amr_wind_lib_name} ${CMAKE_SOURCE_DIR})
 
 # Generate AMR-Wind version header
-configure_file("${CMAKE_SOURCE_DIR}/cmake/AMRWindVersion.H.in"
+configure_file("${PROJECT_SOURCE_DIR}/cmake/AMRWindVersion.H.in"
   "${CMAKE_CURRENT_BINARY_DIR}/AMRWindVersion.H" @ONLY)
 
 target_link_libraries_system(${amr_wind_lib_name} PUBLIC AMReX::amrex AMReX-Hydro::amrex_hydro_api)

amr-wind/incflo.H

@@ -13,6 +13,10 @@
 #include "amr-wind/core/FieldRepo.H"
 #include "amr-wind/overset/OversetOps.H"
 
+#ifdef ERF_AW_MULTIBLOCK
+#include "amr-wind/wind_energy/ABLReadERFFcn.H"
+class MultiBlockContainer;
+#endif
 namespace amr_wind {
 namespace pde {
 class PDEBase;
@@ -58,6 +62,11 @@ public:
     // Evolve solution to final time through repeated calls to Advance()
     void Evolve();
 
+#ifdef ERF_AW_MULTIBLOCK
+    // Evolve solution in a multiblock framework by desired number of steps
+    void Evolve_MB(int MBstep, int max_block_step);
+#endif
+
     // Tag cells for refinement
     void
     ErrorEst(int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow)
@@ -150,6 +159,11 @@ public:
 
     void ReadCheckpointFile();
 
+#ifdef ERF_AW_MULTIBLOCK
+    void SetMultiBlockPointer(MultiBlockContainer *mbc) { m_sim.set_mbc(mbc); }
+    void set_read_erf(ReadERFFcn fcn) {m_sim.set_read_erf(fcn);}
+#endif
+
 private:
     //
     // member variables

amr-wind/incflo.cpp

@@ -23,6 +23,10 @@ incflo::incflo()
     // constructor. No valid BoxArray and DistributionMapping have been defined.
     // But the arrays for them have been resized.
 
+#ifdef ERF_AW_MULTIBLOCK
+    amrex::Print() << std::endl << "AMR-Wind is initializing.. "<< std::endl << std::endl;
+#endif
+
     // Check if dry run is requested and set up if so
     CheckAndSetUpDryRun();
 

amr-wind/utilities/index_operations.H

@@ -19,6 +19,17 @@ closest_index(const amrex::Vector<amrex::Real>& vec, const amrex::Real value)
     return std::max(idx - 1, 0);
 }
 
+//! This version is required for the ERF coupling
+AMREX_FORCE_INLINE int
+closest_index_lbound(const amrex::Vector<amrex::Real>& vec, const amrex::Real value)
+{
+    auto const it = std::lower_bound(vec.begin(), vec.end(), value);
+    AMREX_ALWAYS_ASSERT(it != vec.end());
+
+    const int idx = static_cast<int>(std::distance(vec.begin(), it));
+    return std::max(idx - 1, 0);
+}
+
 //! Return indices perpendicular to normal
 template <typename T = amrex::GpuArray<int, 2>>
 AMREX_FORCE_INLINE T perpendicular_idx(const int normal)

amr-wind/wind_energy/ABLBoundaryPlane.H

@@ -7,6 +7,10 @@
 #include "amr-wind/utilities/ncutils/nc_interface.H"
 #include <AMReX_BndryRegister.H>
 
+#ifdef ERF_AW_MULTIBLOCK
+#include "amr-wind/wind_energy/ABLReadERFFcn.H"
+class MultiBlockContainer;
+#endif
 namespace amr_wind {
 
 enum struct io_mode { output, input, undefined };
@@ -57,7 +61,8 @@ public:
         const int lev,
         const Field* /*fld*/,
         const amrex::Real time,
-        const amrex::Vector<amrex::Real>& /*times*/);
+        const amrex::Vector<amrex::Real>& /*times*/,
+        const bool erf_multiblock = false);
 
     void interpolate(const amrex::Real /*time*/);
     bool is_populated(amrex::Orientation /*ori*/) const;
@@ -158,11 +163,24 @@ public:
 
     io_mode mode() const { return m_io_mode; }
 
+#ifdef ERF_AW_MULTIBLOCK
+    MultiBlockContainer* mbc() {return *m_mbc;}
+#endif
+
 private:
     const amr_wind::SimTime& m_time;
     const FieldRepo& m_repo;
     const amrex::AmrCore& m_mesh;
 
+#ifdef ERF_AW_MULTIBLOCK
+    // pointer to pointer : when ABL boundary plane gets intiialized, amr-wind
+    // CFDsim does not yet have the up to date pointer to the mbc or the read_erf function
+    // be careful: if this is used before amr-wind gets the mbc pointer set
+    // bad things will happen
+    MultiBlockContainer** m_mbc;
+    ReadERFFcn* m_read_erf {nullptr};
+#endif
+
 #ifdef AMR_WIND_USE_NETCDF
     void write_data(
         const ncutils::NCGroup& grp,

amr-wind/wind_energy/ABLBoundaryPlane.cpp

@@ -144,13 +144,22 @@ void InletData::read_data_native(
     const int lev,
     const Field* fld,
     const amrex::Real time,
-    const amrex::Vector<amrex::Real>& times)
+    const amrex::Vector<amrex::Real>& times,
+    const bool erf_multiblock)
 {
     const size_t nc = fld->num_comp();
     const int nstart =
         static_cast<int>(m_components[static_cast<int>(fld->id())]);
 
-    const int idx = utils::closest_index(times, time);
+    int idx;
+    if (erf_multiblock) {
+        // the time is often equal to the new time for multiblock
+        // hence the need for a lower-bound version of this function
+        // otherwise it will abort
+        idx = utils::closest_index_lbound(times, time);
+    } else {
+        idx = utils::closest_index(times, time);
+    }
     const int idxp1 = idx + 1;
 
     m_tn = times[idx];
@@ -245,6 +254,9 @@ bool InletData::is_populated(amrex::Orientation ori) const
 
 ABLBoundaryPlane::ABLBoundaryPlane(CFDSim& sim)
     : m_time(sim.time()), m_repo(sim.repo()), m_mesh(sim.mesh())
+#ifdef ERF_AW_MULTIBLOCK
+    , m_mbc(sim.mbc()), m_read_erf(sim.get_read_erf())
+#endif
 {
     amrex::ParmParse pp("ABL");
     int pp_io_mode = -1;
@@ -281,7 +293,7 @@ ABLBoundaryPlane::ABLBoundaryPlane(CFDSim& sim)
     }
 #endif
 
-    if (!(m_out_fmt == "native" || m_out_fmt == "netcdf")) {
+    if (!(m_out_fmt == "native" || m_out_fmt == "netcdf" || m_out_fmt == "erf-multiblock")) {
         amrex::Print() << "Warning: boundary output format not recognized, "
                           "changing to native format"
                        << std::endl;
@@ -336,6 +348,9 @@ void ABLBoundaryPlane::initialize_data()
                 "ABLBoundaryPlane: invalid variable requested: " + fname);
         }
     }
+    if (m_io_mode == io_mode::output and m_out_fmt == "erf-multiblock") {
+      amrex::Abort("ABLBoundaryPlane: can't output data in erf-multiblock mode");
+    }
 }
 
 void ABLBoundaryPlane::write_header()
@@ -748,6 +763,34 @@ void ABLBoundaryPlane::read_header()
                 m_in_data.define_level_data(ori, pbx, nc);
             }
         }
+    } else if (m_out_fmt == "erf-multiblock") {
+
+        m_in_times.push_back(-1.0e13); // create space for storing time at erf old and new timestep
+        m_in_times.push_back(-1.0e13);
+
+        int nc = 0;
+        for (auto* fld : m_fields) {
+            m_in_data.component(static_cast<int>(fld->id())) = nc;
+            nc += fld->num_comp();
+        }
+
+        // FIXME: need to generalize to lev > 0 somehow
+        const int lev = 0;
+        for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
+            auto ori = oit();
+            // FIXME: would be safer and less storage to not allocate all of
+            // these but we do not use m_planes for input and need to detect
+            // mass inflow from field bcs same for define level data below
+            m_in_data.define_plane(ori);
+            const amrex::Box& minBox = m_mesh.boxArray(lev).minimalBox();
+            amrex::IntVect plo(minBox.loVect());
+            amrex::IntVect phi(minBox.hiVect());
+            const int normal = ori.coordDir();
+            plo[normal] = ori.isHigh() ? minBox.hiVect()[normal] + 1 : -1;
+            phi[normal] = ori.isHigh() ? minBox.hiVect()[normal] + 1 : -1;
+            const amrex::Box pbx(plo, phi);
+            m_in_data.define_level_data(ori, pbx, nc);
+        }
     }
 }
 
@@ -758,6 +801,20 @@ void ABLBoundaryPlane::read_file()
         return;
     }
 
+#ifdef ERF_AW_MULTIBLOCK
+    if (m_out_fmt == "erf-multiblock") {
+        //m_read_erf = sim.get_read_erf();
+        ReadERFFcn read_erf = *m_read_erf;
+        if (read_erf)
+        {
+            read_erf(m_time, m_in_times, m_in_data, m_fields, mbc());
+        } else {
+            amrex::Abort("read_erf function is undefined.");
+        }
+        return;
+    }
+#endif
+
     // populate planes and interpolate
     const amrex::Real time = m_time.new_time();
     AMREX_ALWAYS_ASSERT((m_in_times[0] <= time) && (time < m_in_times.back()));

amr-wind/wind_energy/ABLFieldInit.H

@@ -47,6 +47,20 @@ public:
     void init_tke(const amrex::Geometry& geom, amrex::MultiFab& tke) const;
 
 private:
+    // ABL-with-bubble
+
+    //! Initial bubble location
+    amrex::Vector<amrex::Real> m_bubble_loc{{0.25, 0.25, 0.5}};
+
+    //! Initial bubble radius
+    amrex::Real m_bubble_radius{0.1};
+
+    //! Initial bubble
+    amrex::Real m_bubble_temp_ratio{1.1};
+
+    //! Whether or not to have temperature bubble
+    bool m_use_bubble{0};
+
     //! Initial velocity components
     amrex::Vector<amrex::Real> m_vel;
 

amr-wind/wind_energy/ABLFieldInit.cpp

@@ -24,6 +24,17 @@ void ABLFieldInit::initialize_from_inputfile()
 {
     amrex::ParmParse pp_abl("ABL");
 
+    // ABL-with-bubble
+    pp_abl.query("use_bubble", m_use_bubble);
+    if (m_use_bubble) {
+      amrex::Print() << "Initializing with bubble" << std::endl;
+      pp_abl.getarr("bubble_loc", m_bubble_loc);
+      pp_abl.query("bubble_radius", m_bubble_radius);
+      pp_abl.query("bubble_temp_ratio", m_bubble_temp_ratio);
+    } else {
+      amrex::Print() << "Not initializing with bubble" << std::endl;
+    }
+
     // Temperature variation as a function of height
     pp_abl.getarr("temperature_heights", m_theta_heights);
     pp_abl.getarr("temperature_values", m_theta_values);
@@ -219,8 +230,14 @@ void ABLFieldInit::operator()(
         const amrex::Real bottom_v_vel = m_bottom_vel[1];
         const amrex::Real bottom_w_vel = m_bottom_vel[2];
 
+        const amrex::Real bcx = m_bubble_loc[0];
+        const amrex::Real bcy = m_bubble_loc[1];
+        const amrex::Real bcz = m_bubble_loc[2];
+
         amrex::ParallelFor(
             vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real x = problo[0] + (i + 0.5) * dx[0];
+                const amrex::Real y = problo[1] + (j + 0.5) * dx[1];
                 const amrex::Real z = problo[2] + (k + 0.5) * dx[2];
 
                 density(i, j, k) = rho_init;
@@ -240,6 +257,13 @@ void ABLFieldInit::operator()(
 
                 temperature(i, j, k, 0) += theta;
 
+                amrex::Real ratio = 1.0;
+                if (m_use_bubble) {
+                  amrex::Real radius = std::sqrt((x-bcx)*(x-bcx) + (y-bcy)*(y-bcy) + (z-bcz)*(z-bcz));
+                  ratio = 1.0 + (m_bubble_temp_ratio - 1.0) * exp(-0.5 * radius* radius / (m_bubble_radius * m_bubble_radius));
+                }
+                temperature(i, j, k, 0) *= ratio;
+
                 if (linear_profile) {
                     velocity(i, j, k, 0) =
                         bottom_u_vel + z * (top_u_vel - bottom_u_vel) /

amr-wind/wind_energy/ABLReadERFFcn.H

@@ -0,0 +1,17 @@
+#ifndef ABLREADERFFCN_H
+#define ABLREADERFFCN_H
+#include <functional>
+
+class MultiBlockContainer;
+namespace amr_wind {
+    class InletData;
+}
+
+using ReadERFFcn = std::function <void (
+                    const amr_wind::SimTime&,
+                    amrex::Vector<amrex::Real>&,
+                    amr_wind::InletData&,
+                    const amrex::Vector<amr_wind::Field*>&,
+                    MultiBlockContainer* )>;
+
+#endif /* ABLREADERFFCN_H */