Skip to content

Commit

Permalink
Merge branch 'refactory/linear_correction_viscous_force' into fix/rea…
Browse files Browse the repository at this point in the history 
…ction_dynamics_naming
  • Loading branch information
@Xiangyu-Hu
Xiangyu-Hu committed Jul 23, 2024
2 parents c344c09 + 1317a70 commit f6f70dd
Show file tree
Hide file tree
Showing 12 changed files with 325 additions and 295 deletions.
5 changes: 2 additions & 3 deletions src/shared/particle_dynamics/fluid_dynamics/transport_velocity_correction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,18 +120,17 @@ using TransportVelocityCorrectionComplex =
BaseTransportVelocityCorrectionComplex<SingleResolution, NoLimiter, NoKernelCorrection, ParticleScope>;

template <class ParticleScope>
using TransportVelocityCorrectionCorrectedComplex =
using TransportVelocityCorrectionCorrectedComplex =
BaseTransportVelocityCorrectionComplex<SingleResolution, NoLimiter, LinearGradientCorrection, ParticleScope>;

template <class ParticleScope>
using TransportVelocityLimitedCorrectionComplex =
BaseTransportVelocityCorrectionComplex<SingleResolution, TruncatedLinear, NoKernelCorrection, ParticleScope>;

template <class ParticleScope>
using TransportVelocityKimitedCorrectionCorrectedComplex =
using TransportVelocityLimitedCorrectionCorrectedComplex =
BaseTransportVelocityCorrectionComplex<SingleResolution, TruncatedLinear, LinearGradientCorrection, ParticleScope>;


template <class ParticleScope>
using TransportVelocityCorrectionComplexAdaptive =
BaseTransportVelocityCorrectionComplex<AdaptiveResolution, NoLimiter, NoKernelCorrection, ParticleScope>;
Expand Down
48 changes: 26 additions & 22 deletions src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ class ViscousForce<DataDelegationType>
Real smoothing_length_;
};

template <typename ViscosityType>
class ViscousForce<Inner<>, ViscosityType>
template <typename ViscosityType, class KernelCorrectionType>
class ViscousForce<Inner<>, ViscosityType, KernelCorrectionType>
: public ViscousForce<DataDelegateInner>, public ForcePrior
{
public:
Expand All @@ -91,43 +91,43 @@ class ViscousForce<Inner<>, ViscosityType>

protected:
ViscosityType mu_;
KernelCorrectionType kernel_correction_;
};
using ViscousForceInner = ViscousForce<Inner<>, FixedViscosity>;
using ViscousForceInner = ViscousForce<Inner<>, FixedViscosity, NoKernelCorrection>;

template <typename ViscosityType>
class ViscousForce<Inner<AngularConservative>, ViscosityType>
template <typename ViscosityType, class KernelCorrectionType>
class ViscousForce<Inner<AngularConservative>, ViscosityType, KernelCorrectionType>
: public ViscousForce<DataDelegateInner>, public ForcePrior
{
public:
explicit ViscousForce(BaseInnerRelation &inner_relation)
: ViscousForce<DataDelegateInner>(inner_relation),
ForcePrior(particles_, "ViscousForce"),
mu_(particles_){};
explicit ViscousForce(BaseInnerRelation &inner_relation);
virtual ~ViscousForce(){};
void interaction(size_t index_i, Real dt = 0.0);

protected:
ViscosityType mu_;
KernelCorrectionType kernel_correction_;
};

using BaseViscousForceWithWall = InteractionWithWall<ViscousForce>;

template <typename ViscosityType>
class ViscousForce<Contact<Wall>, ViscosityType> : public BaseViscousForceWithWall
template <typename ViscosityType, class KernelCorrectionType>
class ViscousForce<Contact<Wall>, ViscosityType, KernelCorrectionType>
: public BaseViscousForceWithWall
{
public:
explicit ViscousForce(BaseContactRelation &wall_contact_relation)
: BaseViscousForceWithWall(wall_contact_relation),
mu_(particles_){};
explicit ViscousForce(BaseContactRelation &wall_contact_relation);
virtual ~ViscousForce(){};
void interaction(size_t index_i, Real dt = 0.0);

protected:
ViscosityType mu_;
KernelCorrectionType kernel_correction_;
};

template <typename ViscosityType>
class ViscousForce<Contact<Wall, AngularConservative>, ViscosityType> : public BaseViscousForceWithWall
template <typename ViscosityType, class KernelCorrectionType>
class ViscousForce<Contact<Wall, AngularConservative>, ViscosityType, KernelCorrectionType>
: public BaseViscousForceWithWall
{
public:
explicit ViscousForce(BaseContactRelation &wall_contact_relation);
Expand All @@ -136,11 +136,12 @@ class ViscousForce<Contact<Wall, AngularConservative>, ViscosityType> : public B

protected:
ViscosityType mu_;
StdLargeVec<Vecd> &distance_from_wall_;
KernelCorrectionType kernel_correction_;
};

template <typename ViscosityType>
class ViscousForce<Contact<>, ViscosityType> : public ViscousForce<DataDelegateContact>
template <typename ViscosityType, class KernelCorrectionType>
class ViscousForce<Contact<>, ViscosityType, KernelCorrectionType>
: public ViscousForce<DataDelegateContact>
{
public:
explicit ViscousForce(BaseContactRelation &contact_relation);
Expand All @@ -149,16 +150,19 @@ class ViscousForce<Contact<>, ViscosityType> : public ViscousForce<DataDelegateC

protected:
StdVec<ViscosityType> contact_mu_;
KernelCorrectionType kernel_correction_;
StdVec<KernelCorrectionType> contact_kernel_corrections_;
StdVec<StdLargeVec<Vecd> *> contact_vel_;
StdVec<StdLargeVec<Real> *> wall_Vol_;
};

using ViscousForceWithWall = ComplexInteraction<ViscousForce<Inner<>, Contact<Wall>>, FixedViscosity>;
using MultiPhaseViscousForceWithWall = ComplexInteraction<ViscousForce<Inner<>, Contact<>, Contact<Wall>>, FixedViscosity>;
using ViscousForceWithWall = ComplexInteraction<ViscousForce<Inner<>, Contact<Wall>>, FixedViscosity, NoKernelCorrection>;
using ViscousForceWithWallCorrection = ComplexInteraction<ViscousForce<Inner<>, Contact<Wall>>, FixedViscosity, LinearGradientCorrection>;
using MultiPhaseViscousForceWithWall = ComplexInteraction<ViscousForce<Inner<>, Contact<>, Contact<Wall>>, FixedViscosity, NoKernelCorrection>;

template <typename... FormulationType>
using NonNewtonianViscousForceWithWall =
ComplexInteraction<ViscousForce<Inner<FormulationType...>, Contact<Wall, FormulationType...>>, VariableViscosity>;
ComplexInteraction<ViscousForce<Inner<FormulationType...>, Contact<Wall, FormulationType...>>, VariableViscosity, NoKernelCorrection>;

/**
* @class VorticityInner
Expand Down
86 changes: 54 additions & 32 deletions src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,19 @@ ViscousForce<DataDelegationType>::ViscousForce(BaseRelationType &base_relation)
viscous_force_(*this->particles_->template registerSharedVariable<Vecd>("ViscousForce")),
smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength()) {}
//=================================================================================================//
template <class ViscosityType>
ViscousForce<Inner<>, ViscosityType>::ViscousForce(BaseInnerRelation &inner_relation)
template <typename ViscosityType, class KernelCorrectionType>
ViscousForce<Inner<>, ViscosityType, KernelCorrectionType>::
ViscousForce(BaseInnerRelation &inner_relation)
: ViscousForce<DataDelegateInner>(inner_relation),
ForcePrior(particles_, "ViscousForce"), mu_(particles_)
ForcePrior(particles_, "ViscousForce"), mu_(particles_),
kernel_correction_(particles_)
{
static_assert(std::is_base_of<ParticleAverage, ViscosityType>::value,
"ParticleAverage is not the base of ViscosityType!");
}
//=================================================================================================//
template <class ViscosityType>
void ViscousForce<Inner<>, ViscosityType>::interaction(size_t index_i, Real dt)
template <typename ViscosityType, class KernelCorrectionType>
void ViscousForce<Inner<>, ViscosityType, KernelCorrectionType>::interaction(size_t index_i, Real dt)
{
Vecd force = Vecd::Zero();
const Neighborhood &inner_neighborhood = inner_configuration_[index_i];
Expand All @@ -39,15 +41,23 @@ void ViscousForce<Inner<>, ViscosityType>::interaction(size_t index_i, Real dt)
// viscous force
Vecd vel_derivative = (vel_[index_i] - vel_[index_j]) /
(inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_);
force += 2.0 * mu_(index_i, index_j) * vel_derivative *
force += (kernel_correction_(index_i) + kernel_correction_(index_j)) *
mu_(index_i, index_j) * vel_derivative *
inner_neighborhood.dW_ij_[n] * Vol_[index_j];
}

viscous_force_[index_i] = force * Vol_[index_i];
}
//=================================================================================================//
template <typename ViscosityType>
void ViscousForce<Inner<AngularConservative>, ViscosityType>::interaction(size_t index_i, Real dt)
template <typename ViscosityType, class KernelCorrectionType>
ViscousForce<Inner<AngularConservative>, ViscosityType, KernelCorrectionType>::
ViscousForce(BaseInnerRelation &inner_relation)
: ViscousForce<DataDelegateInner>(inner_relation), ForcePrior(particles_, "ViscousForce"),
mu_(particles_), kernel_correction_(particles_) {}
//=================================================================================================//
template <typename ViscosityType, class KernelCorrectionType>
void ViscousForce<Inner<AngularConservative>, ViscosityType, KernelCorrectionType>::
interaction(size_t index_i, Real dt)
{
Vecd force = Vecd::Zero();
Neighborhood &inner_neighborhood = inner_configuration_[index_i];
Expand All @@ -60,16 +70,24 @@ void ViscousForce<Inner<AngularConservative>, ViscosityType>::interaction(size_t
/** The following viscous force is given in Monaghan 2005 (Rep. Prog. Phys.), it seems that
* this formulation is more accurate than the previous one for Taylor-Green-Vortex flow. */
Real v_r_ij = (vel_[index_i] - vel_[index_j]).dot(e_ij);
Real eta_ij = 2.0 * Real(Dimensions + 2) * mu_(index_i, index_j) * v_r_ij /
Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_j) * v_r_ij /
(r_ij + 0.01 * smoothing_length_);
force += eta_ij * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * e_ij;
force += (kernel_correction_(index_i) + kernel_correction_(index_j)) *
eta_ij * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * e_ij;
}

viscous_force_[index_i] = force * Vol_[index_i];
}
//=================================================================================================//
template <typename ViscosityType>
void ViscousForce<Contact<Wall>, ViscosityType>::interaction(size_t index_i, Real dt)
template <typename ViscosityType, class KernelCorrectionType>
ViscousForce<Contact<Wall>, ViscosityType, KernelCorrectionType>::
ViscousForce(BaseContactRelation &wall_contact_relation)
: BaseViscousForceWithWall(wall_contact_relation),
mu_(particles_), kernel_correction_(particles_) {}
//=================================================================================================//
template <typename ViscosityType, class KernelCorrectionType>
void ViscousForce<Contact<Wall>, ViscosityType, KernelCorrectionType>::
interaction(size_t index_i, Real dt)
{
Vecd force = Vecd::Zero();
for (size_t k = 0; k < contact_configuration_.size(); ++k)
Expand All @@ -84,26 +102,25 @@ void ViscousForce<Contact<Wall>, ViscosityType>::interaction(size_t index_i, Rea

Vecd vel_derivative = 2.0 * (vel_[index_i] - vel_ave_k[index_j]) /
(r_ij + 0.01 * smoothing_length_);
force += 2.0 * mu_(index_i, index_i) * vel_derivative *
contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j];
force += 2.0 * kernel_correction_(index_i) * mu_(index_i, index_i) *
vel_derivative * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j];
}
}

viscous_force_[index_i] += force * Vol_[index_i];
}
//=================================================================================================//
template <typename ViscosityType>
ViscousForce<Contact<Wall, AngularConservative>, ViscosityType>::
template <typename ViscosityType, class KernelCorrectionType>
ViscousForce<Contact<Wall, AngularConservative>, ViscosityType, KernelCorrectionType>::
ViscousForce(BaseContactRelation &wall_contact_relation)
: BaseViscousForceWithWall(wall_contact_relation),
mu_(particles_),
distance_from_wall_(*this->particles_->template registerSharedVariable<Vecd>("DistanceFromWall")){};
mu_(particles_), kernel_correction_(particles_){};
//=================================================================================================//
template <typename ViscosityType>
void ViscousForce<Contact<Wall, AngularConservative>, ViscosityType>::interaction(size_t index_i, Real dt)
template <typename ViscosityType, class KernelCorrectionType>
void ViscousForce<Contact<Wall, AngularConservative>, ViscosityType, KernelCorrectionType>::
interaction(size_t index_i, Real dt)
{
Vecd force = Vecd::Zero();
const Vecd &distance_from_wall = distance_from_wall_[index_i];
for (size_t k = 0; k < contact_configuration_.size(); ++k)
{
StdLargeVec<Vecd> &vel_ave_k = *(wall_vel_ave_[k]);
Expand All @@ -115,46 +132,51 @@ void ViscousForce<Contact<Wall, AngularConservative>, ViscosityType>::interactio
Vecd &e_ij = contact_neighborhood.e_ij_[n];
Real r_ij = contact_neighborhood.r_ij_[n];

Vecd distance_diff = distance_from_wall - r_ij * e_ij;
Real factor = 1.0 - distance_from_wall.dot(distance_diff) / distance_from_wall.squaredNorm();
Real v_r_ij = factor * (vel_[index_i] - vel_ave_k[index_j]).dot(e_ij);
Real eta_ij = 2.0 * Real(Dimensions + 2) * mu_(index_i, index_i) * v_r_ij /
Real v_r_ij = 2.0 * (vel_[index_i] - vel_ave_k[index_j]).dot(e_ij);
Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_i) * v_r_ij /
(r_ij + 0.01 * smoothing_length_);
force += eta_ij * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij;
force += 2.0 * kernel_correction_(index_i) * eta_ij *
contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij;
}
}

viscous_force_[index_i] += force * Vol_[index_i];
}
//=================================================================================================//
template <typename ViscosityType>
ViscousForce<Contact<>, ViscosityType>::ViscousForce(BaseContactRelation &contact_relation)
: ViscousForce<DataDelegateContact>(contact_relation)
template <typename ViscosityType, class KernelCorrectionType>
ViscousForce<Contact<>, ViscosityType, KernelCorrectionType>::
ViscousForce(BaseContactRelation &contact_relation)
: ViscousForce<DataDelegateContact>(contact_relation),
kernel_correction_(particles_)
{
for (size_t k = 0; k != contact_particles_.size(); ++k)
{
contact_mu_.emplace_back(ViscosityType(particles_, contact_particles_[k]));
contact_kernel_corrections_.emplace_back(KernelCorrectionType(contact_particles_[k]));
contact_vel_.push_back(contact_particles_[k]->template getVariableDataByName<Vecd>("Velocity"));
wall_Vol_.push_back(contact_particles_[k]->template getVariableDataByName<Real>("VolumetricMeasure"));
}
}
//=================================================================================================//
template <typename ViscosityType>
void ViscousForce<Contact<>, ViscosityType>::interaction(size_t index_i, Real dt)
template <typename ViscosityType, class KernelCorrectionType>
void ViscousForce<Contact<>, ViscosityType, KernelCorrectionType>::
interaction(size_t index_i, Real dt)
{
Vecd force = Vecd::Zero();
for (size_t k = 0; k < contact_configuration_.size(); ++k)
{
auto &contact_mu_k = contact_mu_[k];
StdLargeVec<Vecd> &vel_k = *(contact_vel_[k]);
StdLargeVec<Real> &wall_Vol_k = *(wall_Vol_[k]);
KernelCorrectionType &kernel_correction_k = contact_kernel_corrections_[k];
Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
{
size_t index_j = contact_neighborhood.j_[n];
Vecd vel_derivative = (vel_[index_i] - vel_k[index_j]) /
(contact_neighborhood.r_ij_[n] + 0.01 * smoothing_length_);
force += 2.0 * contact_mu_k(index_i, index_j) * vel_derivative *
force += (kernel_correction_(index_i) + kernel_correction_k(index_j)) *
contact_mu_k(index_i, index_j) * vel_derivative *
contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j];
}
}
Expand Down
Loading

0 comments on commit f6f70dd

Please sign in to comment.