LBM: Q criterion, curl and vorticity calculation and vtk output

src/lbm/field/MacroscopicValueCalculation.h

@@ -27,6 +27,8 @@
 #include "Equilibrium.h"
 #include "ShearRate.h"
 #include "PressureTensor.h"
+#include "QCriterion.h"
+#include "Vorticity.h"
 
 #include "core/DataTypes.h"
 #include "core/math/Vector3.h"
@@ -433,6 +435,21 @@ inline void getPressureTensor( Matrix3< real_t > & pressureTensor, const Lattice
 }
 
 
+template< typename VelocityField_T, typename Filter_T >
+inline real_t getQCriterion(const VelocityField_T &velocityField, const Filter_T &filter,
+                            const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+                            real_t dx = real_t(1), real_t dy = real_t(1), real_t dz = real_t(1)) {
+   return QCriterion::get(velocityField, filter, x, y, z, dx, dy, dz);
+}
+
+template< typename VelocityField_T, typename Filter_T >
+inline Vector3<real_t> getVorticity(const VelocityField_T &velocityField, const Filter_T &filter,
+                                    const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+                                    real_t dx = real_t(1), real_t dy = real_t(1), real_t dz = real_t(1)) {
+   return Vorticity::get(velocityField, filter, x, y, z, dx, dy, dz);
+}
+
+
 
 } // namespace lbm
 } // namespace walberla

src/lbm/field/QCriterion.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file QCriterion.h
+//! \ingroup lbm
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+
+#include "lbm/lattice_model/CollisionModel.h"
+#include "lbm/lattice_model/EquilibriumDistribution.h"
+
+#include <type_traits>
+
+
+
+// Back-end for calculating macroscopic values
+// You should never use these functions directly, always refer to the member functions
+// of PdfField or the free functions that can be found in MacroscopicValueCalculation.h
+
+namespace walberla {
+namespace lbm {
+
+struct QCriterion {
+   template< typename VelocityField_T, typename Filter_T >
+   static inline real_t get(const VelocityField_T & velocityField, const Filter_T & filter,
+           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+           real_t dx = real_t(1), real_t dy = real_t(1), real_t dz = real_t(1)) {
+      const auto one = cell_idx_t(1);
+
+      if(filter(x,y,z) && filter(x+one,y,z) && filter(x-one,y,z) && filter(x,y+one,z)
+         && filter(x,y-one,z) && filter(x,y,z+one) && filter(x,y,z-one)) {
+         const Vector3<real_t> xa = velocityField.get(x+one,y,z);
+         const Vector3<real_t> xb = velocityField.get(x-one,y,z);
+         const Vector3<real_t> ya = velocityField.get(x,y+one,z);
+         const Vector3<real_t> yb = velocityField.get(x,y-one,z);
+         const Vector3<real_t> za = velocityField.get(x,y,z+one);
+         const Vector3<real_t> zb = velocityField.get(x,y,z-one);
+
+         return calculate(xa, xb, ya, yb, za, zb, dx, dy, dz);
+      }
+
+      return real_t(0);
+   }
+
+   static inline real_t calculate(const Vector3<real_t> xa, const Vector3<real_t> xb,
+           const Vector3<real_t> ya, const Vector3<real_t> yb,
+           const Vector3<real_t> za, const Vector3<real_t> zb,
+           const real_t dx, const real_t dy, const real_t dz) {
+
+      const auto halfInvDx = real_t(0.5) / dx;
+      const auto halfInvDy = real_t(0.5) / dy;
+      const auto halfInvDz = real_t(0.5) / dz;
+
+      const real_t duxdx = (xa[0] - xb[0]) * halfInvDx;
+      const real_t duxdy = (ya[0] - yb[0]) * halfInvDy;
+      const real_t duxdz = (za[0] - zb[0]) * halfInvDz;
+
+      const real_t duydx = (xa[1] - xb[1]) * halfInvDx;
+      const real_t duydy = (ya[1] - yb[1]) * halfInvDy;
+      const real_t duydz = (za[1] - zb[1]) * halfInvDz;
+
+      const real_t duzdx = (xa[2] - xb[2]) * halfInvDx;
+      const real_t duzdy = (ya[2] - yb[2]) * halfInvDy;
+      const real_t duzdz = (za[2] - zb[2]) * halfInvDz;
+
+      // Q = 1/2 * (||W||² - ||S||²)
+      real_t sNormSq = duxdx*duxdx + duydy*duydy + duzdz*duzdz +
+              real_t(0.5)*(duxdy+duydx)*(duxdy+duydx) + real_t(0.5)*(duydz+duzdy)*(duydz+duzdy) +
+              real_t(0.5)*(duxdz+duzdx)*(duxdz+duzdx);
+
+      real_t omegaNormSq = real_t(0.5)*(duxdz-duzdx)*(duxdz-duzdx) +
+              real_t(0.5)*(duxdy-duydx)*(duxdy-duydx) +
+              real_t(0.5)*(duydz-duzdy)*(duydz-duzdy);
+
+      return real_t(0.5) * (omegaNormSq - sNormSq);
+   }
+};
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/field/QCriterionFieldWriter.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file VelocityFieldWriter.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+#include "domain_decomposition/BlockDataID.h"
+#include "domain_decomposition/IBlock.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "field/iterators/IteratorMacros.h"
+#include "lbm/field/QCriterion.h"
+
+namespace walberla {
+namespace lbm {
+
+template< typename VelocityField_T, typename QCriterionField_T, typename Filter_T >
+class QCriterionFieldWriter {
+public:
+   QCriterionFieldWriter(const shared_ptr<StructuredBlockStorage> & blockStorage, const ConstBlockDataID & velocityFieldId,
+           const BlockDataID & qCriterionFieldId, Filter_T & filter) :
+      blockStorage_(blockStorage), velocityFieldId_(velocityFieldId), qCriterionFieldId_(qCriterionFieldId),
+      filter_(filter)
+   {}
+
+   void operator()(IBlock * const block) {
+      const VelocityField_T* velocityField = block->template getData<VelocityField_T>(velocityFieldId_);
+      QCriterionField_T* qCriterionField = block->template getData<QCriterionField_T>(qCriterionFieldId_);
+
+      WALBERLA_ASSERT_EQUAL(velocityField->xyzSize(), qCriterionField->xyzSize());
+
+      const real_t dx = blockStorage_->dx(blockStorage_->getLevel(*block));
+      const real_t dy = blockStorage_->dy(blockStorage_->getLevel(*block));
+      const real_t dz = blockStorage_->dz(blockStorage_->getLevel(*block));
+
+      filter_(*block);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(velocityField,
+         store(velocityField, qCriterionField, x, y, z, dx, dy, dz);
+      )
+   }
+
+private:
+
+   void store(const VelocityField_T* velocityField, QCriterionField_T* qCriterionField,
+           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+           const real_t dx, const real_t dy, const real_t dz) {
+      qCriterionField->get(x,y,z,0) = QCriterion::get(*velocityField, filter_, x, y, z, dx, dy, dz);
+   }
+
+   const shared_ptr<StructuredBlockStorage> blockStorage_;
+   ConstBlockDataID velocityFieldId_;
+   BlockDataID qCriterionFieldId_;
+   Filter_T & filter_;
+};
+
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/field/Vorticity.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file QCriterion.h
+//! \ingroup lbm
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+
+#include "lbm/lattice_model/CollisionModel.h"
+#include "lbm/lattice_model/EquilibriumDistribution.h"
+
+#include <type_traits>
+
+
+
+// Back-end for calculating macroscopic values
+// You should never use these functions directly, always refer to the member functions
+// of PdfField or the free functions that can be found in MacroscopicValueCalculation.h
+
+namespace walberla {
+namespace lbm {
+
+struct Vorticity {
+   template< typename VelocityField_T, typename Filter_T >
+   static inline Vector3<real_t> get(const VelocityField_T & velocityField, const Filter_T & filter,
+           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+           real_t dx = real_t(1), real_t dy = real_t(1), real_t dz = real_t(1)) {
+      const auto one = cell_idx_t(1);
+
+      if(filter(x,y,z) && filter(x+one,y,z) && filter(x-one,y,z) && filter(x,y+one,z)
+         && filter(x,y-one,z) && filter(x,y,z+one) && filter(x,y,z-one)) {
+         const Vector3<real_t> xa = velocityField.get(x+one,y,z);
+         const Vector3<real_t> xb = velocityField.get(x-one,y,z);
+         const Vector3<real_t> ya = velocityField.get(x,y+one,z);
+         const Vector3<real_t> yb = velocityField.get(x,y-one,z);
+         const Vector3<real_t> za = velocityField.get(x,y,z+one);
+         const Vector3<real_t> zb = velocityField.get(x,y,z-one);
+
+         return calculate(xa, xb, ya, yb, za, zb, dx, dy, dz);
+      }
+
+      return Vector3<real_t>(real_t(0));
+   }
+
+   static inline Vector3<real_t> calculate(const Vector3<real_t> xa, const Vector3<real_t> xb,
+           const Vector3<real_t> ya, const Vector3<real_t> yb,
+           const Vector3<real_t> za, const Vector3<real_t> zb,
+           const real_t dx, const real_t dy, const real_t dz) {
+
+      const auto halfInvDx = real_t(0.5) / dx;
+      const auto halfInvDy = real_t(0.5) / dy;
+      const auto halfInvDz = real_t(0.5) / dz;
+
+      const real_t duxdy = (ya[0] - yb[0]) * halfInvDy;
+      const real_t duxdz = (za[0] - zb[0]) * halfInvDz;
+
+      const real_t duydx = (xa[1] - xb[1]) * halfInvDx;
+      const real_t duydz = (za[1] - zb[1]) * halfInvDz;
+
+      const real_t duzdx = (xa[2] - xb[2]) * halfInvDx;
+      const real_t duzdy = (ya[2] - yb[2]) * halfInvDy;
+
+      const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
+
+      return curl;
+   }
+};
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/vtk/CurlMagnitude.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Velocity.h
+//! \ingroup lbm
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+#include "core/debug/Debug.h"
+#include "vtk/BlockCellDataWriter.h"
+
+
+namespace walberla {
+namespace lbm {
+
+/**
+ * lengthScaleWeight is only used to be able to derive quantities for the adaptive mesh refinement criterium of curl.
+ *
+ * @tparam VelocityField_T
+ * @tparam Filter_T
+ * @tparam OutputType
+ */
+template< typename VelocityField_T, typename Filter_T, typename OutputType = float >
+class CurlMagnitudeVTKWriter : public vtk::BlockCellDataWriter< OutputType, 1 >
+{
+public:
+   CurlMagnitudeVTKWriter(const shared_ptr<StructuredBlockStorage> blockStorage, Filter_T & filter,
+         const ConstBlockDataID & velocityFieldId, const std::string & id, const real_t lengthScaleWeight = real_t(-1)) :
+         vtk::BlockCellDataWriter< OutputType, 1 >(id), blockStorage_(blockStorage), filter_(filter),
+         velocityFieldId_(velocityFieldId), velocityField_(NULL), lengthScaleWeight_(lengthScaleWeight) {}
+
+protected:
+
+   void configure() {
+      WALBERLA_ASSERT_NOT_NULLPTR( this->block_ );
+      velocityField_ = this->block_->template getData< VelocityField_T >(velocityFieldId_ );
+   }
+
+   OutputType evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t /*f*/ )
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR(velocityField_ );
+
+      const real_t dx = blockStorage_->dx(blockStorage_->getLevel(*this->block_));
+      const real_t dy = blockStorage_->dy(blockStorage_->getLevel(*this->block_));
+      const real_t dz = blockStorage_->dz(blockStorage_->getLevel(*this->block_));
+
+      filter_(*this->block_);
+
+      const auto one = cell_idx_t(1);
+
+      const auto halfInvDx = real_t(0.5) * real_t(1) / dx;
+      const auto halfInvDy = real_t(0.5) * real_t(1) / dy;
+      const auto halfInvDz = real_t(0.5) * real_t(1) / dz;
+
+      const real_t lengthScale = std::cbrt(dx*dy*dz);
+      const real_t weightedLengthScale = !realIsIdentical(lengthScaleWeight_, real_t(-1)) ?
+            std::pow(lengthScale, (lengthScaleWeight_+1)/lengthScaleWeight_) : real_t(1);
+
+      if(filter_(x,y,z) && filter_(x+one,y,z) && filter_(x-one,y,z) && filter_(x,y+one,z)
+         && filter_(x,y-one,z) && filter_(x,y,z+one) && filter_(x,y,z-one)) {
+         const Vector3< real_t > xa = velocityField_->get(x+one,y,z);
+         const Vector3< real_t > xb = velocityField_->get(x-one,y,z);
+         const Vector3< real_t > ya = velocityField_->get(x,y+one,z);
+         const Vector3< real_t > yb = velocityField_->get(x,y-one,z);
+         const Vector3< real_t > za = velocityField_->get(x,y,z+one);
+         const Vector3< real_t > zb = velocityField_->get(x,y,z-one);
+
+         const real_t duxdy = (ya[0] - yb[0]) * halfInvDy;
+         const real_t duxdz = (za[0] - zb[0]) * halfInvDz;
+
+         const real_t duydx = (xa[1] - xb[1]) * halfInvDx;
+         const real_t duydz = (za[1] - zb[1]) * halfInvDz;
+
+         const real_t duzdx = (xa[2] - xb[2]) * halfInvDx;
+         const real_t duzdy = (ya[2] - yb[2]) * halfInvDy;
+
+         const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
+         const auto curlMag = curl.length();
+
+         const auto curlSensor = curlMag * weightedLengthScale;
+
+         return numeric_cast< OutputType >(curlSensor);
+      }
+
+      return numeric_cast< OutputType >(0);
+   }
+
+   const shared_ptr<StructuredBlockStorage> blockStorage_;
+   Filter_T filter_;
+   const ConstBlockDataID velocityFieldId_;
+   const VelocityField_T* velocityField_;
+   const real_t lengthScaleWeight_;
+};
+
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/vtk/QCriterion.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file QCriterion.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+#include "lbm/field/MacroscopicValueCalculation.h"
+#include "core/debug/Debug.h"
+#include "vtk/BlockCellDataWriter.h"
+
+
+namespace walberla {
+namespace lbm {
+
+
+template< typename VelocityField_T, typename Filter_T, typename OutputType = float >
+class QCriterionVTKWriter : public vtk::BlockCellDataWriter< OutputType, 1 >
+{
+public:
+   QCriterionVTKWriter(const shared_ptr<StructuredBlockStorage> blockStorage, Filter_T & filter,
+                   const ConstBlockDataID & velocityFieldId, const std::string & id ) :
+         vtk::BlockCellDataWriter< OutputType, 1 >(id), blockStorage_(blockStorage), filter_(filter), velocityFieldId_(velocityFieldId), velocityField_(NULL) {}
+
+protected:
+
+   void configure() {
+      WALBERLA_ASSERT_NOT_NULLPTR( this->block_ );
+      velocityField_ = this->block_->template getData< VelocityField_T >(velocityFieldId_ );
+   }
+
+   OutputType evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t /*f*/ )
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR(velocityField_ );
+
+      const real_t dx = blockStorage_->dx(blockStorage_->getLevel(*this->block_));
+      const real_t dy = blockStorage_->dy(blockStorage_->getLevel(*this->block_));
+      const real_t dz = blockStorage_->dz(blockStorage_->getLevel(*this->block_));
+
+      filter_(*this->block_);
+
+      return numeric_cast<OutputType>(getQCriterion(*velocityField_, filter_, x, y, z, dx, dy, dz));
+   }
+
+   const shared_ptr<StructuredBlockStorage> blockStorage_;
+   Filter_T filter_;
+   const ConstBlockDataID velocityFieldId_;
+   const VelocityField_T* velocityField_;
+};
+
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/vtk/QCriterionCellFilter.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file QValueCellFilter.h
+//! \ingroup field
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "field/FlagField.h"
+#include "core/cell/CellSet.h"
+#include "core/debug/Debug.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/field/MacroscopicValueCalculation.h"
+
+
+namespace walberla {
+namespace field {
+
+template< typename VelocityField_T, typename Filter_T >
+class QCriterionCellFilter {
+public:
+   QCriterionCellFilter(const ConstBlockDataID velocityFieldId, Filter_T & filter, real_t lowerLimit,
+         real_t upperLimit = std::numeric_limits<real_t>::infinity())
+         : velocityFieldId_(velocityFieldId), filter_(filter), lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
+
+   void operator()( CellSet& filteredCells, const IBlock& block, const StructuredBlockStorage& storage,
+         const uint_t ghostLayers = uint_t(0) ) {
+      const VelocityField_T* velocityField = block.getData< VelocityField_T >( velocityFieldId_ );
+      WALBERLA_CHECK_NOT_NULLPTR(velocityField);
+
+      const real_t dx = storage.dx(storage.getLevel(block));
+      const real_t dy = storage.dy(storage.getLevel(block));
+      const real_t dz = storage.dz(storage.getLevel(block));
+
+      filter_(block);
+
+      const cell_idx_t gl = cell_idx_c(ghostLayers);
+      const cell_idx_t begin = cell_idx_c(-1) * gl;
+
+      for (cell_idx_t z = begin; z < cell_idx_c(storage.getNumberOfZCells(block))+gl; ++z) {
+         for (cell_idx_t y = begin; y < cell_idx_c(storage.getNumberOfYCells(block))+gl; ++y) {
+            for (cell_idx_t x = begin; x < cell_idx_c(storage.getNumberOfXCells(block))+gl; ++x) {
+               real_t q = lbm::getQCriterion(*velocityField, filter_, x, y, z, dx, dy, dz);
+
+               if (q >= lowerLimit_ && q <= upperLimit_) {
+                  filteredCells.insert(x,y,z);
+               }
+            }
+         }
+      }
+   }
+
+private:
+   const ConstBlockDataID velocityFieldId_;
+   Filter_T filter_;
+   real_t lowerLimit_;
+   real_t upperLimit_;
+};
+
+} // namespace field
+} // namespace walberla

src/lbm/vtk/VelocityCellFilter.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file VectorFieldCellFilter.h
+//! \ingroup field
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "field/FlagField.h"
+#include "core/cell/CellSet.h"
+#include "core/debug/Debug.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "lbm/field/PdfField.h"
+
+
+namespace walberla {
+namespace field {
+
+template< typename PdfField_T, typename Filter_T >
+class VelocityCellFilter {
+public:
+   VelocityCellFilter(const ConstBlockDataID pdfFieldId, Filter_T & filter, real_t lowerLimit,
+         real_t upperLimit = std::numeric_limits<real_t>::infinity())
+         : pdfFieldId_(pdfFieldId), filter_(filter), lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
+
+   void operator()( CellSet& filteredCells, const IBlock& block, const StructuredBlockStorage& storage,
+         const uint_t ghostLayers = uint_t(0) ) {
+      const PdfField_T* pdf = block.getData< PdfField_T >( pdfFieldId_ );
+      WALBERLA_CHECK_NOT_NULLPTR(pdf);
+
+      filter_(block);
+
+      const cell_idx_t gl = cell_idx_c(ghostLayers);
+      const cell_idx_t begin = cell_idx_c(-1) * gl;
+
+      for (cell_idx_t z = begin; z < cell_idx_c(storage.getNumberOfZCells(block))+gl; ++z) {
+         for (cell_idx_t y = begin; y < cell_idx_c(storage.getNumberOfYCells(block))+gl; ++y) {
+            for (cell_idx_t x = begin; x < cell_idx_c(storage.getNumberOfXCells(block))+gl; ++x) {
+               if (filter_(x,y,z)) {
+                  const Vector3<real_t> value = pdf->getVelocity(x,y,z);
+                  const real_t sqrLength = value.sqrLength();
+                  if (sqrLength > lowerLimit_*lowerLimit_
+                     && (realIsIdentical(upperLimit_, std::numeric_limits<real_t>::infinity())
+                        || sqrLength < upperLimit_*upperLimit_)) {
+                     filteredCells.insert(x, y, z);
+                  }
+               }
+            }
+         }
+      }
+   }
+
+private:
+   const ConstBlockDataID pdfFieldId_;
+   Filter_T filter_;
+   real_t lowerLimit_;
+   real_t upperLimit_;
+};
+
+} // namespace field
+} // namespace walberla

src/lbm/vtk/Vorticity.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Vorticity.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/PdfField.h"
+#include "lbm/field/MacroscopicValueCalculation.h"
+#include "core/debug/Debug.h"
+#include "vtk/BlockCellDataWriter.h"
+
+
+namespace walberla {
+namespace lbm {
+
+
+template< typename VelocityField_T, typename Filter_T, typename OutputType = float >
+class VorticityComponentVTKWriter : public vtk::BlockCellDataWriter< OutputType, 1 >
+{
+public:
+   VorticityComponentVTKWriter(const shared_ptr<StructuredBlockStorage> blockStorage, Filter_T & filter,
+         const ConstBlockDataID & velocityFieldId, const uint_t componentIdx, const std::string & id,
+         const real_t normalizationConstant = real_t(1)) :
+         vtk::BlockCellDataWriter< OutputType, 1 >(id), blockStorage_(blockStorage), filter_(filter),
+         velocityFieldId_(velocityFieldId), componentIdx_(componentIdx), velocityField_(NULL),
+         normalizationConstant_(normalizationConstant) {
+      WALBERLA_ASSERT(componentIdx < uint_t(3),
+            "The vorticity vector only has three components, i.e. the highest possible component index is 2.");
+   }
+
+protected:
+
+   void configure() {
+      WALBERLA_ASSERT_NOT_NULLPTR( this->block_ );
+      velocityField_ = this->block_->template getData< VelocityField_T >(velocityFieldId_ );
+   }
+
+   OutputType evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t /*f*/ ) {
+      WALBERLA_ASSERT_NOT_NULLPTR(velocityField_ );
+
+      const real_t dx = blockStorage_->dx(blockStorage_->getLevel(*this->block_));
+      const real_t dy = blockStorage_->dy(blockStorage_->getLevel(*this->block_));
+      const real_t dz = blockStorage_->dz(blockStorage_->getLevel(*this->block_));
+
+      filter_(*this->block_);
+
+      Vector3<real_t> curl = getVorticity(*velocityField_, filter_, x, y, z, dx, dy, dz);
+
+      return numeric_cast<OutputType>(curl[componentIdx_]/normalizationConstant_);
+   }
+
+   const shared_ptr<StructuredBlockStorage> blockStorage_;
+   Filter_T filter_;
+   const ConstBlockDataID velocityFieldId_;
+   const uint_t componentIdx_;
+   const VelocityField_T* velocityField_;
+   const real_t normalizationConstant_;
+};
+
+
+
+} // namespace lbm
+} // namespace walberla

src/lbm/vtk/all.h

@@ -27,5 +27,9 @@
 #include "Velocity.h"
 #include "VTKOutput.h"
 #include "PressureTensor.h"
-
+#include "QCriterion.h"
+#include "Vorticity.h"
+#include "CurlMagnitude.h"
+#include "QCriterionCellFilter.h"
+#include "VelocityCellFilter.h"
 

tests/lbm/CMakeLists.txt

@@ -71,6 +71,9 @@ waLBerla_compile_test( FILES SuViscoelasticityTest.cpp DEPENDS field blockforest
 waLBerla_execute_test( NAME  SuViscoelasticityLongTest COMMAND $<TARGET_FILE:SuViscoelasticityTest> ${CMAKE_CURRENT_SOURCE_DIR}/Su.prm LABELS longrun)
 waLBerla_execute_test( NAME  SuViscoelasticityShortTest COMMAND $<TARGET_FILE:SuViscoelasticityTest> ${CMAKE_CURRENT_SOURCE_DIR}/Su.prm -Parameters.shortrun=true )
 
+waLBerla_compile_test( FILES field/QCriterionTest.cpp DEPENDS field blockforest )
+waLBerla_execute_test( NAME QCriterionTest )
+
 # Code Generation
 if( WALBERLA_BUILD_WITH_CODEGEN )
 

tests/lbm/field/QCriterionTest.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file QCriterionTest.cpp
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+
+#include "blockforest/all.h"
+#include "core/all.h"
+#include "field/all.h"
+#include "lbm/all.h"
+
+namespace walberla {
+
+class Filter {
+public:
+   explicit Filter(uint_t numberOfCells) : numberOfCells_(numberOfCells) {}
+
+   void operator()( const IBlock & /*block*/ ){
+
+   }
+
+   bool operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const {
+      return x >= -1 && x <= cell_idx_t(numberOfCells_) &&
+             y >= -1 && y <= cell_idx_t(numberOfCells_) &&
+             z >= -1 && z <= cell_idx_t(numberOfCells_);
+   }
+
+private:
+   uint_t numberOfCells_;
+};
+
+using VelocityField_T = GhostLayerField<Vector3<real_t>, 1>;
+using FluidFilter_T = Filter;
+
+int main( int argc, char ** argv )
+{
+   debug::enterTestMode();
+   Environment env( argc, argv );
+
+   auto numberOfCells = uint_t(40);
+
+   VelocityField_T velocityField(numberOfCells, numberOfCells, numberOfCells, uint_t(1));
+
+   FluidFilter_T filter(numberOfCells);
+
+   //// initialize field
+
+   real_t xv = 0, yv = 0, zv = 0;
+
+   for (auto cellIt = velocityField.beginWithGhostLayerXYZ(); cellIt != velocityField.end(); ++cellIt) {
+      xv = math::realRandom<real_t>(real_t(-1), real_t(1));
+      yv = math::realRandom<real_t>(real_t(-1), real_t(1));
+      zv = math::realRandom<real_t>(real_t(-1), real_t(1));
+
+      *cellIt = Vector3<real_t>(real_t(xv), real_t(yv), real_t(zv));
+   }
+
+   //// evaluate field
+
+   const auto one = cell_idx_t(1);
+   for (auto cellIt = velocityField.beginXYZ(); cellIt != velocityField.end(); ++cellIt) {
+      cell_idx_t x = cellIt.x();
+      cell_idx_t y = cellIt.y();
+      cell_idx_t z = cellIt.z();
+
+      // ParaView's formula serves as comparison.
+      // From: https://github.com/Kitware/VTK/blob/8b08cf7a0523d88a5602a98c0456f7e397560698/Filters/General/vtkGradientFilter.cxx
+
+      const Vector3<real_t> xa = velocityField.get(x+one,y,z);
+      const Vector3<real_t> xb = velocityField.get(x-one,y,z);
+      const Vector3<real_t> ya = velocityField.get(x,y+one,z);
+      const Vector3<real_t> yb = velocityField.get(x,y-one,z);
+      const Vector3<real_t> za = velocityField.get(x,y,z+one);
+      const Vector3<real_t> zb = velocityField.get(x,y,z-one);
+
+      const real_t duxdx = (xa[0] - xb[0]) * real_t(0.5);
+      const real_t duxdy = (ya[0] - yb[0]) * real_t(0.5);
+      const real_t duxdz = (za[0] - zb[0]) * real_t(0.5);
+
+      const real_t duydx = (xa[1] - xb[1]) * real_t(0.5);
+      const real_t duydy = (ya[1] - yb[1]) * real_t(0.5);
+      const real_t duydz = (za[1] - zb[1]) * real_t(0.5);
+
+      const real_t duzdx = (xa[2] - xb[2]) * real_t(0.5);
+      const real_t duzdy = (ya[2] - yb[2]) * real_t(0.5);
+      const real_t duzdz = (za[2] - zb[2]) * real_t(0.5);
+
+      real_t q_paraview = -real_t(0.5)*(duxdx*duxdx + duydy*duydy + duzdz*duzdz)
+                          -(duxdy*duydx + duxdz*duzdx + duydz*duzdy);
+
+      WALBERLA_CHECK_FLOAT_EQUAL(q_paraview, lbm::getQCriterion(velocityField, filter, x, y, z));
+   }
+
+   return 0;
+}
+}
+
+int main( int argc, char ** argv )
+{
+   return walberla::main(argc, argv);
+}
\ No newline at end of file