Merge pull request #2 from mbkuhn/newcases_whfork

Improvements for multiphase cases: sloshing tank pressure, water level density

include/user_functions/SloshingTankPressureAuxFunction.h

@@ -0,0 +1,48 @@
+// 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 SloshingTankPressureAuxFunction_h
+#define SloshingTankPressureAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class SloshingTankPressureAuxFunction : public AuxFunction
+{
+public:
+  SloshingTankPressureAuxFunction(const std::vector<double>& params);
+
+  virtual ~SloshingTankPressureAuxFunction() {}
+
+  using AuxFunction::do_evaluate;
+  virtual void do_evaluate(
+    const double* coords,
+    const double time,
+    const unsigned spatialDimension,
+    const unsigned numPoints,
+    double* fieldPtr,
+    const unsigned fieldSize,
+    const unsigned beginPos,
+    const unsigned endPos) const;
+
+  double water_level_;
+  double amplitude_;
+  double kappa_;
+  double interface_thickness_;
+  double domain_height_;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

include/user_functions/WaterLevelDensityAuxFunction.h

@@ -0,0 +1,47 @@
+// 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 WaterLevelDensityAuxFunction_h
+#define WaterLevelDensityAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class WaterLevelDensityAuxFunction : public AuxFunction
+{
+public:
+  WaterLevelDensityAuxFunction(const std::vector<double>& params);
+
+  virtual ~WaterLevelDensityAuxFunction() {}
+
+  using AuxFunction::do_evaluate;
+  virtual void do_evaluate(
+    const double* coords,
+    const double time,
+    const unsigned spatialDimension,
+    const unsigned numPoints,
+    double* fieldPtr,
+    const unsigned fieldSize,
+    const unsigned beginPos,
+    const unsigned endPos) const;
+
+  double water_level_;
+  double interface_thickness_;
+  double liquid_density_;
+  double gas_density_;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

src/LowMachEquationSystem.C

@@ -178,6 +178,9 @@
 
 #include <user_functions/DropletVelocityAuxFunction.h>
 
+#include <user_functions/SloshingTankPressureAuxFunction.h>
+#include <user_functions/WaterLevelDensityAuxFunction.h>
+
 // deprecated
 
 // stk_util
@@ -708,6 +711,13 @@ LowMachEquationSystem::register_initial_condition_fcn(
     AuxFunctionAlgorithm* auxAlg = NULL;
     if (fcnName == "flat_interface") {
       theAuxFunc = new FlatDensityAuxFunction();
+    } else if (fcnName == "water_level") {
+      std::map<std::string, std::vector<double>>::const_iterator iterParams =
+        theParams.find(dofName_init_dens);
+      std::vector<double> fcnParams = (iterParams != theParams.end())
+                                        ? (*iterParams).second
+                                        : std::vector<double>();
+      theAuxFunc = new WaterLevelDensityAuxFunction(fcnParams);
     } else {
       throw std::runtime_error(
         "InitialCondFunction::non-supported initial_density IC");
@@ -3547,6 +3557,13 @@ ContinuityEquationSystem::register_initial_condition_fcn(
       theAuxFunc = new TaylorGreenPressureAuxFunction();
     } else if (fcnName == "kovasznay") {
       theAuxFunc = new KovasznayPressureAuxFunction();
+    } else if (fcnName == "sloshing_tank") {
+      std::map<std::string, std::vector<double>>::const_iterator iterParams =
+        theParams.find(dofName);
+      std::vector<double> fcnParams = (iterParams != theParams.end())
+                                        ? (*iterParams).second
+                                        : std::vector<double>();
+      theAuxFunc = new SloshingTankPressureAuxFunction(fcnParams);
     } else {
       throw std::runtime_error("ContinuityEquationSystem::register_initial_"
                                "condition_fcn: limited functions supported");

src/user_functions/CMakeLists.txt

@@ -43,4 +43,6 @@ target_sources(nalu PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/DropletVelocityAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/FlatDensityAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankVOFAuxFunction.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankPressureAuxFunction.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/WaterLevelDensityAuxFunction.C
 )

src/user_functions/SloshingTankPressureAuxFunction.C

@@ -0,0 +1,98 @@
+// 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 <user_functions/SloshingTankPressureAuxFunction.h>
+#include <algorithm>
+
+// basic c++
+#include <cmath>
+#include <vector>
+#include <stdexcept>
+#include <iostream>
+
+namespace sierra {
+namespace nalu {
+
+SloshingTankPressureAuxFunction::SloshingTankPressureAuxFunction(
+  const std::vector<double>& params)
+  : AuxFunction(0, 1),
+    water_level_(0.),
+    amplitude_(0.1),
+    kappa_(0.25),
+    interface_thickness_(0.1),
+    domain_height_(0.5)
+{
+  // check size and populate
+  if (params.size() != 5 && !params.empty())
+    throw std::runtime_error(
+      "Realm::setup_initial_conditions: sloshing_tank requires 5 params: water "
+      "level, amplitude, kappa, interface thickness, and domain height");
+  if (!params.empty()) {
+    water_level_ = params[0];
+    amplitude_ = params[1];
+    kappa_ = params[2];
+    interface_thickness_ = params[3];
+    domain_height_ = params[4];
+  }
+}
+
+void
+SloshingTankPressureAuxFunction::do_evaluate(
+  const double* coords,
+  const double /*time*/,
+  const unsigned spatialDimension,
+  const unsigned numPoints,
+  double* fieldPtr,
+  const unsigned fieldSize,
+  const unsigned /*beginPos*/,
+  const unsigned /*endPos*/) const
+{
+  for (unsigned p = 0; p < numPoints; ++p) {
+
+    const double x = coords[0];
+    const double y = coords[1];
+    const double z0 = coords[2];
+    const double z1 = domain_height_;
+
+    // These need to come from elsewhere
+    const double g = 9.81;
+
+    const double z_init =
+      water_level_ + amplitude_ * std::exp(-kappa_ * (x * x + y * y));
+
+    // Position 1 is top, 0 is current
+    // Density is rho_l * VOF + rho_g * (1-VOF)
+    // Integral of error function is x * erf(x) + e^(-x^2)/sqrt(pi)
+    // Integrate by substitution, where x = (z - z0) / interface_thickness_
+    const double x0 = (z0 - z_init) / interface_thickness_;
+    const double x1 = (z1 - z_init) / interface_thickness_;
+    const double int_factor = interface_thickness_;
+    const double int_erf0 =
+      (x0 * std::erf(x0) + std::exp(-std::pow(x0, 2)) / std::sqrt(M_PI)) *
+      int_factor;
+    const double int_erf1 =
+      (x1 * std::erf(x1) + std::exp(-std::pow(x1, 2)) / std::sqrt(M_PI)) *
+      int_factor;
+
+    const double int_vof0 = 0.5 * z0 - 0.5 * int_erf0;
+    const double int_vof1 = 0.5 * z1 - 0.5 * int_erf1;
+    const double int_rho0 = 1000.0 * int_vof0 + 1.0 * (z0 - int_vof0);
+    const double int_rho1 = 1000.0 * int_vof1 + 1.0 * (z1 - int_vof1);
+    // vof_local = -0.5 * (std::erf((z - z_init) / interface_thickness_) + 1.0) + 1.0;
+
+    // g * integral(rho)dz
+    fieldPtr[0] = g * (int_rho1 - int_rho0);
+
+    fieldPtr += fieldSize;
+    coords += spatialDimension;
+  }
+}
+
+} // namespace nalu
+} // namespace sierra

src/user_functions/WaterLevelDensityAuxFunction.C

@@ -0,0 +1,69 @@
+// 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 <user_functions/WaterLevelDensityAuxFunction.h>
+#include <algorithm>
+
+// basic c++
+#include <cmath>
+#include <vector>
+#include <stdexcept>
+#include <iostream>
+
+namespace sierra {
+namespace nalu {
+
+WaterLevelDensityAuxFunction::WaterLevelDensityAuxFunction(
+  const std::vector<double>& params)
+  : AuxFunction(0, 1),
+    water_level_(0.),
+    interface_thickness_(0.1),
+    liquid_density_(1000.0),
+    gas_density_(1.0)
+{
+  // check size and populate
+  if (params.size() != 4 && !params.empty())
+    throw std::runtime_error("Realm::setup_initial_conditions: "
+                             "water_level requires 4 params: water level, "
+                             "amplitude, kappa, and interface thickness");
+  if (!params.empty()) {
+    water_level_ = params[0];
+    interface_thickness_ = params[1];
+    liquid_density_ = params[2];
+    gas_density_ = params[3];
+  }
+}
+
+void
+WaterLevelDensityAuxFunction::do_evaluate(
+  const double* coords,
+  const double /*time*/,
+  const unsigned spatialDimension,
+  const unsigned numPoints,
+  double* fieldPtr,
+  const unsigned fieldSize,
+  const unsigned /*beginPos*/,
+  const unsigned /*endPos*/) const
+{
+  for (unsigned p = 0; p < numPoints; ++p) {
+
+    const double z = coords[2];
+    const double z0 = water_level_;
+
+    const double local_vof =
+      -0.5 * (std::erf((z - z0) / interface_thickness_) + 1.0) + 1.0;
+    fieldPtr[0] = liquid_density_ * local_vof + gas_density_ * (1.0 - local_vof);
+
+    fieldPtr += fieldSize;
+    coords += spatialDimension;
+  }
+}
+
+} // namespace nalu
+} // namespace sierra