Subtract out normal pressure gradient when mdot correction BC is acti…

…ve (#4)
Exawind · May 30, 2018 · 8f5162d · 8f5162d
1 parent a529c52
commit 8f5162d
Show file tree

Hide file tree

Showing 13 changed files with 432 additions and 38 deletions.
diff --git a/include/AssembleNodalGradPAlgorithmDriver.h b/include/AssembleNodalGradPAlgorithmDriver.h
@@ -0,0 +1,36 @@
+/*------------------------------------------------------------------------*/
+/*  Copyright 2014 Sandia Corporation.                                    */
+/*  This software is released under the license detailed                  */
+/*  in the file, LICENSE, which is located in the top-level Nalu          */
+/*  directory structure                                                   */
+/*------------------------------------------------------------------------*/
+
+
+#ifndef AssembleNodalGradPAlgorithmDriver_h
+#define AssembleNodalGradPAlgorithmDriver_h
+
+#include<AlgorithmDriver.h>
+#include <FieldTypeDef.h>
+
+namespace sierra{
+namespace nalu{
+
+class Realm;
+
+class AssembleNodalGradPAlgorithmDriver : public AlgorithmDriver
+{
+public:
+
+  AssembleNodalGradPAlgorithmDriver(
+    Realm &realm);
+  virtual ~AssembleNodalGradPAlgorithmDriver() {}
+
+  virtual void pre_work() override;
+  virtual void post_work() override;
+
+};
+
+} // namespace nalu
+} // namespace Sierra
+
+#endif
diff --git a/include/LowMachEquationSystem.h b/include/LowMachEquationSystem.h
@@ -23,6 +23,7 @@ namespace nalu{
 class AlgorithmDriver;
 class Realm;
 class AssembleNodalGradAlgorithmDriver;
+class AssembleNodalGradPAlgorithmDriver;
 class AssembleNodalGradUAlgorithmDriver;
 class MomentumEquationSystem;
 class ContinuityEquationSystem;
@@ -259,7 +260,7 @@ class ContinuityEquationSystem : public EquationSystem {
 
   ScalarFieldType *pTmp_;
 
-  AssembleNodalGradAlgorithmDriver *assembleNodalGradAlgDriver_;
+  AssembleNodalGradPAlgorithmDriver *assembleNodalGradPAlgDriver_;
   ComputeMdotAlgorithmDriver *computeMdotAlgDriver_;
   ProjectedNodalGradientEquationSystem *projectedNodalGradEqs_;
 };

diff --git a/include/SolutionOptions.h b/include/SolutionOptions.h
@@ -159,6 +159,7 @@ class SolutionOptions
   double mdotAlgOpenCorrection_;
   size_t mdotAlgOpenIpCount_;
   double mdotAlgOpenPost_;
+  bool explicitlyZeroOpenPressureGradient_;
 
   // turbulence model coeffs
   std::map<TurbulenceModelConstant, double> turbModelConstantMap_;

diff --git a/include/user_functions/CappingInversionTemperatureAuxFunction.h b/include/user_functions/CappingInversionTemperatureAuxFunction.h
@@ -0,0 +1,40 @@
+/*------------------------------------------------------------------------*/
+/*  Copyright 2014 Sandia Corporation.                                    */
+/*  This software is released under the license detailed                  */
+/*  in the file, LICENSE, which is located in the top-level Nalu          */
+/*  directory structure                                                   */
+/*------------------------------------------------------------------------*/
+
+#ifndef CappingInversionTemperatureAuxFunction_h
+#define CappingInversionTemperatureAuxFunction_h
+
+#include <AuxFunction.h>
+
+#include <vector>
+
+namespace sierra{
+namespace nalu{
+
+class CappingInversionTemperatureAuxFunction : public AuxFunction
+{
+public:
+
+  CappingInversionTemperatureAuxFunction();
+
+  virtual ~CappingInversionTemperatureAuxFunction() {}
+
+  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
diff --git a/reg_tests/test_files/ablUnstableEdge/ablUnstableEdge_rst.i b/reg_tests/test_files/ablUnstableEdge/ablUnstableEdge_rst.i
@@ -147,6 +147,7 @@ realms:
       turbulence_model: wale
       interp_rhou_together_for_mdot: yes
       activate_open_mdot_correction: yes
+      explicitly_zero_open_pressure_gradient: yes 
 
       options:
 
@@ -231,6 +232,7 @@ realms:
       name: myOptions
       input_variables_interpolate_in_time: yes
       input_variables_from_file_restoration_time: 0.0
+      explicitly_zero_open_pressure_gradient: yes
 
       options:    
         - input_variables_from_file:

diff --git a/reg_tests/test_files/ablUnstableEdge/ablUnstableEdge_rst.norm.gold b/reg_tests/test_files/ablUnstableEdge/ablUnstableEdge_rst.norm.gold
@@ -1,10 +1,10 @@
-13.79759875908612 1 1.25
-0.8061222546965402 2 2.8125
-0.562023679352072 3 4.76562
-11.90085972620966 4 7.20703
-6.266546368434017 5 9.68282
-9.635424897738176 6 12.1237
-17.95292177378264 7 14.5323
-4.938092546499202 8 16.9087
-10.69998797582451 9 19.2429
-32.6178905646968 10 21.5424
+13.79737695917944 1 1.25
+0.8063289075114444 2 2.8125
+0.5629518334281949 3 4.76562
+11.90121212930885 4 7.20703
+6.266459299114715 5 9.68282
+9.635180348658634 6 12.1237
+17.95413017303376 7 14.5323
+4.939213131231322 8 16.9087
+10.70056020635592 9 19.2429
+32.63860650384363 10 21.5424
diff --git a/src/AssembleNodalGradPAlgorithmDriver.C b/src/AssembleNodalGradPAlgorithmDriver.C
@@ -0,0 +1,125 @@
+/*------------------------------------------------------------------------*/
+/*  Copyright 2014 Sandia Corporation.                                    */
+/*  This software is released under the license detailed                  */
+/*  in the file, LICENSE, which is located in the top-level Nalu          */
+/*  directory structure                                                   */
+/*------------------------------------------------------------------------*/
+
+
+#include <AssembleNodalGradPAlgorithmDriver.h>
+#include <FieldTypeDef.h>
+#include <Realm.h>
+#include <SolutionOptions.h>
+
+
+// stk_mesh/base/fem
+#include <stk_mesh/base/BulkData.hpp>
+#include <stk_mesh/base/Field.hpp>
+#include <stk_mesh/base/FieldParallel.hpp>
+#include <stk_mesh/base/GetBuckets.hpp>
+#include <stk_mesh/base/GetEntities.hpp>
+#include <stk_mesh/base/MetaData.hpp>
+#include <stk_mesh/base/Part.hpp>
+#include <stk_mesh/base/FieldBLAS.hpp>
+
+namespace sierra{
+namespace nalu{
+
+class Realm;
+
+
+void subtract_normal_component(int dim, const double* normal, double* grad_p)
+{
+  double grad_p_dot_n = 0;
+  for (int j = 0; j < dim; ++j) {
+    grad_p_dot_n += normal[j] * grad_p[j];
+  }
+
+  for (int j = 0; j < dim; ++j) {
+    grad_p[j] -= grad_p_dot_n * normal[j];
+  }
+}
+
+void subtract_normal_pressure_gradient(const stk::mesh::BulkData& bulk, VectorFieldType& dqdxField)
+{
+  const auto& meta = bulk.mesh_meta_data();
+  ThrowRequireMsg(meta.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal") != nullptr,
+    "average_open_normal field required");
+  VectorFieldType& meanNormalField = *meta.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal");
+
+  const int nDim = meta.spatial_dimension();
+
+  const stk::mesh::Selector& owned_or_shared_open = (meta.locally_owned_part() | meta.globally_shared_part())
+      & stk::mesh::selectField(meanNormalField);
+  const stk::mesh::BucketVector& node_buckets = bulk.get_buckets(stk::topology::NODE_RANK, owned_or_shared_open);
+  for (const auto* ib : node_buckets) {
+    const auto& b = *ib;
+    const size_t length = b.size();
+    for (size_t k = 0u; k < length; ++k) {
+      const stk::mesh::Entity node = b[k];
+      const double* meanNormal = stk::mesh::field_data(meanNormalField, node);
+      double* dqdx = stk::mesh::field_data(dqdxField, node);
+      subtract_normal_component(nDim, meanNormal, dqdx);
+    }
+  }
+}
+
+//==========================================================================
+// Class Definition
+//==========================================================================
+// AssembleNodalGradPAlgorithmDriver - Drives nodal grad algorithms
+//==========================================================================
+//--------------------------------------------------------------------------
+//-------- constructor -----------------------------------------------------
+//--------------------------------------------------------------------------
+AssembleNodalGradPAlgorithmDriver::AssembleNodalGradPAlgorithmDriver(Realm &realm)
+  : AlgorithmDriver(realm)
+{
+}
+
+//--------------------------------------------------------------------------
+//-------- pre_work --------------------------------------------------------
+//--------------------------------------------------------------------------
+void
+AssembleNodalGradPAlgorithmDriver::pre_work()
+{
+  ThrowRequire(realm_.meta_data().get_field<VectorFieldType>(stk::topology::NODE_RANK, "dpdx") != nullptr);
+  VectorFieldType& dpdxField = *realm_.meta_data().get_field<VectorFieldType>(stk::topology::NODE_RANK, "dpdx");
+  stk::mesh::field_fill(0.0, dpdxField);
+}
+
+//--------------------------------------------------------------------------
+//-------- post_work -------------------------------------------------------
+//--------------------------------------------------------------------------
+void
+AssembleNodalGradPAlgorithmDriver::post_work()
+{
+
+  stk::mesh::BulkData & bulk_data = realm_.bulk_data();
+  stk::mesh::MetaData & meta_data = realm_.meta_data();
+
+  ThrowRequire(realm_.meta_data().get_field<VectorFieldType>(stk::topology::NODE_RANK, "dpdx") != nullptr);
+  VectorFieldType& dpdxField = *realm_.meta_data().get_field<VectorFieldType>(stk::topology::NODE_RANK, "dpdx");
+
+  // extract fields
+  stk::mesh::parallel_sum(bulk_data, {&dpdxField});
+
+  if ( realm_.hasPeriodic_) {
+    const unsigned nDim = meta_data.spatial_dimension();
+    realm_.periodic_field_update(&dpdxField, nDim);
+  }
+
+  if ( realm_.hasOverset_ ) {
+    // this is a tensor
+    const unsigned nDim = meta_data.spatial_dimension();
+    realm_.overset_orphan_node_field_update(&dpdxField, 1, nDim);
+  }
+
+  if (realm_.solutionOptions_->explicitlyZeroOpenPressureGradient_) {
+    subtract_normal_pressure_gradient(bulk_data, dpdxField);
+  }
+
+}
+
+} // namespace nalu
+} // namespace Sierra
diff --git a/src/ComputeGeometryAlgorithmDriver.C b/src/ComputeGeometryAlgorithmDriver.C
@@ -13,11 +13,14 @@
 #include <FieldTypeDef.h>
 #include <master_element/MasterElement.h>
 #include <Realm.h>
+#include <SolutionOptions.h>
+#include <master_element/MasterElement.h>
 
 // stk_mesh/base/fem
 #include <stk_mesh/base/BulkData.hpp>
 #include <stk_mesh/base/Field.hpp>
 #include <stk_mesh/base/FieldParallel.hpp>
+#include <stk_mesh/base/FieldBLAS.hpp>
 #include <stk_mesh/base/GetBuckets.hpp>
 #include <stk_mesh/base/GetEntities.hpp>
 #include <stk_mesh/base/MetaData.hpp>
@@ -28,6 +31,79 @@ namespace nalu{
 
 class Realm;
 
+void compute_open_mean_normal(const stk::mesh::BulkData& bulk)
+{
+  const auto& meta = bulk.mesh_meta_data();
+  GenericFieldType& areavField = *meta.get_field<GenericFieldType>(meta.side_rank(), "exposed_area_vector");
+
+  ThrowRequireMsg(meta.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal") != nullptr,
+    "average_open_normal field required");
+
+  VectorFieldType& meanNormalField = *meta.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal");
+
+  ThrowRequireMsg(meta.get_field<ScalarIntFieldType>(stk::topology::NODE_RANK, "open_face_connection_count") != nullptr,
+    "open_face_connection_count field required");
+  ScalarIntFieldType& faceConnectionCount = *meta.get_field<ScalarIntFieldType>(stk::topology::NODE_RANK, "open_face_connection_count");
+
+  const stk::mesh::Selector& locally_owned_open = meta.locally_owned_part() & stk::mesh::selectField(meanNormalField);
+  const stk::mesh::BucketVector& face_buckets = bulk.get_buckets(meta.side_rank(), locally_owned_open);
+
+  const int nDim = meta.spatial_dimension();
+
+  for (const auto* ib : face_buckets) {
+    const auto& b = *ib;
+    const size_t length = b.size();
+
+    MasterElement& meFC = *MasterElementRepo::get_surface_master_element(b.topology());
+    const int numIp = meFC.numIntPoints_;
+    const int* ipNodeMap = meFC.ipNodeMap();
+
+    for (size_t k = 0u; k < length; ++k) {
+      const double* areaVecs = stk::mesh::field_data(areavField, b, k);
+      const auto* face_node_rels = bulk.begin_nodes(b[k]);
+
+      for (int ip = 0; ip < numIp; ++ip) {
+        const double* const areav = &areaVecs[ip*nDim];
+
+        double aMag = 0.0;
+        for ( int j = 0; j < nDim; ++j ) {
+          aMag += areav[j] * areav[j];
+        }
+        ThrowAssert(aMag > std::numeric_limits<double>::min());
+        const double inv_aMag = 1.0 / std::sqrt(aMag);
+
+        const stk::mesh::Entity nearestNode = face_node_rels[ipNodeMap[ip]];
+        double* meanNormal = stk::mesh::field_data(meanNormalField, nearestNode);
+        ThrowRequire(meanNormal != nullptr);
+
+        for (int j = 0; j < nDim; ++j) {
+          meanNormal[j] += areav[j] * inv_aMag;
+        }
+        ThrowRequire(stk::mesh::field_data(faceConnectionCount, nearestNode) != nullptr);
+        *stk::mesh::field_data(faceConnectionCount, nearestNode) += 1;
+      }
+    }
+  }
+
+  stk::mesh::parallel_sum(bulk, {&meanNormalField});
+  stk::mesh::parallel_sum(bulk, {&faceConnectionCount});
+
+  const stk::mesh::Selector& owned_or_shared_open = (meta.locally_owned_part() | meta.globally_shared_part()) & stk::mesh::selectField(meanNormalField);
+  const stk::mesh::BucketVector& node_buckets = bulk.get_buckets(stk::topology::NODE_RANK, owned_or_shared_open);
+  for (const auto* ib : node_buckets) {
+    const auto& b = *ib;
+    const size_t length = b.size();
+    for (size_t k = 0u; k < length; ++k) {
+      double* meanNormal = stk::mesh::field_data(meanNormalField, b, k);
+      const double avgFactor = *stk::mesh::field_data(faceConnectionCount, b, k);
+      ThrowAssert(avgFactor >= 1);
+      for (int j = 0; j < nDim; ++j) {
+        meanNormal[j] /= avgFactor;
+      }
+    }
+  }
+}
+
 //==========================================================================
 // Class Definition
 //==========================================================================
@@ -104,6 +180,19 @@ ComputeGeometryAlgorithmDriver::pre_work()
     }
   }
 
+  if (realm_.solutionOptions_->explicitlyZeroOpenPressureGradient_) {
+    ThrowRequireMsg(meta_data.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal") != nullptr,
+      "average_open_normal field required");
+
+    VectorFieldType& meanNormalField = *meta_data.get_field<VectorFieldType>(stk::topology::NODE_RANK, "average_open_normal");
+    stk::mesh::field_fill(0.0, meanNormalField );
+
+    ThrowRequireMsg(meta_data.get_field<ScalarIntFieldType>(stk::topology::NODE_RANK, "open_face_connection_count") != nullptr,
+      "open_face_connection_count field required");
+    ScalarIntFieldType& faceConnectionCount = *meta_data.get_field<ScalarIntFieldType>(stk::topology::NODE_RANK, "open_face_connection_count");
+    stk::mesh::field_fill(0, faceConnectionCount );
+  }
+
 }
 
 //--------------------------------------------------------------------------
@@ -139,6 +228,10 @@ ComputeGeometryAlgorithmDriver::post_work()
   if ( realm_.checkJacobians_ ) {
     check_jacobians();
   }
+
+  if (realm_.solutionOptions_->explicitlyZeroOpenPressureGradient_) {
+    compute_open_mean_normal(realm_.bulk_data());
+  }
 }
 
 //--------------------------------------------------------------------------