Balanced buoyancy

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<EffDiffFluxCoeffAlg> diffFluxCoeffAlg_{nullptr};

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_;

include/matrix_free/PolynomialOrders.h

@@ -56,27 +56,27 @@ enum {
 #define P_INVOKEABLE(func)                                                     \
   template <int p, typename... Args>                                           \
   auto func(Args&&... args)                                                    \
-    -> decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))    \
+    ->decltype(IMPLNAME(func) < p > ::invoke(std::forward<Args>(args)...))     \
   {                                                                            \
     return IMPLNAME(func)<p>::invoke(std::forward<Args>(args)...);             \
   }
 
 // can't return a value dependent on template parameter
-#define SWITCH_INVOKEABLE(func)                                                    \
-  template <typename... Args>                                                      \
-  auto func(int p, Args&&... args)                                                 \
-    -> decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
-  {                                                                                \
-    switch (p) {                                                                   \
-    case inst::P2:                                                                 \
-      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);        \
-    case inst::P3:                                                                 \
-      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);        \
-    case inst::P4:                                                                 \
-      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);        \
-    default:                                                                       \
-      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);        \
-    }                                                                              \
+#define SWITCH_INVOKEABLE(func)                                                   \
+  template <typename... Args>                                                     \
+  auto func(int p, Args&&... args)                                                \
+    ->decltype(IMPLNAME(func) < inst::P1 > ::invoke(std::forward<Args>(args)...)) \
+  {                                                                               \
+    switch (p) {                                                                  \
+    case inst::P2:                                                                \
+      return IMPLNAME(func)<inst::P2>::invoke(std::forward<Args>(args)...);       \
+    case inst::P3:                                                                \
+      return IMPLNAME(func)<inst::P3>::invoke(std::forward<Args>(args)...);       \
+    case inst::P4:                                                                \
+      return IMPLNAME(func)<inst::P4>::invoke(std::forward<Args>(args)...);       \
+    default:                                                                      \
+      return IMPLNAME(func)<inst::P1>::invoke(std::forward<Args>(args)...);       \
+    }                                                                             \
   }
 
 } // namespace matrix_free

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 */

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 */

include/node_kernels/MomentumBuoyancyNodeKernel.h

@@ -45,11 +45,15 @@ class MomentumBuoyancyNodeKernel
 private:
   stk::mesh::NgpField<double> dualNodalVolume_;
   stk::mesh::NgpField<double> densityNp1_;
+  stk::mesh::NgpField<double> 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];
 };

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<double>(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<double>(
+      stk::topology::NODE_RANK, "buoyancy_source"));
+    stk::mesh::put_field_on_mesh(*buoyancy_source_, selector, nDim, nullptr);
+    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);
+  }
+
   if (realm_.is_turbulent()) {
     tvisc_ = &(meta_data.declare_field<double>(
       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<BuoyancySourceAlg>(
+      algType, part, "momentum_buoyancy_source", buoyancy_source_,
+      buoyancy_source_weight_);
+  }
+
   const auto theTurbModel = realm_.solutionOptions_->turbulenceModel_;
 
   // solver; interior contribution (advection + diffusion) [possible CMM time]
@@ -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

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",

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<double>(coordinates_);
@@ -65,6 +77,11 @@ ContinuityEdgeSolverAlg::execute()
   auto udiag = fieldMgr.get_field<double>(Udiag_);
   auto edgeAreaVec = fieldMgr.get_field<double>(edgeAreaVec_);
 
+  auto source = !add_balanced_forcing
+                  ? fieldMgr.get_field<double>(Gpdx_)
+                  : fieldMgr.get_field<double>(
+                      get_field_ordinal(realm_.meta_data(), "buoyancy_source"));
+
   stk::mesh::NgpField<double> 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] -