Change initial density to a specified function and introduce wall bcs…

… for VOF

include/user_functions/FlatDensityAuxFunction.h

@@ -0,0 +1,42 @@
+// 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 FlatDensityAuxFunction_h
+#define FlatDensityAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra {
+namespace nalu {
+
+class FlatDensityAuxFunction : public AuxFunction
+{
+public:
+  FlatDensityAuxFunction();
+
+  virtual ~FlatDensityAuxFunction() {}
+
+  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;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif

src/LowMachEquationSystem.C

@@ -166,6 +166,7 @@
 #include <user_functions/KovasznayPressureAuxFunction.h>
 
 #include <user_functions/DropletVelocityAuxFunction.h>
+#include <user_functions/FlatDensityAuxFunction.h>
 
 #include <overset/UpdateOversetFringeAlgorithmDriver.h>
 #include <overset/AssembleOversetPressureAlgorithm.h>
@@ -309,7 +310,7 @@ LowMachEquationSystem::register_nodal_fields(
   density_ = &(meta_data.declare_field<ScalarFieldType>(
     stk::topology::NODE_RANK, "density", numStates));
   initial_density_ = &(meta_data.declare_field<ScalarFieldType>(
-    stk::topology::NODE_RANK, "initial_density", numStates));
+    stk::topology::NODE_RANK, "initial_density"));
 
   stk::mesh::put_field_on_mesh(*initial_density_, selector, nullptr);
   realm_.augment_restart_variable_list("initial_density");
@@ -334,18 +335,18 @@ LowMachEquationSystem::register_nodal_fields(
   if (numVolStates > 1)
     realm_.augment_restart_variable_list("dual_nodal_volume");
 
-  // make sure all states are properly populated (restart can handle this)
-
-  if (!realm_.restarted_simulation() || realm_.support_inconsistent_restart()) {
+ /*if (
+    numStates > 2 &&
+    (!realm_.restarted_simulation() || realm_.support_inconsistent_restart())) {
+    ScalarFieldType& densityN = initial_density_->field_of_state(stk::mesh::StateN);
+    ScalarFieldType& densityNp1 = initial_density_->field_of_state(stk::mesh::StateNP1);
 
-    ScalarFieldType& densityNp1 = density_->field_of_state(stk::mesh::StateNP1);
-    ScalarFieldType& initDensityNp1 = initial_density_->field_of_state(stk::mesh::StateNP1);
-
     CopyFieldAlgorithm* theCopyAlgDens = new CopyFieldAlgorithm(
-      realm_, part_vec, &densityNp1, &initDensityNp1, 0, 1, stk::topology::NODE_RANK);
+      realm_, part_vec, &densityNp1, &densityN, 0, 1, stk::topology::NODE_RANK);
     copyStateAlg_.push_back(theCopyAlgDens);
 
-  }
+
+  }*/
 
   if (
     numStates > 2 &&
@@ -690,6 +691,35 @@ LowMachEquationSystem::register_initial_condition_fcn(
     // push to ic
     realm_.initCondAlg_.push_back(auxAlg);
   }
+
+  // iterate map and check for name
+  const std::string dofName_init_dens = "initial_density";
+  std::map<std::string, std::string>::const_iterator iterName_init_dens =
+    theNames.find(dofName_init_dens);
+  if (iterName_init_dens != theNames.end()) {
+    std::string fcnName = (*iterName_init_dens).second;
+    // save off the field (np1 state)
+    ScalarFieldType* initDensNp1 = meta_data.get_field<ScalarFieldType>(
+      stk::topology::NODE_RANK, "initial_density");
+
+
+    // create a few Aux things
+    AuxFunction* theAuxFunc = NULL;
+    AuxFunctionAlgorithm* auxAlg = NULL;
+    if (fcnName == "flat_interface") {
+      theAuxFunc = new FlatDensityAuxFunction();
+    } else {
+      throw std::runtime_error(
+        "InitialCondFunction::non-supported initial_density IC");
+    }
+    // create the algorithm
+    auxAlg = new AuxFunctionAlgorithm(
+      realm_, part, initDensNp1, theAuxFunc, stk::topology::NODE_RANK);
+
+    // push to ic
+    realm_.initCondAlg_.push_back(auxAlg);
+  }
+
 }
 
 void

src/MomentumBuoyancySrcNodeSuppAlg.C

@@ -45,8 +45,8 @@ MomentumBuoyancySrcNodeSuppAlg::MomentumBuoyancySrcNodeSuppAlg(Realm& realm)
   ScalarFieldType* initDensity =
     meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "initial_density");
 
-  densityNp1_ = &(density->field_of_state(stk::mesh::StateNP1));
-  initDensityNp1_ = &(initDensity->field_of_state(stk::mesh::StateNP1));
+  densityNp1_ = &(density->field_of_state(stk::mesh::StateN));
+  initDensityNp1_ = &(initDensity->field_of_state(stk::mesh::StateNone));
 
   dualNodalVolume_ = meta_data.get_field<ScalarFieldType>(
     stk::topology::NODE_RANK, "dual_nodal_volume");
@@ -58,6 +58,7 @@ MomentumBuoyancySrcNodeSuppAlg::MomentumBuoyancySrcNodeSuppAlg(Realm& realm)
   rhoRef_ = realm_.solutionOptions_->referenceDensity_;
 
   rhoRefIsInitDens_ = realm_.solutionOptions_->rho_ref_to_initial_rho_;
+
 }
 
 //--------------------------------------------------------------------------
@@ -81,6 +82,7 @@ MomentumBuoyancySrcNodeSuppAlg::node_execute(
   const double rhoNp1 = *stk::mesh::field_data(*densityNp1_, node);
   const double initRhoNp1 = *stk::mesh::field_data(*initDensityNp1_, node);
   const double dualVolume = *stk::mesh::field_data(*dualNodalVolume_, node);
+
   const double fac = !rhoRefIsInitDens_ ? (rhoNp1 - rhoRef_) * dualVolume :
     (rhoNp1 - initRhoNp1) * dualVolume;
   const int nDim = nDim_;

src/VolumeOfFluidEquationSystem.C

@@ -272,7 +272,7 @@ VolumeOfFluidEquationSystem::register_inflow_bc(
 
   stk::mesh::MetaData& meta_data = realm_.meta_data();
 
-  // register boundary data; gamma_bc
+  // register boundary data; vof_bc
   ScalarFieldType* theBcField = &(meta_data.declare_field<ScalarFieldType>(
     stk::topology::NODE_RANK, "vof_bc"));
   stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr);
@@ -311,12 +311,12 @@ VolumeOfFluidEquationSystem::register_inflow_bc(
     bcDataAlg_.push_back(auxAlg);
   }
 
-  // copy vof_bc to gamma_transition np1...
+  // copy vof_bc to vof_transition np1...
   CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm(
     realm_, part, theBcField, &vofNp1, 0, 1, stk::topology::NODE_RANK);
   bcDataMapAlg_.push_back(theCopyAlg);
 
-  // non-solver; dgamdx; allow for element-based shifted
+  // non-solver; dvofdx; allow for element-based shifted
   nodalGradAlgDriver_.register_face_algorithm<ScalarNodalGradBndryElemAlg>(
     algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_);
 
@@ -349,7 +349,7 @@ VolumeOfFluidEquationSystem::register_open_bc(
 
   // non-solver; dvofdx; allow for element-based shifted
   nodalGradAlgDriver_.register_face_algorithm<ScalarNodalGradBndryElemAlg>(
-    algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone,
+    algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone,
     edgeNodalGradient_);
 }
 
@@ -360,7 +360,7 @@ void
 VolumeOfFluidEquationSystem::register_wall_bc(
   stk::mesh::Part* part,
   const stk::topology& /*theTopo*/,
-  const WallBoundaryConditionData& /* wallBCData */)
+  const WallBoundaryConditionData& wallBCData )
 {
   // algorithm type
   const AlgorithmType algType = WALL;
@@ -369,10 +369,67 @@ VolumeOfFluidEquationSystem::register_wall_bc(
   VectorFieldType& dvofdxNone =
     dvolumeOfFluiddx_->field_of_state(stk::mesh::StateNone);
 
+  stk::mesh::MetaData& meta_data = realm_.meta_data();
+
+  // register boundary data; vof_bc
+  ScalarFieldType* theBcField = &(meta_data.declare_field<ScalarFieldType>(
+    stk::topology::NODE_RANK, "vof_bc"));
+  stk::mesh::put_field_on_mesh(*theBcField, *part, nullptr);
+
+  // extract the value for user specified tke and save off the AuxFunction
+  WallUserData userData = wallBCData.userData_;
+  std::string vofName = "volume_of_fluid";
+  UserDataType theDataType = get_bc_data_type(userData, vofName);
+
+  AuxFunction* theAuxFunc = NULL;
+
+  if (CONSTANT_UD == theDataType) {
+    VolumeOfFluid volumeOfFluid = userData.volumeOfFluid_;
+    std::vector<double> userSpec(1);
+    userSpec[0] = volumeOfFluid.volumeOfFluid_;
+    theAuxFunc = new ConstantAuxFunction(0, 1, userSpec);
+
+  } else if (FUNCTION_UD == theDataType) {
+    throw std::runtime_error("VolumeOfFluidEquationSystem::register_wall_bc: "
+                             "limited functions supported");
+  } else {
+    throw std::runtime_error("VolumeOfFluidEquationSystem::register_wall_bc: "
+                             "only constant functions supported");
+  }
+
+  // bc data alg
+  AuxFunctionAlgorithm* auxAlg = new AuxFunctionAlgorithm(
+    realm_, part, theBcField, theAuxFunc, stk::topology::NODE_RANK);
+
+  // how to populate the field?
+  if (userData.externalData_) {
+    // xfer will handle population; only need to populate the initial value
+    realm_.initCondAlg_.push_back(auxAlg);
+  } else {
+    // put it on bcData
+    bcDataAlg_.push_back(auxAlg);
+  }
+
+  // copy vof_bc to vof_transition np1...
+  CopyFieldAlgorithm* theCopyAlg = new CopyFieldAlgorithm(
+    realm_, part, theBcField, &vofNp1, 0, 1, stk::topology::NODE_RANK);
+  bcDataMapAlg_.push_back(theCopyAlg);
+
   // non-solver; dvofdx; allow for element-based shifted
   nodalGradAlgDriver_.register_face_algorithm<ScalarNodalGradBndryElemAlg>(
-    algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone,
-    edgeNodalGradient_);
+    algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone, edgeNodalGradient_);
+
+  // Dirichlet bc
+  std::map<AlgorithmType, SolverAlgorithm*>::iterator itd =
+    solverAlgDriver_->solverDirichAlgMap_.find(algType);
+  if (itd == solverAlgDriver_->solverDirichAlgMap_.end()) {
+    DirichletBC* theAlg =
+      new DirichletBC(realm_, this, part, &vofNp1, theBcField, 0, 1);
+    solverAlgDriver_->solverDirichAlgMap_[algType] = theAlg;
+  } else {
+    itd->second->partVec_.push_back(part);
+  }
+
 }
 
 //--------------------------------------------------------------------------
@@ -391,9 +448,14 @@ VolumeOfFluidEquationSystem::register_symmetry_bc(
   VectorFieldType& dvofdxNone =
     dvolumeOfFluiddx_->field_of_state(stk::mesh::StateNone);
 
+  // extract the value for user specified volume_of_fluid and save off the
+  // AuxFunction
+  AuxFunction* theAuxFunc = NULL;
+  Algorithm* auxAlg = NULL;
+
   // non-solver; dvofdx; allow for element-based shifted
   nodalGradAlgDriver_.register_face_algorithm<ScalarNodalGradBndryElemAlg>(
-    algType, part, "gamma_nodal_grad", &vofNp1, &dvofdxNone,
+    algType, part, "vof_nodal_grad", &vofNp1, &dvofdxNone,
     edgeNodalGradient_);
 }
 

src/edge_kernels/MomentumEdgeSolverAlg.C

@@ -221,8 +221,10 @@ MomentumEdgeSolverAlg::execute()
           diff_flux += duidxj[j][j];
         diff_flux *= 2.0 / 3.0 * viscIp * av[i] * includeDivU;
 
-        for (int j = 0; j < ndim; ++j)
+        for (int j = 0; j < ndim; ++j) {
           diff_flux += -viscIp * (duidxj[i][j] + duidxj[j][i]) * av[j];
+        }
+
 
         const DblType maskNode = stk::math::min(
           maskNodeField.get(nodeL, 0), maskNodeField.get(nodeR, 0));

src/edge_kernels/VOFAdvectionEdgeAlg.C

@@ -77,7 +77,7 @@ VOFAdvectionEdgeAlg::VOFAdvectionEdgeAlg(
 void
 VOFAdvectionEdgeAlg::execute()
 {
-  const double eps = 1.0e-16;
+  const double eps = 1.0e-11;
   const double gradient_eps = 1.0e-9;
 
   const int ndim = realm_.meta_data().spatial_dimension();
@@ -144,7 +144,6 @@ VOFAdvectionEdgeAlg::execute()
 
       NALU_ALIGNED DblType densityL = density.get(nodeL, 0);
       NALU_ALIGNED DblType densityR = density.get(nodeR, 0);
-
       const DblType rhoIp = 2.0 / ( 1.0 / densityL + 1.0 / densityR);
 
       const DblType mdot = massFlowRate.get(edge, 0) / rhoIp;

src/user_functions/CMakeLists.txt

@@ -41,5 +41,6 @@ target_sources(nalu PRIVATE
    ${CMAKE_CURRENT_SOURCE_DIR}/ZalesakSphereMassFlowRateKernel.C
    ${CMAKE_CURRENT_SOURCE_DIR}/DropletVOFAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/DropletVelocityAuxFunction.C
+   ${CMAKE_CURRENT_SOURCE_DIR}/FlatDensityAuxFunction.C
    ${CMAKE_CURRENT_SOURCE_DIR}/SloshingTankVOFAuxFunction.C
 )

src/user_functions/DropletVOFAuxFunction.C

@@ -42,10 +42,10 @@ DropletVOFAuxFunction::do_evaluate(
     const double interface_thickness = 0.015;
 
     fieldPtr[0] = 0.0;
-    //fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
+    fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
 
     auto radius = std::sqrt(x * x + y * y + z * z) - 0.075;
-    fieldPtr[0] += -0.5 * (std::erf(radius / interface_thickness) + 1.0) + 1.0;
+    //fieldPtr[0] += -0.5 * (std::erf(radius / interface_thickness) + 1.0) + 1.0;
 
     fieldPtr += fieldSize;
     coords += spatialDimension;

src/user_functions/FlatDensityAuxFunction.C

@@ -0,0 +1,57 @@
+// 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/FlatDensityAuxFunction.h>
+#include <algorithm>
+
+// basic c++
+#include <cmath>
+#include <vector>
+#include <stdexcept>
+#include <iostream>
+
+namespace sierra {
+namespace nalu {
+
+FlatDensityAuxFunction::FlatDensityAuxFunction() : AuxFunction(0, 1)
+{
+  // does nothing
+}
+
+void
+FlatDensityAuxFunction::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 z = coords[2];
+    const double interface_thickness = 0.015;
+
+    fieldPtr[0] = 0.0;
+    fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
+
+    // air-water
+    fieldPtr[0] = 1000.0 * fieldPtr[0] + (1.0 - fieldPtr[0]);
+
+    fieldPtr += fieldSize;
+    coords += spatialDimension;
+  }
+}
+
+} // namespace nalu
+} // namespace sierra