Merge branch 'main' into vof_fillpatch_fixes

CMakeLists.txt

@@ -44,6 +44,7 @@ option(AMR_WIND_ENABLE_OPENFAST "Enable OpenFAST integration" OFF)
 option(AMR_WIND_ENABLE_ASCENT "Enable Ascent visualization library" OFF)
 option(AMR_WIND_ENABLE_UMPIRE "Enable Umpire GPU memory pools" OFF)
 option(AMR_WIND_ENABLE_W2A "Enable Waves2AMR library" OFF)
+option(AMR_WIND_ENABLE_FFT "Enable FFT solver for MAC projection" OFF)
 
 #Options for C++
 set(CMAKE_CXX_STANDARD 17)
@@ -194,6 +195,10 @@ if(AMR_WIND_ENABLE_HELICS)
   target_compile_definitions(${amr_wind_lib_name} PUBLIC AMR_WIND_USE_HELICS)
 endif()
 
+if(AMR_WIND_ENABLE_FFT)
+  target_compile_definitions(${amr_wind_lib_name} PUBLIC AMR_WIND_USE_FFT)
+endif()
+
 if(AMR_WIND_TEST_WITH_PYTHON)
   set(CMAKE_PREFIX_PATH ${PYTHON_DIR} ${CMAKE_PREFIX_PATH})
   find_package(Python3 REQUIRED COMPONENTS Interpreter)

README.md

@@ -70,6 +70,10 @@ through pull-requests on GitHub. Please refer to the
 [coding guidelines](https://exawind.github.io/amr-wind/developer/coding_guidelines.html) as
 a reference for the best practices currently used to develop AMR-Wind.
 
+Please acknowledge as a publication co-author any developer that has
+significantly contributed to implementing or improving specific
+capability that was used for that publication.
+
 ### User discussion, feedback, and community support
 
 The development team manages a mailing list for AMR-Wind users. Invites for quarterly user meetings,

amr-wind/CFDSim.H

@@ -20,6 +20,10 @@ class OversetManager;
 class ExtSolverMgr;
 class HelicsStorage;
 
+namespace transport {
+class TransportModel;
+}
+
 namespace turbulence {
 class TurbulenceModel;
 }
@@ -74,6 +78,12 @@ public:
         return m_physics_mgr.objects();
     }
 
+    transport::TransportModel& transport_model() { return *m_transport; }
+    const transport::TransportModel& transport_model() const
+    {
+        return *m_transport;
+    }
+
     turbulence::TurbulenceModel& turbulence_model() { return *m_turbulence; }
     const turbulence::TurbulenceModel& turbulence_model() const
     {
@@ -103,6 +113,9 @@ public:
 
     bool has_overset() const;
 
+    //! Instantiate the transport model based on user inputs
+    void create_transport_model();
+
     //! Instantiate the turbulence model based on user inputs
     void create_turbulence_model();
 
@@ -137,6 +150,8 @@ private:
 
     PhysicsMgr m_physics_mgr;
 
+    std::unique_ptr<transport::TransportModel> m_transport;
+
     std::unique_ptr<turbulence::TurbulenceModel> m_turbulence;
 
     std::unique_ptr<IOManager> m_io_mgr;

amr-wind/CFDSim.cpp

@@ -1,4 +1,5 @@
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/turbulence/TurbulenceModel.H"
 #include "amr-wind/utilities/IOManager.H"
 #include "amr-wind/utilities/PostProcessing.H"
@@ -22,6 +23,17 @@ CFDSim::CFDSim(amrex::AmrCore& mesh)
 
 CFDSim::~CFDSim() = default;
 
+void CFDSim::create_transport_model()
+{
+    std::string transport_model_name = "ConstTransport";
+    {
+        amrex::ParmParse pp("transport");
+        pp.query("model", transport_model_name);
+    }
+    m_transport =
+        transport::TransportModel::create(transport_model_name, *this);
+}
+
 void CFDSim::create_turbulence_model()
 {
     std::string transport_model_name = "ConstTransport";

amr-wind/equation_systems/icns/icns_advection.H

@@ -56,6 +56,10 @@ private:
     FieldRepo& m_repo;
     PhysicsMgr& m_phy_mgr;
     std::unique_ptr<Hydro::MacProjector> m_mac_proj;
+#ifdef AMR_WIND_USE_FFT
+    std::unique_ptr<Hydro::FFTMacProjector> m_fft_mac_proj;
+    bool m_use_fft{true}; // use fft if possible
+#endif
     MLMGOptions m_options;
     bool m_has_overset{false};
     bool m_need_init{true};

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -1,3 +1,4 @@
+#include <AMReX.H>
 #include <memory>
 
 #include "amr-wind/equation_systems/icns/icns_advection.H"
@@ -25,15 +26,11 @@ amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> get_projection_bc(
             r[dir] = amrex::LinOpBCType::Periodic;
         } else {
             auto bc = bctype[amrex::Orientation(dir, side)];
-            switch (bc) {
-            case BC::pressure_outflow: {
+            if (bc == BC::pressure_outflow) {
                 r[dir] = amrex::LinOpBCType::Dirichlet;
-                break;
-            }
-            default:
+            } else {
                 r[dir] = amrex::LinOpBCType::Neumann;
-                break;
-            };
+            }
         }
     }
     return r;
@@ -58,6 +55,12 @@ MacProjOp::MacProjOp(
 {
     amrex::ParmParse pp("incflo");
     pp.query("density", m_rho_0);
+#ifdef AMR_WIND_USE_FFT
+    {
+        amrex::ParmParse pp("mac_proj");
+        pp.query("use_fft", m_use_fft);
+    }
+#endif
 }
 
 void MacProjOp::enforce_inout_solvability(
@@ -108,6 +111,32 @@ void MacProjOp::init_projector(const amrex::Real beta) noexcept
     }
 
     // Prepare projector
+    m_need_init = false;
+
+    auto& pressure = m_repo.get_field("p");
+    const auto& bctype = pressure.bc_type();
+
+    auto const& lobc = get_projection_bc(
+        amrex::Orientation::low, bctype, m_repo.mesh().Geom());
+    auto const& hibc = get_projection_bc(
+        amrex::Orientation::high, bctype, m_repo.mesh().Geom());
+
+#ifdef AMR_WIND_USE_FFT
+    if (m_use_fft) {
+        if (m_repo.num_active_levels() == 1 && m_has_overset == false) {
+            m_fft_mac_proj = std::make_unique<Hydro::FFTMacProjector>(
+                m_repo.mesh().Geom(0), lobc, hibc);
+            return;
+        } else {
+            amrex::ParmParse pp("mac_proj");
+            if (pp.contains("use_fft")) {
+                amrex::Print() << "WARNING: FFT MAC projection disabled due to "
+                                  "multiple levels/overset\n";
+            }
+        }
+    }
+#endif
+
     m_mac_proj = std::make_unique<Hydro::MacProjector>(
         m_repo.mesh().Geom(0, m_repo.num_active_levels() - 1));
     m_mac_proj->initProjector(
@@ -119,16 +148,7 @@ void MacProjOp::init_projector(const amrex::Real beta) noexcept
 
     m_options(*m_mac_proj);
 
-    auto& pressure = m_repo.get_field("p");
-    const auto& bctype = pressure.bc_type();
-
-    m_mac_proj->setDomainBC(
-        get_projection_bc(
-            amrex::Orientation::low, bctype, m_repo.mesh().Geom()),
-        get_projection_bc(
-            amrex::Orientation::high, bctype, m_repo.mesh().Geom()));
-
-    m_need_init = false;
+    m_mac_proj->setDomainBC(lobc, hibc);
 }
 
 void MacProjOp::set_inflow_velocity(amrex::Real time)
@@ -268,7 +288,12 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
         if (m_need_init) {
             init_projector(factor / m_rho_0);
         } else {
-            m_mac_proj->updateBeta(factor / m_rho_0);
+#ifdef AMR_WIND_USE_FFT
+            if (!m_fft_mac_proj)
+#endif
+            {
+                m_mac_proj->updateBeta(factor / m_rho_0);
+            }
         }
     }
 
@@ -291,7 +316,16 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
         enforce_inout_solvability(mac_vec);
     }
 
-    m_mac_proj->setUMAC(mac_vec);
+#ifdef AMR_WIND_USE_FFT
+    if (m_fft_mac_proj) {
+        // This is set on mac_vec[0] since FFT based projection is restricted to
+        // a single level
+        m_fft_mac_proj->setUMAC(mac_vec[0]);
+    } else
+#endif
+    {
+        m_mac_proj->setUMAC(mac_vec);
+    }
 
     if (m_has_overset) {
         auto phif = m_repo.create_scratch_field(1, 1, amr_wind::FieldLoc::CELL);
@@ -304,7 +338,14 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
             phif->vec_ptrs(), m_options.rel_tol, m_options.abs_tol);
 
     } else {
-        m_mac_proj->project(m_options.rel_tol, m_options.abs_tol);
+#ifdef AMR_WIND_USE_FFT
+        if (m_fft_mac_proj) {
+            m_fft_mac_proj->project();
+        } else
+#endif
+        {
+            m_mac_proj->project(m_options.rel_tol, m_options.abs_tol);
+        }
     }
 
     if (m_is_anelastic) {
@@ -322,7 +363,9 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
         amr_wind::overset_ops::revert_mask_cell_after_mac(m_repo);
     }
 
-    io::print_mlmg_info("MAC_projection", m_mac_proj->getMLMG());
+    if (m_mac_proj) {
+        io::print_mlmg_info("MAC_projection", m_mac_proj->getMLMG());
+    }
 }
 
 void MacProjOp::mac_proj_to_uniform_space(

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

@@ -2,6 +2,7 @@
 #define ABLMEANBOUSSINESQ_H
 
 #include "amr-wind/core/FieldRepo.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/equation_systems/icns/MomentumSource.H"
 #include "amr-wind/utilities/FieldPlaneAveraging.H"
 
@@ -42,11 +43,11 @@ private:
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_ht;
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_vals;
 
-    //! Reference temperature (Kelvin)
-    amrex::Real m_ref_theta{300.0};
+    //! Transport model
+    const transport::TransportModel& m_transport;
 
-    //! Thermal expansion coefficient
-    amrex::Real m_beta{0.0};
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 
     int m_axis{2};