Skip to content

Commit

Permalink
Merge pull request #384 from Xiangyu-Hu/fix/dynamics_level_diffusion
Browse files Browse the repository at this point in the history 
Fix/dynamics level diffusion
  • Loading branch information
@Xiangyu-Hu
Xiangyu-Hu authored Jul 28, 2023
2 parents 172bbbf + c200d58 commit 747627b
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 53 deletions.
16 changes: 8 additions & 8 deletions src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,14 @@ class BaseDiffusionRelaxation
explicit BaseDiffusionRelaxation(SPHBody &sph_body);
virtual ~BaseDiffusionRelaxation(){};
StdVec<BaseDiffusion *> &AllDiffusions() { return material_.AllDiffusions(); };
void initialization(size_t index_i, Real dt = 0.0);
void update(size_t index_i, Real dt = 0.0);
};

class KernelGradientInner
{
public:
KernelGradientInner(BaseParticles *inner_particles){};
explicit KernelGradientInner(BaseParticles *inner_particles){};
Vecd operator()(size_t index_i, size_t index_j, Real dW_ijV_j, const Vecd &e_ij)
{
return dW_ijV_j * e_ij;
Expand All @@ -91,7 +93,7 @@ class CorrectedKernelGradientInner
StdLargeVec<Matd> &B_;

public:
CorrectedKernelGradientInner(BaseParticles *inner_particles)
explicit CorrectedKernelGradientInner(BaseParticles *inner_particles)
: B_(*inner_particles->getVariableByName<Matd>("CorrectionMatrix")){};
Vecd operator()(size_t index_i, size_t index_j, Real dW_ijV_j, const Vecd &e_ij)
{
Expand All @@ -110,16 +112,13 @@ class DiffusionRelaxationInner
{
protected:
KernelGradientType kernel_gradient_;
void initializeDiffusionChangeRate(size_t particle_i);
void getDiffusionChangeRate(size_t particle_i, size_t particle_j, Vecd &e_ij, Real surface_area_ij);
virtual void updateSpeciesDiffusion(size_t particle_i, Real dt);

public:
typedef BaseInnerRelation BodyRelationType;
explicit DiffusionRelaxationInner(BaseInnerRelation &inner_relation);
virtual ~DiffusionRelaxationInner(){};
inline void interaction(size_t index_i, Real dt = 0.0);
void update(size_t index_i, Real dt = 0.0);
};

class KernelGradientContact
Expand Down Expand Up @@ -258,7 +257,6 @@ class SecondStageRK2 : public FirstStageType
{
protected:
StdVec<StdLargeVec<Real>> &diffusion_species_s_;
virtual void updateSpeciesDiffusion(size_t particle_i, Real dt) override;

public:
template <typename... ContactArgsType>
Expand All @@ -267,6 +265,8 @@ class SecondStageRK2 : public FirstStageType
: FirstStageType(body_relation, std::forward<ContactArgsType>(contact_args)...),
diffusion_species_s_(diffusion_species_s){};
virtual ~SecondStageRK2(){};

void update(size_t index_i, Real dt = 0.0);
};

/**
Expand All @@ -280,8 +280,8 @@ class DiffusionRelaxationRK2 : public BaseDynamics<void>
protected:
StdVec<StdLargeVec<Real>> diffusion_species_s_; /**< Intermediate state */
SimpleDynamics<InitializationRK<typename FirstStageType::InnerParticlesType>> rk2_initialization_;
InteractionWithUpdate<FirstStageType> rk2_1st_stage_;
InteractionWithUpdate<SecondStageRK2<FirstStageType>> rk2_2nd_stage_;
Dynamics1Level<FirstStageType> rk2_1st_stage_;
Dynamics1Level<SecondStageRK2<FirstStageType>> rk2_2nd_stage_;
StdVec<BaseDiffusion *> all_diffusions_;

public:
Expand Down
59 changes: 24 additions & 35 deletions src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,45 +42,44 @@ BaseDiffusionRelaxation<ParticlesType>::
std::string &diffusion_species_name = all_species_names[diffusion_species_indexes[i]];
diffusion_dt_[i] = this->particles_->template registerSharedVariable<Real>(diffusion_species_name + "ChangeRate");
}
} //=================================================================================================//
template <class ParticlesType, class KernelGradientType>
DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
DiffusionRelaxationInner(BaseInnerRelation &inner_relation)
: BaseDiffusionRelaxation<ParticlesType>(inner_relation.getSPHBody()),
DataDelegateInner<ParticlesType, DataDelegateEmptyBase>(inner_relation),
kernel_gradient_(this->particles_) {}
}
//=================================================================================================//
template <class ParticlesType, class KernelGradientType>
void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
initializeDiffusionChangeRate(size_t particle_i)
template <class ParticlesType>
void BaseDiffusionRelaxation<ParticlesType>::initialization(size_t index_i, Real dt)
{
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
{
(*this->diffusion_dt_[m])[particle_i] = 0;
(*this->diffusion_dt_[m])[index_i] = 0;
}
}
//=================================================================================================//
template <class ParticlesType, class KernelGradientType>
void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
getDiffusionChangeRate(size_t particle_i, size_t particle_j, Vecd &e_ij, Real surface_area_ij)
template <class ParticlesType>
void BaseDiffusionRelaxation<ParticlesType>::update(size_t index_i, Real dt)
{
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
{
Real diff_coff_ij =
this->all_diffusions_[m]->getInterParticleDiffusionCoff(particle_i, particle_j, e_ij);
StdLargeVec<Real> &gradient_species = *this->gradient_species_[m];
Real phi_ij = gradient_species[particle_i] - gradient_species[particle_j];
(*this->diffusion_dt_[m])[particle_i] += diff_coff_ij * phi_ij * surface_area_ij;
(*this->diffusion_species_[m])[index_i] += dt * (*this->diffusion_dt_[m])[index_i];
}
}
//=================================================================================================//
template <class ParticlesType, class KernelGradientType>
DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
DiffusionRelaxationInner(BaseInnerRelation &inner_relation)
: BaseDiffusionRelaxation<ParticlesType>(inner_relation.getSPHBody()),
DataDelegateInner<ParticlesType, DataDelegateEmptyBase>(inner_relation),
kernel_gradient_(this->particles_) {}
//=================================================================================================//
template <class ParticlesType, class KernelGradientType>
void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
updateSpeciesDiffusion(size_t particle_i, Real dt)
getDiffusionChangeRate(size_t particle_i, size_t particle_j, Vecd &e_ij, Real surface_area_ij)
{
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
{
(*this->diffusion_species_[m])[particle_i] += dt * (*this->diffusion_dt_[m])[particle_i];
Real diff_coff_ij =
this->all_diffusions_[m]->getInterParticleDiffusionCoff(particle_i, particle_j, e_ij);
StdLargeVec<Real> &gradient_species = *this->gradient_species_[m];
Real phi_ij = gradient_species[particle_i] - gradient_species[particle_j];
(*this->diffusion_dt_[m])[particle_i] += diff_coff_ij * phi_ij * surface_area_ij;
}
}
//=================================================================================================//
Expand All @@ -89,8 +88,6 @@ void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
interaction(size_t index_i, Real dt)
{
Neighborhood &inner_neighborhood = this->inner_configuration_[index_i];

initializeDiffusionChangeRate(index_i);
for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
{
size_t index_j = inner_neighborhood.j_[n];
Expand All @@ -104,13 +101,6 @@ void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
}
}
//=================================================================================================//
template <class ParticlesType, class KernelGradientType>
void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
update(size_t index_i, Real dt)
{
updateSpeciesDiffusion(index_i, dt);
}
//=================================================================================================//
template <class ParticlesType, class ContactParticlesType, class KernelGradientType>
BaseDiffusionRelaxationContact<ParticlesType, ContactParticlesType, KernelGradientType>::
BaseDiffusionRelaxationContact(BaseContactRelation &contact_relation)
Expand Down Expand Up @@ -331,14 +321,13 @@ void InitializationRK<ParticlesType>::
}
//=================================================================================================//
template <class FirstStageType>
void SecondStageRK2<FirstStageType>::
updateSpeciesDiffusion(size_t particle_i, Real dt)
void SecondStageRK2<FirstStageType>::update(size_t index_i, Real dt)
{
FirstStageType::update(index_i, dt);
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
{
(*this->diffusion_species_[m])[particle_i] =
0.5 * diffusion_species_s_[m][particle_i] +
0.5 * ((*this->diffusion_species_[m])[particle_i] + dt * (*this->diffusion_dt_[m])[particle_i]);
(*this->diffusion_species_[m])[index_i] = 0.5 * diffusion_species_s_[m][index_i] +
0.5 * (*this->diffusion_species_[m])[index_i];
}
}
//=================================================================================================//
Expand Down
6 changes: 4 additions & 2 deletions tests/2d_examples/test_2d_diffusion_NeumannBC/diffusion_NeumannBC.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,8 @@ class DirichletWallBoundaryInitialCondition
size_t phi_;

public:
DirichletWallBoundaryInitialCondition(SolidBody &diffusion_body) : DiffusionReactionInitialCondition<WallParticles>(diffusion_body)
explicit DirichletWallBoundaryInitialCondition(SolidBody &diffusion_body)
: DiffusionReactionInitialCondition<WallParticles>(diffusion_body)
{
phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"];
}
Expand Down Expand Up @@ -157,7 +158,8 @@ class NeumannWallBoundaryInitialCondition
StdLargeVec<Real> &heat_flux_;

public:
NeumannWallBoundaryInitialCondition(SolidBody &diffusion_body) : DiffusionReactionInitialCondition<WallParticles>(diffusion_body),
explicit NeumannWallBoundaryInitialCondition(SolidBody &diffusion_body)
: DiffusionReactionInitialCondition<WallParticles>(diffusion_body),
heat_flux_(*(particles_->getVariableByName<Real>("HeatFlux")))
{
phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"];
Expand Down
17 changes: 9 additions & 8 deletions tests/2d_examples/test_2d_diffusion_RobinBC/diffusion_RobinBC.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,8 @@ class DirichletWallBoundaryInitialCondition
size_t phi_;

public:
DirichletWallBoundaryInitialCondition(SolidBody &diffusion_body) : DiffusionReactionInitialCondition<WallParticles>(diffusion_body)
explicit DirichletWallBoundaryInitialCondition(SolidBody &diffusion_body)
: DiffusionReactionInitialCondition<WallParticles>(diffusion_body)
{
phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"];
}
Expand Down Expand Up @@ -158,9 +159,10 @@ class RobinWallBoundaryInitialCondition
Real &T_infinity_;

public:
RobinWallBoundaryInitialCondition(SolidBody &diffusion_body) : DiffusionReactionInitialCondition<WallParticles>(diffusion_body),
convection_(*(this->particles_->template getVariableByName<Real>("Convection"))),
T_infinity_(*(this->particles_->template getGlobalVariableByName<Real>("T_infinity")))
explicit RobinWallBoundaryInitialCondition(SolidBody &diffusion_body)
: DiffusionReactionInitialCondition<WallParticles>(diffusion_body),
convection_(*(this->particles_->template getVariableByName<Real>("Convection"))),
T_infinity_(*(this->particles_->template getGlobalVariableByName<Real>("T_infinity")))
{
phi_ = particles_->diffusion_reaction_material_.AllSpeciesIndexMap()["Phi"];
}
Expand Down Expand Up @@ -200,7 +202,7 @@ class DiffusionBodyRelaxation
class TemperatureObserverParticleGenerator : public ObserverParticleGenerator
{
public:
TemperatureObserverParticleGenerator(SPHBody &sph_body) : ObserverParticleGenerator(sph_body)
explicit TemperatureObserverParticleGenerator(SPHBody &sph_body) : ObserverParticleGenerator(sph_body)
{
/** A line of measuring points at the middle line. */
size_t number_of_observation_points = 5;
Expand All @@ -209,9 +211,8 @@ class TemperatureObserverParticleGenerator : public ObserverParticleGenerator

for (size_t i = 0; i < number_of_observation_points; ++i)
{
Vec2d point_coordinate(0.5 * L, range_of_measure * Real(i) /
Real(number_of_observation_points - 1) +
start_of_measure);
Vec2d point_coordinate(
0.5 * L, range_of_measure * Real(i) / Real(number_of_observation_points - 1) + start_of_measure);
positions_.push_back(point_coordinate);
}
}
Expand Down

0 comments on commit 747627b

Please sign in to comment.