Mesh map (Exawind#545)

* make nodal projection consistent with mesh mapping

* framework to take in mesh mapping based on user implement mesh mapping model

* remove alpha scaling for nodal projection

* minor modification to workflow for initializing mesh map

* fix field string identifier

* Fix where mapping happens in nodal projection. also scale back after nodal projection

* make different scaling arrays for cell centers, and nodes. Also change tgv_mol input file to use mesh mapping

* introduce mesh mapping to tgv physics

* Reorganize mesh mapping call to happen before setting any physics

* move mesh mapping manager too CFDSim

* fix bug in mesh mapping of tgv

* account for mesh mapping while computing time step size

* create field methods for mapping mesh to and from mapped spaces

* introduce mesh mapping to mac projection

* add additional checks for mesh mapping of fields

* add comments

* incorporate mesh mapping for computation of source term

* accomadate mesh mapping in diffusion solve

* add comments and general cleanup

* Output netcdf file from 3D field

* output components separately for paraview read

* bug fix

* Bug fix in coordinates

* Output netcdf files along with plot files

* provision for mapping velocity from previous time step

* changes for mesh mapping consistency

* changing mesh space for computation of convective fluxes

* fix bugs in scaling of RHS in presence of mesh mapping

* perform mesh mapping for velocity after prediction of face velocities

* change name of unmapped dimension to physical dimension

* add channel flow scaling

* fix bug with channel flow scaling

* accommodate mesh mapping in convecting taylor vortex

* modify BoussinesqBubble to allow for mesh mapping

* fix bug with handling of diffusion term when performing and implicit solve in presence of mesh mapping

* comment out writing of netcdf file because this does not work in parallel

* refactor velocity solve to account for space-dependent mapping

* Laminar Channel set up

* mesh mapping metric terms for multi-component diffusion solve

* remove obsolete line

* fix formatting

* change channel flow mesh mapping function

* fixes to mesh mapping framework

* revert mesh mapping changes to BoussinesqBubble

* remove mesh mapping from Taylor Green Vortex

* fix bug in channel flow scaling

* Introduce non uniform mesh coordinates as a field

* add framework to compute error norm for laminar channel flow

* fix channel flow analytical solution

* fix for loop iterator

* fix channel flow analytical solution

* scale source term to account for mesh mapping

* add mesh mapping regression test cases

* remove routine to write mesh mapped file to netcdf

* add TODO regarding unclear looping over size of tile boxes

* fix bug with source term and add explicit time stepping regression test

* introduce mesh mapping metric on cell faces

* refactor for a cleaner code

* add missing bracket

* fix formatting

* move to parreduce

* fix bug from rebase

* remove netcdf leftover

* changes to fix compilation on GPUs

* add blueprint for gradient computation in physical space

* update regression test settings

* Revert "add blueprint for gradient computation in physical space"

This reverts commit b116e02.

* fix unit tests

* fix channel flow linear solver

* remove left over mesh mapping settings

* add missing test

* clang-tidy cleanup

* remove left over obsolete input

* refactor mesh mapping

* remove unnecessary comments

* reintroduce constant coefficient pathway for projection in absence of mesh mapping

* change mesh mapping terminology for more clarity. also change access of mapping variables for cleaner code.

* clang-tidy

* store detJ in fields to avoid recomputation

* refactor laminar channel flow to allow flow in any direction

* add documentation

* code refactor to minimize for loops when not using mesh mapping

* update amrex module

* clean up

* make variable density go through single-component sigma pathway for nodal projection

* clean up all mesh mapping code leaks when turned off

* remove left over comments and fix clang-tidy errors

Co-authored-by: shashankNREL <[email protected]>

amr-wind/CFDSim.H

@@ -6,6 +6,7 @@
 #include "amr-wind/core/FieldRepo.H"
 #include "amr-wind/equation_systems/PDEBase.H"
 #include "amr-wind/core/Physics.H"
+#include "amr-wind/core/MeshMap.H"
 
 /** AMR-Wind
  *
@@ -89,6 +90,9 @@ public:
         return m_overset_mgr.get();
     }
 
+    MeshMap* mesh_mapping() { return m_mesh_map.get(); }
+    const MeshMap* mesh_mapping() const { return m_mesh_map.get(); }
+
     ExtSolverMgr& ext_solver_manager() { return *m_ext_solver_mgr; }
     const ExtSolverMgr& ext_solver_manager() const { return *m_ext_solver_mgr; }
 
@@ -103,6 +107,12 @@ public:
     //! Activate overset connectivity
     void activate_overset();
 
+    //! Activate mesh mapping
+    void activate_mesh_map();
+
+    bool has_mesh_mapping() { return m_mesh_mapping; }
+    bool has_mesh_mapping() const { return m_mesh_mapping; }
+
 private:
     amrex::AmrCore& m_mesh;
 
@@ -122,7 +132,11 @@ private:
 
     std::unique_ptr<OversetManager> m_overset_mgr;
 
+    std::unique_ptr<MeshMap> m_mesh_map;
+
     std::unique_ptr<ExtSolverMgr> m_ext_solver_mgr;
+
+    bool m_mesh_mapping{false};
 };
 
 } // namespace amr_wind

amr-wind/CFDSim.cpp

@@ -59,4 +59,15 @@ void CFDSim::activate_overset()
 
 bool CFDSim::has_overset() const { return (static_cast<bool>(m_overset_mgr)); }
 
+void CFDSim::activate_mesh_map()
+{
+    amrex::ParmParse pp("geometry");
+    std::string mesh_map_name; // default
+    m_mesh_mapping = static_cast<bool>(pp.query("mesh_mapping", mesh_map_name));
+    if (m_mesh_mapping) {
+        m_mesh_map = MeshMap::create(mesh_map_name);
+        m_mesh_map->declare_mapping_fields(*this, m_pde_mgr.num_ghost_state());
+    }
+}
+
 } // namespace amr_wind

amr-wind/CMakeLists.txt

@@ -31,6 +31,7 @@ add_subdirectory(turbulence)
 add_subdirectory(physics)
 add_subdirectory(overset)
 add_subdirectory(immersed_boundary)
+add_subdirectory(mesh_mapping_models)
 
 include(AMReXBuildInfo)
 generate_buildinfo(${amr_wind_lib_name} ${CMAKE_SOURCE_DIR})

amr-wind/core/CMakeLists.txt

@@ -8,4 +8,5 @@ target_sources(${amr_wind_lib_name}
   ScratchField.cpp
   ViewField.cpp
   MLMGOptions.cpp
+  MeshMap.cpp
   )

amr-wind/core/Field.H

@@ -225,6 +225,24 @@ public:
         const amrex::BCType::mathematicalBndryTypes bctype =
             amrex::BCType::hoextrap) noexcept;
 
+    /** Map field to the uniform mesh
+     *
+     *  This method transforms the field to the uniform mesh based on
+     *  the fields mesh_scaling_factor_cc or mesh_scaling_factor_nd
+     *  depending on whether the field is cell-centered or node-centered,
+     *  respectively.
+     */
+    void to_uniform_space() noexcept;
+
+    /** Map field to the stretched mesh
+     *
+     *  This method transforms the field to the stretched mesh based on
+     *  the fields mesh_scaling_factor_cc or mesh_scaling_factor_nd
+     *  depending on whether the field is cell-centered or node-centered,
+     *  respectively.
+     */
+    void to_stretched_space() noexcept;
+
     //! Return field at a different time state
     Field& state(const FieldState fstate);
     const Field& state(const FieldState fstate) const;
@@ -357,6 +375,9 @@ public:
         set_inflow(lev, time, mfab, amrex::IntVect(nghost));
     }
 
+    inline bool& in_uniform_space() { return m_mesh_mapped; }
+    inline bool in_uniform_space() const { return m_mesh_mapped; }
+
 protected:
     Field(
         FieldRepo& repo,
@@ -383,6 +404,9 @@ protected:
     //! Flag indicating whether fill patch operation must be performed for this
     //! field during regrid
     bool m_fillpatch_on_regrid{false};
+
+    //! Flag to track mesh mapping (to uniform space) of field
+    bool m_mesh_mapped{false};
 };
 
 } // namespace amr_wind

amr-wind/core/Field.cpp

@@ -347,4 +347,73 @@ void Field::set_default_fillpatch_bc(
     }
 }
 
+void Field::to_uniform_space() noexcept
+{
+    if (m_info->m_ncomp < AMREX_SPACEDIM) {
+        amrex::Abort("Trying to transform a non-vector field:" + m_name);
+    }
+    if (m_mesh_mapped) {
+        amrex::Print() << "WARNING: Field already in uniform mesh space: "
+                       << m_name << std::endl;
+        return;
+    }
+
+    const auto& mesh_fac = m_repo.get_mesh_mapping_field(m_info->m_floc);
+    const auto& mesh_detJ = m_repo.get_mesh_mapping_detJ(m_info->m_floc);
+
+    // scale velocity to accommodate for mesh mapping -> U^bar = U * J/fac
+    for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+        for (amrex::MFIter mfi(mesh_fac(lev)); mfi.isValid(); ++mfi) {
+
+            amrex::Array4<amrex::Real> const& field = operator()(lev).array(
+                mfi);
+            amrex::Array4<amrex::Real const> const& fac =
+                mesh_fac(lev).const_array(mfi);
+            amrex::Array4<amrex::Real const> const& detJ =
+                mesh_detJ(lev).const_array(mfi);
+
+            amrex::ParallelFor(
+                mfi.growntilebox(), AMREX_SPACEDIM,
+                [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
+                    field(i, j, k, n) *= detJ(i, j, k) / fac(i, j, k, n);
+                });
+        }
+    }
+    m_mesh_mapped = true;
+}
+
+void Field::to_stretched_space() noexcept
+{
+    if (m_info->m_ncomp < AMREX_SPACEDIM) {
+        amrex::Abort("Trying to transform a non-vector field:" + m_name);
+    }
+    if (!m_mesh_mapped) {
+        amrex::Print() << "WARNING: Field already in stretched mesh space: "
+                       << m_name << std::endl;
+        return;
+    }
+
+    const auto& mesh_fac = m_repo.get_mesh_mapping_field(m_info->m_floc);
+    const auto& mesh_detJ = m_repo.get_mesh_mapping_detJ(m_info->m_floc);
+
+    // scale field back to stretched mesh -> U = U^bar * fac/J
+    for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+        for (amrex::MFIter mfi(mesh_fac(lev)); mfi.isValid(); ++mfi) {
+            amrex::Array4<amrex::Real> const& field = operator()(lev).array(
+                mfi);
+            amrex::Array4<amrex::Real const> const& fac =
+                mesh_fac(lev).const_array(mfi);
+            amrex::Array4<amrex::Real const> const& detJ =
+                mesh_detJ(lev).const_array(mfi);
+
+            amrex::ParallelFor(
+                mfi.growntilebox(), AMREX_SPACEDIM,
+                [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
+                    field(i, j, k, n) *= fac(i, j, k, n) / detJ(i, j, k);
+                });
+        }
+    }
+    m_mesh_mapped = false;
+}
+
 } // namespace amr_wind

amr-wind/core/FieldRepo.H

@@ -248,6 +248,14 @@ public:
         return *m_field_vec.at(field_id);
     }
 
+    /** Return a previously created mesh mapping field using field location
+     */
+    Field& get_mesh_mapping_field(FieldLoc floc) const;
+
+    /** Return a previously created mesh mapping detJ using field location
+     */
+    Field& get_mesh_mapping_detJ(FieldLoc floc) const;
+
     //! Query if field uniquely identified by name and time state exists in
     //! repository
     bool field_exists(

amr-wind/core/FieldRepo.cpp

@@ -159,6 +159,56 @@ FieldRepo::get_field(const std::string& name, const FieldState fstate) const
     return *m_field_vec[found->second];
 }
 
+Field& FieldRepo::get_mesh_mapping_field(FieldLoc floc) const
+{
+    Field* fac = nullptr;
+    switch (floc) {
+    case FieldLoc::CELL:
+        fac = &(get_field("mesh_scaling_factor_cc"));
+        break;
+    case FieldLoc::NODE:
+        fac = &(get_field("mesh_scaling_factor_nd"));
+        break;
+    case FieldLoc::XFACE:
+        fac = &(get_field("mesh_scaling_factor_xf"));
+        break;
+    case FieldLoc::YFACE:
+        fac = &(get_field("mesh_scaling_factor_yf"));
+        break;
+    case FieldLoc::ZFACE:
+        fac = &(get_field("mesh_scaling_factor_zf"));
+        break;
+    default:
+        amrex::Abort("Invalid field location");
+    }
+    return *fac;
+}
+
+Field& FieldRepo::get_mesh_mapping_detJ(FieldLoc floc) const
+{
+    Field* detJ = nullptr;
+    switch (floc) {
+    case FieldLoc::CELL:
+        detJ = &(get_field("mesh_scaling_detJ_cc"));
+        break;
+    case FieldLoc::NODE:
+        detJ = &(get_field("mesh_scaling_detJ_nd"));
+        break;
+    case FieldLoc::XFACE:
+        detJ = &(get_field("mesh_scaling_detJ_xf"));
+        break;
+    case FieldLoc::YFACE:
+        detJ = &(get_field("mesh_scaling_detJ_yf"));
+        break;
+    case FieldLoc::ZFACE:
+        detJ = &(get_field("mesh_scaling_detJ_zf"));
+        break;
+    default:
+        amrex::Abort("Invalid field location");
+    }
+    return *detJ;
+}
+
 bool FieldRepo::field_exists(
     const std::string& name, const FieldState fstate) const
 {

amr-wind/core/MeshMap.H

@@ -0,0 +1,62 @@
+#ifndef MESHMAP_H
+#define MESHMAP_H
+
+#include "amr-wind/core/CollMgr.H"
+#include "amr-wind/core/Factory.H"
+#include "amr-wind/core/Field.H"
+
+#include "AMReX_MultiFab.H"
+#include "AMReX_Geometry.H"
+
+namespace amr_wind {
+
+class CFDSim;
+
+/**
+ *  \defgroup mesh_map Mesh mapping models
+ *
+ *  AMR-Wind representation of different mesh mapping models
+ *
+ *  In AMR-Wind, different mesh mappings are implemented using the MeshMap
+ * class.
+ */
+
+/** Abstract representation of different mesh mapping models
+ *
+ *  This class defines an abstract API that represents the notion of some
+ *  mesh mapping that will be used to scale the mesh. The most common use-case
+ *  for this class is to perform RANS simulations.
+ */
+class MeshMap : public Factory<MeshMap>
+{
+public:
+    static const std::string base_identifier() { return "MeshMap"; }
+
+    virtual ~MeshMap() = default;
+
+    //! declare mesh mapping fields
+    void declare_mapping_fields(const CFDSim&, int);
+
+    //! Construct mesh scaling field
+    virtual void create_map(int, const amrex::Geometry&) = 0;
+
+protected:
+    Field* m_mesh_scale_fac_cc{nullptr};
+    Field* m_mesh_scale_fac_nd{nullptr};
+    Field* m_mesh_scale_fac_xf{nullptr};
+    Field* m_mesh_scale_fac_yf{nullptr};
+    Field* m_mesh_scale_fac_zf{nullptr};
+
+    Field* m_mesh_scale_detJ_cc{nullptr};
+    Field* m_mesh_scale_detJ_nd{nullptr};
+    Field* m_mesh_scale_detJ_xf{nullptr};
+    Field* m_mesh_scale_detJ_yf{nullptr};
+    Field* m_mesh_scale_detJ_zf{nullptr};
+
+    Field* m_non_uniform_coord_cc{nullptr};
+    Field* m_non_uniform_coord_nd{nullptr};
+};
+
+} // namespace amr_wind
+
+#endif /* MESHMAP_H */

amr-wind/core/MeshMap.cpp

@@ -0,0 +1,42 @@
+#include "amr-wind/core/MeshMap.H"
+#include "amr-wind/CFDSim.H"
+
+namespace amr_wind {
+
+void MeshMap::declare_mapping_fields(const CFDSim& sim, int nghost)
+{
+
+    // declare nodal, cell-centered, and face-centered mesh mapping array
+    m_mesh_scale_fac_cc = &(sim.repo().declare_cc_field(
+        "mesh_scaling_factor_cc", AMREX_SPACEDIM, nghost, 1));
+    m_mesh_scale_fac_nd = &(sim.repo().declare_nd_field(
+        "mesh_scaling_factor_nd", AMREX_SPACEDIM, nghost, 1));
+    m_mesh_scale_fac_xf = &(sim.repo().declare_xf_field(
+        "mesh_scaling_factor_xf", AMREX_SPACEDIM, nghost, 1));
+    m_mesh_scale_fac_yf = &(sim.repo().declare_yf_field(
+        "mesh_scaling_factor_yf", AMREX_SPACEDIM, nghost, 1));
+    m_mesh_scale_fac_zf = &(sim.repo().declare_zf_field(
+        "mesh_scaling_factor_zf", AMREX_SPACEDIM, nghost, 1));
+
+    // declare nodal, cell-centered, and face-centered mesh mapping detJ array
+    m_mesh_scale_detJ_cc =
+        &(sim.repo().declare_cc_field("mesh_scaling_detJ_cc", 1, nghost, 1));
+    m_mesh_scale_detJ_nd =
+        &(sim.repo().declare_nd_field("mesh_scaling_detJ_nd", 1, nghost, 1));
+    m_mesh_scale_detJ_xf =
+        &(sim.repo().declare_xf_field("mesh_scaling_detJ_xf", 1, nghost, 1));
+    m_mesh_scale_detJ_yf =
+        &(sim.repo().declare_yf_field("mesh_scaling_detJ_yf", 1, nghost, 1));
+    m_mesh_scale_detJ_zf =
+        &(sim.repo().declare_zf_field("mesh_scaling_detJ_zf", 1, nghost, 1));
+
+    // declare nodal and cell-centered non-uniform mesh
+    m_non_uniform_coord_cc = &(sim.repo().declare_cc_field(
+        "non_uniform_coord_cc", AMREX_SPACEDIM, nghost, 1));
+    m_non_uniform_coord_nd = &(sim.repo().declare_nd_field(
+        "non_uniform_coord_nd", AMREX_SPACEDIM, nghost, 1));
+
+    // TODO: Create BCNoOP fill patch operators for mesh scaling fields ?
+}
+
+} // namespace amr_wind

amr-wind/diffusion/diffusion.H

@@ -39,6 +39,11 @@ void fixup_eta_on_domain_faces(
     amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>& fc,
     amrex::MultiFab const& cc);
 
+void viscosity_to_uniform_space(
+    amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>& b,
+    const amr_wind::FieldRepo& repo,
+    int lev);
+
 } // namespace diffusion
 
 #endif /* DIFFUSION_H */