diff --git a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
index 9542fad39b7963d76d4c6a7f1a5280ea0ef2dd04..2acb0a54620bef5d8ab2fce03f23255cabc08c91 100644
--- a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
+++ b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
@@ -608,7 +608,7 @@ protected:
const auto & dy = this->blockStorage_->dy( level );
WALBERLA_ASSERT_FLOAT_EQUAL( dy, this->blockStorage_->dz( level ) );
- const real_t acceleration_L = pdf_->latticeModel().forceModel().force()[0]; // force = acceleration (in lattice units of the current level!)
+ const real_t acceleration_L = pdf_->latticeModel().forceModel().forceDensity()[0]; // force = acceleration (in lattice units of the current level!)
const real_t viscosity_L = pdf_->latticeModel().collisionModel().viscosity(); // in lattice units on the current level
const real_t radius_L = channelRadius / dy; // in lattice units on the current level
diff --git a/src/lbm/field/DensityAndVelocity.h b/src/lbm/field/DensityAndVelocity.h
index 540357d565b132ebe23d16a40345cf0847130fca..91ce577d1cf163d3d96ab26091202482dcf01973 100644
--- a/src/lbm/field/DensityAndVelocity.h
+++ b/src/lbm/field/DensityAndVelocity.h
@@ -53,7 +53,7 @@ struct AdaptVelocityToForce< LatticeModel_T, typename std::enable_if< LatticeMod
{
static Vector3<real_t> get( const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t rho )
{
- return velocity - latticeModel.forceModel().force() * real_t(0.5) / rho;
+ return velocity - latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
}
template< typename FieldPtrOrIterator >
@@ -77,7 +77,7 @@ struct AdaptVelocityToForce< LatticeModel_T, typename std::enable_if< ! LatticeM
{
static Vector3<real_t> get( const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t )
{
- return velocity - latticeModel.forceModel().force() * real_t(0.5);
+ return velocity - latticeModel.forceModel().forceDensity() * real_t(0.5);
}
template< typename FieldPtrOrIterator >
@@ -102,20 +102,20 @@ struct AdaptVelocityToForce< LatticeModel_T, typename std::enable_if< LatticeMod
/*
static Vector3<real_t> get( const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t rho )
{
- return velocity - latticeModel.forceModel().force() * real_t(0.5) / rho;
+ return velocity - latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
}
*/
template< typename FieldPtrOrIterator >
static Vector3<real_t> get( FieldPtrOrIterator & it, const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t rho )
{
- return velocity - latticeModel.forceModel().force(it.x(),it.y(),it.z()) * real_t(0.5) / rho;
+ return velocity - latticeModel.forceModel().forceDensity(it.x(),it.y(),it.z()) * real_t(0.5) / rho;
}
static Vector3<real_t> get( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const LatticeModel_T & latticeModel,
const Vector3< real_t > & velocity, const real_t rho )
{
- return velocity - latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+ return velocity - latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
}
};
@@ -128,20 +128,20 @@ struct AdaptVelocityToForce< LatticeModel_T, typename std::enable_if< ! LatticeM
/*
static Vector3<real_t> get( const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t )
{
- return velocity - latticeModel.forceModel().force() * real_t(0.5);
+ return velocity - latticeModel.forceModel().forceDensity() * real_t(0.5);
}
*/
template< typename FieldPtrOrIterator >
static Vector3<real_t> get( FieldPtrOrIterator & it, const LatticeModel_T & latticeModel, const Vector3< real_t > & velocity, const real_t )
{
- return velocity - latticeModel.forceModel().force(it.x(),it.y(),it.z()) * real_t(0.5);
+ return velocity - latticeModel.forceModel().forceDensity(it.x(),it.y(),it.z()) * real_t(0.5);
}
static Vector3<real_t> get( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const LatticeModel_T & latticeModel,
const Vector3< real_t > & velocity, const real_t )
{
- return velocity - latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+ return velocity - latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
}
};
diff --git a/src/lbm/field/DensityVelocityCallback.h b/src/lbm/field/DensityVelocityCallback.h
index d7c76c5ca808b898659e6998dea85a7b94b1aa36..acf99efbe6bb19706bbdda5ca0372fd2229292f2 100644
--- a/src/lbm/field/DensityVelocityCallback.h
+++ b/src/lbm/field/DensityVelocityCallback.h
@@ -75,7 +75,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< Latt
static void apply( Vector3< real_t > & velocity, const cell_idx_t, const cell_idx_t, const cell_idx_t,
const LatticeModel_T & latticeModel, const real_t rho )
{
- velocity += latticeModel.forceModel().force() * real_t(0.5) / rho;
+ velocity += latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
}
};
@@ -89,7 +89,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( Vector3< real_t > & velocity, const cell_idx_t, const cell_idx_t, const cell_idx_t,
const LatticeModel_T & latticeModel, const real_t )
{
- velocity += latticeModel.forceModel().force() * real_t(0.5);
+ velocity += latticeModel.forceModel().forceDensity() * real_t(0.5);
}
};
@@ -103,7 +103,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< Latt
static void apply( Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const LatticeModel_T & latticeModel, const real_t rho )
{
- velocity += latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+ velocity += latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
}
};
@@ -117,7 +117,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const LatticeModel_T & latticeModel, const real_t )
{
- velocity += latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+ velocity += latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
}
};
@@ -131,7 +131,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< Latt
static void apply( Vector3< real_t > & velocity, const cell_idx_t, const cell_idx_t, const cell_idx_t,
const LatticeModel_T & latticeModel, const real_t rho )
{
- velocity -= latticeModel.forceModel().force() * real_t(0.5) / rho;
+ velocity -= latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
}
};
@@ -145,7 +145,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( Vector3< real_t > & velocity, const cell_idx_t, const cell_idx_t, const cell_idx_t,
const LatticeModel_T & latticeModel, const real_t )
{
- velocity -= latticeModel.forceModel().force() * real_t(0.5);
+ velocity -= latticeModel.forceModel().forceDensity() * real_t(0.5);
}
};
@@ -159,7 +159,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< Latt
static void apply( Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const LatticeModel_T & latticeModel, const real_t rho )
{
- velocity -= latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+ velocity -= latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
}
};
@@ -173,7 +173,7 @@ struct VelocityCallbackCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const LatticeModel_T & latticeModel, const real_t )
{
- velocity -= latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+ velocity -= latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
}
};
@@ -495,7 +495,7 @@ private:
// {
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity + storage;
-// storage = latticeModel.forceModel().force() * real_t(0.5) / rho;
+// storage = latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
// }
//};
//
@@ -512,7 +512,7 @@ private:
// {
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity + storage;
-// storage = latticeModel.forceModel().force() * real_t(0.5);
+// storage = latticeModel.forceModel().forceDensity() * real_t(0.5);
// }
//};
//
@@ -529,7 +529,7 @@ private:
// {
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity + storage;
-// storage = latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+// storage = latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
// }
//};
//
@@ -546,7 +546,7 @@ private:
// {
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity + storage;
-// storage = latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+// storage = latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
// }
//};
//
@@ -561,7 +561,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t rho )
// {
-// field->get(x,y,z) = velocity - latticeModel.forceModel().force() * real_t(0.5) / rho;
+// field->get(x,y,z) = velocity - latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
// }
//};
//
@@ -576,7 +576,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t )
// {
-// field->get(x,y,z) = velocity - latticeModel.forceModel().force() * real_t(0.5);
+// field->get(x,y,z) = velocity - latticeModel.forceModel().forceDensity() * real_t(0.5);
// }
//};
//
@@ -591,7 +591,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t rho )
// {
-// field->get(x,y,z) = velocity - latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+// field->get(x,y,z) = velocity - latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
// }
//};
//
@@ -606,7 +606,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t )
// {
-// field->get(x,y,z) = velocity - latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+// field->get(x,y,z) = velocity - latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
// }
//};
//
@@ -621,7 +621,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t rho )
// {
-// const auto current = latticeModel.forceModel().force() * real_t(0.5) / rho;
+// const auto current = latticeModel.forceModel().forceDensity() * real_t(0.5) / rho;
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity - current + storage;
// storage = current;
@@ -639,7 +639,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t )
// {
-// const auto current = latticeModel.forceModel().force() * real_t(0.5);
+// const auto current = latticeModel.forceModel().forceDensity() * real_t(0.5);
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity - current + storage;
// storage = current;
@@ -657,7 +657,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t rho )
// {
-// const auto current = latticeModel.forceModel().force(x,y,z) * real_t(0.5) / rho;
+// const auto current = latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5) / rho;
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity - current + storage;
// storage = current;
@@ -675,7 +675,7 @@ private:
// static void apply( VelocityField_T * const & field, const Vector3< real_t > & velocity, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
// const LatticeModel_T & latticeModel, const real_t )
// {
-// const auto current = latticeModel.forceModel().force(x,y,z) * real_t(0.5);
+// const auto current = latticeModel.forceModel().forceDensity(x,y,z) * real_t(0.5);
// auto & storage = field->get(x,y,z,1);
// field->get(x,y,z) = velocity - current + storage;
// storage = current;
diff --git a/src/lbm/field/MomentumDensity.h b/src/lbm/field/MomentumDensity.h
index 0a79a0fdca59bf03437b4e8f96718588d31dba77..e54602d17d8585654e39ca2c2903c2d903b4073f 100644
--- a/src/lbm/field/MomentumDensity.h
+++ b/src/lbm/field/MomentumDensity.h
@@ -157,7 +157,7 @@ struct MacroscopicForceCorrection< LatticeModel_T, typename std::enable_if< Latt
{
static void apply( const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force();
+ const auto & force = latticeModel.forceModel().forceDensity();
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -185,7 +185,7 @@ struct MacroscopicForceCorrection< LatticeModel_T, typename std::enable_if< ! La
/*
static void apply( const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force();
+ const auto & force = latticeModel.forceModel().forceDensity();
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -197,7 +197,7 @@ struct MacroscopicForceCorrection< LatticeModel_T, typename std::enable_if< ! La
template< typename FieldPtrOrIterator >
static void apply( FieldPtrOrIterator & it, const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force(it.x(), it.y(), it.z());
+ const auto & force = latticeModel.forceModel().forceDensity(it.x(), it.y(), it.z());
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -207,7 +207,7 @@ struct MacroscopicForceCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force(x,y,z);
+ const auto & force = latticeModel.forceModel().forceDensity(x,y,z);
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -246,7 +246,7 @@ struct EquilibriumForceCorrection< LatticeModel_T, typename std::enable_if< Latt
{
static void apply( const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force();
+ const auto & force = latticeModel.forceModel().forceDensity();
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -274,7 +274,7 @@ struct EquilibriumForceCorrection< LatticeModel_T, typename std::enable_if< ! La
/*
static void apply( const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force();
+ const auto & force = latticeModel.forceModel().forceDensity();
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -286,7 +286,7 @@ struct EquilibriumForceCorrection< LatticeModel_T, typename std::enable_if< ! La
template< typename FieldPtrOrIterator >
static void apply( FieldPtrOrIterator & it, const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force(it.x(), it.y(), it.z());
+ const auto & force = latticeModel.forceModel().forceDensity(it.x(), it.y(), it.z());
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
@@ -296,7 +296,7 @@ struct EquilibriumForceCorrection< LatticeModel_T, typename std::enable_if< ! La
static void apply( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const LatticeModel_T & latticeModel, Vector3< real_t > & momentumDensity )
{
- const auto & force = latticeModel.forceModel().force(x,y,z);
+ const auto & force = latticeModel.forceModel().forceDensity(x,y,z);
const real_t dt_2 = real_t(0.5);
momentumDensity[0] += dt_2 * force[0];
diff --git a/src/lbm/lattice_model/ForceModel.h b/src/lbm/lattice_model/ForceModel.h
index c733281ac17544a39e4fac182ee4c67faf280788..527b5f5f395d8164b10cde47e6015e91424b72b4 100644
--- a/src/lbm/lattice_model/ForceModel.h
+++ b/src/lbm/lattice_model/ForceModel.h
@@ -1,15 +1,15 @@
//======================================================================================================================
//
-// This file is part of waLBerla. waLBerla is free software: you can
+// 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
+// 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
+//
+// 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/>.
//
@@ -59,14 +59,14 @@ namespace force_model {
* is a shift by F/2
* 4. static const bool shiftEquVel: specifies whether or not during the evaluation of the equilibrium velocity there
* is a shift by F/2
-* 5. static const bool constant: true, if the body force is constant for all cells throughout the entire simulation
+* 5. static const bool constant: true, if the body force density is constant for all cells throughout the entire simulation
* (false otherwise)
-* 6. void configure( IBlock & block, StructuredBlockStorage & sbs ): this function is called when a force model
+* 6. void configure( IBlock & block, StructuredBlockStorage & sbs ): this function is called when a force density model
* instance is assigned to a block - can be used to
* scale stored values to a different grid level
-* 7. const Vector3<real_t> force() const: must return the body force (this function must only exist if constant == true)
-* 8. const Vector3<real_t> force( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const:
-* -> must return the force for block local cell (x,y,z) - this function must also exist for constant forces!
+* 7. const Vector3<real_t> forceDensity() const: must return the body force density (this function must only exist if constant == true)
+* 8. const Vector3<real_t> forceDensity( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const:
+* -> must return the force density for block local cell (x,y,z) - this function must also exist for constant force densities!
* 9. "template< typename LatticeModel_T >
* DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
* const Vector3<real_t> & velocity, const real_t rho,
@@ -86,11 +86,11 @@ namespace force_model {
* equilibrium distribution, commonTerms = direction-independent terms calculated by function
* "directionIndependentTerms" (see 9), w = lattice model weighting factor, (cx,cy,cz) = lattice
* direction, omega = relaxation parameter that corresponds to the lattice viscosity
-* 11. "bool setConstantBodyForceIfPossible( const Vector3<real_t> & force, const uint_t level = uint_t(0) )":
-* -> Must return true if a constant force can be set. In that case, the force passed in must be set.
-* 'level' is the grid level 'force' corresponds to.
+* 11. "bool setConstantBodyForceIfPossible( const Vector3<real_t> & forceDensity, const uint_t level = uint_t(0) )":
+* -> Must return true if a constant force density can be set. In that case, the force density passed in must be set.
+* 'level' is the grid level 'forceDensity' corresponds to.
*
-* ATTENTION: Both force functions (7+8) and the functions calculating common terms and the force term (9+10) may be
+* ATTENTION: Both forceDensity functions (7+8) and the functions calculating common terms and the force density term (9+10) may be
* called in parallel by multiple threads and therefore must be thread-safe!
*
* For an example see classes SimpleConstant or LuoField.
@@ -101,32 +101,32 @@ namespace force_model {
// a collision and a force model see documentation in LatticeModelBase.h
struct None_tag {};
-struct Simple_tag {}; // Refers to a simple force model which introduces the following additional force term in the
+struct Simple_tag {}; // Refers to a simple force model which introduces the following additional force density term in the
// collision process: 3 * (lattice_weight)_i * (e)_i * (force)_i [often: (force)_i = rho * acceleration]
- // Should only be used with constant forces!
+ // Should only be used with constant force densities!
// Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
struct EDM_tag {}; // Refers to the exact difference method by Kupershtokh
// Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
-struct Luo_tag {}; // Refers to the force model by Luo which introduces a more complex force term in the collision process:
+struct Luo_tag {}; // Refers to the force model by Luo which introduces a more complex force density term in the collision process:
// F_i = w_i * [ (c_i - u) / (c_s^2) + (c_i * (c_i * u)) / (c_s^4) ] * F
// Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
-struct Guo_tag {}; // Refers to the force model by Guo which introduces a more complex force term in the collision process:
+struct Guo_tag {}; // Refers to the force model by Guo which introduces a more complex force density term in the collision process:
// F_i = w_i * [ 1 - 1 / (2 * tau) ] * [ (c_i - u) / (c_s^2) + (c_i * (c_i * u)) / (c_s^4) ] * F
// Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)!
struct Correction_tag{}; // refers to correction terms proposed in Jonas Latt's thesis
-/// Used as return type if there are no common, direction-independent terms for the force term
+/// Used as return type if there are no common, direction-independent terms for the force density term
struct NoDirectionIndependentTerms {};
-/// Returns the body force for level 'targetLevel' given the body force 'force_level' on level 'level'
-inline Vector3< real_t > levelDependentBodyForce( const uint_t targetLevel, const Vector3< real_t > & force_level, const uint_t level )
+/// Returns the body force density for level 'targetLevel' given the body force density 'forceDensity_level' on level 'level'
+inline Vector3< real_t > levelDependentBodyForce( const uint_t targetLevel, const Vector3< real_t > & forceDensity_level, const uint_t level )
{
const real_t powTwoTarget = real_c( uint_t(1) << targetLevel );
const real_t powTwoLevel = real_c( uint_t(1) << level );
- return force_level * ( powTwoLevel / powTwoTarget );
+ return forceDensity_level * ( powTwoLevel / powTwoTarget );
}
/// Returns the acceleration for level 'targetLevel' given the acceleration 'acceleration_level' on level 'level'
@@ -155,10 +155,10 @@ public:
void unpack( mpi::RecvBuffer & ) {}
void configure( IBlock & /*block*/, StructuredBlockStorage & /*sbs*/ ) const {}
-
- const Vector3<real_t> force() const { return Vector3<real_t>(); }
- const Vector3<real_t> force( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return Vector3<real_t>(); }
-
+
+ const Vector3<real_t> forceDensity() const { return Vector3<real_t>(); }
+ const Vector3<real_t> forceDensity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return Vector3<real_t>(); }
+
template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
@@ -174,8 +174,8 @@ public:
-/// For the incompressible LBM, in lattice units, the body force is equal to the acceleration, since
-/// bodyForce_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1
+/// For the incompressible LBM, in lattice units, the body force density is equal to the acceleration, since
+/// bodyForceDensity_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1
class SimpleConstant
{
public:
@@ -188,30 +188,30 @@ public:
static const bool constant = true;
- SimpleConstant( const Vector3< real_t > & bodyForce, const uint_t level = uint_t(0) ) :
- bodyForce_( bodyForce ), level_( level ) {}
+ SimpleConstant( const Vector3< real_t > & bodyForceDensity, const uint_t level = uint_t(0) ) :
+ bodyForceDensity_( bodyForceDensity ), level_( level ) {}
SimpleConstant( const real_t vx, const real_t vy, const real_t vz, const uint_t level = uint_t(0) ) :
- bodyForce_(vx,vy,vz), level_( level ) {}
+ bodyForceDensity_(vx,vy,vz), level_( level ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForce_ << level_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForce_ >> level_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForceDensity_ << level_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForceDensity_ >> level_; }
void configure( IBlock & block, StructuredBlockStorage & sbs )
{
const uint_t level = level_;
level_ = sbs.getLevel( block );
- bodyForce_ = levelDependentBodyForce( level_, bodyForce_, level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity_, level );
}
- /// "force_level" is the level that corresponds to "acceleration"
- void reset( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+ /// "forceDensity_level" is the level that corresponds to "acceleration"
+ void reset( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- bodyForce_ = levelDependentBodyForce( level_, bodyForce, force_level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity, forceDensity_level );
}
- const Vector3< real_t > & force() const { return bodyForce_; }
- const Vector3< real_t > & force( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForce_; }
-
+ const Vector3< real_t > & forceDensity() const { return bodyForceDensity_; }
+ const Vector3< real_t > & forceDensity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForceDensity_; }
+
template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
@@ -222,19 +222,19 @@ public:
const DirectionIndependentTerms_T & /*commonTerms*/, const real_t w,
const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{
- return real_t(3.0) * w * ( cx * bodyForce_[0] + cy * bodyForce_[1] + cz * bodyForce_[2] );
+ return real_t(3.0) * w * ( cx * bodyForceDensity_[0] + cy * bodyForceDensity_[1] + cz * bodyForceDensity_[2] );
}
-
- /// "force_level" is the level that corresponds to "acceleration"
- bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+
+ /// "forceDensity_level" is the level that corresponds to "acceleration"
+ bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- reset( bodyForce, force_level );
+ reset( bodyForceDensity, forceDensity_level );
return true;
}
private:
- Vector3< real_t > bodyForce_;
+ Vector3< real_t > bodyForceDensity_;
uint_t level_;
};
@@ -261,17 +261,17 @@ private:
template< typename LatticeModel_T, class Enable = void >
struct ShiftedVelocity
{
- static Vector3<real_t> get( const Vector3<real_t> & velocity, const Vector3<real_t> & force, const real_t )
+ static Vector3<real_t> get( const Vector3<real_t> & velocity, const Vector3<real_t> & forceDensity, const real_t )
{
- return velocity + force;
+ return velocity + forceDensity;
}
};
template< typename LatticeModel_T >
struct ShiftedVelocity< LatticeModel_T, typename std::enable_if< LatticeModel_T::compressible >::type >
{
- static Vector3<real_t> get( const Vector3<real_t> & velocity, const Vector3<real_t> & force, const real_t rhoInv )
+ static Vector3<real_t> get( const Vector3<real_t> & velocity, const Vector3<real_t> & forceDensity, const real_t rhoInv )
{
- return velocity + force * rhoInv;
+ return velocity + forceDensity * rhoInv;
}
};
@@ -286,20 +286,20 @@ public:
static const bool constant = false;
- EDMField( const BlockDataID & forceFieldId ) :
- forceFieldId_( forceFieldId ), forceField_( NULL ) {}
+ EDMField( const BlockDataID & forceDensityFieldId ) :
+ forceDensityFieldId_( forceDensityFieldId ), forceDensityField_( NULL ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << forceFieldId_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceFieldId_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << forceDensityFieldId_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceDensityFieldId_; }
void configure( IBlock & block, StructuredBlockStorage & /*sbs*/ )
{
- forceField_ = block.getData< ForceField_T >( forceFieldId_ );
+ forceDensityField_ = block.getData< ForceField_T >( forceDensityFieldId_ );
}
- typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector force( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+ typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector forceDensity( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
- return field::VectorFieldAccessor<ForceField_T>::get( forceField_, x,y,z);
+ return field::VectorFieldAccessor<ForceField_T>::get( forceDensityField_, x,y,z);
}
template< typename LatticeModel_T >
@@ -314,7 +314,7 @@ public:
const real_t commonTerm, const real_t w,
const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{
- const auto shiftedVelocity = ShiftedVelocity< LatticeModel_T >::get( velocity, force(x,y,z), commonTerm );
+ const auto shiftedVelocity = ShiftedVelocity< LatticeModel_T >::get( velocity, forceDensity(x,y,z), commonTerm );
return EquilibriumDistribution< LatticeModel_T >::get( cx, cy, cz, w, shiftedVelocity, rho )
- EquilibriumDistribution< LatticeModel_T >::get( cx, cy, cz, w, velocity, rho );
@@ -324,14 +324,14 @@ public:
private:
- BlockDataID forceFieldId_;
- ForceField_T * forceField_;
+ BlockDataID forceDensityFieldId_;
+ ForceField_T * forceDensityField_;
};
/// For the incompressible LBM, in lattice units, the body force is equal to the acceleration, since
-/// bodyForce_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1.
+/// bodyForceDensity_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1.
/// \cite luo1993lattice, \cite shan1994simulation
class LuoConstant
{
@@ -345,30 +345,30 @@ public:
static const bool constant = true;
- LuoConstant( const Vector3< real_t > & bodyForce, const uint_t level = uint_t(0) ) :
- bodyForce_( bodyForce ), level_( level ) {}
+ LuoConstant( const Vector3< real_t > & bodyForceDensity, const uint_t level = uint_t(0) ) :
+ bodyForceDensity_( bodyForceDensity ), level_( level ) {}
LuoConstant( const real_t vx, const real_t vy, const real_t vz, const uint_t level = uint_t(0) ) :
- bodyForce_(vx,vy,vz), level_( level ) {}
+ bodyForceDensity_(vx,vy,vz), level_( level ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForce_ << level_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForce_ >> level_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForceDensity_ << level_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForceDensity_ >> level_; }
void configure( IBlock & block, StructuredBlockStorage & sbs )
{
const uint_t level = level_;
level_ = sbs.getLevel( block );
- bodyForce_ = levelDependentBodyForce( level_, bodyForce_, level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity_, level );
}
- /// "force_level" is the level that corresponds to "bodyForce"
- void reset( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+ /// "forceDensity_level" is the level that corresponds to "bodyForceDensity"
+ void reset( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- bodyForce_ = levelDependentBodyForce( level_, bodyForce, force_level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity, forceDensity_level );
}
- const Vector3< real_t > & force() const { return bodyForce_; }
- const Vector3< real_t > & force( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForce_; }
-
+ const Vector3< real_t > & forceDensity() const { return bodyForceDensity_; }
+ const Vector3< real_t > & forceDensity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForceDensity_; }
+
template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
@@ -380,19 +380,19 @@ public:
const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{
const Vector3<real_t> c( cx, cy, cz );
- return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ );
+ return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForceDensity_ );
}
-
- /// "force_level" is the level that corresponds to "bodyForce"
- bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+
+ /// "forceDensity_level" is the level that corresponds to "bodyForceDensity"
+ bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- reset( bodyForce, force_level );
+ reset( bodyForceDensity, forceDensity_level );
return true;
}
private:
- Vector3< real_t > bodyForce_;
+ Vector3< real_t > bodyForceDensity_;
uint_t level_;
};
@@ -413,22 +413,22 @@ public:
static const bool constant = false;
- LuoField( const BlockDataID & forceFieldId ) :
- forceFieldId_( forceFieldId ), forceField_( NULL ) {}
+ LuoField( const BlockDataID & forceDensityFieldId ) :
+ forceDensityFieldId_( forceDensityFieldId ), forceDensityField_( NULL ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << forceFieldId_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceFieldId_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << forceDensityFieldId_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceDensityFieldId_; }
void configure( IBlock & block, StructuredBlockStorage & /*sbs*/ )
{
- forceField_ = block.getData< ForceField_T >( forceFieldId_ );
+ forceDensityField_ = block.getData< ForceField_T >( forceDensityFieldId_ );
}
- typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector force( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+ typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector forceDensity( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
- return field::VectorFieldAccessor<ForceField_T>::get( forceField_, x,y,z);
+ return field::VectorFieldAccessor<ForceField_T>::get( forceDensityField_, x,y,z);
}
-
+
template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
@@ -440,21 +440,21 @@ public:
const real_t cx, const real_t cy, const real_t cz, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{
const Vector3<real_t> c( cx, cy, cz );
- return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) );
+ return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * forceDensity(x,y,z) );
}
bool setConstantBodyForceIfPossible( const Vector3< real_t > &, const uint_t = uint_t(0) ) { return false; }
private:
- BlockDataID forceFieldId_;
- ForceField_T * forceField_;
+ BlockDataID forceDensityFieldId_;
+ ForceField_T * forceDensityField_;
};
/// For the incompressible LBM, in lattice units, the body force is equal to the acceleration, since
-/// bodyForce_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1.
+/// bodyForceDensity_lattice = density_lattice * acceleration_lattice and density_lattice is implicitly set to 1.
/// \cite guo2002discrete, \cite schiller2008thermal
class GuoConstant
{
@@ -468,30 +468,30 @@ public:
static const bool constant = true;
- GuoConstant( const Vector3< real_t > & bodyForce, const uint_t level = uint_t(0) ) :
- bodyForce_( bodyForce ), level_( level ) {}
+ GuoConstant( const Vector3< real_t > & bodyForceDensity, const uint_t level = uint_t(0) ) :
+ bodyForceDensity_( bodyForceDensity ), level_( level ) {}
GuoConstant( const real_t vx, const real_t vy, const real_t vz, const uint_t level = uint_t(0) ) :
- bodyForce_(vx,vy,vz), level_( level ) {}
+ bodyForceDensity_(vx,vy,vz), level_( level ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForce_ << level_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForce_ >> level_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << bodyForceDensity_ << level_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> bodyForceDensity_ >> level_; }
void configure( IBlock & block, StructuredBlockStorage & sbs )
{
const uint_t level = level_;
level_ = sbs.getLevel( block );
- bodyForce_ = levelDependentBodyForce( level_, bodyForce_, level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity_, level );
}
- /// "force_level" is the level that corresponds to "bodyForce"
- void reset( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+ /// "forceDensity_level" is the level that corresponds to "bodyForceDensity"
+ void reset( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- bodyForce_ = levelDependentBodyForce( level_, bodyForce, force_level );
+ bodyForceDensity_ = levelDependentBodyForce( level_, bodyForceDensity, forceDensity_level );
}
- const Vector3< real_t > & force() const { return bodyForce_; }
- const Vector3< real_t > & force( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForce_; }
-
+ const Vector3< real_t > & forceDensity() const { return bodyForceDensity_; }
+ const Vector3< real_t > & forceDensity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return bodyForceDensity_; }
+
template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk, const real_t /*omega_odd*/ ) const
@@ -499,10 +499,10 @@ public:
if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
{
const real_t one_over_d = real_t(1) / real_t(LatticeModel_T::Stencil::D);
-
- const auto common = Matrix3<real_t>::makeDiagonalMatrix( velocity * bodyForce_ );
- return (tensorProduct( velocity, bodyForce_ ) +
- tensorProduct( bodyForce_, velocity ) -
+
+ const auto common = Matrix3<real_t>::makeDiagonalMatrix( velocity * bodyForceDensity_ );
+ return (tensorProduct( velocity, bodyForceDensity_ ) +
+ tensorProduct( bodyForceDensity_, velocity ) -
common * (real_t(2)*one_over_d) ) * real_t(0.5) * ( real_t(2) - omega )
+ common * ( one_over_d * ( real_t(2) - omega_bulk ) );
}
@@ -522,33 +522,33 @@ public:
if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
{
const real_t one_third = real_t(1) / real_t(3);
-
+
const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace();
- return real_t(3.0) * w * ( bodyForce_ * c + real_t(1.5) * common);
+ return real_t(3.0) * w * ( bodyForceDensity_ * c + real_t(1.5) * common);
}
else if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::TRT_tag >::value)
{
return real_t(3.0) * w * ( ( real_t(1) - real_t(0.5) * omega ) *
- ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ ) +
- ( omega - omega_odd ) * real_t(0.5) * (c * bodyForce_) );
+ ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForceDensity_ ) +
+ ( omega - omega_odd ) * real_t(0.5) * (c * bodyForceDensity_) );
}
else
{
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
- ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ );
+ ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForceDensity_ );
}
}
-
- /// "force_level" is the level that corresponds to "bodyForce"
- bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForce, const uint_t force_level = uint_t(0) )
+
+ /// "forceDensity_level" is the level that corresponds to "bodyForceDensity"
+ bool setConstantBodyForceIfPossible( const Vector3< real_t > & bodyForceDensity, const uint_t forceDensity_level = uint_t(0) )
{
- reset( bodyForce, force_level );
+ reset( bodyForceDensity, forceDensity_level );
return true;
}
private:
- Vector3< real_t > bodyForce_;
+ Vector3< real_t > bodyForceDensity_;
uint_t level_;
};
@@ -569,20 +569,20 @@ public:
static const bool constant = false;
- GuoField( const BlockDataID & forceFieldId ) :
- forceFieldId_( forceFieldId ), forceField_( nullptr ) {}
+ GuoField( const BlockDataID & forceDensityFieldId ) :
+ forceDensityFieldId_( forceDensityFieldId ), forceDensityField_( nullptr ) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << forceFieldId_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceFieldId_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << forceDensityFieldId_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceDensityFieldId_; }
void configure( IBlock & block, StructuredBlockStorage & /*sbs*/ )
{
- forceField_ = block.getData< ForceField_T >( forceFieldId_ );
+ forceDensityField_ = block.getData< ForceField_T >( forceDensityFieldId_ );
}
- typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector force( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+ typename field::VectorFieldAccessor<ForceField_T>::vector_or_constRefVector forceDensity( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
- return field::VectorFieldAccessor<ForceField_T>::get( forceField_, x,y,z);
+ return field::VectorFieldAccessor<ForceField_T>::get( forceDensityField_, x,y,z);
}
template< typename LatticeModel_T >
@@ -592,10 +592,10 @@ public:
if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
{
const real_t one_over_d = real_t(1) / real_t(LatticeModel_T::Stencil::D);
-
- const auto common = Matrix3<real_t>::makeDiagonalMatrix( velocity * force(x,y,z) );
- return (tensorProduct( velocity, force(x,y,z) ) +
- tensorProduct( force(x,y,z), velocity ) -
+
+ const auto common = Matrix3<real_t>::makeDiagonalMatrix( velocity * forceDensity(x,y,z) );
+ return (tensorProduct( velocity, forceDensity(x,y,z) ) +
+ tensorProduct( forceDensity(x,y,z), velocity ) -
common * (real_t(2)*one_over_d) ) * real_t(0.5) * ( real_t(2) - omega )
+ common * ( one_over_d * ( real_t(2) - omega_bulk ) );
}
@@ -615,20 +615,20 @@ public:
if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::MRT_tag >::value)
{
const real_t one_third = real_t(1) / real_t(3);
-
+
const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace();
- return real_t(3.0) * w * ( force(x,y,z) * c + real_t(1.5) * common);
+ return real_t(3.0) * w * ( forceDensity(x,y,z) * c + real_t(1.5) * common);
}
else if (std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::TRT_tag >::value)
{
return real_t(3.0) * w * ( ( real_t(1) - real_t(0.5) * omega ) *
- ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) ) +
- ( omega - omega_odd ) * real_t(0.5) * (c * force(x,y,z)) );
+ ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * forceDensity(x,y,z) ) +
+ ( omega - omega_odd ) * real_t(0.5) * (c * forceDensity(x,y,z)) );
}
else
{
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
- ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) );
+ ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * forceDensity(x,y,z) );
}
}
@@ -636,8 +636,8 @@ public:
private:
- BlockDataID forceFieldId_;
- ForceField_T * forceField_;
+ BlockDataID forceDensityFieldId_;
+ ForceField_T * forceDensityField_;
};
@@ -657,19 +657,19 @@ public:
static const bool constant = true;
Correction( const BlockDataID & previousRhoVelocityId ) :
- force_( real_t(0) ), previousRhoVelocityId_( previousRhoVelocityId ), previousRhoVelocity_(nullptr) {}
+ forceDensity_( real_t(0) ), previousRhoVelocityId_( previousRhoVelocityId ), previousRhoVelocity_(nullptr) {}
- void pack( mpi::SendBuffer & buffer ) const { buffer << force_ << previousRhoVelocityId_; }
- void unpack( mpi::RecvBuffer & buffer ) { buffer >> force_ >> previousRhoVelocityId_; }
+ void pack( mpi::SendBuffer & buffer ) const { buffer << forceDensity_ << previousRhoVelocityId_; }
+ void unpack( mpi::RecvBuffer & buffer ) { buffer >> forceDensity_ >> previousRhoVelocityId_; }
void configure( IBlock & block, StructuredBlockStorage & /*sbs*/ )
{
previousRhoVelocity_ = block.getData< MomentumDensityField_T >( previousRhoVelocityId_ );
}
- const Vector3< real_t > & force() const { return force_; }
- const Vector3< real_t > & force( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return force_; }
-
+ const Vector3< real_t > & forceDensity() const { return forceDensity_; }
+ const Vector3< real_t > & forceDensity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return forceDensity_; }
+
template< typename LatticeModel_T >
Vector3<real_t> directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const Vector3<real_t> & velocity, const real_t rho, const real_t omega, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
@@ -687,12 +687,12 @@ public:
{
return w * ( cx * commonTerm[0] + cy * commonTerm[1] + cz * commonTerm[2] );
}
-
+
bool setConstantBodyForceIfPossible( const Vector3<real_t> &, const uint_t = uint_t(0) ) { return false; }
private:
- Vector3< real_t > force_;
+ Vector3< real_t > forceDensity_;
BlockDataID previousRhoVelocityId_;
MomentumDensityField_T * previousRhoVelocity_;
diff --git a/src/lbm/lattice_model/LatticeModelBase.h b/src/lbm/lattice_model/LatticeModelBase.h
index 3ded2ae7e4b3aecd09547c25265a4c1cfeed9678..568800a4d0434a5e87c6b00078fe8234209c077b 100644
--- a/src/lbm/lattice_model/LatticeModelBase.h
+++ b/src/lbm/lattice_model/LatticeModelBase.h
@@ -24,6 +24,7 @@
#include "CollisionModel.h"
#include "ForceModel.h"
#include "core/DataTypes.h"
+#include "core/logging/Logging.h"
#include "core/mpi/BufferSizeTrait.h"
#include "core/mpi/RecvBuffer.h"
#include "core/mpi/SendBuffer.h"
@@ -80,7 +81,25 @@ public:
LatticeModelBase( const CollisionModel_T & cm, const ForceModel_T & fm ) :
- collisionModel_( cm ), forceModel_( fm ) {}
+ collisionModel_( cm ), forceModel_( fm ) {
+
+ if (Compressible && forceModel_.constant)
+ {
+ WALBERLA_LOG_WARNING_ON_ROOT("WARNING: You are using a compressible lattice model with a constant force "
+ "model. You should consider using a field-based force model, and adjust the body "
+ "force density according to the change in the fluid density in every time step.");
+ }
+
+ if ((std::is_same< typename CollisionModel_T::tag, collision_model::TRT_tag >::value ||
+ std::is_same< typename CollisionModel_T::tag, collision_model::MRT_tag >::value) &&
+ std::is_same< typename ForceModel_T::tag, force_model::Guo_tag >::value)
+ {
+ WALBERLA_LOG_WARNING_ON_ROOT(
+ "WARNING: waLBerla currently does not support using a TRT or MRT collision model in combination with "
+ "the force model from Guo. See https://i10git.cs.fau.de/walberla/walberla/-/issues/176 and "
+ "https://i10git.cs.fau.de/walberla/walberla/-/merge_requests/560 for more information.");
+ }
+ }
virtual ~LatticeModelBase() = default;
diff --git a/src/lbm/srt/CellwiseSweep.impl.h b/src/lbm/srt/CellwiseSweep.impl.h
index ec6d7fa4ba77200c5f8a99e590cd3ef4f4e9de02..62a7bfe925cada300e023c3c870166b5e21d607b 100644
--- a/src/lbm/srt/CellwiseSweep.impl.h
+++ b/src/lbm/srt/CellwiseSweep.impl.h
@@ -599,7 +599,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
dst->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1 * ( vel_trm_E_W + velX ) + three_w1 * force[0];
dst->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1 * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -686,7 +686,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
src->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1 * ( vel_trm_E_W + velX ) + three_w1 * force[0];
src->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1 * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -794,7 +794,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
dst->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1_rho * ( vel_trm_E_W + velX ) + three_w1 * force[0];
dst->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1_rho * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -885,7 +885,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
src->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1_rho * ( vel_trm_E_W + velX ) + three_w1 * force[0];
src->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1_rho * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -1478,7 +1478,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
dst->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1 * ( vel_trm_E_W + velX ) + three_w1 * force[0];
dst->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1 * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -1591,7 +1591,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
src->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1 * ( vel_trm_E_W + velX ) + three_w1 * force[0];
src->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1 * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -1725,7 +1725,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
dst->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1_rho * ( vel_trm_E_W + velX ) + three_w1 * force[0];
dst->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1_rho * ( vel_trm_E_W - velX ) - three_w1 * force[0];
@@ -1844,7 +1844,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
const real_t vel_trm_N_S = dir_indep_trm + real_t(1.5) * velYY;
const real_t vel_trm_T_B = dir_indep_trm + real_t(1.5) * velZZ;
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
src->get(x,y,z,Stencil_T::idx[E]) = omega_trm * vE + omega_w1_rho * ( vel_trm_E_W + velX ) + three_w1 * force[0];
src->get(x,y,z,Stencil_T::idx[W]) = omega_trm * vW + omega_w1_rho * ( vel_trm_E_W - velX ) - three_w1 * force[0];
diff --git a/src/lbm/srt/cell_operations/DefaultCellOperation.impl.h b/src/lbm/srt/cell_operations/DefaultCellOperation.impl.h
index 0cf052e6b5ea7f04dbb5cbec8af475002d3b67d5..0a56869bf06cd4d7cddb09a98b79c6c65e1f82c4 100644
--- a/src/lbm/srt/cell_operations/DefaultCellOperation.impl.h
+++ b/src/lbm/srt/cell_operations/DefaultCellOperation.impl.h
@@ -145,7 +145,7 @@ public:
{
omega_ = latticeModel.collisionModel().omega();
latticeModel_ = &latticeModel;
- force_ = latticeModel.forceModel().force();
+ force_ = latticeModel.forceModel().forceDensity();
}
void operator()( PdfField_T * src, PdfField_T * dst, cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
diff --git a/src/lbm/trt/CellwiseSweep.impl.h b/src/lbm/trt/CellwiseSweep.impl.h
index 4780803d80873db706e808b91452586f303ea41c..c38526e9d8bb7019a42376c61d094ce3df1daf65 100644
--- a/src/lbm/trt/CellwiseSweep.impl.h
+++ b/src/lbm/trt/CellwiseSweep.impl.h
@@ -640,7 +640,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
dst->get( x, y, z, Stencil_T::idx[C] ) = vC * (real_t(1.0) - lambda_e) + lambda_e * t0 * feq_common; // no force term
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
const real_t velXPY = velX + velY;
const real_t sym_NE_SW = lambda_e_scaled * ( vNE + vSW - fac2 * velXPY * velXPY - t2x2 * feq_common );
@@ -735,7 +735,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
src->get( x, y, z, Stencil_T::idx[C] ) = vC * (real_t(1.0) - lambda_e) + lambda_e * t0 * feq_common; // no force term
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
const real_t velXPY = velX + velY;
const real_t sym_NE_SW = lambda_e_scaled * ( vNE + vSW - fac2 * velXPY * velXPY - t2x2 * feq_common );
@@ -1369,7 +1369,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
dst->get( x, y, z, Stencil_T::idx[C] ) = vC * (real_t(1.0) - lambda_e) + lambda_e * t0 * feq_common;// no force term
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
const real_t velXPY = velX + velY;
const real_t sym_NE_SW = lambda_e_scaled * ( vNE + vSW - fac2 * velXPY * velXPY - t2x2 * feq_common );
@@ -1492,7 +1492,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
src->get( x, y, z, Stencil_T::idx[C] ) = vC * (real_t(1.0) - lambda_e) + lambda_e * t0 * feq_common; // no force term
- const Vector3< real_t > & force = src->latticeModel().forceModel().force(x,y,z);
+ const Vector3< real_t > & force = src->latticeModel().forceModel().forceDensity(x,y,z);
const real_t velXPY = velX + velY;
const real_t sym_NE_SW = lambda_e_scaled * ( vNE + vSW - fac2 * velXPY * velXPY - t2x2 * feq_common );
@@ -1626,7 +1626,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
real_t rho = this->densityVelocityIn( velocity, dst, x, y, z );
if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::Guo_tag >::value)
- velocity -= real_c(0.5) * lm.forceModel().force(x,y,z);
+ velocity -= real_c(0.5) * lm.forceModel().forceDensity(x,y,z);
this->densityVelocityOut( x, y, z, lm, velocity, rho );
@@ -1677,7 +1677,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
real_t rho = this->densityVelocityIn( velocity, src, x, y, z );
if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::Guo_tag >::value)
- velocity -= real_c(0.5) * lm.forceModel().force(x,y,z);
+ velocity -= real_c(0.5) * lm.forceModel().forceDensity(x,y,z);
this->densityVelocityOut( x, y, z, lm, velocity, rho );
diff --git a/src/pe_coupling/discrete_particle_methods/gns_lbm/GNSSweep.h b/src/pe_coupling/discrete_particle_methods/gns_lbm/GNSSweep.h
index a36e475ef06c216441d20894c6db17fd0670b22e..58f35a75acf50c9474ec9e76d50da66626c71994 100644
--- a/src/pe_coupling/discrete_particle_methods/gns_lbm/GNSSweep.h
+++ b/src/pe_coupling/discrete_particle_methods/gns_lbm/GNSSweep.h
@@ -156,7 +156,7 @@ void GNSSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
- const Vector3<real_t> extForce = lm.forceModel().force(x,y,z);
+ const Vector3<real_t> extForce = lm.forceModel().forceDensity(x,y,z);
// collide
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
@@ -224,7 +224,7 @@ void GNSSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
- const Vector3<real_t> extForce = lm.forceModel().force(x,y,z);
+ const Vector3<real_t> extForce = lm.forceModel().forceDensity(x,y,z);
// collide
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )