From 763116e1f0c064c6dc2b1d149abe6f69f2075a7f Mon Sep 17 00:00:00 2001 From: whorne Date: Tue, 6 Aug 2024 13:43:13 -0600 Subject: [PATCH 1/9] Add balanced buoyancy --- include/LowMachEquationSystem.h | 5 + include/SolutionOptions.h | 1 + include/ngp_algorithms/BuoyancySourceAlg.h | 52 ++++++++ .../ngp_algorithms/NodalBuoyancyAlgDriver.h | 41 ++++++ .../node_kernels/MomentumBuoyancyNodeKernel.h | 4 + src/LowMachEquationSystem.C | 27 +++- src/SolutionOptions.C | 4 + src/edge_kernels/ContinuityEdgeSolverAlg.C | 30 ++++- src/ngp_algorithms/BuoyancySourceAlg.C | 122 +++++++++++++++++ src/ngp_algorithms/CMakeLists.txt | 2 + src/ngp_algorithms/MdotEdgeAlg.C | 32 ++++- src/ngp_algorithms/NodalBuoyancyAlgDriver.C | 124 ++++++++++++++++++ src/node_kernels/MomentumBuoyancyNodeKernel.C | 14 +- 13 files changed, 451 insertions(+), 7 deletions(-) create mode 100644 include/ngp_algorithms/BuoyancySourceAlg.h create mode 100644 include/ngp_algorithms/NodalBuoyancyAlgDriver.h create mode 100644 src/ngp_algorithms/BuoyancySourceAlg.C create mode 100644 src/ngp_algorithms/NodalBuoyancyAlgDriver.C diff --git a/include/LowMachEquationSystem.h b/include/LowMachEquationSystem.h index 29c9ffb7c6..d0327ee260 100644 --- a/include/LowMachEquationSystem.h +++ b/include/LowMachEquationSystem.h @@ -18,6 +18,7 @@ #include "AMSAlgDriver.h" #include "ngp_algorithms/NodalGradAlgDriver.h" +#include "ngp_algorithms/NodalBuoyancyAlgDriver.h" #include "ngp_algorithms/WallFricVelAlgDriver.h" #include "ngp_algorithms/EffDiffFluxCoeffAlg.h" #include "ngp_algorithms/CourantReAlgDriver.h" @@ -204,7 +205,11 @@ class MomentumEquationSystem : public EquationSystem ScalarFieldType* evisc_; ScalarFieldType* iddesRansIndicator_; + VectorFieldType* buoyancy_source_; + ScalarFieldType* buoyancy_source_weight_; + TensorNodalGradAlgDriver nodalGradAlgDriver_; + NodalBuoyancyAlgDriver nodalBuoyancyAlgDriver_; WallFricVelAlgDriver wallFuncAlgDriver_; NgpAlgDriver dynPressAlgDriver_; std::unique_ptr diffFluxCoeffAlg_{nullptr}; diff --git a/include/SolutionOptions.h b/include/SolutionOptions.h index 4282073162..5257abd9be 100644 --- a/include/SolutionOptions.h +++ b/include/SolutionOptions.h @@ -100,6 +100,7 @@ class SolutionOptions bool has_set_boussinesq_time_scale(); + bool use_balanced_buoyancy_force_{false}; bool realm_has_vof_{false}; double hybridDefault_; diff --git a/include/ngp_algorithms/BuoyancySourceAlg.h b/include/ngp_algorithms/BuoyancySourceAlg.h new file mode 100644 index 0000000000..9ea50f6a6d --- /dev/null +++ b/include/ngp_algorithms/BuoyancySourceAlg.h @@ -0,0 +1,52 @@ +// 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 BUOYANCYSOURCEALG_H +#define BUOYANCYSOURCEALG_H + +#include "Algorithm.h" +#include "FieldTypeDef.h" + +#include "stk_mesh/base/Types.hpp" + +namespace sierra { +namespace nalu { + +class BuoyancySourceAlg : public Algorithm +{ + +public: + using DblType = double; + + BuoyancySourceAlg( + Realm&, + stk::mesh::Part*, + VectorFieldType* source, + ScalarFieldType* source_weight); + + virtual ~BuoyancySourceAlg() = default; + + virtual void execute() override; + +private: + unsigned source_{stk::mesh::InvalidOrdinal}; + unsigned source_weight_{stk::mesh::InvalidOrdinal}; + unsigned edgeAreaVec_{stk::mesh::InvalidOrdinal}; + unsigned dualNodalVol_{stk::mesh::InvalidOrdinal}; + unsigned coordinates_{stk::mesh::InvalidOrdinal}; + unsigned density_{stk::mesh::InvalidOrdinal}; + + //! Maximum size for static arrays used within device loops + static constexpr int NDimMax = 3; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* NODALGRADEDGEALG_H */ diff --git a/include/ngp_algorithms/NodalBuoyancyAlgDriver.h b/include/ngp_algorithms/NodalBuoyancyAlgDriver.h new file mode 100644 index 0000000000..e167babd83 --- /dev/null +++ b/include/ngp_algorithms/NodalBuoyancyAlgDriver.h @@ -0,0 +1,41 @@ +// 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 NODALBUOYANCYALGDRIVER_H +#define NODALBUOYANCYALGDRIVER_H + +#include "ngp_algorithms/NgpAlgDriver.h" +#include "FieldTypeDef.h" + +namespace sierra { +namespace nalu { + +class NodalBuoyancyAlgDriver : public NgpAlgDriver +{ +public: + NodalBuoyancyAlgDriver(Realm&, const std::string&, const std::string&); + + virtual ~NodalBuoyancyAlgDriver() = default; + + //! Reset fields before calling algorithms + virtual void pre_work() override; + + //! Synchronize fields after algorithms have done their work + virtual void post_work() override; + +private: + //! Field that is synchronized pre/post updates + const std::string sourceName_; + const std::string sourceweightName_; +}; + +} // namespace nalu +} // namespace sierra + +#endif /* NODALBUOYANCYALGDRIVER_H */ diff --git a/include/node_kernels/MomentumBuoyancyNodeKernel.h b/include/node_kernels/MomentumBuoyancyNodeKernel.h index f938db064e..83a809d6b5 100644 --- a/include/node_kernels/MomentumBuoyancyNodeKernel.h +++ b/include/node_kernels/MomentumBuoyancyNodeKernel.h @@ -45,11 +45,15 @@ class MomentumBuoyancyNodeKernel private: stk::mesh::NgpField dualNodalVolume_; stk::mesh::NgpField densityNp1_; + stk::mesh::NgpField source_; + const int nDim_; + const bool use_balanced_buoyancy_; NodeKernelTraits::DblType rhoRef_; unsigned dualNodalVolumeID_{stk::mesh::InvalidOrdinal}; unsigned densityNp1ID_{stk::mesh::InvalidOrdinal}; + unsigned sourceID_{stk::mesh::InvalidOrdinal}; NALU_ALIGNED NodeKernelTraits::DblType gravity_[NodeKernelTraits::NDimMax]; }; diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C index 5797d6762d..92f90328f3 100644 --- a/src/LowMachEquationSystem.C +++ b/src/LowMachEquationSystem.C @@ -101,6 +101,7 @@ // ngp #include "ngp_algorithms/ABLWallFrictionVelAlg.h" #include "ngp_algorithms/ABLWallFluxesAlg.h" +#include "ngp_algorithms/BuoyancySourceAlg.h" #include "ngp_algorithms/CourantReAlg.h" #include "ngp_algorithms/GeometryAlgDriver.h" #include "ngp_algorithms/MdotEdgeAlg.h" @@ -1050,6 +1051,8 @@ MomentumEquationSystem::MomentumEquationSystem(EquationSystems& eqSystems) tvisc_(NULL), evisc_(NULL), nodalGradAlgDriver_(realm_, "velocity", "dudx"), + nodalBuoyancyAlgDriver_( + realm_, "buoyancy_source", "buoyancy_source_weight"), wallFuncAlgDriver_(realm_), dynPressAlgDriver_(realm_), cflReAlgDriver_(realm_), @@ -1187,6 +1190,15 @@ MomentumEquationSystem::register_nodal_fields( &(meta_data.declare_field(stk::topology::NODE_RANK, "viscosity")); stk::mesh::put_field_on_mesh(*visc_, selector, nullptr); + if (realm_.solutionOptions_->use_balanced_buoyancy_force_) { + buoyancy_source_ = &(meta_data.declare_field( + stk::topology::NODE_RANK, "buoyancy_source")); + stk::mesh::put_field_on_mesh(*buoyancy_source_, selector, nDim, nullptr); + 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); + } + if (realm_.is_turbulent()) { tvisc_ = &(meta_data.declare_field( stk::topology::NODE_RANK, "turbulent_viscosity")); @@ -1255,7 +1267,7 @@ MomentumEquationSystem::register_nodal_fields( stk::topology::NODE_RANK, "abl_wall_no_slip_wall_func_node_mask"); double one = 1; stk::mesh::put_field_on_mesh(node_mask, selector, &one); -} +} // namespace nalu //-------------------------------------------------------------------------- //-------- register_element_fields ----------------------------------------- @@ -1313,6 +1325,12 @@ MomentumEquationSystem::register_interior_algorithm(stk::mesh::Part* part) edgeNodalGradient_); } + if (realm_.solutionOptions_->use_balanced_buoyancy_force_) { + nodalBuoyancyAlgDriver_.register_edge_algorithm( + algType, part, "momentum_buoyancy_source", buoyancy_source_, + buoyancy_source_weight_); + } + const auto theTurbModel = realm_.solutionOptions_->turbulenceModel_; // solver; interior contribution (advection + diffusion) [possible CMM time] @@ -1371,8 +1389,8 @@ MomentumEquationSystem::register_interior_algorithm(stk::mesh::Part* part) // Check if the user has requested CMM or LMM algorithms; if so, do not // include Nodal Mass algorithms - std::vector checkAlgNames = { - "momentum_time_derivative", "lumped_momentum_time_derivative"}; + std::vector checkAlgNames = {"momentum_time_derivative", + "lumped_momentum_time_derivative"}; bool elementMassAlg = supp_alg_is_requested(checkAlgNames); // solver; time contribution (lumped mass matrix) if (!elementMassAlg || nodal_src_is_requested()) { @@ -2517,6 +2535,9 @@ MomentumEquationSystem::compute_projected_nodal_gradient() if (!managePNG_) { const double timeA = -NaluEnv::self().nalu_time(); nodalGradAlgDriver_.execute(); + if (realm_.solutionOptions_->use_balanced_buoyancy_force_) { + nodalBuoyancyAlgDriver_.execute(); + } timerMisc_ += (NaluEnv::self().nalu_time() + timeA); } else { // this option is more complex... Rather than solving a nDim*nDim system, we diff --git a/src/SolutionOptions.C b/src/SolutionOptions.C index b108024148..092b1b0083 100644 --- a/src/SolutionOptions.C +++ b/src/SolutionOptions.C @@ -132,6 +132,10 @@ SolutionOptions::load(const YAML::Node& y_node) y_solution_options, "interp_rhou_together_for_mdot", mdotInterpRhoUTogether_, mdotInterpRhoUTogether_); + get_if_present( + y_solution_options, "use_balanced_buoyancy_force", + use_balanced_buoyancy_force_, use_balanced_buoyancy_force_); + // Solve for incompressible continuity get_if_present( y_solution_options, "solve_incompressible_continuity", diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C index 11b7778356..61c7610db1 100644 --- a/src/edge_kernels/ContinuityEdgeSolverAlg.C +++ b/src/edge_kernels/ContinuityEdgeSolverAlg.C @@ -55,6 +55,18 @@ ContinuityEdgeSolverAlg::execute() const DblType solveIncompressibleEqn = realm_.get_incompressible_solve(); const DblType om_solveIncompressibleEqn = 1.0 - solveIncompressibleEqn; + const bool add_balanced_forcing = + realm_.solutionOptions_->use_balanced_buoyancy_force_; + + const auto& solnOptsGravity = + realm_.solutionOptions_->get_gravity_vector(ndim); + DblType gravity[3] = {}; + if (add_balanced_forcing) { + for (int idim = 0; idim < ndim; ++idim) { + gravity[idim] = solnOptsGravity[idim]; + } + } + // STK stk::mesh::NgpField instances for capture by lambda const auto& fieldMgr = realm_.ngp_field_manager(); auto coordinates = fieldMgr.get_field(coordinates_); @@ -65,6 +77,11 @@ ContinuityEdgeSolverAlg::execute() auto udiag = fieldMgr.get_field(Udiag_); auto edgeAreaVec = fieldMgr.get_field(edgeAreaVec_); + auto source = !add_balanced_forcing + ? fieldMgr.get_field(Gpdx_) + : fieldMgr.get_field( + get_field_ordinal(realm_.meta_data(), "buoyancy_source")); + stk::mesh::NgpField edgeFaceVelMag; bool needs_gcl = false; if (realm_.has_mesh_deformation()) { @@ -81,6 +98,7 @@ ContinuityEdgeSolverAlg::execute() pressure.sync_to_device(); udiag.sync_to_device(); edgeAreaVec.sync_to_device(); + source.sync_to_device(); run_algorithm( realm_.bulk_data(), @@ -119,6 +137,12 @@ ContinuityEdgeSolverAlg::execute() const DblType inv_axdx = 1.0 / axdx; DblType tmdot = -projTimeScale * (pressureR - pressureL) * asq * inv_axdx; + + if (add_balanced_forcing) { + for (int d = 0; d < ndim; ++d) { + tmdot += projTimeScale * av[d] * gravity[d] * rhoIp; + } + } if (needs_gcl) { tmdot -= rhoIp * edgeFaceVelMag.get(edge, 0); } @@ -132,8 +156,12 @@ ContinuityEdgeSolverAlg::execute() densityL * velocity.get(nodeL, d)); const DblType ujIp = 0.5 * (velocity.get(nodeR, d) + velocity.get(nodeL, d)); - const DblType GjIp = + 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)); + } tmdot += (interpTogether * rhoUjIp + om_interpTogether * rhoIp * ujIp + GjIp) * av[d] - diff --git a/src/ngp_algorithms/BuoyancySourceAlg.C b/src/ngp_algorithms/BuoyancySourceAlg.C new file mode 100644 index 0000000000..3d787e77c8 --- /dev/null +++ b/src/ngp_algorithms/BuoyancySourceAlg.C @@ -0,0 +1,122 @@ +// 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. +// + +#include "ngp_algorithms/BuoyancySourceAlg.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "ngp_utils/NgpFieldOps.h" +#include "ngp_utils/NgpFieldManager.h" +#include "Realm.h" +#include "utils/StkHelpers.h" +#include "stk_mesh/base/NgpMesh.hpp" +#include "SolutionOptions.h" + +namespace sierra { +namespace nalu { + +BuoyancySourceAlg::BuoyancySourceAlg( + Realm& realm, + stk::mesh::Part* part, + VectorFieldType* source, + ScalarFieldType* source_weight) + : Algorithm(realm, part), + source_(source->mesh_meta_data_ordinal()), + source_weight_(source_weight->mesh_meta_data_ordinal()), + edgeAreaVec_(get_field_ordinal( + realm_.meta_data(), "edge_area_vector", stk::topology::EDGE_RANK)), + dualNodalVol_(get_field_ordinal(realm_.meta_data(), "dual_nodal_volume")), + coordinates_( + get_field_ordinal(realm_.meta_data(), realm.get_coordinates_name())), + density_( + get_field_ordinal(realm_.meta_data(), "density", stk::mesh::StateNP1)) +{ +} + +void +BuoyancySourceAlg::execute() +{ + using EntityInfoType = nalu_ngp::EntityInfo; + const auto& meshInfo = realm_.mesh_info(); + const auto& meta = meshInfo.meta(); + const int ndim = meta.spatial_dimension(); + const auto ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + const auto coordinates = fieldMgr.template get_field(coordinates_); + const auto density = fieldMgr.template get_field(density_); + const auto edgeAreaVec = fieldMgr.template get_field(edgeAreaVec_); + const auto dualVol = fieldMgr.template get_field(dualNodalVol_); + auto source = fieldMgr.template get_field(source_); + auto sourceweight = fieldMgr.template get_field(source_weight_); + const auto sourceOps = nalu_ngp::edge_nodal_field_updater(ngpMesh, source); + const auto sourceweightOps = + nalu_ngp::edge_nodal_field_updater(ngpMesh, sourceweight); + + const stk::mesh::Selector sel = + meta.locally_owned_part() & stk::mesh::selectUnion(partVec_); // & + //!(realm_.get_inactive_selector()); + + double gravity[3] = {0.0, 0.0, 0.0}; + + if (realm_.solutionOptions_->gravity_.size() >= ndim) + for (int idim = 0; idim < ndim; ++idim) + gravity[idim] = realm_.solutionOptions_->gravity_[idim]; + + source.sync_to_device(); + + const std::string algName = meta.get_fields()[source_]->name() + "_edge"; + nalu_ngp::run_edge_algorithm( + algName, ngpMesh, sel, KOKKOS_LAMBDA(const EntityInfoType& einfo) { + NALU_ALIGNED DblType av[NDimMax]; + + for (int d = 0; d < ndim; ++d) + av[d] = edgeAreaVec.get(einfo.meshIdx, d); + + const auto nodeL = ngpMesh.fast_mesh_index(einfo.entityNodes[0]); + const auto nodeR = ngpMesh.fast_mesh_index(einfo.entityNodes[1]); + + const DblType invVolL = 1.0 / dualVol.get(nodeL, 0); + const DblType invVolR = 1.0 / dualVol.get(nodeR, 0); + + const DblType rhoIp = + 0.5 * (density.get(nodeL, 0) + density.get(nodeR, 0)); + + DblType cc_face[3] = {0.0, 0.0, 0.0}; + for (int i = 0; i < ndim; ++i) + cc_face[i] = + 0.5 * (coordinates.get(nodeL, i) + coordinates.get(nodeR, i)); + + DblType weight_l = 0.0; + DblType weight_r = 0.0; + + for (int i = 0; i < ndim; ++i) { + weight_l += stk::math::pow( + stk::math::abs((cc_face[i] - coordinates.get(nodeL, i)) * av[i]) * + invVolL, + 2); + weight_r += stk::math::pow( + stk::math::abs((cc_face[i] - coordinates.get(nodeR, i)) * av[i]) * + invVolR, + 2); + } + + weight_l = 1.0; // stk::math::sqrt(weight_l); + weight_r = 1.0; // stk::math::sqrt(weight_r); + + for (int i = 0; i < ndim; ++i) { + sourceOps(einfo, 0, i) += weight_l * rhoIp * gravity[i]; + sourceOps(einfo, 1, i) += weight_r * rhoIp * gravity[i]; + } + sourceweightOps(einfo, 0, 0) += weight_l; + sourceweightOps(einfo, 1, 0) += weight_r; + }); + source.modify_on_device(); + sourceweight.modify_on_device(); +} + +} // namespace nalu +} // namespace sierra diff --git a/src/ngp_algorithms/CMakeLists.txt b/src/ngp_algorithms/CMakeLists.txt index 87ba9cec7f..2679affb3a 100644 --- a/src/ngp_algorithms/CMakeLists.txt +++ b/src/ngp_algorithms/CMakeLists.txt @@ -1,5 +1,6 @@ target_sources(nalu PRIVATE # Algorithms + ${CMAKE_CURRENT_SOURCE_DIR}/BuoyancySourceAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/CourantReAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/DynamicPressureOpenAlg.C ${CMAKE_CURRENT_SOURCE_DIR}/NodalGradEdgeAlg.C @@ -38,6 +39,7 @@ target_sources(nalu PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/NgpAlgDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/MdotAlgDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/NodalGradAlgDriver.C + ${CMAKE_CURRENT_SOURCE_DIR}/NodalBuoyancyAlgDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/TKEWallFuncAlgDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/GeometryAlgDriver.C ${CMAKE_CURRENT_SOURCE_DIR}/WallFricVelAlgDriver.C diff --git a/src/ngp_algorithms/MdotEdgeAlg.C b/src/ngp_algorithms/MdotEdgeAlg.C index 02a454a2e1..52253915fe 100644 --- a/src/ngp_algorithms/MdotEdgeAlg.C +++ b/src/ngp_algorithms/MdotEdgeAlg.C @@ -78,6 +78,23 @@ MdotEdgeAlg::execute() } auto mdot = fieldMgr.get_field(massFlowRate_); + const bool add_balanced_forcing = + realm_.solutionOptions_->use_balanced_buoyancy_force_; + + const auto& solnOptsGravity = + realm_.solutionOptions_->get_gravity_vector(ndim); + DblType gravity[3] = {}; + if (add_balanced_forcing) { + for (int idim = 0; idim < ndim; ++idim) { + gravity[idim] = solnOptsGravity[idim]; + } + } + + auto source = !add_balanced_forcing + ? fieldMgr.get_field(Gpdx_) + : fieldMgr.get_field( + get_field_ordinal(realm_.meta_data(), "buoyancy_source")); + mdot.clear_sync_state(); coordinates.sync_to_device(); velocity.sync_to_device(); @@ -86,6 +103,7 @@ MdotEdgeAlg::execute() pressure.sync_to_device(); udiag.sync_to_device(); edgeAreaVec.sync_to_device(); + source.sync_to_device(); const stk::mesh::Selector sel = meta.locally_owned_part() & stk::mesh::selectUnion(partVec_) & @@ -126,6 +144,13 @@ MdotEdgeAlg::execute() const DblType inv_axdx = 1.0 / axdx; DblType tmdot = -projTimeScale * (pressureR - pressureL) * asq * inv_axdx; + + if (add_balanced_forcing) { + for (int d = 0; d < ndim; ++d) { + tmdot += projTimeScale * av[d] * gravity[d] * rhoIp; + } + } + if (needs_gcl) { tmdot -= rhoIp * edgeFaceVelMag.get(einfo.meshIdx, 0); } @@ -138,8 +163,13 @@ MdotEdgeAlg::execute() densityL * velocity.get(nodeL, d)); const DblType ujIp = 0.5 * (velocity.get(nodeR, d) + velocity.get(nodeL, d)); - const DblType GjIp = + 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)); + } + tmdot += (interpTogether * rhoUjIp + om_interpTogether * rhoIp * ujIp + GjIp) * av[d] - diff --git a/src/ngp_algorithms/NodalBuoyancyAlgDriver.C b/src/ngp_algorithms/NodalBuoyancyAlgDriver.C new file mode 100644 index 0000000000..cf3e82941c --- /dev/null +++ b/src/ngp_algorithms/NodalBuoyancyAlgDriver.C @@ -0,0 +1,124 @@ +// 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. +// + +#include "ngp_algorithms/NodalBuoyancyAlgDriver.h" +#include "ngp_utils/NgpFieldUtils.h" +#include "ngp_utils/NgpLoopUtils.h" +#include "Realm.h" + +#include "stk_mesh/base/Field.hpp" +#include "stk_mesh/base/FieldParallel.hpp" +#include "stk_mesh/base/FieldBLAS.hpp" +#include "stk_mesh/base/MetaData.hpp" +#include "stk_mesh/base/NgpFieldParallel.hpp" + +namespace sierra { +namespace nalu { + +NodalBuoyancyAlgDriver::NodalBuoyancyAlgDriver( + Realm& realm, + const std::string& sourceName, + const std::string& sourceweightName) + : NgpAlgDriver(realm), + sourceName_(sourceName), + sourceweightName_(sourceweightName) +{ +} + +void +NodalBuoyancyAlgDriver::pre_work() +{ + const auto& meta = realm_.meta_data(); + + auto* source = + meta.template get_field(stk::topology::NODE_RANK, sourceName_); + + auto* sourceweight = meta.template get_field( + stk::topology::NODE_RANK, sourceweightName_); + + stk::mesh::field_fill(0.0, *source); + stk::mesh::field_fill(0.0, *sourceweight); + + const auto& meshInfo = realm_.mesh_info(); + const auto ngpMesh = meshInfo.ngp_mesh(); + const auto& fieldMgr = meshInfo.ngp_field_manager(); + auto ngpsource = + fieldMgr.template get_field(source->mesh_meta_data_ordinal()); + auto ngpsourceweight = + fieldMgr.template get_field(sourceweight->mesh_meta_data_ordinal()); + + ngpsource.set_all(ngpMesh, 0.0); + ngpsource.clear_sync_state(); + + ngpsourceweight.set_all(ngpMesh, 0.0); + ngpsourceweight.clear_sync_state(); +} + +void +NodalBuoyancyAlgDriver::post_work() +{ + // TODO: Revisit logic after STK updates to ngp parallel updates + const auto& meta = realm_.meta_data(); + const auto& bulk = realm_.bulk_data(); + const auto& meshInfo = realm_.mesh_info(); + + auto* sourceweight = meta.template get_field( + stk::topology::NODE_RANK, sourceweightName_); + auto& ngpsourceweight = nalu_ngp::get_ngp_field(meshInfo, sourceweightName_); + ngpsourceweight.sync_to_host(); + + auto* source = + meta.template get_field(stk::topology::NODE_RANK, sourceName_); + auto& ngpsource = nalu_ngp::get_ngp_field(meshInfo, sourceName_); + ngpsource.sync_to_host(); + + const std::vector fVec{&ngpsource, &ngpsourceweight}; + bool doFinalSyncToDevice = false; + stk::mesh::parallel_sum(bulk, fVec, doFinalSyncToDevice); + + const int dim2 = meta.spatial_dimension(); + + if (realm_.hasPeriodic_) { + realm_.periodic_field_update(source, dim2); + realm_.periodic_field_update(sourceweight, 1); + } + + ngpsource.modify_on_host(); + ngpsource.sync_to_device(); + + ngpsourceweight.modify_on_host(); + ngpsourceweight.sync_to_device(); + + // Divide by weight here + + using Traits = nalu_ngp::NGPMeshTraits<>; + using EntityInfoType = nalu_ngp::EntityInfo; + + const auto& ngpMesh = meshInfo.ngp_mesh(); + + stk::mesh::Selector sel = + (meta.locally_owned_part() | meta.globally_shared_part()) & + stk::mesh::selectField(*source); + + nalu_ngp::run_entity_algorithm( + "apply_weight_to_buoyancy", ngpMesh, stk::topology::NODE_RANK, sel, + KOKKOS_LAMBDA(const Traits::MeshIndex& mi) { + if (ngpsourceweight.get(mi, 0) > 1e-12) { + for (int idim = 0; idim < dim2; ++idim) { + ngpsource.get(mi, idim) = + ngpsource.get(mi, idim) / ngpsourceweight.get(mi, 0); + } + } + }); + + ngpsource.modify_on_device(); +} + +} // namespace nalu +} // namespace sierra diff --git a/src/node_kernels/MomentumBuoyancyNodeKernel.C b/src/node_kernels/MomentumBuoyancyNodeKernel.C index a45d7591db..19c7a4c723 100644 --- a/src/node_kernels/MomentumBuoyancyNodeKernel.C +++ b/src/node_kernels/MomentumBuoyancyNodeKernel.C @@ -23,12 +23,15 @@ MomentumBuoyancyNodeKernel::MomentumBuoyancyNodeKernel( const stk::mesh::BulkData& bulk, const SolutionOptions& solnOpts) : NGPNodeKernel(), nDim_(bulk.mesh_meta_data().spatial_dimension()), + use_balanced_buoyancy_(solnOpts.use_balanced_buoyancy_force_), rhoRef_(solnOpts.referenceDensity_) { const auto& meta = bulk.mesh_meta_data(); dualNodalVolumeID_ = get_field_ordinal(meta, "dual_nodal_volume"); densityNp1ID_ = get_field_ordinal(meta, "density"); + if (use_balanced_buoyancy_) + sourceID_ = get_field_ordinal(meta, "buoyancy_source"); const std::vector& solnOptsGravity = solnOpts.get_gravity_vector(nDim_); @@ -42,6 +45,7 @@ MomentumBuoyancyNodeKernel::setup(Realm& realm) const auto& fieldMgr = realm.ngp_field_manager(); dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); densityNp1_ = fieldMgr.get_field(densityNp1ID_); + source_ = fieldMgr.get_field(sourceID_); } KOKKOS_FUNCTION @@ -55,8 +59,14 @@ MomentumBuoyancyNodeKernel::execute( const NodeKernelTraits::DblType dualVolume = dualNodalVolume_.get(node, 0); const double fac = (rhoNp1 - rhoRef_) * dualVolume; - for (int i = 0; i < nDim_; ++i) { - rhs(i) += fac * gravity_[i]; + if (use_balanced_buoyancy_) { + for (int i = 0; i < nDim_; ++i) { + rhs(i) += source_.get(node, i) * dualVolume; + } + } else { + for (int i = 0; i < nDim_; ++i) { + rhs(i) += fac * gravity_[i]; + } } } From 97310fa2aa883fed84c521780f9a129099c7ebdd Mon Sep 17 00:00:00 2001 From: whorne Date: Tue, 6 Aug 2024 13:44:55 -0600 Subject: [PATCH 2/9] Fix weighting --- src/ngp_algorithms/BuoyancySourceAlg.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ngp_algorithms/BuoyancySourceAlg.C b/src/ngp_algorithms/BuoyancySourceAlg.C index 3d787e77c8..f7a277f93d 100644 --- a/src/ngp_algorithms/BuoyancySourceAlg.C +++ b/src/ngp_algorithms/BuoyancySourceAlg.C @@ -104,8 +104,8 @@ BuoyancySourceAlg::execute() 2); } - weight_l = 1.0; // stk::math::sqrt(weight_l); - weight_r = 1.0; // stk::math::sqrt(weight_r); + weight_l = stk::math::sqrt(weight_l); + weight_r = stk::math::sqrt(weight_r); for (int i = 0; i < ndim; ++i) { sourceOps(einfo, 0, i) += weight_l * rhoIp * gravity[i]; From 8eb1e55677f337131a156f1547840757ab8b11aa Mon Sep 17 00:00:00 2001 From: whorne Date: Tue, 6 Aug 2024 14:03:56 -0600 Subject: [PATCH 3/9] Formatting --- include/matrix_free/PolynomialOrders.h | 32 +++++++++++++------------- src/LowMachEquationSystem.C | 4 ++-- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/include/matrix_free/PolynomialOrders.h b/include/matrix_free/PolynomialOrders.h index 9abc85e513..140f8a6e89 100644 --- a/include/matrix_free/PolynomialOrders.h +++ b/include/matrix_free/PolynomialOrders.h @@ -56,27 +56,27 @@ enum { #define P_INVOKEABLE(func) \ template \ auto func(Args&&... args) \ - -> decltype(IMPLNAME(func) < p > ::invoke(std::forward(args)...)) \ + ->decltype(IMPLNAME(func) < p > ::invoke(std::forward(args)...)) \ { \ return IMPLNAME(func)

::invoke(std::forward(args)...); \ } // can't return a value dependent on template parameter -#define SWITCH_INVOKEABLE(func) \ - template \ - auto func(int p, Args&&... args) \ - -> decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward(args)...)) \ - { \ - switch (p) { \ - case inst::P2: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - case inst::P3: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - case inst::P4: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - default: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - } \ +#define SWITCH_INVOKEABLE(func) \ + template \ + auto func(int p, Args&&... args) \ + ->decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward(args)...)) \ + { \ + switch (p) { \ + case inst::P2: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + case inst::P3: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + case inst::P4: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + default: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + } \ } } // namespace matrix_free diff --git a/src/LowMachEquationSystem.C b/src/LowMachEquationSystem.C index 92f90328f3..4096224a04 100644 --- a/src/LowMachEquationSystem.C +++ b/src/LowMachEquationSystem.C @@ -1389,8 +1389,8 @@ MomentumEquationSystem::register_interior_algorithm(stk::mesh::Part* part) // Check if the user has requested CMM or LMM algorithms; if so, do not // include Nodal Mass algorithms - std::vector checkAlgNames = {"momentum_time_derivative", - "lumped_momentum_time_derivative"}; + std::vector checkAlgNames = { + "momentum_time_derivative", "lumped_momentum_time_derivative"}; bool elementMassAlg = supp_alg_is_requested(checkAlgNames); // solver; time contribution (lumped mass matrix) if (!elementMassAlg || nodal_src_is_requested()) { From d8aa4e209cc755f02d00801f1f579920dbaafc82 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Mon, 12 Aug 2024 09:00:19 -0600 Subject: [PATCH 4/9] git clang-format applied to PR commits --- include/matrix_free/PolynomialOrders.h | 32 +++++++++++++------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/include/matrix_free/PolynomialOrders.h b/include/matrix_free/PolynomialOrders.h index 140f8a6e89..9abc85e513 100644 --- a/include/matrix_free/PolynomialOrders.h +++ b/include/matrix_free/PolynomialOrders.h @@ -56,27 +56,27 @@ enum { #define P_INVOKEABLE(func) \ template \ auto func(Args&&... args) \ - ->decltype(IMPLNAME(func) < p > ::invoke(std::forward(args)...)) \ + -> decltype(IMPLNAME(func) < p > ::invoke(std::forward(args)...)) \ { \ return IMPLNAME(func)

::invoke(std::forward(args)...); \ } // can't return a value dependent on template parameter -#define SWITCH_INVOKEABLE(func) \ - template \ - auto func(int p, Args&&... args) \ - ->decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward(args)...)) \ - { \ - switch (p) { \ - case inst::P2: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - case inst::P3: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - case inst::P4: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - default: \ - return IMPLNAME(func)::invoke(std::forward(args)...); \ - } \ +#define SWITCH_INVOKEABLE(func) \ + template \ + auto func(int p, Args&&... args) \ + -> decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward(args)...)) \ + { \ + switch (p) { \ + case inst::P2: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + case inst::P3: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + case inst::P4: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + default: \ + return IMPLNAME(func)::invoke(std::forward(args)...); \ + } \ } } // namespace matrix_free From 1dc70705e916d2d93796a3e51ca8cff0f05b37c9 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Mon, 12 Aug 2024 09:04:15 -0600 Subject: [PATCH 5/9] remove comment --- src/ngp_algorithms/BuoyancySourceAlg.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/ngp_algorithms/BuoyancySourceAlg.C b/src/ngp_algorithms/BuoyancySourceAlg.C index f7a277f93d..9936817f3d 100644 --- a/src/ngp_algorithms/BuoyancySourceAlg.C +++ b/src/ngp_algorithms/BuoyancySourceAlg.C @@ -57,8 +57,7 @@ BuoyancySourceAlg::execute() nalu_ngp::edge_nodal_field_updater(ngpMesh, sourceweight); const stk::mesh::Selector sel = - meta.locally_owned_part() & stk::mesh::selectUnion(partVec_); // & - //!(realm_.get_inactive_selector()); + meta.locally_owned_part() & stk::mesh::selectUnion(partVec_); double gravity[3] = {0.0, 0.0, 0.0}; From dfcfb238809cc1ad2f05b8aeb9962aac1180a2ff Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Tue, 13 Aug 2024 10:50:06 -0600 Subject: [PATCH 6/9] putting things inside conditionals to avoid segfaults --- src/edge_kernels/ContinuityEdgeSolverAlg.C | 4 ++-- src/ngp_algorithms/MdotEdgeAlg.C | 4 ++-- src/node_kernels/MomentumBuoyancyNodeKernel.C | 3 ++- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C index 61c7610db1..c607dd1dbf 100644 --- a/src/edge_kernels/ContinuityEdgeSolverAlg.C +++ b/src/edge_kernels/ContinuityEdgeSolverAlg.C @@ -58,10 +58,10 @@ ContinuityEdgeSolverAlg::execute() const bool add_balanced_forcing = realm_.solutionOptions_->use_balanced_buoyancy_force_; - const auto& solnOptsGravity = - realm_.solutionOptions_->get_gravity_vector(ndim); DblType gravity[3] = {}; if (add_balanced_forcing) { + const auto& solnOptsGravity = + realm_.solutionOptions_->get_gravity_vector(ndim); for (int idim = 0; idim < ndim; ++idim) { gravity[idim] = solnOptsGravity[idim]; } diff --git a/src/ngp_algorithms/MdotEdgeAlg.C b/src/ngp_algorithms/MdotEdgeAlg.C index 52253915fe..9565c3e30c 100644 --- a/src/ngp_algorithms/MdotEdgeAlg.C +++ b/src/ngp_algorithms/MdotEdgeAlg.C @@ -81,10 +81,10 @@ MdotEdgeAlg::execute() const bool add_balanced_forcing = realm_.solutionOptions_->use_balanced_buoyancy_force_; - const auto& solnOptsGravity = - realm_.solutionOptions_->get_gravity_vector(ndim); DblType gravity[3] = {}; if (add_balanced_forcing) { + const auto& solnOptsGravity = + realm_.solutionOptions_->get_gravity_vector(ndim); for (int idim = 0; idim < ndim; ++idim) { gravity[idim] = solnOptsGravity[idim]; } diff --git a/src/node_kernels/MomentumBuoyancyNodeKernel.C b/src/node_kernels/MomentumBuoyancyNodeKernel.C index 19c7a4c723..b8787db600 100644 --- a/src/node_kernels/MomentumBuoyancyNodeKernel.C +++ b/src/node_kernels/MomentumBuoyancyNodeKernel.C @@ -45,7 +45,8 @@ MomentumBuoyancyNodeKernel::setup(Realm& realm) const auto& fieldMgr = realm.ngp_field_manager(); dualNodalVolume_ = fieldMgr.get_field(dualNodalVolumeID_); densityNp1_ = fieldMgr.get_field(densityNp1ID_); - source_ = fieldMgr.get_field(sourceID_); + if (use_balanced_buoyancy_) + source_ = fieldMgr.get_field(sourceID_); } KOKKOS_FUNCTION From a7b7d975b7997e65a59fead46827c6585ef7348b Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Tue, 13 Aug 2024 13:52:32 -0600 Subject: [PATCH 7/9] add balanced buoyancy to a reg test --- reg_tests/test_files/VOFDroplet/VOFDroplet.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml index 84905c5ced..0d11a3c8ff 100755 --- a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml +++ b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml @@ -106,6 +106,7 @@ realms: solution_options: name: myOptions + use_balanced_buoyancy_force: yes options: - hybrid_factor: From bc2e7ffc0b84730a76da79d012f5a12892c8b125 Mon Sep 17 00:00:00 2001 From: whorne Date: Thu, 15 Aug 2024 10:47:31 -0600 Subject: [PATCH 8/9] Add gravity-direction weighting instead of distance weighting --- src/ngp_algorithms/BuoyancySourceAlg.C | 33 ++++++-------------------- 1 file changed, 7 insertions(+), 26 deletions(-) diff --git a/src/ngp_algorithms/BuoyancySourceAlg.C b/src/ngp_algorithms/BuoyancySourceAlg.C index 9936817f3d..981a8c8086 100644 --- a/src/ngp_algorithms/BuoyancySourceAlg.C +++ b/src/ngp_algorithms/BuoyancySourceAlg.C @@ -46,10 +46,8 @@ BuoyancySourceAlg::execute() const int ndim = meta.spatial_dimension(); const auto ngpMesh = meshInfo.ngp_mesh(); const auto& fieldMgr = meshInfo.ngp_field_manager(); - const auto coordinates = fieldMgr.template get_field(coordinates_); const auto density = fieldMgr.template get_field(density_); const auto edgeAreaVec = fieldMgr.template get_field(edgeAreaVec_); - const auto dualVol = fieldMgr.template get_field(dualNodalVol_); auto source = fieldMgr.template get_field(source_); auto sourceweight = fieldMgr.template get_field(source_weight_); const auto sourceOps = nalu_ngp::edge_nodal_field_updater(ngpMesh, source); @@ -78,40 +76,23 @@ BuoyancySourceAlg::execute() const auto nodeL = ngpMesh.fast_mesh_index(einfo.entityNodes[0]); const auto nodeR = ngpMesh.fast_mesh_index(einfo.entityNodes[1]); - const DblType invVolL = 1.0 / dualVol.get(nodeL, 0); - const DblType invVolR = 1.0 / dualVol.get(nodeR, 0); - const DblType rhoIp = 0.5 * (density.get(nodeL, 0) + density.get(nodeR, 0)); - DblType cc_face[3] = {0.0, 0.0, 0.0}; - for (int i = 0; i < ndim; ++i) - cc_face[i] = - 0.5 * (coordinates.get(nodeL, i) + coordinates.get(nodeR, i)); - - DblType weight_l = 0.0; - DblType weight_r = 0.0; + DblType weight = 0.0; for (int i = 0; i < ndim; ++i) { - weight_l += stk::math::pow( - stk::math::abs((cc_face[i] - coordinates.get(nodeL, i)) * av[i]) * - invVolL, - 2); - weight_r += stk::math::pow( - stk::math::abs((cc_face[i] - coordinates.get(nodeR, i)) * av[i]) * - invVolR, - 2); + weight += stk::math::pow(gravity[i] * av[i], 2); } - weight_l = stk::math::sqrt(weight_l); - weight_r = stk::math::sqrt(weight_r); + weight = stk::math::sqrt(weight); for (int i = 0; i < ndim; ++i) { - sourceOps(einfo, 0, i) += weight_l * rhoIp * gravity[i]; - sourceOps(einfo, 1, i) += weight_r * rhoIp * gravity[i]; + sourceOps(einfo, 0, i) += weight * rhoIp * gravity[i]; + sourceOps(einfo, 1, i) += weight * rhoIp * gravity[i]; } - sourceweightOps(einfo, 0, 0) += weight_l; - sourceweightOps(einfo, 1, 0) += weight_r; + sourceweightOps(einfo, 0, 0) += weight; + sourceweightOps(einfo, 1, 0) += weight; }); source.modify_on_device(); sourceweight.modify_on_device(); From 8260c9120f330d5949c66ca838e189fbf910cbac Mon Sep 17 00:00:00 2001 From: whorne Date: Thu, 15 Aug 2024 12:10:34 -0600 Subject: [PATCH 9/9] 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 d0327ee260..254c8994cb 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 0000000000..f1c2bdc055 --- /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 4096224a04..f728d91f68 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 c607dd1dbf..da8b4ac673 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 9565c3e30c..a95305505c 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 +=