diff --git a/apps/tutorials/lbm/01_BasicLBM.cpp b/apps/tutorials/lbm/01_BasicLBM.cpp
index 1380a1f10e052d8f09032cda1e943698a193d668..ada14422b7ea0d2ac4b73e93ec206621ad97b19b 100644
--- a/apps/tutorials/lbm/01_BasicLBM.cpp
+++ b/apps/tutorials/lbm/01_BasicLBM.cpp
@@ -28,12 +28,20 @@
#include "gui/all.h"
#include "lbm/all.h"
#include "timeloop/all.h"
+#include "field/interpolators/TrilinearFieldInterpolator.h"
+#include "field/interpolators/FieldInterpolatorCreators.h"
+#include "field/distributors/DistributorCreators.h"
+#include "field/distributors/KernelDistributor.h"
+#include "CustomLatticeModel.h"
+#include "TurbineModel.h"
namespace walberla {
-using LatticeModel_T = lbm::D2Q9<lbm::collision_model::SRT>;
+//using LatticeModel_T = lbm::D3Q19<lbm::collision_model::SRT>;
+using LatticeModel_T = lbm::CustomLatticeModel;
+
using Stencil_T = LatticeModel_T::Stencil;
using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
@@ -42,7 +50,10 @@ using PdfField_T = lbm::PdfField<LatticeModel_T>;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField<flag_t>;
+using VectorField_T = GhostLayerField<real_t, 3>;
+using Vec3Interpolator_T = field::TrilinearFieldInterpolator<VectorField_T, FlagField_T>;
+using Vec3Distributor_T = field::KernelDistributor<VectorField_T, FlagField_T>;
int main( int argc, char ** argv )
{
@@ -52,6 +63,7 @@ int main( int argc, char ** argv )
// read parameters
auto parameters = walberlaEnv.config()->getOneBlock( "Parameters" );
+ auto turbineCfg = walberlaEnv.config()->getOneBlock( "Turbine" );
const real_t omega = parameters.getParameter< real_t > ( "omega", real_c( 1.4 ) );
const Vector3<real_t> initialVelocity = parameters.getParameter< Vector3<real_t> >( "initialVelocity", Vector3<real_t>() );
@@ -60,10 +72,15 @@ int main( int argc, char ** argv )
const double remainingTimeLoggerFrequency = parameters.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); // in seconds
// create fields
- LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::SRT( omega ) );
+ //LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::SRT( omega ) );
+ BlockDataID velocityFieldId = field::addToStorage<VectorField_T>( blocks, "vel field");
+ BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "force field");
+
+ LatticeModel_T latticeModel( forceFieldId, velocityFieldId, omega );
BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_t(1) );
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+
// create and initialize boundary handling
const FlagUID fluidFlagUID( "Fluid" );
@@ -80,7 +97,25 @@ int main( int argc, char ** argv )
geometry::initBoundaryHandling<BHFactory::BoundaryHandling>( *blocks, boundaryHandlingId, boundariesConfig );
geometry::setNonBoundaryCellsToDomain<BHFactory::BoundaryHandling> ( *blocks, boundaryHandlingId );
- // create time loop
+
+ BlockDataID velocityFieldInterpolator = field::addFieldInterpolator<Vec3Interpolator_T, FlagField_T>(
+ blocks, velocityFieldId, flagFieldId, fluidFlagUID, "Velocity Interpolator" );
+ BlockDataID forceFieldDistributor = field::addDistributor<Vec3Distributor_T, FlagField_T>(
+ blocks, forceFieldId, flagFieldId, fluidFlagUID, "Force distributor"
+ );
+
+
+ Vector3<real_t> midPoint = turbineCfg.getParameter< Vector3<real_t> >("midPoint");
+ const real_t innerRadius = turbineCfg.getParameter<real_t>("innerRadius");
+ const real_t outerRadius = turbineCfg.getParameter<real_t>("outerRadius");
+ const real_t spacing = turbineCfg.getParameter<real_t>("spacing");
+ const real_t drag = turbineCfg.getParameter<real_t>("drag", 1e-2);
+ const real_t lift = turbineCfg.getParameter<real_t>("lift", 1.0);
+ TurbineModelLine<Vec3Interpolator_T, Vec3Distributor_T> turbineModel(blocks, midPoint,
+ innerRadius, outerRadius, spacing,
+ drag, lift);
+
+ // create time loop
SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
// create communication for PdfField
@@ -90,7 +125,17 @@ int main( int argc, char ** argv )
// add LBM sweep and communication to time loop
timeloop.add() << BeforeFunction( communication, "communication" )
<< Sweep( BHFactory::BoundaryHandling::getBlockSweep( boundaryHandlingId ), "boundary handling" );
- timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID ) ), "LB stream & collide" );
+ timeloop.add() << Sweep( LatticeModel_T::Sweep(pdfFieldId ) )
+ << AfterFunction( [&]() {
+
+ for( auto & block : *blocks ) {
+ auto forceField = block.getData<VectorField_T>( forceFieldId );
+ forceField->set( 0.0 );
+ }
+ turbineModel.rotate( 2 * math::M_PI / 360.0 );
+ turbineModel.evaluate(velocityFieldInterpolator, forceFieldDistributor);
+ } );
+ //timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID ) ), "LB stream & collide" );
// LBM stability check
timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >( walberlaEnv.config(), blocks, pdfFieldId,
diff --git a/apps/tutorials/lbm/01_BasicLBM.prm b/apps/tutorials/lbm/01_BasicLBM.prm
index 15efce2094b505f5063d99fe84e62847161c9db7..7f9d0094569b72b81823cd123dd95a2de76b4808 100644
--- a/apps/tutorials/lbm/01_BasicLBM.prm
+++ b/apps/tutorials/lbm/01_BasicLBM.prm
@@ -1,24 +1,37 @@
+DomainSetup
+{
+ blocks < 1, 1, 1 >;
+ cellsPerBlock < 120, 30, 30 >;
+ periodic < 0, 0, 0 >;
+}
+
+
+Turbine
+{
+ midpoint < 20, 15, 15>;
+ innerRadius 1;
+ outerRadius 8;
+ spacing 1;
+ drag 1e-2;
+ lift 10;
+}
+
Parameters
{
omega 1.8;
- initialVelocity < 0.1, 0, 0 >;
+ initialVelocity < 0.05, 0, 0 >;
timesteps 10000;
- useGui 1;
- remainingTimeLoggerFrequency 3; // in seconds
+ useGui 0;
+ remainingTimeLoggerFrequency 5; // in seconds
}
-DomainSetup
-{
- blocks < 1, 1, 1 >;
- cellsPerBlock < 300, 80, 1 >;
- periodic < 0, 0, 1 >;
-}
+
StabilityChecker
{
- checkFrequency 1;
+ checkFrequency 100;
streamOutput false;
vtkOutput true;
}
@@ -39,23 +52,10 @@ Boundaries
- GrayScaleImage BoundaryFromImage.h
- RGBAImage BoundaryFromImage.h
*/
-
- Border { direction W; walldistance -1; Velocity0 {} }
- Border { direction E; walldistance -1; Pressure1 {} }
- Border { direction S,N; walldistance -1; NoSlip {} }
-
- // load obstacle geometry from file
- GrayScaleImage
- {
- file waLBerla.png;
- extrusionCoordinate 2;
- rescaleToDomain true;
-
- BoundaryValueMapping {
- NoSlip {}
- value 0;
- }
- }
+
+ Border { direction S,N,T,B; walldistance -1; FreeSlip {} }
+ Border { direction W; walldistance -1; Velocity0 {} }
+ Border { direction E; walldistance -1; Pressure1 {} }
}
diff --git a/apps/tutorials/lbm/CMakeLists.txt b/apps/tutorials/lbm/CMakeLists.txt
index 6a36280444bc9db6ba2e96606a4d594f12fe52f2..78ea6b299ce998eff260556699fde21f3a6241da 100644
--- a/apps/tutorials/lbm/CMakeLists.txt
+++ b/apps/tutorials/lbm/CMakeLists.txt
@@ -1,15 +1,19 @@
waLBerla_link_files_to_builddir( *.prm )
waLBerla_link_files_to_builddir( *.png )
+waLBerla_python_file_generates(CustomLatticeModel.py
+ CustomLatticeModel.cpp CustomLatticeModel.h)
+
+
waLBerla_add_executable ( NAME 01_BasicLBM
- FILES 01_BasicLBM.cpp
+ FILES 01_BasicLBM.cpp CustomLatticeModel.py
DEPENDS blockforest core field lbm geometry timeloop gui )
waLBerla_add_executable ( NAME 02_BasicLBM_ExemplaryExtensions
FILES 02_BasicLBM_ExemplaryExtensions.cpp
DEPENDS blockforest core field lbm geometry timeloop gui )
-
-waLBerla_add_executable ( NAME 03_LBLidDrivenCavity
+
+waLBerla_add_executable ( NAME 03_LBLidDrivenCavity
FILES 03_LBLidDrivenCavity.cpp
DEPENDS blockforest boundary core field lbm stencil timeloop vtk )
diff --git a/apps/tutorials/lbm/CustomLatticeModel.py b/apps/tutorials/lbm/CustomLatticeModel.py
new file mode 100644
index 0000000000000000000000000000000000000000..01a858c8c9efa13aadf7e7f34dca700d657e2fa3
--- /dev/null
+++ b/apps/tutorials/lbm/CustomLatticeModel.py
@@ -0,0 +1,25 @@
+import sympy as sp
+import pystencils as ps
+from lbmpy.creationfunctions import create_lb_method
+from pystencils_walberla import CodeGeneration
+from lbmpy_walberla import generate_lattice_model, RefinementScaling, generate_boundary
+
+
+with CodeGeneration() as ctx:
+ omega = sp.Symbol("omega")
+ velocity_field = ps.fields("velocity(3) : [3D]", layout='fzyx')
+ force_field = ps.fields("force(3): [3D]", layout='fzyx')
+
+ # lattice Boltzmann method
+ params = {
+ 'stencil': 'D3Q19',
+ 'method': 'srt',
+ 'relaxation_rates': [omega],
+ 'smagorinsky': True,
+ 'output': {'velocity': velocity_field},
+ 'force_model': 'guo',
+ 'force': force_field.center_vector,
+ }
+ lb_method = create_lb_method(**params)
+ generate_lattice_model(ctx, 'CustomLatticeModel', lb_method,
+ update_rule_params=params)
diff --git a/apps/tutorials/lbm/TurbineModel.h b/apps/tutorials/lbm/TurbineModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..a70b1c6ccade7efa96ace5cca69f46dad117e9b4
--- /dev/null
+++ b/apps/tutorials/lbm/TurbineModel.h
@@ -0,0 +1,113 @@
+//======================================================================================================================
+//
+// 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 TurbineModel.h
+//
+//======================================================================================================================
+
+#pragma once
+#include "core/math/Matrix3.h"
+#include "core/math/Vector3.h"
+#include "core/math/Rot3.h"
+
+
+namespace walberla {
+
+template<typename VelocityInterpolator_T, typename ForceDistributor_T>
+class TurbineModelLine
+{
+public:
+ TurbineModelLine(const shared_ptr <StructuredBlockStorage> &blocks,
+ const Vector3 <real_t> midPoint, real_t innerRadius, real_t outerRadius,
+ real_t pointSpacing, real_t drag, real_t lift)
+ : blocks_(blocks), midPoint_(midPoint), pointSpacing_(pointSpacing),
+ drag_(drag), lift_(lift) {
+ // initialize in z-direction
+ for (real_t curZ = -outerRadius; curZ <= -innerRadius; curZ += pointSpacing) {
+ points_.push_back(Vector3<real_t>(real_t(0), real_t(0), curZ) + midPoint_);
+ }
+ for (real_t curZ = innerRadius; curZ <= outerRadius; curZ += pointSpacing) {
+ points_.push_back(Vector3<real_t>(real_t(0), real_t(0), curZ) + midPoint_);
+ }
+ }
+
+ // compute forces on fluid, sets them into forceField
+ void evaluate(ConstBlockDataID velocityInterpolatorId, BlockDataID forceDistributorId)
+ {
+ for(auto & block : *blocks_)
+ {
+ auto velocity = block.getData<VelocityInterpolator_T>(velocityInterpolatorId);
+ auto force = block.getData<ForceDistributor_T>(forceDistributorId);
+ auto blockBB = block.getAABB();
+ for( const auto & point : points_ )
+ {
+ if( blockBB.contains(point) )
+ {
+ Vector3<real_t> freeStreamVelocity;
+ velocity->get( point, freeStreamVelocity.data() );
+ auto velMagnitude = freeStreamVelocity.length();
+ auto angleOfAttack = real_t(0.0);
+
+ auto drag = real_t(0);
+ auto lift = real_t(0);
+ airfoilTable(velMagnitude, angleOfAttack, drag, lift);
+
+ auto wingVec = (point - midPoint_).getNormalized();
+
+ auto area = pointSpacing_;
+
+ auto forceToDistribute = Vector3<real_t>(
+ - 0.5 * drag * area * freeStreamVelocity[0] * freeStreamVelocity[0],
+ - 0.5 * (lift * wingVec[2]) * area * freeStreamVelocity[0] * freeStreamVelocity[0],
+ - 0.5 * (-lift * wingVec[1]) * area * freeStreamVelocity[0] * freeStreamVelocity[0]
+ );
+ force->distribute( point, forceToDistribute.data() );
+ }
+ //TODO communicate force field with "+=" communication
+ }
+ }
+ }
+
+ void airfoilTable(real_t velocityMagnitude, real_t angleOfAttack,
+ real_t &drag, real_t &lift)
+ {
+ drag = drag_;
+ lift = lift_;
+ }
+
+ void rotate(real_t angle)
+ {
+ const auto rotationAxis = Vector3<real_t>(real_t(1.0), real_t(0.0), real_t(0.0));
+ auto rotationMatrix = Matrix3<real_t>(rotationAxis, angle);
+ for (auto &point : points_) {
+ point -= midPoint_;
+ point = rotationMatrix * point;
+ point += midPoint_;
+ }
+ absoluteAngle_ += angle;
+ }
+
+private:
+ std::vector<Vector3 < real_t> > points_;
+ shared_ptr <StructuredBlockStorage> blocks_;
+ Vector3 <real_t> midPoint_;
+ real_t absoluteAngle_;
+ real_t pointSpacing_;
+ real_t drag_;
+ real_t lift_;
+};
+
+
+} // namespace walberla
\ No newline at end of file