Skip to content

GitLab

Commit 6e75c8a8 authored by Markus Holzer's avatar Markus Holzer
Browse files

Merge branch 'guo' into 'master' 

Fix TRT Guo

Closes #155

See merge request !476
parents 7be319e9 7a293a3f
Hide whitespace changes
Inline Side-by-side
Showing with 54 additions and 38 deletions
src/lbm/lattice_model/CollisionModel.h
View file @ 6e75c8a8
......@@ -109,6 +109,7 @@ public:
real_t omega() 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 omega( const cell_idx_t /*x*/, const cell_idx_t /*y*/, const cell_idx_t /*z*/,
......@@ -274,6 +275,7 @@ public:
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_odd() const { return lambda_d(); }
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 )
......@@ -484,6 +486,7 @@ public:
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_odd() const { return s_[4]; }
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_; }
......
src/lbm/lattice_model/ForceModel.h
View file @ 6e75c8a8
......@@ -161,13 +161,13 @@ public:
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
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(); }
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*/,
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; }
};
......@@ -214,13 +214,13 @@ public:
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
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(); }
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*/,
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] );
}
......@@ -304,7 +304,7 @@ public:
template< typename LatticeModel_T >
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 );
}
......@@ -312,7 +312,7 @@ public:
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,
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 );
......@@ -371,13 +371,13 @@ public:
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
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(); }
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*/,
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 );
return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * bodyForce_ );
......@@ -431,13 +431,13 @@ public:
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
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(); }
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*/,
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 );
return real_t(3) * w * ( ( c - velocity + ( real_t(3) * ( c * velocity ) * c ) ) * force(x,y,z) );
......@@ -494,9 +494,9 @@ public:
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
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);
......@@ -508,6 +508,7 @@ public:
}
else
{
WALBERLA_ASSERT_FLOAT_EQUAL( omega, omega_bulk );
return DirectionIndependentTerms_T();
}
}
......@@ -515,16 +516,22 @@ public:
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*/,
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 );
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 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);
}
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
{
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
......@@ -580,9 +587,9 @@ public:
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
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);
......@@ -602,16 +609,22 @@ public:
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*/,
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 );
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 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);
}
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
{
return real_t(3.0) * w * ( real_t(1) - real_t(0.5) * omega ) *
......@@ -659,7 +672,7 @@ public:
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
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;
Vector3<real_t> & previousRhoVelocity = previousRhoVelocity_->get(x,y,z);
......@@ -670,7 +683,7 @@ public:
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*/,
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] );
}
......
src/lbm/mrt/CellwiseSweep.impl.h
View file @ 6e75c8a8
......@@ -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)
{
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 )
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
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 )
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() );
}
}
......
src/lbm/srt/CellwiseSweep.impl.h
View file @ 6e75c8a8
......@@ -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 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
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() ],
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() ) +
omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
......@@ -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 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 )
{
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() ) +
omega * EquilibriumDistribution< LatticeModel_T >::get( *d, velocity, rho ) +
......
src/lbm/trt/CellwiseSweep.impl.h
View file @ 6e75c8a8
......@@ -1630,13 +1630,13 @@ WALBERLA_LBM_CELLWISE_SWEEP_STREAM_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPE
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
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() ],
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 fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
......@@ -1681,7 +1681,7 @@ WALBERLA_LBM_CELLWISE_SWEEP_COLLIDE_HEAD( WALBERLA_LBM_CELLWISE_SWEEP_SPECIALIZA
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 )
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
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() ],
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 fasym = EquilibriumDistribution< LatticeModel_T >::getAsymmetricPart( *d, velocity, rho );
......
src/pe_coupling/partially_saturated_cells_method/PSMSweep.h
View file @ 6e75c8a8
......@@ -187,7 +187,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
// 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
if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
......@@ -252,7 +252,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{
// external forcing
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
+ omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
......@@ -265,7 +265,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{
// external forcing
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;
}
......@@ -340,7 +340,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
auto pdfs_equ = lbm::EquilibriumDistribution< LatticeModel_T >::get( velocity, rho );
// 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
if( bodyAndVolumeFractionField->get(x,y,z).size() != size_t(0) )
......@@ -406,7 +406,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{
// external forcing
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
+ omega_n[d.toIdx()] + ( real_c(1) - Bn ) * forceTerm;
......@@ -419,7 +419,7 @@ void PSMSweep< LatticeModel_T, Filter_T, DensityVelocityIn_T, DensityVelocityOut
{
// external forcing
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;
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment