Commit 7a293a3f authored by Michael Kuron's avatar Michael Kuron
Browse files

Fix TRT Guo

parent 3455bf3e
...@@ -109,6 +109,7 @@ public: ...@@ -109,6 +109,7 @@ public:
real_t omega() const { return omega_; } real_t omega() const { return omega_; }
inline real_t omega_bulk() const { return omega(); } inline real_t omega_bulk() const { return omega(); }
inline real_t omega_odd() const { return omega(); }
real_t viscosity() const { return viscosity_; } real_t viscosity() const { return viscosity_; }
real_t omega( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, real_t omega( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
...@@ -274,6 +275,7 @@ public: ...@@ -274,6 +275,7 @@ public:
const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); } const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); }
inline real_t omega_bulk() const { return omega(); } inline real_t omega_bulk() const { return omega(); }
inline real_t omega_odd() const { return lambda_d(); }
static real_t lambda_e( const real_t _omega ) { return _omega; } static real_t lambda_e( const real_t _omega ) { return _omega; }
static real_t lambda_d( const real_t _omega, const real_t _magicNumber = threeSixteenth ) static real_t lambda_d( const real_t _omega, const real_t _magicNumber = threeSixteenth )
...@@ -484,6 +486,7 @@ public: ...@@ -484,6 +486,7 @@ public:
const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); } const Vector3<real_t> & /*velocity*/ = Vector3<real_t>(), const real_t /*rho*/ = real_t(1) ) const { return omega(); }
real_t omega_bulk() const { return s_[1]; } real_t omega_bulk() const { return s_[1]; }
real_t omega_odd() const { return s_[4]; }
real_t viscosity() const { return viscosity_; } real_t viscosity() const { return viscosity_; }
real_t viscosity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return viscosity_; } real_t viscosity( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/ ) const { return viscosity_; }
......
...@@ -161,13 +161,13 @@ public: ...@@ -161,13 +161,13 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ return DirectionIndependentTerms_T(); } { return DirectionIndependentTerms_T(); }
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & /*commonTerms*/, const real_t /*w*/, 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 { return real_t(0); } 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(0); }
bool setConstantBodyForceIfPossible( const Vector3<real_t> &, const uint_t = uint_t(0) ) { return false; } bool setConstantBodyForceIfPossible( const Vector3<real_t> &, const uint_t = uint_t(0) ) { return false; }
}; };
...@@ -214,13 +214,13 @@ public: ...@@ -214,13 +214,13 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ return DirectionIndependentTerms_T(); } { return DirectionIndependentTerms_T(); }
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & /*commonTerms*/, const real_t w, 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 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 * bodyForce_[0] + cy * bodyForce_[1] + cz * bodyForce_[2] );
} }
...@@ -304,7 +304,7 @@ public: ...@@ -304,7 +304,7 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & /*velocity*/, const real_t rho, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ {
return DirectionIndependentTerm< LatticeModel_T >::get( rho ); return DirectionIndependentTerm< LatticeModel_T >::get( rho );
} }
...@@ -312,7 +312,7 @@ public: ...@@ -312,7 +312,7 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t rho, real_t forceTerm( 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 commonTerm, const real_t w, 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 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, force(x,y,z), commonTerm );
...@@ -371,13 +371,13 @@ public: ...@@ -371,13 +371,13 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ return DirectionIndependentTerms_T(); } { return DirectionIndependentTerms_T(); }
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & velocity, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & /*commonTerms*/, const real_t w, 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 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 ); 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 ) ) * bodyForce_ );
...@@ -431,13 +431,13 @@ public: ...@@ -431,13 +431,13 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, const real_t /*omega*/, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ return DirectionIndependentTerms_T(); } { return DirectionIndependentTerms_T(); }
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & /*commonTerms*/, const real_t w, 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 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 ); 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 ) ) * force(x,y,z) );
...@@ -494,9 +494,9 @@ public: ...@@ -494,9 +494,9 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, 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 const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk, const real_t /*omega_odd*/ ) const
{ {
if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value) 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 real_t one_over_d = real_t(1) / real_t(LatticeModel_T::Stencil::D);
...@@ -508,6 +508,7 @@ public: ...@@ -508,6 +508,7 @@ public:
} }
else else
{ {
WALBERLA_ASSERT_FLOAT_EQUAL( omega, omega_bulk );
return DirectionIndependentTerms_T(); return DirectionIndependentTerms_T();
} }
} }
...@@ -515,16 +516,22 @@ public: ...@@ -515,16 +516,22 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & velocity, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & commonTerms, const real_t w, 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 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 ); const Vector3<real_t> c( cx, cy, cz );
if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value) 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 one_third = real_t(1) / real_t(3);
const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace(); 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 * ( bodyForce_ * 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_) );
}
else else
{ {
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) * return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
...@@ -580,9 +587,9 @@ public: ...@@ -580,9 +587,9 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
DirectionIndependentTerms_T directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, 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 const Vector3<real_t> & velocity, const real_t /*rho*/, const real_t omega, const real_t omega_bulk, const real_t /*omega_odd*/ ) const
{ {
if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value) 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 real_t one_over_d = real_t(1) / real_t(LatticeModel_T::Stencil::D);
...@@ -602,16 +609,22 @@ public: ...@@ -602,16 +609,22 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const Vector3<real_t> & velocity, const real_t /*rho*/, real_t forceTerm( 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 DirectionIndependentTerms_T & commonTerms, const real_t w, 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 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 ); const Vector3<real_t> c( cx, cy, cz );
if (!std::is_same< typename LatticeModel_T::CollisionModel::tag, collision_model::SRT_tag >::value) 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 one_third = real_t(1) / real_t(3);
const real_t common = (commonTerms * ( tensorProduct(c,c) - Matrix3<real_t>::makeDiagonalMatrix(one_third) )).trace(); 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 * ( force(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)) );
}
else else
{ {
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) * return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
...@@ -659,7 +672,7 @@ public: ...@@ -659,7 +672,7 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
Vector3<real_t> directionIndependentTerms( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, 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 const Vector3<real_t> & velocity, const real_t rho, const real_t omega, const real_t /*omega_bulk*/, const real_t /*omega_odd*/ ) const
{ {
const auto rhoVelocity = rho * velocity; const auto rhoVelocity = rho * velocity;
Vector3<real_t> & previousRhoVelocity = previousRhoVelocity_->get(x,y,z); Vector3<real_t> & previousRhoVelocity = previousRhoVelocity_->get(x,y,z);
...@@ -670,7 +683,7 @@ public: ...@@ -670,7 +683,7 @@ public:
template< typename LatticeModel_T > template< typename LatticeModel_T >
real_t forceTerm( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/, const Vector3<real_t> & /*velocity*/, const real_t /*rho*/, real_t forceTerm( 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 Vector3<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 const Vector3<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
{ {
return w * ( cx * commonTerm[0] + cy * commonTerm[1] + cz * commonTerm[2] ); return w * ( cx * commonTerm[0] + cy * commonTerm[1] + cz * commonTerm[2] );
} }
......
...@@ -185,9 +185,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE ...@@ -185,9 +185,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false) if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false)
{ {
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk() ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
dst->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk() ); dst->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
} }
} }
...@@ -324,9 +324,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA ...@@ -324,9 +324,9 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false) if (std::is_same< typename LatticeModel_T::ForceModel::tag, force_model::None_tag >::value == false)
{ {
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk() ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
src->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk() ); src->get( x, y, z, d.toIdx() ) += lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho + real_t(1.0), commonForceTerms, LatticeModel_T::w[ d.toIdx() ], real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), collisionModel.omega(), collisionModel.omega_bulk(), collisionModel.omega_odd() );
} }
} }
......
...@@ -1965,13 +1965,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE ...@@ -1965,13 +1965,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho ); const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega, omega );
// collide // collide
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
{ {
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega, omega );
dst->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * dst->get( x, y, z, d.toIdx() ) + dst->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * dst->get( x, y, z, d.toIdx() ) +
omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) + omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
...@@ -1996,12 +1996,12 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA ...@@ -1996,12 +1996,12 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho ); const real_t omega = lm.collisionModel().omega( x, y, z, velocity, rho );
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, omega, omega );
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
{ {
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, omega, omega );
src->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * src->get( x, y, z, d.toIdx() ) + src->get( x, y, z, d.toIdx() ) = ( real_t(1.0) - omega ) * src->get( x, y, z, d.toIdx() ) +
omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) + omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
......
...@@ -1630,13 +1630,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE ...@@ -1630,13 +1630,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
this->densityVelocityOut( x, y, z, lm, velocity, rho ); this->densityVelocityOut( x, y, z, lm, velocity, rho );
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e, lambda_d );
// collide // collide
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
{ {
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e, lambda_d );
const real_t fsym = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho ); const real_t fsym = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho );
const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho ); const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
...@@ -1681,7 +1681,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA ...@@ -1681,7 +1681,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
this->densityVelocityOut( x, y, z, lm, velocity, rho ); this->densityVelocityOut( x, y, z, lm, velocity, rho );
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, lambda_e, lambda_e, lambda_d );
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
pdfs[ d.toIdx() ] = src->get( x, y, z, d.toIdx() ); pdfs[ d.toIdx() ] = src->get( x, y, z, d.toIdx() );
...@@ -1689,7 +1689,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA ...@@ -1689,7 +1689,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d ) for( auto d = Stencil_T::begin(); d != Stencil_T::end(); ++d )
{ {
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), lambda_e, lambda_e, lambda_d );
const real_t fsym = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho ); const real_t fsym = EquilibriumDistribution< LatticeModel_T >::getSymmetricPart ( *d, velocity, rho );
const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho ); const real_t fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
......
...@@ -187,7 +187,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -187,7 +187,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho ); auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
// possible external forcing on fluid // possible external forcing on fluid
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk() ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
// check if body is present // check if body is present
if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) ) if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
...@@ -252,7 +252,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -252,7 +252,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{ {
// external forcing // external forcing
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT
+ omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm; + omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
...@@ -265,7 +265,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -265,7 +265,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{ {
// external forcing // external forcing
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) + forceTerm; dst->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) + forceTerm;
} }
...@@ -340,7 +340,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -340,7 +340,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho ); auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
// possible external forcing on fluid // possible external forcing on fluid
const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk() ); const auto commonForceTerms = lm.forceModel().template directionIndependentTerms< LatticeModel_T >( x, y, z, velocity, rho, omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
// check if a body is present // check if a body is present
if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) ) if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
...@@ -406,7 +406,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -406,7 +406,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{ {
// external forcing // external forcing
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( real_c(1) - Bn ) * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) //SRT
+ omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm; + omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
...@@ -419,7 +419,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut ...@@ -419,7 +419,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{ {
// external forcing // external forcing
const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ], const real_t forceTerm = lm.forceModel().template forceTerm< LatticeModel_T >( x, y, z, velocity, rho, commonForceTerms, LatticeModel_T::w[ d.toIdx() ],
real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk() ); real_c(d.cx()), real_c(d.cy()), real_c(d.cz()), omega, collisionModel.omega_bulk(), collisionModel.omega_odd() );
src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) + forceTerm; src->get( x, y, z, d.toIdx() ) = pdfs[d.toIdx()] - omega * ( pdfs[d.toIdx()] - pdfs_equ[d.toIdx()] ) + forceTerm;
} }
......