From 584df175f4e42507224ed5d8723ae0802b5ed45f Mon Sep 17 00:00:00 2001 From: ashesh2512 <36968394+ashesh2512@users.noreply.github.com> Date: Thu, 17 Feb 2022 18:03:07 -0700 Subject: [PATCH] Mesh map (#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 b116e029286df5a0839a6e56e18af7bbc234d3bd. * 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 --- amr-wind/CFDSim.H | 14 + amr-wind/CFDSim.cpp | 11 + amr-wind/CMakeLists.txt | 1 + amr-wind/core/CMakeLists.txt | 1 + amr-wind/core/Field.H | 24 ++ amr-wind/core/Field.cpp | 69 ++++ amr-wind/core/FieldRepo.H | 8 + amr-wind/core/FieldRepo.cpp | 50 +++ amr-wind/core/MeshMap.H | 62 ++++ amr-wind/core/MeshMap.cpp | 42 +++ amr-wind/diffusion/diffusion.H | 5 + amr-wind/diffusion/incflo_diffusion.cpp | 63 ++++ amr-wind/equation_systems/AdvOp_Godunov.H | 3 +- amr-wind/equation_systems/AdvOp_MOL.H | 3 +- amr-wind/equation_systems/CompRHSOps.H | 84 ++++- amr-wind/equation_systems/DiffusionOps.H | 38 +- amr-wind/equation_systems/DiffusionOps.cpp | 28 +- amr-wind/equation_systems/PDE.H | 22 +- amr-wind/equation_systems/PDEOps.H | 2 +- .../equation_systems/density/density_ops.H | 6 +- .../equation_systems/icns/icns_advection.H | 38 +- .../equation_systems/icns/icns_advection.cpp | 111 +++++- .../equation_systems/icns/icns_diffusion.H | 75 +++- amr-wind/equation_systems/icns/icns_ops.H | 25 +- .../equation_systems/levelset/levelset_ops.H | 6 +- amr-wind/equation_systems/sdr/sdr_ops.H | 6 +- amr-wind/equation_systems/tke/tke_ops.H | 6 +- amr-wind/equation_systems/vof/vof_advection.H | 2 +- amr-wind/equation_systems/vof/vof_ops.H | 4 +- amr-wind/incflo.cpp | 27 +- amr-wind/incflo_advance.cpp | 28 +- amr-wind/incflo_compute_dt.cpp | 103 ++++-- amr-wind/mesh_mapping_models/CMakeLists.txt | 6 + amr-wind/mesh_mapping_models/ChannelFlowMap.H | 43 +++ .../mesh_mapping_models/ChannelFlowMap.cpp | 338 ++++++++++++++++++ amr-wind/mesh_mapping_models/ConstantMap.H | 41 +++ amr-wind/mesh_mapping_models/ConstantMap.cpp | 169 +++++++++ amr-wind/physics/ChannelFlow.H | 26 ++ amr-wind/physics/ChannelFlow.cpp | 262 +++++++++++--- amr-wind/physics/ConvectingTaylorVortex.H | 3 + amr-wind/physics/ConvectingTaylorVortex.cpp | 98 +++-- .../incflo_apply_nodal_projection.cpp | 87 ++++- amr-wind/setup/init.cpp | 4 + amr-wind/setup/set_background_pressure.cpp | 18 +- amr-wind/utilities/IOManager.cpp | 4 + docs/sphinx/theory/mapping.rst | 130 +++++++ docs/sphinx/theory/theory.rst | 11 + docs/sphinx/user/inputs_geometry.rst | 15 + submods/amrex | 2 +- test/CMakeLists.txt | 5 + .../channel_mol_mesh_map_x.i | 87 +++++ .../channel_mol_mesh_map_y.i | 91 +++++ .../channel_mol_mesh_map_z.i | 91 +++++ .../ctv_mol_mesh_map/ctv_mol_mesh_map.i | 53 +++ .../ctv_mol_mesh_map_explicit.i | 53 +++ unit_tests/equation_systems/test_pde.cpp | 2 + 56 files changed, 2384 insertions(+), 222 deletions(-) create mode 100644 amr-wind/core/MeshMap.H create mode 100644 amr-wind/core/MeshMap.cpp create mode 100644 amr-wind/mesh_mapping_models/CMakeLists.txt create mode 100644 amr-wind/mesh_mapping_models/ChannelFlowMap.H create mode 100644 amr-wind/mesh_mapping_models/ChannelFlowMap.cpp create mode 100644 amr-wind/mesh_mapping_models/ConstantMap.H create mode 100644 amr-wind/mesh_mapping_models/ConstantMap.cpp create mode 100644 docs/sphinx/theory/mapping.rst create mode 100644 test/test_files/channel_mol_mesh_map_x/channel_mol_mesh_map_x.i create mode 100644 test/test_files/channel_mol_mesh_map_y/channel_mol_mesh_map_y.i create mode 100644 test/test_files/channel_mol_mesh_map_z/channel_mol_mesh_map_z.i create mode 100644 test/test_files/ctv_mol_mesh_map/ctv_mol_mesh_map.i create mode 100644 test/test_files/ctv_mol_mesh_map_explicit/ctv_mol_mesh_map_explicit.i diff --git a/amr-wind/CFDSim.H b/amr-wind/CFDSim.H index 41b67c2aed..60e33e4978 100644 --- a/amr-wind/CFDSim.H +++ b/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 m_overset_mgr; + std::unique_ptr m_mesh_map; + std::unique_ptr m_ext_solver_mgr; + + bool m_mesh_mapping{false}; }; } // namespace amr_wind diff --git a/amr-wind/CFDSim.cpp b/amr-wind/CFDSim.cpp index 69055ca471..380d6283c6 100644 --- a/amr-wind/CFDSim.cpp +++ b/amr-wind/CFDSim.cpp @@ -59,4 +59,15 @@ void CFDSim::activate_overset() bool CFDSim::has_overset() const { return (static_cast(m_overset_mgr)); } +void CFDSim::activate_mesh_map() +{ + amrex::ParmParse pp("geometry"); + std::string mesh_map_name; // default + m_mesh_mapping = static_cast(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 diff --git a/amr-wind/CMakeLists.txt b/amr-wind/CMakeLists.txt index a4b1a50956..b68753a0c2 100644 --- a/amr-wind/CMakeLists.txt +++ b/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}) diff --git a/amr-wind/core/CMakeLists.txt b/amr-wind/core/CMakeLists.txt index f5c926f574..b84b8f6901 100644 --- a/amr-wind/core/CMakeLists.txt +++ b/amr-wind/core/CMakeLists.txt @@ -8,4 +8,5 @@ target_sources(${amr_wind_lib_name} ScratchField.cpp ViewField.cpp MLMGOptions.cpp + MeshMap.cpp ) diff --git a/amr-wind/core/Field.H b/amr-wind/core/Field.H index 90fdd11581..20bfe1445f 100644 --- a/amr-wind/core/Field.H +++ b/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 diff --git a/amr-wind/core/Field.cpp b/amr-wind/core/Field.cpp index 4cf77946c0..9732c66d46 100644 --- a/amr-wind/core/Field.cpp +++ b/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 const& field = operator()(lev).array( + mfi); + amrex::Array4 const& fac = + mesh_fac(lev).const_array(mfi); + amrex::Array4 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 const& field = operator()(lev).array( + mfi); + amrex::Array4 const& fac = + mesh_fac(lev).const_array(mfi); + amrex::Array4 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 diff --git a/amr-wind/core/FieldRepo.H b/amr-wind/core/FieldRepo.H index b4d79efda6..3d1cbc9216 100644 --- a/amr-wind/core/FieldRepo.H +++ b/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( diff --git a/amr-wind/core/FieldRepo.cpp b/amr-wind/core/FieldRepo.cpp index bdaad3f6f7..6104bcce6d 100644 --- a/amr-wind/core/FieldRepo.cpp +++ b/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 { diff --git a/amr-wind/core/MeshMap.H b/amr-wind/core/MeshMap.H new file mode 100644 index 0000000000..a4a5f6932f --- /dev/null +++ b/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 +{ +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 */ diff --git a/amr-wind/core/MeshMap.cpp b/amr-wind/core/MeshMap.cpp new file mode 100644 index 0000000000..f096b13596 --- /dev/null +++ b/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 diff --git a/amr-wind/diffusion/diffusion.H b/amr-wind/diffusion/diffusion.H index 18687c5af1..6b78f84400 100644 --- a/amr-wind/diffusion/diffusion.H +++ b/amr-wind/diffusion/diffusion.H @@ -39,6 +39,11 @@ void fixup_eta_on_domain_faces( amrex::Array& fc, amrex::MultiFab const& cc); +void viscosity_to_uniform_space( + amrex::Array& b, + const amr_wind::FieldRepo& repo, + int lev); + } // namespace diffusion #endif /* DIFFUSION_H */ diff --git a/amr-wind/diffusion/incflo_diffusion.cpp b/amr-wind/diffusion/incflo_diffusion.cpp index 9a3b1b3123..c09080b044 100644 --- a/amr-wind/diffusion/incflo_diffusion.cpp +++ b/amr-wind/diffusion/incflo_diffusion.cpp @@ -197,4 +197,67 @@ void fixup_eta_on_domain_faces( } } } + +void viscosity_to_uniform_space( + amrex::Array& b, + const amr_wind::FieldRepo& repo, + int lev) +{ + const auto& mesh_fac_xf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::XFACE); + const auto& mesh_fac_yf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::YFACE); + const auto& mesh_fac_zf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::ZFACE); + const auto& mesh_detJ_xf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::XFACE); + const auto& mesh_detJ_yf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::YFACE); + const auto& mesh_detJ_zf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::ZFACE); + + // beta accounted for mesh mapping (x-face) = J/fac^2 * mu + for (amrex::MFIter mfi(b[0]); mfi.isValid(); ++mfi) { + amrex::Array4 const& mu = b[0].array(mfi); + amrex::Array4 const& fac = + mesh_fac_xf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_xf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + mu(i, j, k) = + mu(i, j, k) * detJ(i, j, k) / std::pow(fac(i, j, k, 0), 2); + }); + } + // beta accounted for mesh mapping (y-face) = J/fac^2 * mu + for (amrex::MFIter mfi(b[1]); mfi.isValid(); ++mfi) { + amrex::Array4 const& mu = b[1].array(mfi); + amrex::Array4 const& fac = + mesh_fac_yf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_yf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + mu(i, j, k) = + mu(i, j, k) * detJ(i, j, k) / std::pow(fac(i, j, k, 1), 2); + }); + } + // beta accounted for mesh mapping (z-face) = J/fac^2 * mu + for (amrex::MFIter mfi(b[2]); mfi.isValid(); ++mfi) { + amrex::Array4 const& mu = b[2].array(mfi); + amrex::Array4 const& fac = + mesh_fac_zf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_zf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + mu(i, j, k) = + mu(i, j, k) * detJ(i, j, k) / std::pow(fac(i, j, k, 2), 2); + }); + } +} + } // namespace diffusion diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index d5008e02e5..ae14e822d8 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -27,7 +27,8 @@ struct AdvectionOp< AdvectionOp( PDEFields& fields_in, bool /* has_overset */, - bool /* variable density */) + bool /* variable density */, + bool /* mesh mapping */) : fields(fields_in) , density(fields_in.repo.get_field("density")) , u_mac(fields_in.repo.get_field("u_mac")) diff --git a/amr-wind/equation_systems/AdvOp_MOL.H b/amr-wind/equation_systems/AdvOp_MOL.H index 89d3316d45..ec95acca3e 100644 --- a/amr-wind/equation_systems/AdvOp_MOL.H +++ b/amr-wind/equation_systems/AdvOp_MOL.H @@ -27,7 +27,8 @@ struct AdvectionOp< AdvectionOp( PDEFields& fields_in, bool /* has_overset */, - bool /* variable density */) + bool /* variable density */, + bool /* mesh mapping */) : fields(fields_in) , density(fields_in.repo.get_field("density")) , u_mac(fields_in.repo.get_field("u_mac")) diff --git a/amr-wind/equation_systems/CompRHSOps.H b/amr-wind/equation_systems/CompRHSOps.H index 08d0d01166..bfc0b996ad 100644 --- a/amr-wind/equation_systems/CompRHSOps.H +++ b/amr-wind/equation_systems/CompRHSOps.H @@ -26,7 +26,8 @@ struct ComputeRHSOp * \param difftype Indicating whether time-integration is explicit/implicit * \param dt time step size */ - void predictor_rhs(const DiffusionType difftype, const amrex::Real dt) + void predictor_rhs( + const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping) { amrex::Real factor = 0.0; switch (difftype) { @@ -53,14 +54,26 @@ struct ComputeRHSOp : FieldState::Old; const int nlevels = fields.repo.num_active_levels(); + + // for RHS evaluation velocity field should be in stretched space auto& field = fields.field; + if (field.in_uniform_space() && mesh_mapping) { + field.to_stretched_space(); + } auto& field_old = field.state(FieldState::Old); + if (field_old.in_uniform_space() && mesh_mapping) { + field_old.to_stretched_space(); + } + auto& den_new = density.state(FieldState::New); auto& den_old = density.state(FieldState::Old); auto& src_term = fields.src_term; auto& diff_term = fields.diff_term.state(fstate); auto& conv_term = fields.conv_term.state(fstate); auto& mask_cell = fields.repo.get_int_field("mask_cell"); + Field const* mesh_detJ = + mesh_mapping ? &(fields.repo.get_mesh_mapping_detJ(FieldLoc::CELL)) + : nullptr; for (int lev = 0; lev < nlevels; ++lev) { #ifdef _OPENMP @@ -76,6 +89,9 @@ struct ComputeRHSOp const auto diff = diff_term(lev).const_array(mfi); const auto ddt_o = conv_term(lev).const_array(mfi); const auto imask = mask_cell(lev).const_array(mfi); + amrex::Array4 detJ = + mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) + : amrex::Array4(); if (PDE::multiply_rho) { // Remove multiplication by density as it will be added back @@ -84,24 +100,40 @@ struct ComputeRHSOp bx, PDE::ndim, [=] AMREX_GPU_DEVICE( int i, int j, int k, int n) noexcept { + amrex::Real det_j = + mesh_mapping ? (detJ(i, j, k)) : 1.0; + fld(i, j, k, n) = - rho_o(i, j, k) * fld_o(i, j, k, n) + + rho_o(i, j, k) * det_j * fld_o(i, j, k, n) + static_cast(imask(i, j, k)) * dt * - (ddt_o(i, j, k, n) + src(i, j, k, n) + + (ddt_o(i, j, k, n) + + det_j * src(i, j, k, n) + factor * diff(i, j, k, n)); fld(i, j, k, n) /= rho(i, j, k); + + if (difftype == DiffusionType::Explicit) { + fld(i, j, k, n) /= det_j; + } }); } else { amrex::ParallelFor( bx, PDE::ndim, [=] AMREX_GPU_DEVICE( int i, int j, int k, int n) noexcept { + amrex::Real det_j = + mesh_mapping ? (detJ(i, j, k)) : 1.0; + fld(i, j, k, n) = - fld_o(i, j, k, n) + + det_j * fld_o(i, j, k, n) + static_cast(imask(i, j, k)) * dt * - (ddt_o(i, j, k, n) + src(i, j, k, n) + + (ddt_o(i, j, k, n) + + det_j * src(i, j, k, n) + factor * diff(i, j, k, n)); + + if (difftype == DiffusionType::Explicit) { + fld(i, j, k, n) /= det_j; + } }); } } @@ -113,7 +145,8 @@ struct ComputeRHSOp * \param difftype Indicating whether time-integration is explicit/implicit * \param dt time step size */ - void corrector_rhs(const DiffusionType difftype, const amrex::Real dt) + void corrector_rhs( + const DiffusionType difftype, const amrex::Real dt, bool mesh_mapping) { amrex::Real ofac = 0.0; amrex::Real nfac = 0.0; @@ -138,8 +171,17 @@ struct ComputeRHSOp } const int nlevels = fields.repo.num_active_levels(); + + // for RHS evaluation velocity field should be in stretched space auto& field = fields.field; + if (field.in_uniform_space() && mesh_mapping) { + field.to_stretched_space(); + } auto& field_old = field.state(FieldState::Old); + if (field_old.in_uniform_space() && mesh_mapping) { + field_old.to_stretched_space(); + } + auto& den_new = density.state(FieldState::New); auto& den_old = density.state(FieldState::Old); auto& src_term = fields.src_term; @@ -148,6 +190,9 @@ struct ComputeRHSOp auto& diff_term_old = fields.diff_term.state(FieldState::Old); auto& conv_term_old = fields.conv_term.state(FieldState::Old); auto& mask_cell = fields.repo.get_int_field("mask_cell"); + Field const* mesh_detJ = + mesh_mapping ? &(fields.repo.get_mesh_mapping_detJ(FieldLoc::CELL)) + : nullptr; for (int lev = 0; lev < nlevels; ++lev) { #ifdef _OPENMP @@ -165,6 +210,9 @@ struct ComputeRHSOp const auto diff_o = diff_term_old(lev).const_array(mfi); const auto ddt_o = conv_term_old(lev).const_array(mfi); const auto imask = mask_cell(lev).const_array(mfi); + amrex::Array4 detJ = + mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) + : amrex::Array4(); if (PDE::multiply_rho) { // Remove multiplication by density as it will be added back @@ -173,28 +221,44 @@ struct ComputeRHSOp bx, PDE::ndim, [=] AMREX_GPU_DEVICE( int i, int j, int k, int n) noexcept { + amrex::Real det_j = + mesh_mapping ? (detJ(i, j, k)) : 1.0; + fld(i, j, k, n) = - rho_o(i, j, k) * fld_o(i, j, k, n) + + rho_o(i, j, k) * det_j * fld_o(i, j, k, n) + static_cast(imask(i, j, k)) * dt * (0.5 * (ddt_o(i, j, k, n) + ddt(i, j, k, n)) + ofac * diff_o(i, j, k, n) + - nfac * diff(i, j, k, n) + src(i, j, k, n)); + nfac * diff(i, j, k, n) + + det_j * src(i, j, k, n)); fld(i, j, k, n) /= rho(i, j, k); + + if (difftype == DiffusionType::Explicit) { + fld(i, j, k, n) /= det_j; + } }); } else { amrex::ParallelFor( bx, PDE::ndim, [=] AMREX_GPU_DEVICE( int i, int j, int k, int n) noexcept { + amrex::Real det_j = + mesh_mapping ? (detJ(i, j, k)) : 1.0; + fld(i, j, k, n) = - fld_o(i, j, k, n) + + det_j * fld_o(i, j, k, n) + static_cast(imask(i, j, k)) * dt * (0.5 * (ddt_o(i, j, k, n) + ddt(i, j, k, n)) + ofac * diff_o(i, j, k, n) + - nfac * diff(i, j, k, n) + src(i, j, k, n)); + nfac * diff(i, j, k, n) + + det_j * src(i, j, k, n)); + + if (difftype == DiffusionType::Explicit) { + fld(i, j, k, n) /= det_j; + } }); } } diff --git a/amr-wind/equation_systems/DiffusionOps.H b/amr-wind/equation_systems/DiffusionOps.H index 42e83f3ca8..fd29bce1f8 100644 --- a/amr-wind/equation_systems/DiffusionOps.H +++ b/amr-wind/equation_systems/DiffusionOps.H @@ -12,6 +12,7 @@ #include "AMReX_MLTensorOp.H" #include "AMReX_MLMG.H" #include "AMReX_AmrCore.H" +#include "AMReX_MultiFabUtil.H" namespace amr_wind { namespace pde { @@ -32,6 +33,7 @@ public: DiffSolverIface( PDEFields& fields, const bool has_overset, + const bool mesh_mapping, const std::string& prefix = "diffusion"); virtual ~DiffSolverIface() = default; @@ -44,16 +46,6 @@ public: virtual void linsys_solve_impl(); -protected: - //! Sets up the linear operator (e.g., setup BCs, etc.) - virtual void setup_operator( - LinOp& linop, - const amrex::Real alpha, - const amrex::Real beta, - const FieldState fstate); - - virtual void setup_solver(amrex::MLMG& mlmg); - virtual void set_acoeffs(LinOp& linop, const FieldState fstate); template @@ -65,9 +57,13 @@ protected: const int nlevels = m_pdefields.repo.num_active_levels(); const auto& viscosity = m_pdefields.mueff; const auto& geom = m_pdefields.repo.mesh().Geom(); + for (int lev = 0; lev < nlevels; ++lev) { auto b = diffusion::average_velocity_eta_to_faces( geom[lev], viscosity(lev)); + if (m_mesh_mapping) { + diffusion::viscosity_to_uniform_space(b, m_pdefields.repo, lev); + } linop.setShearViscosity(lev, amrex::GetArrOfConstPtrs(b)); } } @@ -81,18 +77,34 @@ protected: const int nlevels = m_pdefields.repo.num_active_levels(); const auto& viscosity = m_pdefields.mueff; const auto& geom = m_pdefields.repo.mesh().Geom(); + for (int lev = 0; lev < nlevels; ++lev) { auto b = diffusion::average_velocity_eta_to_faces( geom[lev], viscosity(lev)); + if (m_mesh_mapping) { + diffusion::viscosity_to_uniform_space(b, m_pdefields.repo, lev); + } linop.setBCoeffs(lev, amrex::GetArrOfConstPtrs(b)); } } +protected: + //! Sets up the linear operator (e.g., setup BCs, etc.) + virtual void setup_operator( + LinOp& linop, + const amrex::Real alpha, + const amrex::Real beta, + const FieldState fstate); + + virtual void setup_solver(amrex::MLMG& mlmg); + PDEFields& m_pdefields; Field& m_density; MLMGOptions m_options; + bool m_mesh_mapping{false}; + std::unique_ptr m_solver; std::unique_ptr m_applier; }; @@ -113,8 +125,10 @@ struct DiffusionOp< std::is_same::value, "Invalid linear operator for scalar diffusion operator"); - DiffusionOp(PDEFields& fields, const bool has_overset) - : DiffSolverIface(fields, has_overset) + DiffusionOp( + PDEFields& fields, const bool has_overset, const bool mesh_mapping) + : DiffSolverIface( + fields, has_overset, mesh_mapping) { this->m_solver->setDomainBC( diffusion::get_diffuse_scalar_bc( diff --git a/amr-wind/equation_systems/DiffusionOps.cpp b/amr-wind/equation_systems/DiffusionOps.cpp index 49ea5235c5..57d6071a36 100644 --- a/amr-wind/equation_systems/DiffusionOps.cpp +++ b/amr-wind/equation_systems/DiffusionOps.cpp @@ -8,10 +8,14 @@ namespace pde { template DiffSolverIface::DiffSolverIface( - PDEFields& fields, const bool has_overset, const std::string& prefix) + PDEFields& fields, + const bool has_overset, + const bool mesh_mapping, + const std::string& prefix) : m_pdefields(fields) , m_density(fields.repo.get_field("density")) , m_options(prefix, m_pdefields.field.name() + "_" + prefix) + , m_mesh_mapping(mesh_mapping) { amrex::LPInfo isolve = m_options.lpinfo(); amrex::LPInfo iapply; @@ -73,9 +77,24 @@ void DiffSolverIface::set_acoeffs(LinOp& linop, const FieldState fstate) auto& repo = m_pdefields.repo; const int nlevels = repo.num_active_levels(); const auto& density = m_density.state(fstate); + Field const* mesh_detJ = m_mesh_mapping + ? &(repo.get_mesh_mapping_detJ(FieldLoc::CELL)) + : nullptr; + std::unique_ptr rho_times_detJ = + m_mesh_mapping ? repo.create_scratch_field( + 1, m_density.num_grow()[0], FieldLoc::CELL) + : nullptr; for (int lev = 0; lev < nlevels; ++lev) { - linop.setACoeffs(lev, density(lev)); + if (m_mesh_mapping) { + (*rho_times_detJ)(lev).setVal(0.0); + amrex::MultiFab::AddProduct( + (*rho_times_detJ)(lev), density(lev), 0, (*mesh_detJ)(lev), 0, + 0, 1, m_density.num_grow()[0]); + linop.setACoeffs(lev, (*rho_times_detJ)(lev)); + } else { + linop.setACoeffs(lev, density(lev)); + } } } @@ -93,6 +112,11 @@ void DiffSolverIface::linsys_solve_impl() FieldState fstate = FieldState::New; auto& repo = this->m_pdefields.repo; auto& field = this->m_pdefields.field; + if (field.in_uniform_space()) { + amrex::Abort( + "For diffusion solve, velocity should not be in uniform mesh " + "space."); + } const auto& density = m_density.state(fstate); const int nlevels = repo.num_active_levels(); const int ndim = field.num_comp(); diff --git a/amr-wind/equation_systems/PDE.H b/amr-wind/equation_systems/PDE.H index 2b87b8caf7..f900083814 100644 --- a/amr-wind/equation_systems/PDE.H +++ b/amr-wind/equation_systems/PDE.H @@ -59,8 +59,8 @@ public: { if (PDE::has_diffusion) { BL_PROFILE("amr-wind::" + this->identifier() + "::initialize"); - m_diff_op.reset( - new DiffusionOp(m_fields, m_sim.has_overset())); + m_diff_op.reset(new DiffusionOp( + m_fields, m_sim.has_overset(), m_sim.has_mesh_mapping())); m_turb_op.reset( new TurbulenceOp(m_sim.turbulence_model(), m_fields)); } @@ -70,7 +70,8 @@ public: m_sim.physics_manager().contains("MultiPhase")); m_adv_op.reset(new AdvectionOp( - m_fields, m_sim.has_overset(), variable_density)); + m_fields, m_sim.has_overset(), variable_density, + m_sim.has_mesh_mapping())); m_src_op.init_source_terms(m_sim); } @@ -80,8 +81,8 @@ public: if (PDE::has_diffusion) { BL_PROFILE( "amr-wind::" + this->identifier() + "::post_regrid_actions"); - m_diff_op.reset( - new DiffusionOp(m_fields, m_sim.has_overset())); + m_diff_op.reset(new DiffusionOp( + m_fields, m_sim.has_overset(), m_sim.has_mesh_mapping())); } bool variable_density = @@ -89,7 +90,8 @@ public: m_sim.physics_manager().contains("MultiPhase")); m_adv_op.reset(new AdvectionOp( - m_fields, m_sim.has_overset(), variable_density)); + m_fields, m_sim.has_overset(), variable_density, + m_sim.has_mesh_mapping())); } //! Return the object holding the fields necessary for solving this PDE @@ -99,7 +101,7 @@ public: void compute_source_term(const FieldState fstate) override { BL_PROFILE("amr-wind::" + this->identifier() + "::compute_source_term"); - m_src_op(fstate); + m_src_op(fstate, m_sim.has_mesh_mapping()); } void compute_mueff(const FieldState) override @@ -138,14 +140,16 @@ public: { BL_PROFILE( "amr-wind::" + this->identifier() + "::compute_predictor_rhs"); - m_rhs_op.predictor_rhs(difftype, m_time.deltaT()); + m_rhs_op.predictor_rhs( + difftype, m_time.deltaT(), m_sim.has_mesh_mapping()); } virtual void compute_corrector_rhs(const DiffusionType difftype) override { BL_PROFILE( "amr-wind::" + this->identifier() + "::compute_corrector_rhs"); - m_rhs_op.corrector_rhs(difftype, m_time.deltaT()); + m_rhs_op.corrector_rhs( + difftype, m_time.deltaT(), m_sim.has_mesh_mapping()); } void solve(const amrex::Real dt) override diff --git a/amr-wind/equation_systems/PDEOps.H b/amr-wind/equation_systems/PDEOps.H index 4e83d9331b..1ee1b584cc 100644 --- a/amr-wind/equation_systems/PDEOps.H +++ b/amr-wind/equation_systems/PDEOps.H @@ -119,7 +119,7 @@ struct SrcTermOpBase } //! Update source terms during time-integration procedure - void operator()(const FieldState fstate) + void operator()(const FieldState fstate, const bool /* mesh_mapping */) { // Zero out source term this->fields.src_term.setVal(0.0); diff --git a/amr-wind/equation_systems/density/density_ops.H b/amr-wind/equation_systems/density/density_ops.H index c16bab1dea..a6d05e8803 100644 --- a/amr-wind/equation_systems/density/density_ops.H +++ b/amr-wind/equation_systems/density/density_ops.H @@ -11,7 +11,8 @@ struct ComputeRHSOp { explicit ComputeRHSOp(PDEFields& fields_in) : fields(fields_in) {} - void predictor_rhs(const DiffusionType, const amrex::Real dt) + void predictor_rhs( + const DiffusionType, const amrex::Real dt, bool /*mesh_mapping*/) { // Field states for diffusion and advection terms. In Godunov scheme // these terms only have one state. @@ -45,7 +46,8 @@ struct ComputeRHSOp } } - void corrector_rhs(const DiffusionType, const amrex::Real dt) + void corrector_rhs( + const DiffusionType, const amrex::Real dt, bool /*mesh_mapping*/) { const int nlevels = fields.repo.num_active_levels(); auto& field = fields.field; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index f6ea6cec87..2f208ec964 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -17,10 +17,19 @@ public: using FaceFabPtrVec = amrex::Vector>; - MacProjOp(FieldRepo&, bool, bool); + MacProjOp(FieldRepo&, bool, bool, bool); void operator()(const FieldState fstate, const amrex::Real dt); + static void mac_proj_to_uniform_space( + const amr_wind::FieldRepo&, + amr_wind::Field&, + amr_wind::Field&, + amr_wind::Field&, + amrex::Array&, + amrex::Real, + int) noexcept; + private: void init_projector(const FaceFabPtrVec&) noexcept; void init_projector(const amrex::Real) noexcept; @@ -31,6 +40,7 @@ private: bool m_has_overset{false}; bool m_need_init{true}; bool m_variable_density{false}; + bool m_mesh_mapping{false}; amrex::Real m_rho_0{1.0}; }; @@ -40,12 +50,16 @@ private: template <> struct AdvectionOp { - AdvectionOp(PDEFields& fields_in, bool has_overset, bool variable_density) + AdvectionOp( + PDEFields& fields_in, + bool has_overset, + bool variable_density, + bool mesh_mapping) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) , w_mac(fields_in.repo.get_field("w_mac")) - , m_macproj_op(fields.repo, has_overset, variable_density) + , m_macproj_op(fields.repo, has_overset, variable_density, mesh_mapping) { amrex::ParmParse pp("incflo"); @@ -337,12 +351,18 @@ struct AdvectionOp template <> struct AdvectionOp { - AdvectionOp(PDEFields& fields_in, bool has_overset, bool variable_density) + AdvectionOp( + PDEFields& fields_in, + bool has_overset, + bool variable_density, + bool mesh_mapping) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) , w_mac(fields_in.repo.get_field("w_mac")) - , m_macproj_op(fields.repo, has_overset, variable_density) + , m_mesh_mapping(mesh_mapping) + , m_macproj_op( + fields.repo, has_overset, variable_density, m_mesh_mapping) {} void preadvect(const FieldState fstate, const amrex::Real dt) @@ -352,6 +372,12 @@ struct AdvectionOp const auto& geom = repo.mesh().Geom(); auto& dof_field = fields.field.state(fstate); + // computation of velocity on faces requires + // dof field to be in uniform mesh space + if (dof_field.in_uniform_space() && m_mesh_mapping) { + dof_field.to_stretched_space(); + } + // // Predict velocities // @@ -439,6 +465,8 @@ struct AdvectionOp Field& v_mac; Field& w_mac; + bool m_mesh_mapping; + MacProjOp m_macproj_op; }; diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp index c6b180d32f..3137f4db6d 100644 --- a/amr-wind/equation_systems/icns/icns_advection.cpp +++ b/amr-wind/equation_systems/icns/icns_advection.cpp @@ -41,11 +41,13 @@ amrex::Array get_projection_bc( } // namespace -MacProjOp::MacProjOp(FieldRepo& repo, bool has_overset, bool variable_density) +MacProjOp::MacProjOp( + FieldRepo& repo, bool has_overset, bool variable_density, bool mesh_mapping) : m_repo(repo) , m_options("mac_proj") , m_has_overset(has_overset) , m_variable_density(variable_density) + , m_mesh_mapping(mesh_mapping) { amrex::ParmParse pp("incflo"); pp.query("rho_0", m_rho_0); @@ -130,32 +132,27 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) auto& w_mac = m_repo.get_field("w_mac"); const auto& density = m_repo.get_field("density", fstate); - // This will hold (1/rho) on faces - std::unique_ptr rho_xf, rho_yf, rho_zf; - - amrex::Vector> rho_face( - m_repo.num_active_levels()); amrex::Vector> mac_vec( m_repo.num_active_levels()); - // fixme todo clean this up, this was done to replace - // GetVecOfArrOfConstPtrs() below - amrex::Vector> - rho_face_const; - rho_face_const.reserve(m_repo.num_active_levels()); - amrex::Real factor = m_has_overset ? 0.5 * dt : 1.0; // TODO: remove the or in the if statement for m_has_overset // For now assume variable viscosity for overset // this can be removed once the nsolve overset // masking is implemented in cell based AMReX poisson solvers + if (m_variable_density || m_has_overset || m_mesh_mapping) { + amrex::Vector> + rho_face_const; + rho_face_const.reserve(m_repo.num_active_levels()); - if (m_variable_density || m_has_overset) { - + // This will hold density on faces + std::unique_ptr rho_xf, rho_yf, rho_zf; rho_xf = m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::XFACE); rho_yf = m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::YFACE); rho_zf = m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::ZFACE); + amrex::Vector> rho_face( + m_repo.num_active_levels()); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { rho_face[lev][0] = &(*rho_xf)(lev); @@ -164,8 +161,14 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) amrex::average_cellcenter_to_face( rho_face[lev], density(lev), geom[lev]); - for (int idim = 0; idim < ICNS::ndim; ++idim) { - rho_face[lev][idim]->invert(factor, 0); + + if (m_mesh_mapping) { + mac_proj_to_uniform_space( + m_repo, u_mac, v_mac, w_mac, rho_face[lev], factor, lev); + } else { + for (int idim = 0; idim < ICNS::ndim; ++idim) { + rho_face[lev][idim]->invert(factor, 0); + } } rho_face_const.push_back(GetArrOfConstPtrs(rho_face[lev])); @@ -176,9 +179,7 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) } else { m_mac_proj->updateBeta(rho_face_const); } - } else { - if (m_need_init) { init_projector(factor / m_rho_0); } else { @@ -212,5 +213,79 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) io::print_mlmg_info("MAC_projection", m_mac_proj->getMLMG()); } +void MacProjOp::mac_proj_to_uniform_space( + const amr_wind::FieldRepo& repo, + amr_wind::Field& u_mac, + amr_wind::Field& v_mac, + amr_wind::Field& w_mac, + amrex::Array& rho_face, + amrex::Real ovst_fac, + int lev) noexcept +{ + const auto& mesh_fac_xf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::XFACE); + const auto& mesh_fac_yf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::YFACE); + const auto& mesh_fac_zf = + repo.get_mesh_mapping_field(amr_wind::FieldLoc::ZFACE); + const auto& mesh_detJ_xf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::XFACE); + const auto& mesh_detJ_yf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::YFACE); + const auto& mesh_detJ_zf = + repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::ZFACE); + + // scale U^mac to accommodate for mesh mapping -> U^bar = J/fac * + // U^mac beta accounted for mesh mapping = J/fac^2 * 1/rho construct + // rho and mesh map u_mac on x-face + for (amrex::MFIter mfi(*(rho_face[0])); mfi.isValid(); ++mfi) { + amrex::Array4 const& u = u_mac(lev).array(mfi); + amrex::Array4 const& rho = rho_face[0]->array(mfi); + amrex::Array4 const& fac = + mesh_fac_xf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_xf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + u(i, j, k) *= detJ(i, j, k) / fac(i, j, k, 0); + rho(i, j, k) = ovst_fac * detJ(i, j, k) / + std::pow(fac(i, j, k, 0), 2) / rho(i, j, k); + }); + } + // construct rho on y-face + for (amrex::MFIter mfi(*(rho_face[1])); mfi.isValid(); ++mfi) { + amrex::Array4 const& v = v_mac(lev).array(mfi); + amrex::Array4 const& rho = rho_face[1]->array(mfi); + amrex::Array4 const& fac = + mesh_fac_yf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_yf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + v(i, j, k) *= detJ(i, j, k) / fac(i, j, k, 1); + rho(i, j, k) = ovst_fac * detJ(i, j, k) / + std::pow(fac(i, j, k, 1), 2) / rho(i, j, k); + }); + } + // construct rho on z-face + for (amrex::MFIter mfi(*(rho_face[2])); mfi.isValid(); ++mfi) { + amrex::Array4 const& w = w_mac(lev).array(mfi); + amrex::Array4 const& rho = rho_face[2]->array(mfi); + amrex::Array4 const& fac = + mesh_fac_zf(lev).array(mfi); + amrex::Array4 const& detJ = + mesh_detJ_zf(lev).const_array(mfi); + + amrex::ParallelFor( + mfi.tilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + w(i, j, k) *= detJ(i, j, k) / fac(i, j, k, 2); + rho(i, j, k) = ovst_fac * detJ(i, j, k) / + std::pow(fac(i, j, k, 2), 2) / rho(i, j, k); + }); + } +} + } // namespace pde } // namespace amr_wind diff --git a/amr-wind/equation_systems/icns/icns_diffusion.H b/amr-wind/equation_systems/icns/icns_diffusion.H index d4e4a14b78..bcfddb2f92 100644 --- a/amr-wind/equation_systems/icns/icns_diffusion.H +++ b/amr-wind/equation_systems/icns/icns_diffusion.H @@ -14,8 +14,9 @@ namespace pde { class ICNSDiffTensorOp : public DiffSolverIface { public: - ICNSDiffTensorOp(PDEFields& fields, const bool has_overset) - : DiffSolverIface(fields, has_overset) + ICNSDiffTensorOp( + PDEFields& fields, const bool has_overset, const bool mesh_mapping) + : DiffSolverIface(fields, has_overset, mesh_mapping) { this->m_solver->setDomainBC( diffusion::get_diffuse_tensor_bc( @@ -74,10 +75,12 @@ public: ICNSDiffScalarOp( PDEFields& fields, const bool has_overset, + const bool mesh_mapping, const std::string& prefix = "diffusion") : m_pdefields(fields) , m_density(fields.repo.get_field("density")) , m_options(prefix, m_pdefields.field.name() + "_" + prefix) + , m_mesh_mapping(mesh_mapping) { amrex::LPInfo isolve = m_options.lpinfo(); amrex::LPInfo iapply; @@ -136,15 +139,38 @@ public: auto& divtau = m_pdefields.diff_term.state(tau_state); const auto& density = m_density.state(fstate); const auto& viscosity = m_pdefields.mueff; + Field const* mesh_detJ = + m_mesh_mapping ? &(repo.get_mesh_mapping_detJ(FieldLoc::CELL)) + : nullptr; + std::unique_ptr rho_times_detJ = + m_mesh_mapping ? repo.create_scratch_field( + 1, m_density.num_grow()[0], FieldLoc::CELL) + : nullptr; + const amrex::Real alpha = 0.0; const amrex::Real beta = -1.0; - m_applier_scalar->setScalars(alpha, beta); + for (int lev = 0; lev < nlevels; ++lev) { m_applier_scalar->setLevelBC(lev, &m_pdefields.field(lev)); - m_applier_scalar->setACoeffs(lev, density(lev)); + + // A coeffs + if (m_mesh_mapping) { + (*rho_times_detJ)(lev).setVal(0.0); + amrex::MultiFab::AddProduct( + (*rho_times_detJ)(lev), density(lev), 0, (*mesh_detJ)(lev), + 0, 0, 1, m_density.num_grow()[0]); + m_applier_scalar->setACoeffs(lev, (*rho_times_detJ)(lev)); + } else { + m_applier_scalar->setACoeffs(lev, density(lev)); + } + + // B coeffs auto b = diffusion::average_velocity_eta_to_faces( geom[lev], viscosity(lev)); + if (m_mesh_mapping) { + diffusion::viscosity_to_uniform_space(b, repo, lev); + } m_applier_scalar->setBCoeffs(lev, amrex::GetArrOfConstPtrs(b)); } @@ -176,22 +202,45 @@ public: { const FieldState fstate = FieldState::New; auto& repo = m_pdefields.repo; + const auto& geom = repo.mesh().Geom(); auto& field = m_pdefields.field; const auto& density = m_density.state(fstate); const int nlevels = repo.num_active_levels(); const int ndim = field.num_comp(); auto rhs_ptr = repo.create_scratch_field("rhs", field.num_comp(), 0); const auto& viscosity = m_pdefields.mueff; - const auto& geom = repo.mesh().Geom(); + Field const* mesh_detJ = + m_mesh_mapping ? &(repo.get_mesh_mapping_detJ(FieldLoc::CELL)) + : nullptr; + std::unique_ptr rho_times_detJ = + m_mesh_mapping ? repo.create_scratch_field( + 1, m_density.num_grow()[0], FieldLoc::CELL) + : nullptr; + const amrex::Real alpha = 1.0; const amrex::Real beta = dt; - m_solver_scalar->setScalars(alpha, beta); + for (int lev = 0; lev < nlevels; ++lev) { m_solver_scalar->setLevelBC(lev, &m_pdefields.field(lev)); - m_solver_scalar->setACoeffs(lev, density(lev)); + + // A coeffs + if (m_mesh_mapping) { + (*rho_times_detJ)(lev).setVal(0.0); + amrex::MultiFab::AddProduct( + (*rho_times_detJ)(lev), density(lev), 0, (*mesh_detJ)(lev), + 0, 0, 1, m_density.num_grow()[0]); + m_solver_scalar->setACoeffs(lev, (*rho_times_detJ)(lev)); + } else { + m_solver_scalar->setACoeffs(lev, density(lev)); + } + + // B coeffs auto b = diffusion::average_velocity_eta_to_faces( geom[lev], viscosity(lev)); + if (m_mesh_mapping) { + diffusion::viscosity_to_uniform_space(b, repo, lev); + } m_solver_scalar->setBCoeffs(lev, amrex::GetArrOfConstPtrs(b)); } @@ -230,6 +279,7 @@ protected: PDEFields& m_pdefields; Field& m_density; MLMGOptions m_options; + bool m_mesh_mapping{false}; std::unique_ptr m_solver_scalar; std::unique_ptr m_applier_scalar; @@ -248,7 +298,8 @@ struct DiffusionOp ICNS::ndim == AMREX_SPACEDIM, "DiffusionOp invoked for scalar PDE type"); - DiffusionOp(PDEFields& fields, const bool has_overset) + DiffusionOp( + PDEFields& fields, const bool has_overset, const bool mesh_mapping) { bool use_tensor_op = true; @@ -256,11 +307,11 @@ struct DiffusionOp pp.query("use_tensor_operator", use_tensor_op); if (use_tensor_op) { - m_tensor_op = - std::make_unique(fields, has_overset); + m_tensor_op = std::make_unique( + fields, has_overset, mesh_mapping); } else { - m_scalar_op = - std::make_unique(fields, has_overset); + m_scalar_op = std::make_unique( + fields, has_overset, mesh_mapping); } } diff --git a/amr-wind/equation_systems/icns/icns_ops.H b/amr-wind/equation_systems/icns/icns_ops.H index 5dc3706065..d6fbfa2fe8 100644 --- a/amr-wind/equation_systems/icns/icns_ops.H +++ b/amr-wind/equation_systems/icns/icns_ops.H @@ -80,10 +80,14 @@ struct SrcTermOp : SrcTermOpBase : SrcTermOpBase(fields_in), grad_p(fields_in.repo.get_field("gp")) {} - void operator()(const FieldState fstate) + void operator()(const FieldState fstate, const bool mesh_mapping) { const auto rhostate = field_impl::phi_state(fstate); const auto& density = m_density.state(rhostate); + Field const* mesh_fac = + mesh_mapping + ? &(this->fields.repo.get_mesh_mapping_field(FieldLoc::CELL)) + : nullptr; const int nlevels = this->fields.repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { @@ -97,13 +101,26 @@ struct SrcTermOp : SrcTermOpBase const auto& vf = src_term.array(mfi); const auto& rho = density(lev).const_array(mfi); const auto& gp = grad_p(lev).const_array(mfi); + amrex::Array4 fac = + mesh_mapping ? ((*mesh_fac)(lev).const_array(mfi)) + : amrex::Array4(); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { amrex::Real rhoinv = 1.0 / rho(i, j, k); - vf(i, j, k, 0) = -(gp(i, j, k, 0)) * rhoinv; - vf(i, j, k, 1) = -(gp(i, j, k, 1)) * rhoinv; - vf(i, j, k, 2) = -(gp(i, j, k, 2)) * rhoinv; + amrex::Real fac_x = + mesh_mapping ? (fac(i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac(i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac(i, j, k, 2)) : 1.0; + + vf(i, j, k, 0) = + -(1.0 / fac_x * gp(i, j, k, 0)) * rhoinv; + vf(i, j, k, 1) = + -(1.0 / fac_y * gp(i, j, k, 1)) * rhoinv; + vf(i, j, k, 2) = + -(1.0 / fac_z * gp(i, j, k, 2)) * rhoinv; }); for (auto& src : this->sources) { diff --git a/amr-wind/equation_systems/levelset/levelset_ops.H b/amr-wind/equation_systems/levelset/levelset_ops.H index 52ddd3d1e2..59ae73728d 100644 --- a/amr-wind/equation_systems/levelset/levelset_ops.H +++ b/amr-wind/equation_systems/levelset/levelset_ops.H @@ -51,7 +51,8 @@ struct ComputeRHSOp { explicit ComputeRHSOp(PDEFields& fields_in) : fields(fields_in) {} - void predictor_rhs(const DiffusionType, const amrex::Real dt) + void predictor_rhs( + const DiffusionType, const amrex::Real dt, bool /*mesh_mapping*/) { // Field states for diffusion and advection terms. In Godunov scheme // these terms only have one state. @@ -85,7 +86,8 @@ struct ComputeRHSOp } } - void corrector_rhs(const DiffusionType, const amrex::Real dt) + void corrector_rhs( + const DiffusionType, const amrex::Real dt, bool /*mesh_mapping*/) { const int nlevels = fields.repo.num_active_levels(); auto& field = fields.field; diff --git a/amr-wind/equation_systems/sdr/sdr_ops.H b/amr-wind/equation_systems/sdr/sdr_ops.H index 1ec469fa3f..ceb374cbe5 100644 --- a/amr-wind/equation_systems/sdr/sdr_ops.H +++ b/amr-wind/equation_systems/sdr/sdr_ops.H @@ -78,8 +78,10 @@ struct DiffusionOp : public DiffSolverIface std::is_same::value, "Invalid linear operator for scalar diffusion operator"); - DiffusionOp(PDEFields& fields, const bool has_overset) - : DiffSolverIface(fields, has_overset) + DiffusionOp( + PDEFields& fields, const bool has_overset, const bool mesh_mapping) + : DiffSolverIface( + fields, has_overset, mesh_mapping) , m_lhs_src_term( fields.repo.get_field(SDR::var_name() + "_lhs_src_term")) { diff --git a/amr-wind/equation_systems/tke/tke_ops.H b/amr-wind/equation_systems/tke/tke_ops.H index 638b147092..5f19ffe27b 100644 --- a/amr-wind/equation_systems/tke/tke_ops.H +++ b/amr-wind/equation_systems/tke/tke_ops.H @@ -113,8 +113,10 @@ struct DiffusionOp : public DiffSolverIface std::is_same::value, "Invalid linear operator for scalar diffusion operator"); - DiffusionOp(PDEFields& fields, const bool has_overset) - : DiffSolverIface(fields, has_overset) + DiffusionOp( + PDEFields& fields, const bool has_overset, const bool mesh_mapping) + : DiffSolverIface( + fields, has_overset, mesh_mapping) , m_lhs_src_term( fields.repo.get_field(TKE::var_name() + "_lhs_src_term")) { diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H index f6e746f363..4bddeb0b0a 100644 --- a/amr-wind/equation_systems/vof/vof_advection.H +++ b/amr-wind/equation_systems/vof/vof_advection.H @@ -14,7 +14,7 @@ namespace pde { template <> struct AdvectionOp { - AdvectionOp(PDEFields& fields_in, bool, bool) + AdvectionOp(PDEFields& fields_in, bool, bool, bool) : fields(fields_in) , u_mac(fields_in.repo.get_field("u_mac")) , v_mac(fields_in.repo.get_field("v_mac")) diff --git a/amr-wind/equation_systems/vof/vof_ops.H b/amr-wind/equation_systems/vof/vof_ops.H index cb44a8ead8..af889383fd 100644 --- a/amr-wind/equation_systems/vof/vof_ops.H +++ b/amr-wind/equation_systems/vof/vof_ops.H @@ -45,9 +45,9 @@ struct ComputeRHSOp { explicit ComputeRHSOp(PDEFields& fields_in) : fields(fields_in) {} - void predictor_rhs(const DiffusionType, const amrex::Real) {} + void predictor_rhs(const DiffusionType, const amrex::Real, bool) {} - void corrector_rhs(const DiffusionType, const amrex::Real) {} + void corrector_rhs(const DiffusionType, const amrex::Real, bool) {} // data members PDEFields& fields; diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 4cb00874ad..0cac9f24e7 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -4,10 +4,13 @@ #include "amr-wind/utilities/tagging/RefinementCriteria.H" #include "amr-wind/equation_systems/PDEBase.H" #include "amr-wind/turbulence/TurbulenceModel.H" +#include "amr-wind/equation_systems/SchemeTraits.H" #include "amr-wind/utilities/IOManager.H" #include "amr-wind/utilities/PostProcessing.H" #include "amr-wind/overset/OversetManager.H" +#include "AMReX_ParmParse.H" + using namespace amrex; incflo::incflo() @@ -178,6 +181,20 @@ bool incflo::regrid_and_update() printGridSummary(amrex::OutStream(), 0, finest_level); } + // update mesh map + { + if (m_sim.has_mesh_mapping()) { + // TODO: Is this the only change required in presence of regrid + // ? + amrex::Print() << "Creating mesh mapping after regrid ... "; + + for (int lev = 0; lev <= finest_level; lev++) { + m_sim.mesh_mapping()->create_map(lev, Geom(lev)); + } + amrex::Print() << "done" << std::endl; + } + } + if (m_sim.has_overset()) { m_sim.overset_manager()->post_regrid_actions(); } else { @@ -250,6 +267,7 @@ void incflo::Evolve() amrex::Real time0 = amrex::ParallelDescriptor::second(); regrid_and_update(); + pre_advance_stage1(); pre_advance_stage2(); @@ -282,7 +300,6 @@ void incflo::Evolve() if (m_time.write_last_plot_file()) { m_sim.io_manager().write_plot_file(); } - if (m_time.write_last_checkpoint()) { m_sim.io_manager().write_checkpoint_file(); } @@ -312,6 +329,11 @@ void incflo::MakeNewLevelFromScratch( m_repo.make_new_level_from_scratch(lev, time, new_grids, new_dmap); + // initialize the mesh map before initializing physics + if (m_sim.has_mesh_mapping()) { + m_sim.mesh_mapping()->create_map(lev, Geom(lev)); + } + for (auto& pp : m_sim.physics()) { pp->initialize_fields(lev, Geom(lev)); } @@ -319,6 +341,9 @@ void incflo::MakeNewLevelFromScratch( void incflo::init_physics_and_pde() { + // Check for mesh mapping + m_sim.activate_mesh_map(); + { // Query and activate overset before initializing PDEs and physics amrex::ParmParse pp("incflo"); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 0e6d798caa..3caff0e716 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -192,6 +192,8 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* + // TODO: This sub-section has not been adjusted for mesh mapping - adjust in + // corrector too m_sim.turbulence_model().update_turbulent_viscosity( amr_wind::FieldState::Old); icns().compute_mueff(amr_wind::FieldState::Old); @@ -202,6 +204,8 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* // Define the forcing terms to use in the Godunov prediction // ************************************************************************************* + // TODO: Godunov has not been adjusted for mesh mapping - adjust in + // corrector too if (m_use_godunov) { icns().compute_source_term(amr_wind::FieldState::Old); for (auto& seqn : scalar_eqns()) { @@ -209,6 +213,8 @@ void incflo::ApplyPredictor(bool incremental_projection) } } + // TODO: This sub-section has not been adjusted for mesh mapping - adjust in + // corrector too if (need_divtau()) { // ************************************************************************************* // Compute explicit viscous term @@ -265,6 +271,8 @@ void incflo::ApplyPredictor(bool incremental_projection) // if (!m_use_godunov) Compute the explicit advective terms // R_u^n , R_s^n and R_t^n // ************************************************************************************* + // TODO: Advection computation for scalar equations have not been adjusted + // for mesh mapping for (auto& seqn : scalar_eqns()) { seqn->compute_advection_term(amr_wind::FieldState::Old); } @@ -275,9 +283,12 @@ void incflo::ApplyPredictor(bool incremental_projection) if (m_constant_density) { amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1); } - // Perform scalar update one at a time. This is to allow an updated - // density at `n+1/2` to be computed before other scalars use it when - // computing their source terms. + + // TODO: This sub-section has not been adjusted for mesh mapping - adjust in + // corrector too. + // Perform scalar update one at a time. This is to allow an + // updated density at `n+1/2` to be computed before other scalars use it + // when computing their source terms. for (auto& eqn : scalar_eqns()) { // Compute (recompute for Godunov) the scalar forcing terms eqn->compute_source_term(amr_wind::FieldState::NPH); @@ -316,7 +327,7 @@ void incflo::ApplyPredictor(bool incremental_projection) icns().compute_source_term(amr_wind::FieldState::New); // ************************************************************************************* - // Update the velocity + // Evaluate right hand side and store in velocity // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); @@ -489,9 +500,10 @@ void incflo::ApplyCorrector() amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1); } - // Perform scalar update one at a time. This is to allow an updated density - // at `n+1/2` to be computed before other scalars use it when computing - // their source terms. + // TODO: This sub-section has not been adjusted for mesh mapping - adjust in + // corrector too Perform scalar update one at a time. This is to allow an + // updated density at `n+1/2` to be computed before other scalars use it + // when computing their source terms. for (auto& eqn : scalar_eqns()) { // Compute (recompute for Godunov) the scalar forcing terms // Note this is (rho * scalar) and not just scalar @@ -527,7 +539,7 @@ void incflo::ApplyCorrector() icns().compute_source_term(amr_wind::FieldState::New); // ************************************************************************************* - // Update velocity + // Evaluate right hand side and store in velocity // ************************************************************************************* icns().compute_corrector_rhs(m_diff_type); diff --git a/amr-wind/incflo_compute_dt.cpp b/amr-wind/incflo_compute_dt.cpp index 394b1e8e71..1b9a5ca174 100644 --- a/amr-wind/incflo_compute_dt.cpp +++ b/amr-wind/incflo_compute_dt.cpp @@ -38,8 +38,13 @@ void incflo::ComputeDt(bool explicit_diffusion) Real conv_cfl = 0.0; Real diff_cfl = 0.0; Real force_cfl = 0.0; + const bool mesh_mapping = m_sim.has_mesh_mapping(); const auto& den = density(); + amr_wind::Field const* mesh_fac = + mesh_mapping + ? &(m_repo.get_mesh_mapping_field(amr_wind::FieldLoc::CELL)) + : nullptr; for (int lev = 0; lev <= finest_level; ++lev) { auto const dxinv = geom[lev].InvCellSizeArray(); @@ -48,57 +53,83 @@ void incflo::ComputeDt(bool explicit_diffusion) MultiFab const& mu = icns().fields().mueff(lev); MultiFab const& rho = den(lev); + auto const& vel_arr = vel.const_arrays(); + MultiArray4 fac_arr = + mesh_mapping ? ((*mesh_fac)(lev).const_arrays()) + : MultiArray4(); + Real conv_lev = 0.0; Real diff_lev = 0.0; Real force_lev = 0.0; - conv_lev = amrex::ReduceMax( - vel, 0, + conv_lev += amrex::ParReduce( + TypeList{}, TypeList{}, vel, IntVect(0), [=] AMREX_GPU_HOST_DEVICE( - Box const& b, Array4 const& v) -> Real { - Real mx = -1.0; - amrex::Loop(b, [=, &mx](int i, int j, int k) noexcept { - mx = amrex::max( - amrex::Math::abs(v(i, j, k, 0)) * dxinv[0], - amrex::Math::abs(v(i, j, k, 1)) * dxinv[1], - amrex::Math::abs(v(i, j, k, 2)) * dxinv[2], mx); - }); - return mx; + int box_no, int i, int j, int k) -> GpuTuple { + auto const& v_bx = vel_arr[box_no]; + + amrex::Real fac_x = + mesh_mapping ? (fac_arr[box_no](i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac_arr[box_no](i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac_arr[box_no](i, j, k, 2)) : 1.0; + + return amrex::max( + amrex::Math::abs(v_bx(i, j, k, 0)) * dxinv[0] / fac_x, + amrex::Math::abs(v_bx(i, j, k, 1)) * dxinv[1] / fac_y, + amrex::Math::abs(v_bx(i, j, k, 2)) * dxinv[2] / fac_z, + -1.0); }); if (explicit_diffusion) { + auto const& mu_arr = mu.const_arrays(); + auto const& rho_arr = rho.const_arrays(); + diff_lev += amrex::ParReduce( + TypeList{}, TypeList{}, rho, IntVect(0), + [=] AMREX_GPU_HOST_DEVICE( + int box_no, int i, int j, int k) -> GpuTuple { + auto const& mu_bx = mu_arr[box_no]; + auto const& rho_bx = rho_arr[box_no]; - const Real dxinv2 = - 2.0 * (dxinv[0] * dxinv[0] + dxinv[1] * dxinv[1] + - dxinv[2] * dxinv[2]); + amrex::Real fac_x = + mesh_mapping ? (fac_arr[box_no](i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac_arr[box_no](i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac_arr[box_no](i, j, k, 2)) : 1.0; - diff_lev = amrex::ReduceMax( - rho, mu, 0, - [=] AMREX_GPU_HOST_DEVICE( - Box const& b, Array4 const& rho_arr, - Array4 const& mu_arr) -> Real { - Real mx = -1.0; - amrex::Loop(b, [=, &mx](int i, int j, int k) noexcept { - mx = amrex::max( - mu_arr(i, j, k) * dxinv2 / rho_arr(i, j, k), mx); - }); - return mx; + const Real dxinv2 = + 2.0 * (dxinv[0] / fac_x * dxinv[0] / fac_x + + dxinv[1] / fac_y * dxinv[1] / fac_y + + dxinv[2] / fac_z * dxinv[2] / fac_z); + + return amrex::max( + mu_bx(i, j, k) * dxinv2 / rho_bx(i, j, k), -1.0); }); } if (m_time.use_force_cfl()) { - force_lev = amrex::ReduceMax( - vel_force, 0, + auto const& vf_arr = vel_force.const_arrays(); + force_lev += amrex::ParReduce( + TypeList{}, TypeList{}, vel_force, + IntVect(0), [=] AMREX_GPU_HOST_DEVICE( - Box const& b, Array4 const& vf) -> Real { - Real mx = -1.0; - amrex::Loop(b, [=, &mx](int i, int j, int k) noexcept { - mx = amrex::max( - amrex::Math::abs(vf(i, j, k, 0)) * dxinv[0], - amrex::Math::abs(vf(i, j, k, 1)) * dxinv[1], - amrex::Math::abs(vf(i, j, k, 2)) * dxinv[2], mx); - }); - return mx; + int box_no, int i, int j, int k) -> GpuTuple { + auto const& vf_bx = vf_arr[box_no]; + + amrex::Real fac_x = + mesh_mapping ? (fac_arr[box_no](i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac_arr[box_no](i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac_arr[box_no](i, j, k, 2)) : 1.0; + + return amrex::max( + amrex::Math::abs(vf_bx(i, j, k, 0)) * dxinv[0] / fac_x, + amrex::Math::abs(vf_bx(i, j, k, 1)) * dxinv[1] / fac_y, + amrex::Math::abs(vf_bx(i, j, k, 2)) * dxinv[2] / fac_z, + -1.0); }); } diff --git a/amr-wind/mesh_mapping_models/CMakeLists.txt b/amr-wind/mesh_mapping_models/CMakeLists.txt new file mode 100644 index 0000000000..6568e14cd9 --- /dev/null +++ b/amr-wind/mesh_mapping_models/CMakeLists.txt @@ -0,0 +1,6 @@ +target_sources(${amr_wind_lib_name} + PRIVATE + + ChannelFlowMap.cpp + ConstantMap.cpp + ) diff --git a/amr-wind/mesh_mapping_models/ChannelFlowMap.H b/amr-wind/mesh_mapping_models/ChannelFlowMap.H new file mode 100644 index 0000000000..6fc6f9ae90 --- /dev/null +++ b/amr-wind/mesh_mapping_models/ChannelFlowMap.H @@ -0,0 +1,43 @@ +#ifndef CHANNELFLOWMAP_H +#define CHANNELFLOWMAP_H + +#include "amr-wind/core/MeshMap.H" + +namespace amr_wind { +namespace channel_map { + +/** Channel flow scaling mesh map + * \ingroup mesh_map + */ +class ChannelFlowMap : public MeshMap::Register +{ +public: + static const std::string identifier() { return "ChannelFlowMap"; } + + explicit ChannelFlowMap(); + + virtual ~ChannelFlowMap() = default; + + //! Construct the mesh scaling field + void create_map(int, const amrex::Geometry&) override; + + //! Construct mesh scaling field on cell centers and nodes + void create_cell_node_map(int, const amrex::Geometry&); + + //! Construct mesh scaling field on cell faces + void create_face_map(int, const amrex::Geometry&); + + //! Construct the non-uniform mesh field + void create_non_uniform_mesh(int, const amrex::Geometry&); + +private: + //! User input parameters + amrex::Vector m_beta{{0.0, 3.0, 0.0}}; + + amrex::Real m_eps{1e-11}; +}; + +} // namespace channel_map +} // namespace amr_wind + +#endif /* CHANNELFLOWMAP_H */ diff --git a/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp b/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp new file mode 100644 index 0000000000..35a7b6b8ad --- /dev/null +++ b/amr-wind/mesh_mapping_models/ChannelFlowMap.cpp @@ -0,0 +1,338 @@ +#include + +#include "amr-wind/mesh_mapping_models/ChannelFlowMap.H" + +#include "AMReX_ParmParse.H" + +namespace amr_wind { +namespace channel_map { + +namespace { + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eval_fac( + const amrex::Real x, + const amrex::Real beta, + const amrex::Real prob_lo, + const amrex::Real len) +{ + return (beta == 0.0) + ? 1.0 + : (beta * + (1 - + std::pow( + std::tanh(beta * (1 - 2 * (x - prob_lo) / len)), 2)) / + std::tanh(beta)); +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eval_coord( + const amrex::Real x, + const amrex::Real beta, + const amrex::Real prob_lo, + const amrex::Real len) +{ + return (beta == 0.0) + ? x + : (prob_lo + + len / 2 * + (1 - std::tanh(beta * (1 - 2 * (x - prob_lo) / len)) / + std::tanh(beta))); +} + +} // namespace + +ChannelFlowMap::ChannelFlowMap() +{ + amrex::ParmParse pp("ChannelFlowMap"); + pp.queryarr("beta", m_beta, 0, AMREX_SPACEDIM); +} + +/** Construct the mesh mapping field + */ +void ChannelFlowMap::create_map(int lev, const amrex::Geometry& geom) +{ + create_cell_node_map(lev, geom); + create_face_map(lev, geom); + create_non_uniform_mesh(lev, geom); +} + +/** Construct the mesh mapping field on cell centers and nodes + */ +void ChannelFlowMap::create_cell_node_map(int lev, const amrex::Geometry& geom) +{ + amrex::GpuArray beta{ + {m_beta[0], m_beta[1], m_beta[2]}}; + const auto eps = m_eps; + + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + + amrex::GpuArray len{ + {prob_hi[0] - prob_lo[0], prob_hi[1] - prob_lo[1], + prob_hi[2] - prob_lo[2]}}; + + for (amrex::MFIter mfi((*m_mesh_scale_fac_cc)(lev)); mfi.isValid(); ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_cc = + (*m_mesh_scale_fac_cc)(lev).array(mfi); + amrex::Array4 const& scale_detJ_cc = + (*m_mesh_scale_detJ_cc)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + amrex::Real fac_x = eval_fac(x, beta[0], prob_lo[0], len[0]); + amrex::Real fac_y = eval_fac(y, beta[1], prob_lo[1], len[1]); + amrex::Real fac_z = eval_fac(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x > prob_lo[0]) && (x < prob_hi[0]) && (y > prob_lo[1]) && + (y < prob_hi[1]) && (z > prob_lo[2]) && (z < prob_hi[2])); + + scale_fac_cc(i, j, k, 0) = in_domain ? fac_x : 1.0; + scale_fac_cc(i, j, k, 1) = in_domain ? fac_y : 1.0; + scale_fac_cc(i, j, k, 2) = in_domain ? fac_z : 1.0; + + scale_detJ_cc(i, j, k) = scale_fac_cc(i, j, k, 0) * + scale_fac_cc(i, j, k, 1) * + scale_fac_cc(i, j, k, 2); + }); + + const auto& nbx = mfi.grownnodaltilebox(); + amrex::Array4 const& scale_fac_nd = + (*m_mesh_scale_fac_nd)(lev).array(mfi); + amrex::Array4 const& scale_detJ_nd = + (*m_mesh_scale_detJ_nd)(lev).array(mfi); + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + i * dx[0]; + amrex::Real y = prob_lo[1] + j * dx[1]; + amrex::Real z = prob_lo[2] + k * dx[2]; + + amrex::Real fac_x = eval_fac(x, beta[0], prob_lo[0], len[0]); + amrex::Real fac_y = eval_fac(y, beta[1], prob_lo[1], len[1]); + amrex::Real fac_z = eval_fac(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x >= prob_lo[0] - eps) && (x <= prob_hi[0] + eps) && + (y >= prob_lo[1] - eps) && (y <= prob_hi[1] + eps) && + (z >= prob_lo[2] - eps) && (z <= prob_hi[2] + eps)); + + scale_fac_nd(i, j, k, 0) = in_domain ? fac_x : 1.0; + scale_fac_nd(i, j, k, 1) = in_domain ? fac_y : 1.0; + scale_fac_nd(i, j, k, 2) = in_domain ? fac_z : 1.0; + + scale_detJ_nd(i, j, k) = scale_fac_nd(i, j, k, 0) * + scale_fac_nd(i, j, k, 1) * + scale_fac_nd(i, j, k, 2); + }); + } + + // TODO: Call fill patch operators ? +} + +/** Construct the mesh mapping field on cell faces + */ +void ChannelFlowMap::create_face_map(int lev, const amrex::Geometry& geom) +{ + amrex::GpuArray beta{ + {m_beta[0], m_beta[1], m_beta[2]}}; + const auto eps = m_eps; + + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + + amrex::GpuArray len{ + {prob_hi[0] - prob_lo[0], prob_hi[1] - prob_lo[1], + prob_hi[2] - prob_lo[2]}}; + + for (amrex::MFIter mfi((*m_mesh_scale_fac_xf)(lev)); mfi.isValid(); ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_xf = + (*m_mesh_scale_fac_xf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_xf = + (*m_mesh_scale_detJ_xf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + i * dx[0]; + amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + amrex::Real fac_x = eval_fac(x, beta[0], prob_lo[0], len[0]); + amrex::Real fac_y = eval_fac(y, beta[1], prob_lo[1], len[1]); + amrex::Real fac_z = eval_fac(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x >= prob_lo[0] - eps) && (x <= prob_hi[0] + eps) && + (y > prob_lo[1]) && (y < prob_hi[1]) && (z > prob_lo[2]) && + (z < prob_hi[2])); + + scale_fac_xf(i, j, k, 0) = in_domain ? fac_x : 1.0; + scale_fac_xf(i, j, k, 1) = in_domain ? fac_y : 1.0; + scale_fac_xf(i, j, k, 2) = in_domain ? fac_z : 1.0; + + scale_detJ_xf(i, j, k) = scale_fac_xf(i, j, k, 0) * + scale_fac_xf(i, j, k, 1) * + scale_fac_xf(i, j, k, 2); + }); + } + + for (amrex::MFIter mfi((*m_mesh_scale_fac_yf)(lev)); mfi.isValid(); ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_yf = + (*m_mesh_scale_fac_yf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_yf = + (*m_mesh_scale_detJ_yf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + amrex::Real y = prob_lo[1] + j * dx[1]; + amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + amrex::Real fac_x = eval_fac(x, beta[0], prob_lo[0], len[0]); + amrex::Real fac_y = eval_fac(y, beta[1], prob_lo[1], len[1]); + amrex::Real fac_z = eval_fac(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x > prob_lo[0]) && (x < prob_hi[0]) && + (y >= prob_lo[1] - eps) && (y <= prob_hi[1] + eps) && + (z > prob_lo[2]) && (z < prob_hi[2])); + + scale_fac_yf(i, j, k, 0) = in_domain ? fac_x : 1.0; + scale_fac_yf(i, j, k, 1) = in_domain ? fac_y : 1.0; + scale_fac_yf(i, j, k, 2) = in_domain ? fac_z : 1.0; + + scale_detJ_yf(i, j, k) = scale_fac_yf(i, j, k, 0) * + scale_fac_yf(i, j, k, 1) * + scale_fac_yf(i, j, k, 2); + }); + } + + for (amrex::MFIter mfi((*m_mesh_scale_fac_zf)(lev)); mfi.isValid(); ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_zf = + (*m_mesh_scale_fac_zf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_zf = + (*m_mesh_scale_detJ_zf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + amrex::Real z = prob_lo[2] + k * dx[2]; + + amrex::Real fac_x = eval_fac(x, beta[0], prob_lo[0], len[0]); + amrex::Real fac_y = eval_fac(y, beta[1], prob_lo[1], len[1]); + amrex::Real fac_z = eval_fac(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x > prob_lo[0]) && (x < prob_hi[0]) && (y > prob_lo[1]) && + (y < prob_hi[1]) && (z >= prob_lo[2] - eps) && + (z <= prob_hi[2] + eps)); + + scale_fac_zf(i, j, k, 0) = in_domain ? fac_x : 1.0; + scale_fac_zf(i, j, k, 1) = in_domain ? fac_y : 1.0; + scale_fac_zf(i, j, k, 2) = in_domain ? fac_z : 1.0; + + scale_detJ_zf(i, j, k) = scale_fac_zf(i, j, k, 0) * + scale_fac_zf(i, j, k, 1) * + scale_fac_zf(i, j, k, 2); + }); + } + + // TODO: Call fill patch operators ? +} + +/** Construct the non-uniform mesh field + */ +void ChannelFlowMap::create_non_uniform_mesh( + int lev, const amrex::Geometry& geom) +{ + amrex::Vector probhi_physical{{0.0, 0.0, 0.0}}; + { + amrex::ParmParse pp("geometry"); + if (pp.contains("prob_hi_physical")) { + pp.getarr("prob_hi_physical", probhi_physical); + } else { + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + probhi_physical[d] = geom.ProbHiArray()[d]; + } + } + } + + amrex::GpuArray beta{ + {m_beta[0], m_beta[1], m_beta[2]}}; + const auto eps = m_eps; + + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + + amrex::GpuArray len{ + {probhi_physical[0] - prob_lo[0], probhi_physical[1] - prob_lo[1], + probhi_physical[2] - prob_lo[2]}}; + + for (amrex::MFIter mfi((*m_non_uniform_coord_cc)(lev)); mfi.isValid(); + ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& nu_coord_cc = + (*m_non_uniform_coord_cc)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + amrex::Real x_non_uni = + eval_coord(x, beta[0], prob_lo[0], len[0]); + amrex::Real y_non_uni = + eval_coord(y, beta[1], prob_lo[1], len[1]); + amrex::Real z_non_uni = + eval_coord(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x > prob_lo[0]) && (x < prob_hi[0]) && (y > prob_lo[1]) && + (y < prob_hi[1]) && (z > prob_lo[2]) && (z < prob_hi[2])); + + nu_coord_cc(i, j, k, 0) = in_domain ? x_non_uni : x; + nu_coord_cc(i, j, k, 1) = in_domain ? y_non_uni : y; + nu_coord_cc(i, j, k, 2) = in_domain ? z_non_uni : z; + }); + + const auto& nbx = mfi.grownnodaltilebox(); + amrex::Array4 const& nu_coord_nd = + (*m_non_uniform_coord_nd)(lev).array(mfi); + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real x = prob_lo[0] + i * dx[0]; + amrex::Real y = prob_lo[1] + j * dx[1]; + amrex::Real z = prob_lo[2] + k * dx[2]; + + amrex::Real x_non_uni = + eval_coord(x, beta[0], prob_lo[0], len[0]); + amrex::Real y_non_uni = + eval_coord(y, beta[1], prob_lo[1], len[1]); + amrex::Real z_non_uni = + eval_coord(z, beta[2], prob_lo[2], len[2]); + + bool in_domain = + ((x >= prob_lo[0] - eps) && (x <= prob_hi[0] + eps) && + (y >= prob_lo[1] - eps) && (y <= prob_hi[1] + eps) && + (z >= prob_lo[2] - eps) && (z <= prob_hi[2] + eps)); + + nu_coord_nd(i, j, k, 0) = in_domain ? x_non_uni : x; + nu_coord_nd(i, j, k, 1) = in_domain ? y_non_uni : y; + nu_coord_nd(i, j, k, 2) = in_domain ? z_non_uni : z; + }); + } +} + +} // namespace channel_map +} // namespace amr_wind diff --git a/amr-wind/mesh_mapping_models/ConstantMap.H b/amr-wind/mesh_mapping_models/ConstantMap.H new file mode 100644 index 0000000000..42f83dd83c --- /dev/null +++ b/amr-wind/mesh_mapping_models/ConstantMap.H @@ -0,0 +1,41 @@ +#ifndef CONSTANTMAP_H +#define CONSTANTMAP_H + +#include "amr-wind/core/MeshMap.H" + +namespace amr_wind { +namespace const_map { + +/** Coonstant scaling mesh map + * \ingroup mesh_map + */ +class ConstantMap : public MeshMap::Register +{ +public: + static const std::string identifier() { return "ConstantMap"; } + + explicit ConstantMap(); + + virtual ~ConstantMap() = default; + + //! Construct the mesh scaling field + void create_map(int, const amrex::Geometry&) override; + + //! Construct mesh scaling field on cell centers and nodes + void create_cell_node_map(int); + + //! Construct mesh scaling field on cell faces + void create_face_map(int); + + //! Construct the non-uniform mesh field + void create_non_uniform_mesh(int, const amrex::Geometry&); + +private: + //! Factor to scale the mesh by + amrex::Vector m_fac{{1.0, 1.0, 1.0}}; +}; + +} // namespace const_map +} // namespace amr_wind + +#endif /* CONSTANTMAP_H */ diff --git a/amr-wind/mesh_mapping_models/ConstantMap.cpp b/amr-wind/mesh_mapping_models/ConstantMap.cpp new file mode 100644 index 0000000000..f047fa8e92 --- /dev/null +++ b/amr-wind/mesh_mapping_models/ConstantMap.cpp @@ -0,0 +1,169 @@ +#include "amr-wind/mesh_mapping_models/ConstantMap.H" + +#include "AMReX_ParmParse.H" + +namespace amr_wind { +namespace const_map { + +ConstantMap::ConstantMap() +{ + amrex::ParmParse pp("ConstantMap"); + pp.queryarr("scaling_factor", m_fac, 0, AMREX_SPACEDIM); +} + +/** Construct the mesh mapping field + */ +void ConstantMap::create_map(int lev, const amrex::Geometry& geom) +{ + create_cell_node_map(lev); + create_face_map(lev); + create_non_uniform_mesh(lev, geom); +} + +/** Construct the mesh mapping field on cell centers and nodes + */ +void ConstantMap::create_cell_node_map(int lev) +{ + amrex::Real fac_x = m_fac[0]; + amrex::Real fac_y = m_fac[1]; + amrex::Real fac_z = m_fac[2]; + + for (amrex::MFIter mfi((*m_mesh_scale_fac_cc)(lev)); mfi.isValid(); ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_cc = + (*m_mesh_scale_fac_cc)(lev).array(mfi); + amrex::Array4 const& scale_detJ_cc = + (*m_mesh_scale_detJ_cc)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + scale_fac_cc(i, j, k, 0) = fac_x; + scale_fac_cc(i, j, k, 1) = fac_y; + scale_fac_cc(i, j, k, 2) = fac_z; + scale_detJ_cc(i, j, k) = scale_fac_cc(i, j, k, 0) * + scale_fac_cc(i, j, k, 1) * + scale_fac_cc(i, j, k, 2); + }); + + const auto& nbx = mfi.grownnodaltilebox(); + amrex::Array4 const& scale_fac_nd = + (*m_mesh_scale_fac_nd)(lev).array(mfi); + amrex::Array4 const& scale_detJ_nd = + (*m_mesh_scale_detJ_nd)(lev).array(mfi); + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + scale_fac_nd(i, j, k, 0) = fac_x; + scale_fac_nd(i, j, k, 1) = fac_y; + scale_fac_nd(i, j, k, 2) = fac_z; + scale_detJ_nd(i, j, k) = scale_fac_nd(i, j, k, 0) * + scale_fac_nd(i, j, k, 1) * + scale_fac_nd(i, j, k, 2); + }); + } +} + +/** Construct the mesh mapping field on cell faces + */ +void ConstantMap::create_face_map(int lev) +{ + amrex::Real fac_x = m_fac[0]; + amrex::Real fac_y = m_fac[1]; + amrex::Real fac_z = m_fac[2]; + + for (amrex::MFIter mfi((*m_mesh_scale_fac_xf)(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_xf = + (*m_mesh_scale_fac_xf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_xf = + (*m_mesh_scale_detJ_xf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + scale_fac_xf(i, j, k, 0) = fac_x; + scale_fac_xf(i, j, k, 1) = fac_y; + scale_fac_xf(i, j, k, 2) = fac_z; + scale_detJ_xf(i, j, k) = scale_fac_xf(i, j, k, 0) * + scale_fac_xf(i, j, k, 1) * + scale_fac_xf(i, j, k, 2); + }); + } + + for (amrex::MFIter mfi((*m_mesh_scale_fac_yf)(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_yf = + (*m_mesh_scale_fac_yf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_yf = + (*m_mesh_scale_detJ_yf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + scale_fac_yf(i, j, k, 0) = fac_x; + scale_fac_yf(i, j, k, 1) = fac_y; + scale_fac_yf(i, j, k, 2) = fac_z; + scale_detJ_yf(i, j, k) = scale_fac_yf(i, j, k, 0) * + scale_fac_yf(i, j, k, 1) * + scale_fac_yf(i, j, k, 2); + }); + } + + for (amrex::MFIter mfi((*m_mesh_scale_fac_zf)(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_zf = + (*m_mesh_scale_fac_zf)(lev).array(mfi); + amrex::Array4 const& scale_detJ_zf = + (*m_mesh_scale_detJ_zf)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + scale_fac_zf(i, j, k, 0) = fac_x; + scale_fac_zf(i, j, k, 1) = fac_y; + scale_fac_zf(i, j, k, 2) = fac_z; + scale_detJ_zf(i, j, k) = scale_fac_zf(i, j, k, 0) * + scale_fac_zf(i, j, k, 1) * + scale_fac_zf(i, j, k, 2); + }); + } +} + +/** Construct the non-uniform mesh field + */ +void ConstantMap::create_non_uniform_mesh(int lev, const amrex::Geometry& geom) +{ + + const auto& problo = geom.ProbLoArray(); + const auto& dx = geom.CellSizeArray(); + + for (amrex::MFIter mfi((*m_non_uniform_coord_cc)(lev)); mfi.isValid(); + ++mfi) { + + const auto& bx = mfi.growntilebox(); + amrex::Array4 const& scale_fac_cc = + (*m_mesh_scale_fac_cc)(lev).array(mfi); + amrex::Array4 const& nu_coord_cc = + (*m_non_uniform_coord_cc)(lev).array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + nu_coord_cc(i, j, k, 0) = + problo[0] + (i + 0.5) * dx[0] * scale_fac_cc(i, j, k, 0); + nu_coord_cc(i, j, k, 1) = + problo[1] + (j + 0.5) * dx[1] * scale_fac_cc(i, j, k, 1); + nu_coord_cc(i, j, k, 2) = + problo[2] + (k + 0.5) * dx[2] * scale_fac_cc(i, j, k, 2); + }); + + const auto& nbx = mfi.grownnodaltilebox(); + amrex::Array4 const& scale_fac_nd = + (*m_mesh_scale_fac_nd)(lev).array(mfi); + amrex::Array4 const& nu_coord_nd = + (*m_non_uniform_coord_nd)(lev).array(mfi); + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + nu_coord_nd(i, j, k, 0) = + problo[0] + i * dx[0] * scale_fac_nd(i, j, k, 0); + nu_coord_nd(i, j, k, 1) = + problo[1] + j * dx[1] * scale_fac_nd(i, j, k, 1); + nu_coord_nd(i, j, k, 2) = + problo[2] + k * dx[2] * scale_fac_nd(i, j, k, 2); + }); + } +} + +} // namespace const_map +} // namespace amr_wind diff --git a/amr-wind/physics/ChannelFlow.H b/amr-wind/physics/ChannelFlow.H index ed7d03f123..931b9b06c3 100644 --- a/amr-wind/physics/ChannelFlow.H +++ b/amr-wind/physics/ChannelFlow.H @@ -4,6 +4,7 @@ #include "amr-wind/core/Physics.H" #include "amr-wind/core/Field.H" #include "amr-wind/CFDSim.H" +#include namespace amr_wind { namespace channel_flow { @@ -30,6 +31,11 @@ public: const IndexSelector& idxOp, const int n_idx); + template + amrex::Real compute_error(const IndexSelector& idxOp); + + void output_error(); + void post_init_actions() override; void post_regrid_actions() override {} @@ -49,6 +55,9 @@ private: //! initial density value amrex::Real m_rho{1.0}; + //! viscosity + amrex::Real m_mu{1.0}; + //! Re_tau amrex::Real m_re_tau{1000.0}; @@ -66,6 +75,23 @@ private: //! Von-Karman constant amrex::Real m_kappa{0.41}; + + bool m_laminar{false}; + + bool m_dns{false}; + + bool m_mesh_mapping{false}; + + amrex::Real m_mean_vel; + + //! direction of mean velocity + int m_mean_vel_dir{0}; + + //! output precision + const int m_w{18}; + + //! error log file + std::string m_output_fname{"channel_flow.log"}; }; } // namespace channel_flow } // namespace amr_wind diff --git a/amr-wind/physics/ChannelFlow.cpp b/amr-wind/physics/ChannelFlow.cpp index d33189b303..efb91ad6be 100644 --- a/amr-wind/physics/ChannelFlow.cpp +++ b/amr-wind/physics/ChannelFlow.cpp @@ -10,25 +10,44 @@ namespace amr_wind { namespace channel_flow { ChannelFlow::ChannelFlow(CFDSim& sim) - : m_time(sim.time()), m_repo(sim.repo()), m_mesh(sim.mesh()) + : m_time(sim.time()) + , m_repo(sim.repo()) + , m_mesh(sim.mesh()) + , m_mesh_mapping(sim.has_mesh_mapping()) { - { amrex::ParmParse pp("ChannelFlow"); + pp.query("flow_direction", m_mean_vel_dir); pp.query("normal_direction", m_norm_dir); + pp.query("Laminar", m_laminar); + pp.query("Turbulent_DNS", m_dns); + pp.query("error_log_file", m_output_fname); + + if (m_laminar) { + pp.query("density", m_rho); + pp.query("Mean_Velocity", m_mean_vel); + } else { + pp.query("density", m_rho); + pp.query("re_tau", m_re_tau); - pp.query("density", m_rho); - pp.query("re_tau", m_re_tau); - pp.query("tke0", m_tke0); - pp.query("sdr0", m_sdr0); + if (!m_dns) { + pp.query("tke0", m_tke0); + pp.query("sdr0", m_sdr0); + } + } } { - amrex::Real mu; amrex::ParmParse pp("transport"); - pp.query("viscosity", mu); + pp.query("viscosity", m_mu); // Assumes a boundary layer height of 1.0 - m_utau = mu * m_re_tau / (m_rho * 1.0); - m_ytau = mu / (m_utau * m_rho); + m_utau = m_mu * m_re_tau / (m_rho * 1.0); + m_ytau = m_mu / (m_utau * m_rho); + } + if (amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream f; + f.open(m_output_fname.c_str()); + f << std::setw(m_w) << "time" << std::setw(m_w) << "L2_u" << std::endl; + f.close(); } } @@ -37,8 +56,10 @@ ChannelFlow::ChannelFlow(CFDSim& sim) */ void ChannelFlow::initialize_fields(int level, const amrex::Geometry& geom) { - switch (m_norm_dir) { + case 0: + initialize_fields(level, geom, XDir(), 0); + break; case 1: initialize_fields(level, geom, YDir(), 1); break; @@ -58,51 +79,200 @@ void ChannelFlow::initialize_fields( const IndexSelector& idxOp, const int n_idx) { - const amrex::Real kappa = m_kappa; const amrex::Real y_tau = m_ytau; const amrex::Real utau = m_utau; - auto& velocity = m_repo.get_field("velocity")(level); + auto& velocity_field = m_repo.get_field("velocity"); + auto& velocity = velocity_field(level); auto& density = m_repo.get_field("density")(level); - auto& tke = m_repo.get_field("tke")(level); - auto& sdr = m_repo.get_field("sdr")(level); - auto& walldist = m_repo.get_field("wall_dist")(level); - density.setVal(m_rho); - tke.setVal(m_tke0); - sdr.setVal(m_sdr0); - - for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) { - const auto& vbx = mfi.validbox(); - - const auto& dx = geom.CellSizeArray(); - const auto& problo = geom.ProbLoArray(); - auto vel = velocity.array(mfi); - auto wd = walldist.array(mfi); - - amrex::ParallelFor( - vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const int n_ind = idxOp(i, j, k); - amrex::Real h = problo[n_idx] + (n_ind + 0.5) * dx[n_idx]; - if (h > 1.0) { - h = 2.0 - h; - } - wd(i, j, k) = h; - const amrex::Real hp = h / y_tau; - vel(i, j, k, 0) = - utau * (1. / kappa * std::log1p(kappa * hp) + - 7.8 * (1.0 - std::exp(-hp / 11.0) - - (hp / 11.0) * std::exp(-hp / 3.0))); - // vel(i,j,k,0) = 22.0; - vel(i, j, k, 1) = 0.0; - vel(i, j, k, 2) = 0.0; + + if (!m_laminar && !m_dns) { + auto& tke = m_repo.get_field("tke")(level); + auto& sdr = m_repo.get_field("sdr")(level); + auto& walldist = m_repo.get_field("wall_dist")(level); + tke.setVal(m_tke0); + sdr.setVal(m_sdr0); + + for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) { + const auto& vbx = mfi.validbox(); + + const auto& dx = geom.CellSizeArray(); + const auto& problo = geom.ProbLoArray(); + auto vel = velocity.array(mfi); + auto wd = walldist.array(mfi); + + amrex::ParallelFor( + vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const int n_ind = idxOp(i, j, k); + amrex::Real h = problo[n_idx] + (n_ind + 0.5) * dx[n_idx]; + if (h > 1.0) { + h = 2.0 - h; + } + wd(i, j, k) = h; + const amrex::Real hp = h / y_tau; + vel(i, j, k, 0) = + utau * (1. / kappa * std::log1p(kappa * hp) + + 7.8 * (1.0 - std::exp(-hp / 11.0) - + (hp / 11.0) * std::exp(-hp / 3.0))); + vel(i, j, k, 1) = 0.0; + vel(i, j, k, 2) = 0.0; + }); + } + } else if (m_laminar) { + velocity.setVal(0.0); + velocity.setVal( + m_mean_vel, m_mean_vel_dir, 1, + velocity_field.num_grow()[m_mean_vel_dir]); + } +} + +template +amrex::Real ChannelFlow::compute_error(const IndexSelector& idxOp) +{ + amrex::Real error = 0.0; + const auto flow_dir = m_mean_vel_dir; + const auto norm_dir = m_norm_dir; + const auto mu = m_mu; + const auto mesh_mapping = m_mesh_mapping; + + amrex::Real dpdx = 0; + { + amrex::Vector body_force{{0.0, 0.0, 0.0}}; + amrex::ParmParse pp("BodyForce"); + pp.queryarr("magnitude", body_force, 0, AMREX_SPACEDIM); + dpdx = body_force[flow_dir]; + } + + const auto& problo = m_mesh.Geom(0).ProbLoArray(); + amrex::Real ht = 0.0; + { + amrex::Vector probhi_physical{{0.0, 0.0, 0.0}}; + amrex::ParmParse pp("geometry"); + if (pp.contains("prob_hi_physical")) { + pp.getarr("prob_hi_physical", probhi_physical); + } else { + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + probhi_physical[d] = m_mesh.Geom(0).ProbHiArray()[d]; + } + } + ht = probhi_physical[norm_dir] - problo[norm_dir]; + } + + const auto& velocity = m_repo.get_field("velocity"); + Field const* nu_coord_cc = + mesh_mapping ? &(m_repo.get_field("non_uniform_coord_cc")) : nullptr; + Field const* mesh_fac_cc = + mesh_mapping + ? &(m_repo.get_mesh_mapping_field(amr_wind::FieldLoc::CELL)) + : nullptr; + + const int nlevels = m_repo.num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + + amrex::iMultiFab level_mask; + if (lev < nlevels - 1) { + level_mask = makeFineMask( + m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), + m_mesh.boxArray(lev + 1), amrex::IntVect(2), 1, 0); + } else { + level_mask.define( + m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& dx = m_mesh.Geom(lev).CellSizeArray(); + const auto& prob_lo = m_mesh.Geom(lev).ProbLoArray(); + + const auto& vel = velocity(lev); + auto const& vel_arr = vel.const_arrays(); + auto const& mask_arr = level_mask.const_arrays(); + amrex::MultiArray4 fac_arr = + mesh_mapping ? ((*mesh_fac_cc)(lev).const_arrays()) + : amrex::MultiArray4(); + amrex::MultiArray4 nu_cc = + mesh_mapping ? ((*nu_coord_cc)(lev).const_arrays()) + : amrex::MultiArray4(); + + error += amrex::ParReduce( + amrex::TypeList{}, + amrex::TypeList{}, vel, amrex::IntVect(0), + [=] AMREX_GPU_HOST_DEVICE(int box_no, int i, int j, int k) + -> amrex::GpuTuple { + auto const& vel_bx = vel_arr[box_no]; + auto const& mask_bx = mask_arr[box_no]; + + amrex::Real y = mesh_mapping + ? (nu_cc[box_no](i, j, k, norm_dir)) + : (prob_lo[norm_dir] + + (idxOp(i, j, k) + 0.5) * dx[norm_dir]); + amrex::Real fac_x = + mesh_mapping ? (fac_arr[box_no](i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac_arr[box_no](i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac_arr[box_no](i, j, k, 2)) : 1.0; + + const amrex::Real u = vel_bx(i, j, k, flow_dir); + const amrex::Real u_exact = + 1 / (2 * mu) * -dpdx * (y * y - y * ht); + + const amrex::Real cell_vol = + dx[0] * fac_x * dx[1] * fac_y * dx[2] * fac_z; + + return cell_vol * mask_bx(i, j, k) * (u - u_exact) * + (u - u_exact); }); } + + amrex::ParallelDescriptor::ReduceRealSum(error); + + const amrex::Real total_vol = m_mesh.Geom(0).ProbDomain().volume(); + return std::sqrt(error / total_vol); } -void ChannelFlow::post_init_actions() {} +void ChannelFlow::output_error() +{ + // analytical solution exists only for flow in streamwise direction + amrex::Real u_err = 0.0; + switch (m_norm_dir) { + case 0: + u_err = compute_error(XDir()); + break; + case 1: + u_err = compute_error(YDir()); + break; + case 2: + u_err = compute_error(ZDir()); + break; + default: + amrex::Abort("axis must be equal to 1 or 2"); + break; + } -void ChannelFlow::post_advance_work() {} + if (amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream f; + f.open(m_output_fname.c_str(), std::ios_base::app); + f << std::setprecision(12) << std::setw(m_w) << m_time.new_time() + << std::setw(m_w) << u_err << std::endl; + f.close(); + } +} + +void ChannelFlow::post_init_actions() +{ + if (m_laminar) { + output_error(); + } +} + +void ChannelFlow::post_advance_work() +{ + if (m_laminar) { + output_error(); + } +} } // namespace channel_flow } // namespace amr_wind diff --git a/amr-wind/physics/ConvectingTaylorVortex.H b/amr-wind/physics/ConvectingTaylorVortex.H index 9e6d9f4633..868929df97 100644 --- a/amr-wind/physics/ConvectingTaylorVortex.H +++ b/amr-wind/physics/ConvectingTaylorVortex.H @@ -119,9 +119,11 @@ private: const amr_wind::CFDSim& m_sim; const FieldRepo& m_repo; const amrex::AmrCore& m_mesh; + Field& m_velocity; Field& m_gradp; Field& m_density; + void output_error(); //! initial density value @@ -143,6 +145,7 @@ private: std::string m_output_fname{"ctv.log"}; bool m_activate_pressure{false}; + bool m_mesh_mapping{false}; }; } // namespace ctv } // namespace amr_wind diff --git a/amr-wind/physics/ConvectingTaylorVortex.cpp b/amr-wind/physics/ConvectingTaylorVortex.cpp index 98acbf4cd5..5be83aa69f 100644 --- a/amr-wind/physics/ConvectingTaylorVortex.cpp +++ b/amr-wind/physics/ConvectingTaylorVortex.cpp @@ -3,6 +3,7 @@ #include "AMReX_iMultiFab.H" #include "AMReX_MultiFabUtil.H" #include "AMReX_ParmParse.H" +#include "AMReX_ParReduce.H" #include "amr-wind/utilities/trig_ops.H" namespace amr_wind { @@ -94,6 +95,7 @@ ConvectingTaylorVortex::ConvectingTaylorVortex(const CFDSim& sim) , m_velocity(sim.repo().get_field("velocity")) , m_gradp(sim.repo().get_field("gp")) , m_density(sim.repo().get_field("density")) + , m_mesh_mapping(sim.has_mesh_mapping()) { { amrex::ParmParse pp("CTV"); @@ -128,15 +130,23 @@ void ConvectingTaylorVortex::initialize_fields( { using namespace utils; + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto u0 = m_u0; const auto v0 = m_v0; const auto omega = m_omega; const bool activate_pressure = m_activate_pressure; + const auto mesh_mapping = m_mesh_mapping; auto& velocity = m_velocity(level); auto& density = m_density(level); auto& pressure = m_repo.get_field("p")(level); auto& gradp = m_repo.get_field("gp")(level); + Field const* nu_coord_cc = + mesh_mapping ? &(m_repo.get_field("non_uniform_coord_cc")) : nullptr; + Field const* nu_coord_nd = + mesh_mapping ? &(m_repo.get_field("non_uniform_coord_nd")) : nullptr; density.setVal(m_rho); @@ -150,15 +160,19 @@ void ConvectingTaylorVortex::initialize_fields( for (amrex::MFIter mfi(velocity); mfi.isValid(); ++mfi) { const auto& vbx = mfi.validbox(); - const auto& dx = geom.CellSizeArray(); - const auto& problo = geom.ProbLoArray(); auto vel = velocity.array(mfi); auto gp = gradp.array(mfi); + amrex::Array4 nu_cc = + mesh_mapping ? ((*nu_coord_cc)(level).const_array(mfi)) + : amrex::Array4(); 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]; + amrex::Real x = mesh_mapping ? (nu_cc(i, j, k, 0)) + : (prob_lo[0] + (i + 0.5) * dx[0]); + amrex::Real y = mesh_mapping ? (nu_cc(i, j, k, 1)) + : (prob_lo[1] + (j + 0.5) * dx[1]); + vel(i, j, k, 0) = u_exact(u0, v0, omega, x, y, 0.0); vel(i, j, k, 1) = v_exact(u0, v0, omega, x, y, 0.0); vel(i, j, k, 2) = w_exact(u0, v0, omega, x, y, 0.0); @@ -173,11 +187,17 @@ void ConvectingTaylorVortex::initialize_fields( if (activate_pressure) { const auto& nbx = mfi.nodaltilebox(); auto pres = pressure.array(mfi); + amrex::Array4 nu_nd = + mesh_mapping ? ((*nu_coord_nd)(level).const_array(mfi)) + : amrex::Array4(); amrex::ParallelFor( nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + i * dx[0]; - const amrex::Real y = problo[1] + j * dx[1]; + amrex::Real x = mesh_mapping ? (nu_nd(i, j, k, 0)) + : (prob_lo[0] + i * dx[0]); + amrex::Real y = mesh_mapping ? (nu_nd(i, j, k, 1)) + : (prob_lo[1] + j * dx[1]); + pres(i, j, k, 0) = -0.25 * (std::cos(2.0 * utils::pi() * x) + std::cos(2.0 * utils::pi() * y)); @@ -189,7 +209,6 @@ void ConvectingTaylorVortex::initialize_fields( template amrex::Real ConvectingTaylorVortex::compute_error(const Field& field) { - amrex::Real error = 0.0; const amrex::Real time = m_time.new_time(); const auto u0 = m_u0; @@ -197,6 +216,14 @@ amrex::Real ConvectingTaylorVortex::compute_error(const Field& field) const auto omega = m_omega; T f_exact; const auto comp = f_exact.m_comp; + const auto mesh_mapping = m_mesh_mapping; + + Field const* nu_coord_cc = + mesh_mapping ? &(m_repo.get_field("non_uniform_coord_cc")) : nullptr; + Field const* mesh_fac_cc = + mesh_mapping + ? &(m_repo.get_mesh_mapping_field(amr_wind::FieldLoc::CELL)) + : nullptr; const int nlevels = m_repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { @@ -230,28 +257,44 @@ amrex::Real ConvectingTaylorVortex::compute_error(const Field& field) } const auto& dx = m_mesh.Geom(lev).CellSizeArray(); - const auto& problo = m_mesh.Geom(lev).ProbLoArray(); - const amrex::Real cell_vol = dx[0] * dx[1] * dx[2]; + const auto& prob_lo = m_mesh.Geom(lev).ProbLoArray(); const auto& fld = field(lev); - error += amrex::ReduceSum( - fld, level_mask, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, - amrex::Array4 const& fld_arr, - amrex::Array4 const& mask_arr) -> amrex::Real { - amrex::Real err_fab = 0.0; - - amrex::Loop(bx, [=, &err_fab](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 u = fld_arr(i, j, k, comp); - const amrex::Real u_exact = - f_exact(u0, v0, omega, x, y, time); - err_fab += cell_vol * mask_arr(i, j, k) * (u - u_exact) * - (u - u_exact); - }); - return err_fab; + auto const& fld_arr = fld.const_arrays(); + auto const& mask_arr = level_mask.const_arrays(); + amrex::MultiArray4 fac_arr = + mesh_mapping ? ((*mesh_fac_cc)(lev).const_arrays()) + : amrex::MultiArray4(); + amrex::MultiArray4 nu_cc = + mesh_mapping ? ((*nu_coord_cc)(lev).const_arrays()) + : amrex::MultiArray4(); + + error += amrex::ParReduce( + amrex::TypeList{}, + amrex::TypeList{}, fld, amrex::IntVect(0), + [=] AMREX_GPU_HOST_DEVICE(int box_no, int i, int j, int k) + -> amrex::GpuTuple { + auto const& fld_bx = fld_arr[box_no]; + auto const& mask_bx = mask_arr[box_no]; + + amrex::Real x = mesh_mapping ? (nu_cc[box_no](i, j, k, 0)) + : (prob_lo[0] + (i + 0.5) * dx[0]); + amrex::Real y = mesh_mapping ? (nu_cc[box_no](i, j, k, 1)) + : (prob_lo[1] + (j + 0.5) * dx[1]); + amrex::Real fac_x = + mesh_mapping ? (fac_arr[box_no](i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac_arr[box_no](i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac_arr[box_no](i, j, k, 2)) : 1.0; + + const amrex::Real u = fld_bx(i, j, k, comp); + const amrex::Real u_exact = f_exact(u0, v0, omega, x, y, time); + const amrex::Real cell_vol = + dx[0] * fac_x * dx[1] * fac_y * dx[2] * fac_z; + + return cell_vol * mask_bx(i, j, k) * (u - u_exact) * + (u - u_exact); }); } @@ -263,6 +306,7 @@ amrex::Real ConvectingTaylorVortex::compute_error(const Field& field) void ConvectingTaylorVortex::output_error() { + // TODO: gradp analytical solution has not been adjusted for mesh mapping const amrex::Real u_err = compute_error(m_velocity); const amrex::Real v_err = compute_error(m_velocity); const amrex::Real w_err = compute_error(m_velocity); diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp index 408caf252f..d4712d9ded 100644 --- a/amr-wind/projection/incflo_apply_nodal_projection.cpp +++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp @@ -125,16 +125,34 @@ void incflo::ApplyProjection( (!m_sim.pde_manager().constant_density() || m_sim.physics_manager().contains("MultiPhase")); + bool mesh_mapping = m_sim.has_mesh_mapping(); + auto& grad_p = m_repo.get_field("gp"); auto& pressure = m_repo.get_field("p"); auto& velocity = icns().fields().field; - + auto& velocity_old = icns().fields().field.state(amr_wind::FieldState::Old); + amr_wind::Field const* mesh_fac = + mesh_mapping + ? &(m_repo.get_mesh_mapping_field(amr_wind::FieldLoc::CELL)) + : nullptr; + amr_wind::Field const* mesh_detJ = + mesh_mapping ? &(m_repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::CELL)) + : nullptr; + + // TODO: Mesh mapping doesn't work with immersed boundaries // Do the pre pressure correction work -- this applies to IB only for (auto& pp : m_sim.physics()) { pp->pre_pressure_correction_work(); } + // ensure velocity is in stretched mesh space + if (velocity.in_uniform_space() && mesh_mapping) { + velocity.to_stretched_space(); + } + // Add the ( grad p /ro ) back to u* (note the +dt) + // Also account for mesh mapping in ( grad p /ro ) -> 1/fac * grad(p) * + // dt/rho if (!incremental) { for (int lev = 0; lev <= finest_level; lev++) { @@ -147,12 +165,23 @@ void incflo::ApplyProjection( Array4 const& u = velocity(lev).array(mfi); Array4 const& rho = density[lev]->const_array(mfi); Array4 const& gp = grad_p(lev).const_array(mfi); + amrex::Array4 fac = + mesh_mapping ? ((*mesh_fac)(lev).const_array(mfi)) + : amrex::Array4(); + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Real soverrho = scaling_factor / rho(i, j, k); - u(i, j, k, 0) += gp(i, j, k, 0) * soverrho; - u(i, j, k, 1) += gp(i, j, k, 1) * soverrho; - u(i, j, k, 2) += gp(i, j, k, 2) * soverrho; + amrex::Real fac_x = + mesh_mapping ? (fac(i, j, k, 0)) : 1.0; + amrex::Real fac_y = + mesh_mapping ? (fac(i, j, k, 1)) : 1.0; + amrex::Real fac_z = + mesh_mapping ? (fac(i, j, k, 2)) : 1.0; + + u(i, j, k, 0) += 1 / fac_x * gp(i, j, k, 0) * soverrho; + u(i, j, k, 1) += 1 / fac_y * gp(i, j, k, 1) * soverrho; + u(i, j, k, 2) += 1 / fac_z * gp(i, j, k, 2) * soverrho; }); } } @@ -192,21 +221,32 @@ void incflo::ApplyProjection( } } + // ensure velocity is in stretched mesh space + if (velocity_old.in_uniform_space() && mesh_mapping) { + velocity_old.to_stretched_space(); + } + // Define "vel" to be U^* - U^n rather than U^* if (proj_for_small_dt || incremental) { for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Subtract( - velocity(lev), velocity.state(amr_wind::FieldState::Old)(lev), - 0, 0, AMREX_SPACEDIM, 0); + velocity(lev), velocity_old(lev), 0, 0, AMREX_SPACEDIM, 0); } } - // Create sigma + // scale U^* to accommodate for mesh mapping -> U^bar = J/fac * U + if (mesh_mapping) { + velocity.to_uniform_space(); + } + + // Create sigma while accounting for mesh mapping + // sigma = 1/(fac^2)*J * dt/rho Vector sigma(finest_level + 1); - if (variable_density) { + if (variable_density || mesh_mapping) { + int ncomp = mesh_mapping ? AMREX_SPACEDIM : 1; for (int lev = 0; lev <= finest_level; ++lev) { sigma[lev].define( - grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev)); + grids[lev], dmap[lev], ncomp, 0, MFInfo(), Factory(lev)); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -215,9 +255,22 @@ void incflo::ApplyProjection( Box const& bx = mfi.tilebox(); Array4 const& sig = sigma[lev].array(mfi); Array4 const& rho = density[lev]->const_array(mfi); + amrex::Array4 fac = + mesh_mapping ? ((*mesh_fac)(lev).const_array(mfi)) + : amrex::Array4(); + amrex::Array4 detJ = + mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) + : amrex::Array4(); + amrex::ParallelFor( - bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - sig(i, j, k) = scaling_factor / rho(i, j, k); + bx, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { + amrex::Real fac_cc = + mesh_mapping ? (fac(i, j, k, n)) : 1.0; + amrex::Real det_j = + mesh_mapping ? (detJ(i, j, k)) : 1.0; + sig(i, j, k, n) = std::pow(fac_cc, -2.) * det_j * + scaling_factor / rho(i, j, k); }); } } @@ -235,12 +288,12 @@ void incflo::ApplyProjection( vel[lev]->setBndry(0.0); if (!proj_for_small_dt and !incremental) { set_inflow_velocity(lev, time, *vel[lev], 1); - // velocity.fillphysbc(lev, time, *vel[lev], 1); } } amr_wind::MLMGOptions options("nodal_proj"); - if (variable_density) { + + if (variable_density || mesh_mapping) { nodal_projector = std::make_unique( vel, GetVecOfConstPtrs(sigma), Geom(0, finest_level), options.lpinfo()); @@ -301,12 +354,16 @@ void incflo::ApplyProjection( amr_wind::io::print_mlmg_info( "Nodal_projection", nodal_projector->getMLMG()); + // scale U^* back to -> U = fac/J * U^bar + if (mesh_mapping) { + velocity.to_stretched_space(); + } + // Define "vel" to be U^{n+1} rather than (U^{n+1}-U^n) if (proj_for_small_dt || incremental) { for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Add( - velocity(lev), velocity.state(amr_wind::FieldState::Old)(lev), - 0, 0, AMREX_SPACEDIM, 0); + velocity(lev), velocity_old(lev), 0, 0, AMREX_SPACEDIM, 0); } } diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp index 2dfc4a657d..26f1213751 100644 --- a/amr-wind/setup/init.cpp +++ b/amr-wind/setup/init.cpp @@ -110,6 +110,10 @@ void incflo::InitialIterations() { auto& vel = icns().fields().field; + // ensure velocity is in stretched mesh space + if (vel.in_uniform_space() && m_sim.has_mesh_mapping()) { + vel.to_stretched_space(); + } vel.copy_state( amr_wind::FieldState::New, amr_wind::FieldState::Old); diff --git a/amr-wind/setup/set_background_pressure.cpp b/amr-wind/setup/set_background_pressure.cpp index c117fe1b64..a814b28a8a 100644 --- a/amr-wind/setup/set_background_pressure.cpp +++ b/amr-wind/setup/set_background_pressure.cpp @@ -1,13 +1,27 @@ #include "amr-wind/incflo.H" +#include "AMReX_ParmParse.H" + using namespace amrex; void incflo::set_background_pressure() { const auto problo = geom[0].ProbLoArray(); - const auto probhi = geom[0].ProbHiArray(); + // determine probhi based on if mesh is stretched + amrex::Vector probhi_physical{{0.0, 0.0, 0.0}}; + { + amrex::ParmParse pp("geometry"); + if (pp.contains("prob_hi_physical")) { + pp.getarr("prob_hi_physical", probhi_physical); + } else { + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + probhi_physical[d] = geom[0].ProbHiArray()[d]; + } + } + } GpuArray problen{ - {probhi[0] - problo[0], probhi[1] - problo[1], probhi[2] - problo[2]}}; + {probhi_physical[0] - problo[0], probhi_physical[1] - problo[1], + probhi_physical[2] - problo[2]}}; amrex::Vector m_gp0{{0.0, 0.0, 0.0}}; diff --git a/amr-wind/utilities/IOManager.cpp b/amr-wind/utilities/IOManager.cpp index 5e6df659ae..36e3a227ca 100644 --- a/amr-wind/utilities/IOManager.cpp +++ b/amr-wind/utilities/IOManager.cpp @@ -1,5 +1,8 @@ +#include +#include #include #include +#include #include "amr-wind/utilities/IOManager.H" #include "amr-wind/CFDSim.H" @@ -7,6 +10,7 @@ #include "amr-wind/utilities/io_utils.H" #include "amr-wind/utilities/DerivedQuantity.H" #include "amr-wind/utilities/DerivedQtyDefs.H" +#include "amr-wind/utilities/ncutils/nc_interface.H" #include "AMReX_ParmParse.H" #include "AMReX_PlotFileUtil.H" diff --git a/docs/sphinx/theory/mapping.rst b/docs/sphinx/theory/mapping.rst new file mode 100644 index 0000000000..73b8ba2fd3 --- /dev/null +++ b/docs/sphinx/theory/mapping.rst @@ -0,0 +1,130 @@ +.. _mapping: + +Mapping definition +================== + +Let the mapping of stretched mesh to uniform mesh be defined by :math:`{\mathbf X} = \Phi (\Xi) = (x(\chi), y(\eta), z(\xi))`. For demonstration purposes, we assume mapping only in the :math:`z`-direction. Note that +the normals to faces continue to be aligned with the grid axes. + +We define the velocity + +.. math:: \overline{U} = T U + +where + +.. math:: T = J [\nabla_\Xi \Phi]^{-1} + +and + +.. math:: J = {\rm det} [\nabla_\Xi \Phi] + +In our case :math:`T` is + +.. math:: + + \begin{bmatrix} + z_\xi & 0 & 0 \\ + 0 & z_\xi & 0\\ + 0 & 0 & 1 \\ + \end{bmatrix}, + \hspace{0.5in} J = ( z_\xi), + \hspace{0.5in} \overline{U} = (z_\xi u, \; z_\xi v, \; w). + +We can write the mesh spacing in :math:`{\mathbf X}`-space as + +.. math:: ( x_\chi dx, y_\eta dy, z_\xi dz) = ( dx, dy, z_\xi dz) + +where :math:`(dx,dy,dz)` are the mesh spacings in +:math:`{\mathbf \Xi}`-space. + +Handy identities +---------------- + +.. math:: \nabla_X = \frac{1}{J} T^T \nabla_\Xi + +.. math:: \nabla_X \cdot A = \frac{1}{J} \nabla_\Xi \cdot (T A) = \frac{1}{J} \nabla_\Xi \cdot (\overline{A}) + +for any vector field :math:`A` + +Governing equations in stretched mesh space +------------------------------------------- + +.. math:: + + JU_t + (\overline{U^{MAC}}\cdot \nabla_\Xi) \; U = + J(F_b -\frac{1}{\rho} \frac{T^{T}}{J} \nabla_\Xi \; p) + \frac{\mu}{\rho} \nabla_\Xi \cdot (\frac{1}{J} T T^T \nabla_\Xi U) , + +.. math:: J \rho_t + \nabla_\Xi \cdot (\rho \overline{U^{MAC}}) = 0, + +.. math:: J (\rho s)_t + \nabla_\Xi \cdot (\rho s \overline{U^{MAC}}) = S, + +.. math:: \nabla_\Xi \cdot \overline{U}= 0 \;\;\;, + +MAC Projection +-------------- + +In the MAC projection we want to solve + +.. math:: \nabla_X \cdot (\frac{1}{\rho} \; \nabla_X p) = \nabla_X \cdot U^{pred} + +then set + +.. math:: U^{MAC} = U^{pred} - \frac{1}{\rho} \nabla_X \; p + +Mapping the above to uniform coordinates yields + +.. math:: \nabla_\Xi \cdot (\frac{1}{J \rho} T T^T \nabla_\Xi \; p) = \nabla_\Xi \cdot \overline{U^{pred}} + +The MAC projection from AMReX will then return + +.. math:: \overline{U^{MAC}} = \overline{U^{pred}} - \frac{1}{\rho} T (\frac{1}{J} T^T \nabla_\Xi p) = \overline{U^{pred}} - \frac{1}{J \rho} T T^T \nabla_\Xi \; p + +Nodal Projection +---------------- + +Analogously to the MAC projection we want to solve + +.. math:: \nabla_\Xi \cdot (\frac{1}{J \rho} T T^T \nabla_\Xi \; p) = T^T \nabla_\Xi \cdot ( U^{n+1,*} + \frac{1}{J} T^T \nabla_\Xi p^{n-1/2} ) = \nabla_\Xi \overline{U}^{n+1,*} + +where here :math:`\overline{U}^{n+1,*} = T U^{n+1,*}` is cell-centered +with +:math:`U^{n+1,*} = U^{n+1,*} + \frac{1}{J} T^T \nabla_\Xi p^{n-1/2}`. +Note, although :math:`p` is defined on the nodes, :math:`\nabla p` is +cell-centered. + +We construct the divergence and gradient operators with a finite element +approach, so the mapping enters in here via integrals over cell volumes. +Thus we will only need :math:`z_\xi` at cell centers; this is the same +metric term used in the MAC projector. + +As with the MAC projection we want to create +:math:`\overline{U}^{n+1,*} = T U^{n+1,*}` +as the velocity field we will send in to the nodal projector. + +The nodal projection will return + +.. math:: \overline{U}^{n+1} = \overline{U}^{n+1,*} - \frac{1}{J \rho} T T^T \nabla_\Xi \; p + +Note that – unlike in the MAC projection – we want :math:`U^{n+1}` +rather than :math:`\overline{U}^{n+1},` so after the projection we must +define + +.. math:: U = T^{-1} \overline{U} + +Initial conditions & computing the timestep +------------------------------------------- + +The initial conditions for any problem as well as any other global computations +such as evaluating the CFL or the adaptive timestep size should be done so in the +stretched coordinates. + +Exceptions to the current implementation +---------------------------------------- + +- As a first pass, the stretched mesh capability is implemented only for the MOL scheme. + +- The stretched mesh capability has been tested and verified only for single-level AMR meshes. + +- The mesh stretching capability requires the use of multi-component velocity solves by setting ``velocity_diffusion.use_tensor_operator = false`` + +- Efforts are underway to extend the capability beyond laminar physics. diff --git a/docs/sphinx/theory/theory.rst b/docs/sphinx/theory/theory.rst index e92a1e412a..d8cb8f2659 100644 --- a/docs/sphinx/theory/theory.rst +++ b/docs/sphinx/theory/theory.rst @@ -99,6 +99,17 @@ and .. math:: {p}^{n+1/2, \ast} = \phi +Solving physics on a stretched AMReX mesh +------------------------------------------ + +Often for simulations involving walls, (e.g., channel flows, complex terrains etc.) it is desirable to have finer mesh near the wall which gradually coarsens only in the wall-normal direction. Consequently, modifications within AMR-Wind are underway to support solving a non-uniformly spaced mesh while still using most of AMReX's machinery directed at uniformly-spaced Cartesian meshes. The governing equations solved on a non-uniform stretched mesh are further explained below - + +.. toctree:: + :glob: + :maxdepth: 2 + + mapping.rst + Multiphase flow modelling ------------------------------------ diff --git a/docs/sphinx/user/inputs_geometry.rst b/docs/sphinx/user/inputs_geometry.rst index 780256a2dc..2ca6157f6b 100644 --- a/docs/sphinx/user/inputs_geometry.rst +++ b/docs/sphinx/user/inputs_geometry.rst @@ -14,6 +14,14 @@ This section deals with inputs related to the problem domain. **type:** List of 3 real numbers, mandatory The coordinates of the *upper corner* of the computational domain bounding box. + +.. input_param:: geometry.prob_hi_physical + + **type:** List of 3 real numbers, optional, default = ``geometry.prob_hi`` + + The coordinates of the *upper corner* of the computational domain bounding box for a + `Stretched Mesh + `_. .. input_param:: geometry.is_periodic @@ -21,3 +29,10 @@ This section deals with inputs related to the problem domain. Flags indicating whether the flow is periodic in the ``x``, ``y``, or ``z`` directions respectively. + +.. input_param:: geometry.mesh_mapping + + **type:** String, optional, default = ``ConstantMap`` + + Define the style of mapping between the stretched coordinates and the uniform coordinates. The default map is a constant scaling map with + ``ConstantMap.scaling_factor = 1.0 1.0 1.0``. diff --git a/submods/amrex b/submods/amrex index fa46bc9eff..0f55de0e2e 160000 --- a/submods/amrex +++ b/submods/amrex @@ -1 +1 @@ -Subproject commit fa46bc9eff1306b8c3eeaf1ae32d070271606bba +Subproject commit 0f55de0e2ef6a66ae2c0e5aabc7979c8cc140e04 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 573e61a8b5..1b221472f5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -201,6 +201,8 @@ add_test_re(abl_bndry_output_native) if (NOT AMR_WIND_ENABLE_CUDA) add_test_re(ctv_godunov_plm) + add_test_re(ctv_mol_mesh_map) + add_test_re(ctv_mol_mesh_map_explicit) endif() if (AMR_WIND_ENABLE_NETCDF) @@ -218,6 +220,9 @@ endif() if (AMR_WIND_ENABLE_HYPRE) add_test_re(abl_godunov_hypre) add_test_re(channel_kwsst_hypre) + add_test_re(channel_mol_mesh_map_x) + add_test_re(channel_mol_mesh_map_y) + add_test_re(channel_mol_mesh_map_z) endif() #============================================================================= diff --git a/test/test_files/channel_mol_mesh_map_x/channel_mol_mesh_map_x.i b/test/test_files/channel_mol_mesh_map_x/channel_mol_mesh_map_x.i new file mode 100644 index 0000000000..19ae1f17e5 --- /dev/null +++ b/test/test_files/channel_mol_mesh_map_x/channel_mol_mesh_map_x.i @@ -0,0 +1,87 @@ +#include +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 100.00 # Max (simulated) time to evolve +time.max_step = -1000 # Max number of time steps +incflo.initial_iterations = 0 +incflo.do_initial_proj = 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.03 # Use this constant dt if > 0 +time.cfl = 0.5 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1000 # Steps between plot files +time.checkpoint_interval = -100 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.density = 1.0 # Reference density +incflo.use_godunov = 0 +transport.viscosity = 0.005 +turbulence.model = Laminar + +ICNS.source_terms = BodyForce +BodyForce.magnitude = 6e-2 0 0 +incflo.physics = ChannelFlow +ChannelFlow.density = 1.0 +ChannelFlow.Laminar = true +ChannelFlow.Turbulent_DNS = false +ChannelFlow.Mean_Velocity = 1.0 + +io.output_default_variables = 1 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 32 32 16 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0.0 0.0 0.0 # Lo corner coordinates +geometry.prob_hi = 6.0 1.0 1.0 # Hi corner coordinates +geometry.is_periodic = 1 0 1 # Periodicity x y z (0/1) + +geometry.mesh_mapping = ChannelFlowMap +ChannelFlowMap.beta = 0 3 0 + +velocity_diffusion.use_tensor_operator = false +incflo.diffusion_type = 2 + +# Boundary conditions +ylo.type = "no_slip_wall" +yhi.type = "no_slip_wall" + +diffusion.max_coarsening_level = 0 + +incflo.verbose = 0 +nodal_proj.mg_atol = 1.0e-09 +nodal_proj.verbose = 0 +nodal_proj.bottom_solver = hypre +nodal_proj.max_coarsening_level = 0 +# +mac_proj.mg_rtol = 1.0e-11 +mac_proj.mg_atol = 1.0e-09 +mac_proj.do_semicoarsening = true +mac_proj.bottom_solver = hypre +mac_proj.bottom_verbose = 0 +mac_proj.max_coarsening_level = 0 +mac_proj.bottom_rtol = 1.0e-12 +mac_proj.bottom_atol = 1.0e-16 +# +hypre.hypre_solver = GMRES +hypre.hypre_preconditioner = BoomerAMG +hypre.verbose = 0 +hypre.bamg_verbose = 0 +hypre.num_krylov = 20 +hypre.bamg_max_levels = 4 +hypre.bamg_num_sweeps = 1 +# diff --git a/test/test_files/channel_mol_mesh_map_y/channel_mol_mesh_map_y.i b/test/test_files/channel_mol_mesh_map_y/channel_mol_mesh_map_y.i new file mode 100644 index 0000000000..44d684242a --- /dev/null +++ b/test/test_files/channel_mol_mesh_map_y/channel_mol_mesh_map_y.i @@ -0,0 +1,91 @@ +#include +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 300.00 # Max (simulated) time to evolve +time.max_step = -1000 # Max number of time steps +incflo.initial_iterations = 0 +incflo.do_initial_proj = 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.03 # Use this constant dt if > 0 +time.cfl = 0.5 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1000 # Steps between plot files +time.checkpoint_interval = -100 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.density = 1.0 # Reference density +incflo.use_godunov = 0 +transport.viscosity = 0.005 +turbulence.model = Laminar + +ICNS.source_terms = BodyForce +BodyForce.magnitude = 0.0 6e-2 0.0 +incflo.physics = ChannelFlow +ChannelFlow.density = 1.0 +ChannelFlow.Laminar = true +ChannelFlow.flow_direction = 1 +ChannelFlow.normal_direction = 2 +ChannelFlow.Turbulent_DNS = false +ChannelFlow.Mean_Velocity = 1.0 + +io.output_default_variables = 1 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 16 32 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0.0 0.0 0.0 # Lo corner coordinates +geometry.prob_hi = 1.0 6.0 1.0 # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) + +geometry.mesh_mapping = ChannelFlowMap +ChannelFlowMap.beta = 0 0 3 + +velocity_diffusion.use_tensor_operator = false +incflo.diffusion_type = 2 + +# Boundary conditions +zlo.type = "no_slip_wall" +zhi.type = "no_slip_wall" + +incflo.verbose = 0 + +diffusion.max_coarsening_level = 0 + +nodal_proj.mg_atol = 1.0e-09 +nodal_proj.verbose = 0 +nodal_proj.bottom_solver = hypre +nodal_proj.max_coarsening_level = 0 + +mac_proj.mg_rtol = 1.0e-11 +mac_proj.mg_atol = 1.0e-09 +mac_proj.do_semicoarsening = true +mac_proj.bottom_solver = hypre +mac_proj.bottom_verbose = 0 +mac_proj.max_coarsening_level = 0 +mac_proj.bottom_rtol = 1.0e-12 +mac_proj.bottom_atol = 1.0e-16 +# +hypre.hypre_solver = GMRES +hypre.hypre_preconditioner = BoomerAMG +hypre.recompute_preconditioner = 0 +hypre.verbose = 0 +hypre.bamg_verbose = 0 +hypre.num_krylov = 20 +hypre.bamg_max_levels = 4 +hypre.bamg_num_sweeps = 1 +# diff --git a/test/test_files/channel_mol_mesh_map_z/channel_mol_mesh_map_z.i b/test/test_files/channel_mol_mesh_map_z/channel_mol_mesh_map_z.i new file mode 100644 index 0000000000..5cee2eb144 --- /dev/null +++ b/test/test_files/channel_mol_mesh_map_z/channel_mol_mesh_map_z.i @@ -0,0 +1,91 @@ +#include +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 300.00 # Max (simulated) time to evolve +time.max_step = -1000 # Max number of time steps +incflo.initial_iterations = 0 +incflo.do_initial_proj = 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.03 # Use this constant dt if > 0 +time.cfl = 0.5 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1000 # Steps between plot files +time.checkpoint_interval = -100 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.density = 1.0 # Reference density +incflo.use_godunov = 0 +transport.viscosity = 0.005 +turbulence.model = Laminar + +ICNS.source_terms = BodyForce +BodyForce.magnitude = 0.0 0.0 6e-2 +incflo.physics = ChannelFlow +ChannelFlow.density = 1.0 +ChannelFlow.Laminar = true +ChannelFlow.flow_direction = 2 +ChannelFlow.normal_direction = 0 +ChannelFlow.Turbulent_DNS = false +ChannelFlow.Mean_Velocity = 1.0 + +io.output_default_variables = 1 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 32 16 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0.0 0.0 0.0 # Lo corner coordinates +geometry.prob_hi = 1.0 1.0 6.0 # Hi corner coordinates +geometry.is_periodic = 0 1 1 # Periodicity x y z (0/1) + +geometry.mesh_mapping = ChannelFlowMap +ChannelFlowMap.beta = 3 0 0 + +velocity_diffusion.use_tensor_operator = false +incflo.diffusion_type = 2 + +# Boundary conditions +xlo.type = "no_slip_wall" +xhi.type = "no_slip_wall" + +incflo.verbose = 0 + +diffusion.max_coarsening_level = 0 + +nodal_proj.mg_atol = 1.0e-09 +nodal_proj.verbose = 0 +nodal_proj.bottom_solver = hypre +nodal_proj.max_coarsening_level = 0 + +mac_proj.mg_rtol = 1.0e-11 +mac_proj.mg_atol = 1.0e-09 +mac_proj.do_semicoarsening = true +mac_proj.bottom_solver = hypre +mac_proj.bottom_verbose = 0 +mac_proj.max_coarsening_level = 0 +mac_proj.bottom_rtol = 1.0e-12 +mac_proj.bottom_atol = 1.0e-16 +# +hypre.hypre_solver = GMRES +hypre.hypre_preconditioner = BoomerAMG +hypre.recompute_preconditioner = 0 +hypre.verbose = 0 +hypre.bamg_verbose = 0 +hypre.num_krylov = 20 +hypre.bamg_max_levels = 4 +hypre.bamg_num_sweeps = 1 +# diff --git a/test/test_files/ctv_mol_mesh_map/ctv_mol_mesh_map.i b/test/test_files/ctv_mol_mesh_map/ctv_mol_mesh_map.i new file mode 100644 index 0000000000..855f23664a --- /dev/null +++ b/test/test_files/ctv_mol_mesh_map/ctv_mol_mesh_map.i @@ -0,0 +1,53 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 0.2 # Max (simulated) time to evolve +time.max_step = 20 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.01 # Use this constant dt if > 0 +time.cfl = 0.5 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 20 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.output_default_variables = 0 +io.outputs = density p +io.derived_outputs = "components(velocity,0,1)" "components(gp,0,1)" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.use_godunov = 0 +incflo.godunov_type = "plm" +transport.viscosity = 1e-3 +transport.laminar_prandtl = 1.0 +turbulence.model = Laminar + + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 32 32 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1. 1. 1. # Hi corner coordinates +geometry.prob_hi_physical = 2. 2. 2. +geometry.is_periodic = 1 1 1 # Periodicity x y z (0/1) + +geometry.mesh_mapping = ConstantMap +ConstantMap.scaling_factor = 2.0 2.0 2.0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INITIAL CONDITIONS # +#.......................................# +incflo.physics = ConvectingTaylorVortex +CTV.activate_pressure = 1 diff --git a/test/test_files/ctv_mol_mesh_map_explicit/ctv_mol_mesh_map_explicit.i b/test/test_files/ctv_mol_mesh_map_explicit/ctv_mol_mesh_map_explicit.i new file mode 100644 index 0000000000..af56c6f06e --- /dev/null +++ b/test/test_files/ctv_mol_mesh_map_explicit/ctv_mol_mesh_map_explicit.i @@ -0,0 +1,53 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 0.2 # Max (simulated) time to evolve +time.max_step = 20 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.01 # Use this constant dt if > 0 +time.cfl = 0.5 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 20 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.output_default_variables = 0 +io.outputs = density p +io.derived_outputs = "components(velocity,0,1)" "components(gp,0,1)" + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.use_godunov = 0 +incflo.godunov_type = "plm" +transport.viscosity = 1e-3 +transport.laminar_prandtl = 1.0 +turbulence.model = Laminar +incflo.diffusion_type = 0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 32 32 32 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1. 1. 1. # Hi corner coordinates +geometry.prob_hi_physical = 2. 2. 2. +geometry.is_periodic = 1 1 1 # Periodicity x y z (0/1) + +geometry.mesh_mapping = ConstantMap +ConstantMap.scaling_factor = 2.0 2.0 2.0 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INITIAL CONDITIONS # +#.......................................# +incflo.physics = ConvectingTaylorVortex +CTV.activate_pressure = 1 diff --git a/unit_tests/equation_systems/test_pde.cpp b/unit_tests/equation_systems/test_pde.cpp index e7bc898d59..2fddb6c6cc 100644 --- a/unit_tests/equation_systems/test_pde.cpp +++ b/unit_tests/equation_systems/test_pde.cpp @@ -13,6 +13,7 @@ TEST_F(PDETest, test_pde_create_godunov) pp.add("use_godunov", 1); initialize_mesh(); + auto& pde_mgr = mesh().sim().pde_manager(); pde_mgr.register_icns(); pde_mgr.register_transport_pde("Temperature"); @@ -29,6 +30,7 @@ TEST_F(PDETest, test_pde_create_mol) pp.add("use_godunov", 0); initialize_mesh(); + auto& pde_mgr = mesh().sim().pde_manager(); pde_mgr.register_icns(); pde_mgr.register_transport_pde("Temperature");