From 081b9e08ab22e66126d294643b7cdae66d2301af Mon Sep 17 00:00:00 2001
From: Lukas Werner <lks.werner@fau.de>
Date: Fri, 16 Oct 2020 09:50:08 +0200
Subject: [PATCH] LBM: Q criterion, curl and vorticity calculation and vtk
output
---
src/lbm/field/MacroscopicValueCalculation.h | 17 +++
src/lbm/field/QCriterion.h | 99 +++++++++++++++++
src/lbm/field/QCriterionFieldWriter.h | 78 +++++++++++++
src/lbm/field/Vorticity.h | 89 +++++++++++++++
src/lbm/vtk/CurlMagnitude.h | 114 +++++++++++++++++++
src/lbm/vtk/QCriterion.h | 72 ++++++++++++
src/lbm/vtk/QCriterionCellFilter.h | 77 +++++++++++++
src/lbm/vtk/VelocityCellFilter.h | 76 +++++++++++++
src/lbm/vtk/Vorticity.h | 81 ++++++++++++++
src/lbm/vtk/all.h | 6 +-
tests/lbm/CMakeLists.txt | 3 +
tests/lbm/field/QCriterionTest.cpp | 116 ++++++++++++++++++++
12 files changed, 827 insertions(+), 1 deletion(-)
create mode 100644 src/lbm/field/QCriterion.h
create mode 100644 src/lbm/field/QCriterionFieldWriter.h
create mode 100644 src/lbm/field/Vorticity.h
create mode 100644 src/lbm/vtk/CurlMagnitude.h
create mode 100644 src/lbm/vtk/QCriterion.h
create mode 100644 src/lbm/vtk/QCriterionCellFilter.h
create mode 100644 src/lbm/vtk/VelocityCellFilter.h
create mode 100644 src/lbm/vtk/Vorticity.h
create mode 100644 tests/lbm/field/QCriterionTest.cpp
diff --git a/src/lbm/field/MacroscopicValueCalculation.h b/src/lbm/field/MacroscopicValueCalculation.h
index dae0a43cb..50059b370 100644
--- a/src/lbm/field/MacroscopicValueCalculation.h
+++ b/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
diff --git a/src/lbm/field/QCriterion.h b/src/lbm/field/QCriterion.h
new file mode 100644
index 000000000..c399218de
--- /dev/null
+++ b/src/lbm/field/QCriterion.h
@@ -0,0 +1,99 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/field/QCriterionFieldWriter.h b/src/lbm/field/QCriterionFieldWriter.h
new file mode 100644
index 000000000..9948ea9ec
--- /dev/null
+++ b/src/lbm/field/QCriterionFieldWriter.h
@@ -0,0 +1,78 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/field/Vorticity.h b/src/lbm/field/Vorticity.h
new file mode 100644
index 000000000..7766ddb3b
--- /dev/null
+++ b/src/lbm/field/Vorticity.h
@@ -0,0 +1,89 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/CurlMagnitude.h b/src/lbm/vtk/CurlMagnitude.h
new file mode 100644
index 000000000..0aa89feab
--- /dev/null
+++ b/src/lbm/vtk/CurlMagnitude.h
@@ -0,0 +1,114 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/QCriterion.h b/src/lbm/vtk/QCriterion.h
new file mode 100644
index 000000000..e3b2748ec
--- /dev/null
+++ b/src/lbm/vtk/QCriterion.h
@@ -0,0 +1,72 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/QCriterionCellFilter.h b/src/lbm/vtk/QCriterionCellFilter.h
new file mode 100644
index 000000000..5310bc172
--- /dev/null
+++ b/src/lbm/vtk/QCriterionCellFilter.h
@@ -0,0 +1,77 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/VelocityCellFilter.h b/src/lbm/vtk/VelocityCellFilter.h
new file mode 100644
index 000000000..0e290db21
--- /dev/null
+++ b/src/lbm/vtk/VelocityCellFilter.h
@@ -0,0 +1,76 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/Vorticity.h b/src/lbm/vtk/Vorticity.h
new file mode 100644
index 000000000..237b711d5
--- /dev/null
+++ b/src/lbm/vtk/Vorticity.h
@@ -0,0 +1,81 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/vtk/all.h b/src/lbm/vtk/all.h
index 07ab8b35a..7f4285777 100644
--- a/src/lbm/vtk/all.h
+++ b/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"
diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt
index b9631921b..b17f65cfa 100644
--- a/tests/lbm/CMakeLists.txt
+++ b/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 )
diff --git a/tests/lbm/field/QCriterionTest.cpp b/tests/lbm/field/QCriterionTest.cpp
new file mode 100644
index 000000000..0e5dc2dd0
--- /dev/null
+++ b/tests/lbm/field/QCriterionTest.cpp
@@ -0,0 +1,116 @@
+//======================================================================================================================
+//
+// 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
--
GitLab