Consistent limiting of Pk for SST DES (#462)

* Missing Pk scaling in SST DES SDR equation; courtesy
  of RCK.

* diff accepted in elemBackStepLRSST.

  Scary.... Coverage for read from file lost in transition...
  Subsequent input file named "_Input" was stripped from the
  suite thereby removing "initialization" transfer objective
  coverage. Added it back as per the "_rst" test attribute.

* remove odd comment in WALE.

* add theory comment on limiting approach for Pk when using
SST/DES.

docs/source/theory/supportedEquationSet.rst

@@ -796,6 +796,11 @@ and :math:`c_{DES}` represents a blended set of DES constants:
 scale, :math:`l_{DES}` is the maximum edge length scale touching a given
 node.
 
+Note that the production term appearing in both the turbulent kinetic
+energy and specific dissipation rate equation is limited to a 
+user-supplied scaling of the above dissipation term. Currently, all
+SST variants use a production to dissipation ratio limiting of ten.
+
 Solid Stress
 ++++++++++++
 

include/SpecificDissipationRateSSTDESNodeSourceSuppAlg.h

@@ -0,0 +1,58 @@
+/*------------------------------------------------------------------------*/
+/*  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 SpecificDissipationRateSSTDESNodeSourceSuppAlg_h
+#define SpecificDissipationRateSSTDESNodeSourceSuppAlg_h
+
+#include <SupplementalAlgorithm.h>
+#include <FieldTypeDef.h>
+
+#include <stk_mesh/base/Entity.hpp>
+
+namespace sierra{
+namespace nalu{
+
+class Realm;
+
+class SpecificDissipationRateSSTDESNodeSourceSuppAlg : public SupplementalAlgorithm
+{
+public:
+  SpecificDissipationRateSSTDESNodeSourceSuppAlg(
+    Realm &realm);
+
+  virtual ~SpecificDissipationRateSSTDESNodeSourceSuppAlg() {}
+
+  virtual void setup();
+
+  virtual void node_execute(
+    double *lhs,
+    double *rhs,
+    stk::mesh::Entity node);
+
+  const double sigmaWTwo_, betaStar_, betaOne_, betaTwo_, gammaOne_, gammaTwo_;
+  ScalarFieldType *sdrNp1_;
+  ScalarFieldType *tkeNp1_;
+  ScalarFieldType *densityNp1_;
+  ScalarFieldType *fOneBlend_;
+  ScalarFieldType *tvisc_;
+  GenericFieldType *dudx_;
+  VectorFieldType *dkdx_;
+  VectorFieldType *dwdx_;
+  ScalarFieldType *dualNodalVolume_;
+  ScalarFieldType *maxLengthScale_;
+  double tkeProdLimitRatio_;
+  int nDim_;
+  double cDESke_;
+  double cDESkw_;
+
+};
+
+} // namespace nalu
+} // namespace Sierra
+
+#endif

reg_tests/CTestList.cmake

@@ -90,7 +90,7 @@ add_test_r(ductElemWedge 2)
 add_test_r(ductWedge 2)
 add_test_r(edgeHybridFluids 8)
 add_test_r(edgePipeCHT 4)
-add_test_r(elemBackStepLRSST 4)
+add_test_r_rst(elemBackStepLRSST 4)
 add_test_r(elemClosedDomain 2)
 add_test_r(elemHybridFluids 8)
 add_test_r(elemHybridFluidsShift 8)

reg_tests/test_files/elemBackStepLRSST/elemBackStepLRSST.norm.gold

@@ -1,16 +1,16 @@
-12.87979495951468 1 0.0012
-7.898634260346332 2 0.00264
-4.866439930207051 3 0.00423449
-3.346588485701302 4 0.00590796
-2.273816685700415 5 0.00756626
-1.6506181242032 6 0.00921964
-1.268443842337867 7 0.0108755
-1.011900042572801 8 0.0125448
-0.8251254018370128 9 0.0142252
-0.6998735424347657 10 0.0159145
-0.5916860656595713 11 0.0176111
-0.5138579971671057 12 0.0193122
-0.4521984490856337 13 0.0210175
-0.4019927100476395 14 0.0227266
-0.3593993848099154 15 0.0244392
-0.3215127817442578 16 0.0261548
+12.8805390875359 1 0.0012
+7.902020224670594 2 0.00264
+4.864988638339074 3 0.00423439
+3.346334565810965 4 0.00590777
+2.273361569265452 5 0.00756599
+1.650618482632744 6 0.0092194
+1.268415324198112 7 0.0108752
+1.012269498614699 8 0.0125446
+0.8260282722975488 9 0.014225
+0.6925275376679251 10 0.0159145
+0.5908504803484264 11 0.0176107
+0.5140319986486012 12 0.0193117
+0.4523338503065832 13 0.021017
+0.4020473166775929 14 0.0227262
+0.3593900935222692 15 0.0244389
+0.3214471710030072 16 0.0261545

reg_tests/test_files/elemBackStepLRSST/elemBackStepLRSST_Input.i → reg_tests/test_files/elemBackStepLRSST/elemBackStepLRSST_rst.i

@@ -166,7 +166,7 @@ realms:
             specific_dissipation_rate: element 
 
     output:
-      output_data_base_name: elemBackStepLRSST_Input.e
+      output_data_base_name: elemBackStepLRSST_rst.e
       output_frequency: 1
       output_start: 0 
       output_node_set: no 

reg_tests/test_files/elemBackStepLRSST/elemBackStepLRSST_rst.norm.gold

@@ -0,0 +1,20 @@
+11.27319387063573 1 0.0171145
+27.950695295086 2 0.0185545
+27.94180500283688 3 0.0202585
+27.93029944851022 4 0.0219658
+0.3726086834841603 5 0.0236753
+27.91121367035497 6 0.0253869
+0.2993386763400038 7 0.0271008
+27.89533280265562 8 0.028817
+0.2441522398792308 9 0.0305356
+27.8832921001239 10 0.0322561
+0.2042377135506065 11 0.0339782
+27.87450793627816 12 0.0357016
+0.1756522151120993 13 0.037426
+27.86783986249212 14 0.0391512
+27.86531919000372 15 0.0408771
+27.8631699884068 16 0.0426037
+27.8613266894294 17 0.0443308
+27.85981392202672 18 0.0460584
+27.85848781032781 19 0.0477863
+27.85734006195818 20 0.0495145

src/SpecificDissipationRateEquationSystem.C

@@ -45,6 +45,7 @@
 #include "SolutionOptions.h"
 #include "TimeIntegrator.h"
 #include "SpecificDissipationRateSSTNodeSourceSuppAlg.h"
+#include "SpecificDissipationRateSSTDESNodeSourceSuppAlg.h"
 #include "SolverAlgorithmDriver.h"
 
 // template for supp algs
@@ -312,11 +313,24 @@ SpecificDissipationRateEquationSystem::register_interior_algorithm(
       }
 
       // now create the src alg for sdr source
-      SpecificDissipationRateSSTNodeSourceSuppAlg *theSrc
-        = new SpecificDissipationRateSSTNodeSourceSuppAlg(realm_);
+      SupplementalAlgorithm *theSrc = NULL;
+      switch(realm_.solutionOptions_->turbulenceModel_) {
+      case SST:
+        {
+          theSrc = new SpecificDissipationRateSSTNodeSourceSuppAlg(realm_);
+        }
+        break;
+      case SST_DES:
+        {
+          theSrc = new SpecificDissipationRateSSTDESNodeSourceSuppAlg(realm_);
+        }
+        break;
+      default:
+        throw std::runtime_error("Unsupported turbulence model in TurbKe: only SST and SST_DES supported");
+      }
       theAlg->supplementalAlg_.push_back(theSrc);
 
-      // Add src term supp alg...; limited number supported
+      // Add nodal src term supp alg...; limited number supported
       std::map<std::string, std::vector<std::string> >::iterator isrc
         = realm_.solutionOptions_->srcTermsMap_.find("specific_dissipation_rate");
       if (isrc != realm_.solutionOptions_->srcTermsMap_.end()) {

src/SpecificDissipationRateSSTDESNodeSourceSuppAlg.C

@@ -0,0 +1,142 @@
+/*------------------------------------------------------------------------*/
+/*  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 <SpecificDissipationRateSSTDESNodeSourceSuppAlg.h>
+#include <FieldTypeDef.h>
+#include <Realm.h>
+#include <SupplementalAlgorithm.h>
+#include <TimeIntegrator.h>
+
+// stk_mesh/base/fem
+#include <stk_mesh/base/Entity.hpp>
+#include <stk_mesh/base/MetaData.hpp>
+#include <stk_mesh/base/BulkData.hpp>
+#include <stk_mesh/base/Field.hpp>
+
+namespace sierra{
+namespace nalu{
+
+//==========================================================================
+// Class Definition
+//==========================================================================
+// SpecificDissipationRateSSTDESNodeSourceSuppAlg - base class for algorithm
+//==========================================================================
+//--------------------------------------------------------------------------
+//-------- constructor -----------------------------------------------------
+//--------------------------------------------------------------------------
+SpecificDissipationRateSSTDESNodeSourceSuppAlg::SpecificDissipationRateSSTDESNodeSourceSuppAlg(
+  Realm &realm)
+  : SupplementalAlgorithm(realm),
+    sigmaWTwo_(realm.get_turb_model_constant(TM_sigmaWTwo)),
+    betaStar_(realm.get_turb_model_constant(TM_betaStar)),
+    betaOne_(realm.get_turb_model_constant(TM_betaOne)),
+    betaTwo_(realm.get_turb_model_constant(TM_betaTwo)),
+    gammaOne_(realm.get_turb_model_constant(TM_gammaOne)),
+    gammaTwo_(realm.get_turb_model_constant(TM_gammaTwo)),
+    sdrNp1_(NULL),
+    tkeNp1_(NULL),
+    densityNp1_(NULL),
+    fOneBlend_(NULL),
+    tvisc_(NULL),
+    dkdx_(NULL),
+    dwdx_(NULL),
+    dualNodalVolume_(NULL),
+    maxLengthScale_(NULL),
+    tkeProdLimitRatio_(realm.get_turb_model_constant(TM_tkeProdLimitRatio)),
+    nDim_(realm_.meta_data().spatial_dimension()),
+    cDESke_(realm.get_turb_model_constant(TM_cDESke)),
+    cDESkw_(realm.get_turb_model_constant(TM_cDESkw))
+{
+  // save off fields
+  stk::mesh::MetaData & meta_data = realm_.meta_data();
+  ScalarFieldType *sdr = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "specific_dissipation_rate");
+  sdrNp1_ = &(sdr->field_of_state(stk::mesh::StateNP1));
+  ScalarFieldType *tke = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "turbulent_ke");
+  tkeNp1_ = &(tke->field_of_state(stk::mesh::StateNP1));
+  ScalarFieldType *density = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "density");
+  densityNp1_ = &(density->field_of_state(stk::mesh::StateNP1));
+  fOneBlend_ = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "sst_f_one_blending");
+  tvisc_ = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "turbulent_viscosity");
+  dudx_ = meta_data.get_field<GenericFieldType>(stk::topology::NODE_RANK, "dudx");
+  dkdx_ = meta_data.get_field<VectorFieldType>(stk::topology::NODE_RANK, "dkdx");
+  dwdx_ = meta_data.get_field<VectorFieldType>(stk::topology::NODE_RANK, "dwdx");
+  maxLengthScale_ = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "sst_max_length_scale");
+  dualNodalVolume_ = meta_data.get_field<ScalarFieldType>(stk::topology::NODE_RANK, "dual_nodal_volume");
+}
+
+//--------------------------------------------------------------------------
+//-------- setup -----------------------------------------------------------
+//--------------------------------------------------------------------------
+void
+SpecificDissipationRateSSTDESNodeSourceSuppAlg::setup()
+{
+  // could extract user-based values
+}
+
+//--------------------------------------------------------------------------
+//-------- node_execute ----------------------------------------------------
+//--------------------------------------------------------------------------
+void
+SpecificDissipationRateSSTDESNodeSourceSuppAlg::node_execute(
+  double *lhs,
+  double *rhs,
+  stk::mesh::Entity node)
+{
+  const double sdr        = *stk::mesh::field_data(*sdrNp1_, node );
+  const double tke        = *stk::mesh::field_data(*tkeNp1_, node );
+  const double rho        = *stk::mesh::field_data(*densityNp1_, node );
+  const double fOneBlend  = *stk::mesh::field_data(*fOneBlend_, node );
+  const double tvisc      = *stk::mesh::field_data(*tvisc_, node );
+  const double *dudx      =  stk::mesh::field_data(*dudx_, node );
+  const double *dkdx      =  stk::mesh::field_data(*dkdx_, node );
+  const double *dwdx      =  stk::mesh::field_data(*dwdx_, node );
+  const double dualVolume = *stk::mesh::field_data(*dualNodalVolume_, node );
+  const double maxLengthScale = *stk::mesh::field_data(*maxLengthScale_, node );
+
+  int nDim = nDim_;
+  double Pk = 0.0;
+  double crossDiff = 0.0;
+  for ( int i = 0; i < nDim; ++i ) {
+    crossDiff += dkdx[i]*dwdx[i];
+    const int offSet = nDim*i;
+    for ( int j = 0; j < nDim; ++j ) {
+      Pk += dudx[offSet+j]*(dudx[offSet+j] + dudx[nDim*j+i]);
+    }
+  }
+  Pk *= tvisc;
+
+  // blend cDES  
+  const double cDES = fOneBlend*cDESkw_ + (1.0-fOneBlend)*cDESke_;
+
+  const double sqrtTke = std::sqrt(tke);
+  const double lSST = sqrtTke/betaStar_/sdr;
+  // prevent divide by zero (possible at resolved tke bcs = 0.0)
+  const double lDES = std::max(1.0e-16, std::min(lSST, cDES*maxLengthScale));
+
+  double Dk = rho*tke*sqrtTke/lDES;
+
+  if ( Pk > tkeProdLimitRatio_*Dk )
+    Pk = tkeProdLimitRatio_*Dk;
+
+  // start the blending and constants
+  const double om_fOneBlend = 1.0 - fOneBlend;
+  const double beta = fOneBlend*betaOne_ + om_fOneBlend*betaTwo_;
+  const double gamma = fOneBlend*gammaOne_ + om_fOneBlend*gammaTwo_;
+  const double sigmaD = 2.0*om_fOneBlend*sigmaWTwo_;
+
+  // Pw includes 1/tvisc scaling; tvisc may be zero at a dirichlet low Re approach (clip)
+  const double Pw = gamma*rho*Pk/std::max(tvisc, 1.0e-16);
+  const double Dw = beta*rho*sdr*sdr;
+  const double Sw = sigmaD*rho*crossDiff/sdr;
+
+  rhs[0] += (Pw - Dw + Sw)*dualVolume;
+  lhs[0] += (2.0*beta*rho*sdr + std::max(Sw/sdr,0.0))*dualVolume;
+}
+
+} // namespace nalu
+} // namespace Sierra

src/TurbViscWaleAlgorithm.C

@@ -122,7 +122,7 @@ TurbViscWaleAlgorithm::execute()
       }
 
       const double filter = std::pow(dualNodalVolume[k], invNdim);
-      const double Ls = Cw_*filter; // min (Cw_*filter, kappa*ndtw)
+      const double Ls = Cw_*filter;
       const double numer = std::pow(SijdSq, threeHalves) + small*small;
       const double demom = std::pow(SijSq, fiveHalves)+std::pow(SijdSq, fiveFourths) + small;
       tvisc[k] = density[k]*Ls*Ls*numer/demom;