From 8260c9120f330d5949c66ca838e189fbf910cbac Mon Sep 17 00:00:00 2001 From: whorne Date: Thu, 15 Aug 2024 12:10:34 -0600 Subject: [PATCH] Mask balanced buoyancy terms from walls --- include/LowMachEquationSystem.h | 1 + .../ngp_algorithms/NodalBuoyancyFuncUtil.h | 65 +++++++++++++++++++ src/LowMachEquationSystem.C | 16 +++++ src/edge_kernels/ContinuityEdgeSolverAlg.C | 16 ++++- src/ngp_algorithms/MdotEdgeAlg.C | 15 ++++- 5 files changed, 107 insertions(+), 6 deletions(-) create mode 100644 include/ngp_algorithms/NodalBuoyancyFuncUtil.h diff --git a/include/LowMachEquationSystem.h b/include/LowMachEquationSystem.h index d0327ee26..254c8994c 100644 --- a/include/LowMachEquationSystem.h +++ b/include/LowMachEquationSystem.h @@ -216,6 +216,7 @@ class MomentumEquationSystem : public EquationSystem std::unique_ptr tviscAlg_{nullptr}; std::unique_ptr pecletAlg_{nullptr}; std::unique_ptr ablWallNodeMask_{nullptr}; + std::unique_ptr buoyancySrcMask_{nullptr}; CourantReAlgDriver cflReAlgDriver_; std::unique_ptr AMSAlgDriver_{nullptr}; diff --git a/include/ngp_algorithms/NodalBuoyancyFuncUtil.h b/include/ngp_algorithms/NodalBuoyancyFuncUtil.h new file mode 100644 index 000000000..f1c2bdc05 --- /dev/null +++ b/include/ngp_algorithms/NodalBuoyancyFuncUtil.h @@ -0,0 +1,65 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#ifndef NODALBUOYANCYFUNCUTIL_H_ +#define NODALBUOYANCYFUNCUTIL_H_ + +#include +#include +#include +#include "stk_mesh/base/NgpField.hpp" +#include +#include +#include +#include +#include + +namespace sierra { +namespace nalu { + +class Realm; + +class NodalBuoyancyFuncUtil : public Algorithm +{ +public: + NodalBuoyancyFuncUtil(Realm& realm, stk::mesh::Part* part) + : Algorithm(realm, part), + maskNodeIndex_(get_field_ordinal( + realm.meta_data(), "buoyancy_source_mask", stk::topology::NODE_RANK)) + { + } + virtual ~NodalBuoyancyFuncUtil() = default; + void execute() override + { + + using Traits = nalu_ngp::NGPMeshTraits; + + const auto& meta = realm_.meta_data(); + const auto ngpMesh = realm_.ngp_mesh(); + const auto& fieldMgr = realm_.ngp_field_manager(); + auto myNodeMask = fieldMgr.get_field(maskNodeIndex_); + + const stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectUnion(partVec_) & !(realm_.get_inactive_selector()); + + nalu_ngp::run_entity_algorithm( + "BuoyancySourceFuncMaskUtil", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const Traits::MeshIndex& meshIdx) { + myNodeMask.get(meshIdx, 0) = 0; + }); + myNodeMask.modify_on_device(); + } + +private: + unsigned maskNodeIndex_{stk::mesh::InvalidOrdinal}; +}; +} // namespace nalu +} // namespace sierra +#endif diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C index 4096224a0..f728d91f6 100644 --- a/src/LowMachEquationSystem.C +++ b/src/LowMachEquationSystem.C @@ -122,6 +122,7 @@ #include "ngp_algorithms/WallFuncGeometryAlg.h" #include "ngp_algorithms/DynamicPressureOpenAlg.h" #include "ngp_algorithms/MomentumABLWallFuncMaskUtil.h" +#include "ngp_algorithms/NodalBuoyancyFuncUtil.h" #include "ngp_utils/NgpLoopUtils.h" #include "ngp_utils/NgpFieldBLAS.h" #include "ngp_utils/NgpFieldUtils.h" @@ -1103,6 +1104,9 @@ MomentumEquationSystem::initial_work() if (ablWallNodeMask_) ablWallNodeMask_->execute(); + if (buoyancySrcMask_) + buoyancySrcMask_->execute(); + // proceed with a bunch of initial work; wrap in timer { const double timeA = NaluEnv::self().nalu_time(); @@ -1197,6 +1201,10 @@ MomentumEquationSystem::register_nodal_fields( buoyancy_source_weight_ = &(meta_data.declare_field( stk::topology::NODE_RANK, "buoyancy_source_weight")); stk::mesh::put_field_on_mesh(*buoyancy_source_weight_, selector, nullptr); + ScalarFieldType& node_mask = realm_.meta_data().declare_field( + stk::topology::NODE_RANK, "buoyancy_source_mask"); + double one = 1; + stk::mesh::put_field_on_mesh(node_mask, selector, &one); } if (realm_.is_turbulent()) { @@ -1948,6 +1956,14 @@ MomentumEquationSystem::register_wall_bc( stk::mesh::put_field_on_mesh( *wallNormalDistanceBip, *part, numScsBip, nullptr); + if (realm_.solutionOptions_->use_balanced_buoyancy_force_) { + if (!buoyancySrcMask_) { + buoyancySrcMask_.reset(new NodalBuoyancyFuncUtil(realm_, part)); + } else { + buoyancySrcMask_->partVec_.push_back(part); + } + } + // need wall friction velocity for TKE boundary condition if (RANSAblBcApproach_) { const AlgorithmType wfAlgType = WALL_FCN; diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C index c607dd1db..da8b4ac67 100644 --- a/src/edge_kernels/ContinuityEdgeSolverAlg.C +++ b/src/edge_kernels/ContinuityEdgeSolverAlg.C @@ -82,6 +82,11 @@ ContinuityEdgeSolverAlg::execute() : fieldMgr.get_field( get_field_ordinal(realm_.meta_data(), "buoyancy_source")); + auto source_mask = !add_balanced_forcing + ? fieldMgr.get_field(densityNp1_) + : fieldMgr.get_field(get_field_ordinal( + realm_.meta_data(), "buoyancy_source_mask")); + stk::mesh::NgpField edgeFaceVelMag; bool needs_gcl = false; if (realm_.has_mesh_deformation()) { @@ -99,6 +104,7 @@ ContinuityEdgeSolverAlg::execute() udiag.sync_to_device(); edgeAreaVec.sync_to_device(); source.sync_to_device(); + source_mask.sync_to_device(); run_algorithm( realm_.bulk_data(), @@ -139,8 +145,10 @@ ContinuityEdgeSolverAlg::execute() DblType tmdot = -projTimeScale * (pressureR - pressureL) * asq * inv_axdx; if (add_balanced_forcing) { + const DblType masked_weights = + 0.5 * (source_mask.get(nodeL, 0) + source_mask.get(nodeR, 0)); for (int d = 0; d < ndim; ++d) { - tmdot += projTimeScale * av[d] * gravity[d] * rhoIp; + tmdot += projTimeScale * av[d] * gravity[d] * rhoIp * masked_weights; } } if (needs_gcl) { @@ -159,8 +167,10 @@ ContinuityEdgeSolverAlg::execute() DblType GjIp = 0.5 * (Gpdx.get(nodeR, d) / (udiagR) + Gpdx.get(nodeL, d) / (udiagL)); if (add_balanced_forcing) { - GjIp -= 0.5 * ((source.get(nodeR, d)) / (udiagR) + - (source.get(nodeL, d)) / (udiagL)); + GjIp -= + 0.5 * + ((source_mask.get(nodeR, 0) * source.get(nodeR, d)) / (udiagR) + + (source_mask.get(nodeL, 0) * source.get(nodeL, d)) / (udiagL)); } tmdot += (interpTogether * rhoUjIp + om_interpTogether * rhoIp * ujIp + GjIp) * diff --git a/src/ngp_algorithms/MdotEdgeAlg.C b/src/ngp_algorithms/MdotEdgeAlg.C index 9565c3e30..a95305505 100644 --- a/src/ngp_algorithms/MdotEdgeAlg.C +++ b/src/ngp_algorithms/MdotEdgeAlg.C @@ -94,6 +94,10 @@ MdotEdgeAlg::execute() ? fieldMgr.get_field(Gpdx_) : fieldMgr.get_field( get_field_ordinal(realm_.meta_data(), "buoyancy_source")); + auto source_mask = !add_balanced_forcing + ? fieldMgr.get_field(densityNp1_) + : fieldMgr.get_field(get_field_ordinal( + realm_.meta_data(), "buoyancy_source_mask")); mdot.clear_sync_state(); coordinates.sync_to_device(); @@ -104,6 +108,7 @@ MdotEdgeAlg::execute() udiag.sync_to_device(); edgeAreaVec.sync_to_device(); source.sync_to_device(); + source_mask.sync_to_device(); const stk::mesh::Selector sel = meta.locally_owned_part() & stk::mesh::selectUnion(partVec_) & @@ -146,8 +151,10 @@ MdotEdgeAlg::execute() DblType tmdot = -projTimeScale * (pressureR - pressureL) * asq * inv_axdx; if (add_balanced_forcing) { + const DblType masked_weights = + 0.5 * (source_mask.get(nodeL, 0) + source_mask.get(nodeR, 0)); for (int d = 0; d < ndim; ++d) { - tmdot += projTimeScale * av[d] * gravity[d] * rhoIp; + tmdot += projTimeScale * av[d] * gravity[d] * rhoIp * masked_weights; } } @@ -166,8 +173,10 @@ MdotEdgeAlg::execute() DblType GjIp = 0.5 * (Gpdx.get(nodeR, d) / udiagR + Gpdx.get(nodeL, d) / udiagL); if (add_balanced_forcing) { - GjIp -= 0.5 * ((source.get(nodeR, d)) / (udiagR) + - (source.get(nodeL, d)) / (udiagL)); + GjIp -= + 0.5 * + ((source_mask.get(nodeR, 0) * source.get(nodeR, d)) / (udiagR) + + (source_mask.get(nodeL, 0) * source.get(nodeL, d)) / (udiagL)); } tmdot +=