Mask balanced buoyancy terms from walls

include/LowMachEquationSystem.h

@@ -216,6 +216,7 @@ class MomentumEquationSystem : public EquationSystem
   std::unique_ptr<Algorithm> tviscAlg_{nullptr};
   std::unique_ptr<Algorithm> pecletAlg_{nullptr};
   std::unique_ptr<Algorithm> ablWallNodeMask_{nullptr};
+  std::unique_ptr<Algorithm> buoyancySrcMask_{nullptr};
 
   CourantReAlgDriver cflReAlgDriver_;
   std::unique_ptr<AMSAlgDriver> AMSAlgDriver_{nullptr};

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 <Realm.h>
+#include <Algorithm.h>
+#include <stk_mesh/base/Types.hpp>
+#include "stk_mesh/base/NgpField.hpp"
+#include <ngp_utils/NgpLoopUtils.h>
+#include <ngp_utils/NgpTypes.h>
+#include <stk_mesh/base/Part.hpp>
+#include <stk_mesh/base/Selector.hpp>
+#include <ngp_utils/NgpFieldUtils.h>
+
+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<stk::mesh::NgpMesh>;
+
+    const auto& meta = realm_.meta_data();
+    const auto ngpMesh = realm_.ngp_mesh();
+    const auto& fieldMgr = realm_.ngp_field_manager();
+    auto myNodeMask = fieldMgr.get_field<double>(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

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<double>(
       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<double>(
+      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;

src/edge_kernels/ContinuityEdgeSolverAlg.C

@@ -82,6 +82,11 @@ ContinuityEdgeSolverAlg::execute()
                   : fieldMgr.get_field<double>(
                       get_field_ordinal(realm_.meta_data(), "buoyancy_source"));
 
+  auto source_mask = !add_balanced_forcing
+                       ? fieldMgr.get_field<double>(densityNp1_)
+                       : fieldMgr.get_field<double>(get_field_ordinal(
+                           realm_.meta_data(), "buoyancy_source_mask"));
+
   stk::mesh::NgpField<double> 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) *

src/ngp_algorithms/MdotEdgeAlg.C

@@ -94,6 +94,10 @@ MdotEdgeAlg::execute()
                   ? fieldMgr.get_field<double>(Gpdx_)
                   : fieldMgr.get_field<double>(
                       get_field_ordinal(realm_.meta_data(), "buoyancy_source"));
+  auto source_mask = !add_balanced_forcing
+                       ? fieldMgr.get_field<double>(densityNp1_)
+                       : fieldMgr.get_field<double>(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 +=