Merge pull request #5 from marchdf/add-init-channel

Adding channel initial condition utility and documentation

doc/manual/source/user/nalu_preprocess.rst

@@ -8,14 +8,15 @@ This utility loads an input mesh and performs various pre-processing *tasks* so
 that the resulting output database can be used in a wind LES simulation.
 Currently, the following *tasks* have been implemented within this utility.
 
-======================  ===========================================================
-Task type               Description
-======================  ===========================================================
-``init_abl_fields``     Initialize ABL velocity and temperature fields
-``generate_planes``     Generate horizontal sampling planes for ``dp/dx`` forcing
-``create_bdy_io_mesh``  Create an I/O transfer mesh for sampling inflow planes
-``rotate_mesh``         Rotate mesh
-======================  ===========================================================
+=======================  ===========================================================
+Task type                Description
+=======================  ===========================================================
+``init_abl_fields``      Initialize ABL velocity and temperature fields
+``init_channel_fields``  Initialize channel velocity fields
+``generate_planes``      Generate horizontal sampling planes for ``dp/dx`` forcing
+``create_bdy_io_mesh``   Create an I/O transfer mesh for sampling inflow planes
+``rotate_mesh``          Rotate mesh
+=======================  ===========================================================
 
 .. warning::
 
@@ -170,6 +171,35 @@ and writes it out to the :confval:`output_db`. It is safe to run
    process if the user intends to initialize uniform velocity throughout the
    domain within Nalu.
 
+``init_channel_fields``
+-----------------------
+
+This task initializes the velocity fields for channel flow simulations
+based on the parameters provided by the user and writes it out to the
+:confval:`output_db`. It is safe to run ``init_channel_fields`` in
+parallel. A sample invocation is shown below
+
+.. code-block:: yaml
+   :linenos:
+
+   init_channel_fields:
+     fluid_parts: [Unspecified-2-HEX]
+
+     velocity:
+       Re_tau : 550
+       viscosity : 0.0000157
+
+.. confval:: fluid_parts
+
+   A list of element block names where the velocity fields are to be
+   initialized.
+
+.. confval:: velocity
+
+   A YAML dictionary containing two values: the friction Reynolds
+   number, ``Re_tau``, and the kinematic ``viscosity``
+   (:math:`m^2/s`).
+
 ``generate_planes``
 -------------------
 

src/core/CFDMesh.h

@@ -85,7 +85,7 @@ class CFDMesh
 
     /** Register a field for output during write
      *
-     * \param Name of the field to be output
+     * \param field Name of the field to be output
      */
     inline void add_output_field(const std::string field)
     {
@@ -126,7 +126,7 @@ class CFDMesh
 
     /** Set automatic mesh decomposition property
      *
-     *  \param The decomposition type
+     *  \param decompType The decomposition type
      *
      *  Valid decomposition types are: rcb, rib, block, linear
      */

src/preprocessing/CMakeLists.txt

@@ -4,6 +4,7 @@ add_executable(nalu_preprocess
   NDTW2D.cpp
   SamplingPlanes.cpp
   ABLFields.cpp
+  ChannelFields.cpp
   RotateMesh.cpp
   BdyIOPlanes.cpp
 

src/preprocessing/ChannelFields.cpp

@@ -0,0 +1,145 @@
+//  Copyright 2016 National Renewable Energy Laboratory
+//
+//  Licensed under the Apache License, Version 2.0 (the "License");
+//  you may not use this file except in compliance with the License.
+//  You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+//  Unless required by applicable law or agreed to in writing, software
+//  distributed under the License is distributed on an "AS IS" BASIS,
+//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//  See the License for the specific language governing permissions and
+//  limitations under the License.
+//
+
+
+#include "ChannelFields.h"
+#include "core/LinearInterpolation.h"
+
+namespace sierra {
+namespace nalu {
+
+REGISTER_DERIVED_CLASS(PreProcessingTask, ChannelFields, "init_channel_fields");
+
+ChannelFields::ChannelFields(
+    CFDMesh& mesh,
+    const YAML::Node& node
+) : PreProcessingTask(mesh),
+    meta_(mesh.meta()),
+    bulk_(mesh.bulk()),
+    fluid_parts_(0),
+    ndim_(meta_.spatial_dimension()),
+    doVelocity_(false),
+    Re_tau_(0.0),
+    viscosity_(0.0)
+{
+    load(node);
+    srand(seed_);
+}
+
+void ChannelFields::load(const YAML::Node& channel)
+{
+    auto fluid_partnames = channel["fluid_parts"].as<std::vector<std::string>>();
+
+    if (channel["velocity"]) {
+        doVelocity_ = true;
+        load_velocity_info(channel["velocity"]);
+    }
+
+    fluid_parts_.resize(fluid_partnames.size());
+
+    for(size_t i=0; i<fluid_partnames.size(); i++) {
+        stk::mesh::Part* part = meta_.get_part(fluid_partnames[i]);
+        if (NULL == part) {
+            throw std::runtime_error("Missing fluid part in mesh database: " +
+                                     fluid_partnames[i]);
+        } else {
+            fluid_parts_[i] = part;
+        }
+    }
+}
+
+void ChannelFields::initialize()
+{
+    if (doVelocity_) {
+        VectorFieldType& velocity = meta_.declare_field<VectorFieldType>(
+            stk::topology::NODE_RANK, "velocity");
+        for(auto part: fluid_parts_) {
+            stk::mesh::put_field(velocity, *part);
+        }
+        mesh_.add_output_field("velocity");
+    }
+}
+
+void ChannelFields::run()
+{
+    if (bulk_.parallel_rank() == 0)
+        std::cout << "Generating channel fields" << std::endl;
+    if (doVelocity_) init_velocity_field();
+
+    mesh_.set_write_flag();
+}
+
+void ChannelFields::load_velocity_info(const YAML::Node& channel)
+{
+  if (channel["Re_tau"])
+    Re_tau_ = channel["Re_tau"].as<double>();
+  else
+    throw std::runtime_error("ChannelFields: missing mandatory Re_tau parameter");
+
+  if (channel["viscosity"])
+    viscosity_ = channel["viscosity"].as<double>();
+  else
+    throw std::runtime_error("ChannelFields: missing mandatory viscosity parameter");
+}
+
+void ChannelFields::init_velocity_field()
+{
+    stk::mesh::Selector fluid_union = stk::mesh::selectUnion(fluid_parts_);
+    auto bbox = mesh_.calc_bounding_box(fluid_union,false);
+    const stk::mesh::BucketVector& fluid_bkts = bulk_.get_buckets(
+        stk::topology::NODE_RANK, fluid_union);
+
+    VectorFieldType* coords = meta_.get_field<VectorFieldType>(
+        stk::topology::NODE_RANK, "coordinates");
+    VectorFieldType* velocity = meta_.get_field<VectorFieldType>(
+        stk::topology::NODE_RANK, "velocity");
+
+    const double length = bbox.get_x_max() - bbox.get_x_min();
+    const double height = bbox.get_y_max() - bbox.get_y_min();
+    const double width = bbox.get_z_max() - bbox.get_z_min();
+    const double delta = 0.5 * height;
+    const double utau = Re_tau_ * viscosity_ / delta;
+    const double yph = height * utau / viscosity_;
+    const double C = (1.0/kappa_ * log(1.0 + kappa_ * yph)) / (1 - exp(- yph / 11.0) - yph / 11.0 * exp(- yph / 3)) + log(kappa_) / kappa_;
+
+    for(size_t ib=0; ib < fluid_bkts.size(); ib++) {
+        stk::mesh::Bucket& fbkt = *fluid_bkts[ib];
+        double* xyz = stk::mesh::field_data(*coords, fbkt);
+        double* vel = stk::mesh::field_data(*velocity, fbkt);
+
+        for (size_t in=0; in < fbkt.size(); in++) {
+	  const double x = xyz[in*ndim_ + 0];
+	  const double y = xyz[in*ndim_ + 1];
+	  const double z = xyz[in*ndim_ + 2];
+
+	  const double yp = std::min(y, height - y) * utau / viscosity_;
+	  const double reichardt = (1.0/kappa_ * log(1.0 + kappa_ * yp)) + (C - log(kappa_) / kappa_) * (1 - exp(- yp / 11.0) - yp / 11.0 * exp(- yp / 3));
+
+	  const double pert_u = sin(k_pert_u_ * M_PI / width * z);
+	  const double pert_w = sin(k_pert_w_ * M_PI / length * x);
+
+	  const double rand_u = a_rand_u_ * 2. * (double)rand() / RAND_MAX - 1;
+	  const double rand_w = a_rand_w_ * 2. * (double)rand() / RAND_MAX - 1;
+
+	  vel[in * ndim_ + 0] = utau * (reichardt + pert_u + rand_u);
+	  vel[in * ndim_ + 1] = 0.0;
+	  vel[in * ndim_ + 2] = utau * (pert_w + rand_w);
+        }
+    }
+}
+
+
+} // nalu
+} // sierra

src/preprocessing/ChannelFields.h

@@ -0,0 +1,117 @@
+//  Copyright 2016 National Renewable Energy Laboratory
+//
+//  Licensed under the Apache License, Version 2.0 (the "License");
+//  you may not use this file except in compliance with the License.
+//  You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+//  Unless required by applicable law or agreed to in writing, software
+//  distributed under the License is distributed on an "AS IS" BASIS,
+//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//  See the License for the specific language governing permissions and
+//  limitations under the License.
+//
+
+#ifndef CHANNELFIELDS_H
+#define CHANNELFIELDS_H
+
+#include "PreProcessingTask.h"
+
+namespace sierra {
+namespace nalu {
+
+/**
+ * Initialize velocity fields for channel flow simulations
+ *
+ * This task is activated by using the `init_channel_fields` task in the
+ * preprocessing input file. It requires a section `init_channel_fields` in the
+ * `nalu_preprocess` section with the following parameters:
+ *
+ *  ```
+ *  init_channel_fields:
+ *    fluid_parts: [Unspecified-2-HEX]
+ *
+ *    velocity:
+ *      Re_tau : 550
+ *      viscosity : 0.0000157
+ *  ```
+ *
+ * The user specified the friction Reynolds number, `Re_tau`, and the
+ * kinematic `viscosity` (in m^2/s). The velocity field is initialized
+ * to a Reichardt function, with an imposed sinusoidal perturbation
+ * and random perturbation in the wall parallel directions.
+ */
+class ChannelFields: public PreProcessingTask
+{
+public:
+    template<typename T>
+    using Array2D = std::vector<std::vector<T>>;
+
+    ChannelFields(CFDMesh&, const YAML::Node&);
+
+    virtual ~ChannelFields() {}
+
+    //! Declare velocity fields and register them for output
+    void initialize();
+
+    //! Initialize the velocity fields by linear interpolation
+    void run();
+
+private:
+    ChannelFields() = delete;
+    ChannelFields(const ChannelFields&) = delete;
+
+    //! Parse the YAML file and initialize parameters
+    void load(const YAML::Node&);
+
+    //! Helper function to parse and initialize velocity inputs
+    void load_velocity_info(const YAML::Node&);
+
+    //! Initialize the velocity field through linear interpolation
+    void init_velocity_field();
+
+    //! STK Metadata object
+    stk::mesh::MetaData& meta_;
+
+    //! STK Bulkdata object
+    stk::mesh::BulkData& bulk_;
+
+    //! Parts of the fluid mesh where velocity is initialized
+    stk::mesh::PartVector fluid_parts_;
+
+    //! Dimensionality of the mesh
+    int ndim_;
+
+    //! Flag indicating whether velocity is initialized
+    bool doVelocity_;
+
+    //! Skin friction velocity
+    double Re_tau_;
+
+    //! Kinematic viscosity
+    double viscosity_;
+
+    //! Von Karman constant
+    const double kappa_ = 0.4;
+
+    //! Seed for RNG
+    const int seed_ = 2864;
+
+    //! Wavenumber of sinusoidal perturbation for u
+    const int k_pert_u_ = 16;
+
+    //! Wavenumber of sinusoidal perturbation for w
+    const int k_pert_w_ = 16;
+
+    //! Amplitude of random perturbation for u
+    const double a_rand_u_ = 0.2;
+
+    //! Amplitude of random perturbation for w
+    const double a_rand_w_ = 0.2;
+};
+
+}
+}
+
+#endif /* CHANNELFIELDS_H */