Skip to content

Commit

Permalink
Merge pull request #389 from ChenxiZhaoTUM/change_diff_coff_name
Browse files Browse the repository at this point in the history
change getInterParticleDiffusionCoeff() of Dirichlet wall boundary
  • Loading branch information
Xiangyu-Hu authored Aug 1, 2023
2 parents f44c31b + b3312e7 commit 5f36441
Show file tree
Hide file tree
Showing 13 changed files with 49 additions and 49 deletions.
16 changes: 8 additions & 8 deletions src/shared/materials/diffusion_reaction.h
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ class BaseDiffusion : public BaseMaterial
size_t gradient_species_index_;

virtual Real getReferenceDiffusivity() = 0;
virtual Real getInterParticleDiffusionCoff(size_t particle_i, size_t particle_j, Vecd &direction_from_j_to_i) = 0;
virtual Real getInterParticleDiffusionCoeff(size_t particle_i, size_t particle_j, const Vecd &direction_from_j_to_i) = 0;
};

/**
Expand All @@ -81,7 +81,7 @@ class IsotropicDiffusion : public BaseDiffusion
virtual ~IsotropicDiffusion(){};

virtual Real getReferenceDiffusivity() override { return diff_cf_; };
virtual Real getInterParticleDiffusionCoff(size_t particle_i, size_t particle_j, Vecd &direction_from_j_to_i) override
virtual Real getInterParticleDiffusionCoeff(size_t particle_i, size_t particle_j, const Vecd &direction_from_j_to_i) override
{
return diff_cf_;
};
Expand Down Expand Up @@ -117,8 +117,8 @@ class DirectionalDiffusion : public IsotropicDiffusion
return SMAX(diff_cf_, diff_cf_ + bias_diff_cf_);
};

virtual Real getInterParticleDiffusionCoff(size_t particle_index_i,
size_t particle_index_j, Vecd &inter_particle_direction) override
virtual Real getInterParticleDiffusionCoeff(size_t particle_index_i,
size_t particle_index_j, const Vecd &inter_particle_direction) override
{
Vecd grad_ij = transformed_diffusivity_ * inter_particle_direction;
return 1.0 / grad_ij.squaredNorm();
Expand Down Expand Up @@ -147,7 +147,7 @@ class LocalDirectionalDiffusion : public DirectionalDiffusion
virtual void registerReloadLocalParameters(BaseParticles *base_particles) override;
virtual void initializeLocalParameters(BaseParticles *base_particles) override;

virtual Real getInterParticleDiffusionCoff(size_t particle_index_i, size_t particle_index_j, Vecd &inter_particle_direction) override
virtual Real getInterParticleDiffusionCoeff(size_t particle_index_i, size_t particle_index_j, const Vecd &inter_particle_direction) override
{
Matd trans_diffusivity = getAverageValue(local_transformed_diffusivity_[particle_index_i], local_transformed_diffusivity_[particle_index_j]);
Vecd grad_ij = trans_diffusivity * inter_particle_direction;
Expand Down Expand Up @@ -293,10 +293,10 @@ class DiffusionReaction : public BaseMaterialType
*/
Real getDiffusionTimeStepSize(Real smoothing_length)
{
Real diff_coff_max = 0.0;
Real diff_coeff_max = 0.0;
for (size_t k = 0; k < all_diffusions_.size(); ++k)
diff_coff_max = SMAX(diff_coff_max, all_diffusions_[k]->getReferenceDiffusivity());
return 0.5 * smoothing_length * smoothing_length / diff_coff_max / Real(Dimensions);
diff_coeff_max = SMAX(diff_coeff_max, all_diffusions_[k]->getReferenceDiffusivity());
return 0.5 * smoothing_length * smoothing_length / diff_coeff_max / Real(Dimensions);
};

/** Initialize a diffusion material. */
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ class DiffusionRelaxationDirichlet
{
protected:
StdVec<StdVec<StdLargeVec<Real> *>> contact_gradient_species_;
void getDiffusionChangeRateDirichletContact(
void getDiffusionChangeRateDirichlet(
size_t particle_i, size_t particle_j, Vecd &e_ij, Real surface_area_ij,
const StdVec<StdLargeVec<Real> *> &gradient_species_k);

Expand All @@ -199,7 +199,7 @@ class DiffusionRelaxationNeumann
StdVec<StdLargeVec<Vecd> *> contact_n_;

protected:
void getDiffusionChangeRateNeumannContact(size_t particle_i, size_t particle_j,
void getDiffusionChangeRateNeumann(size_t particle_i, size_t particle_j,
Real surface_area_ij_Neumann, StdLargeVec<Real> &heat_flux_k);

public:
Expand All @@ -223,7 +223,7 @@ class DiffusionRelaxationRobin
StdVec<StdLargeVec<Vecd> *> contact_n_;

protected:
void getDiffusionChangeRateRobinContact(size_t particle_i, size_t particle_j, Real surface_area_ij_Robin, StdLargeVec<Real> &convection_k, Real &T_infinity_k);
void getDiffusionChangeRateRobin(size_t particle_i, size_t particle_j, Real surface_area_ij_Robin, StdLargeVec<Real> &convection_k, Real &T_infinity_k);

public:
explicit DiffusionRelaxationRobin(BaseContactRelation &contact_relation);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -75,11 +75,11 @@ void DiffusionRelaxationInner<ParticlesType, KernelGradientType>::
{
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);
Real diff_coeff_ij =
this->all_diffusions_[m]->getInterParticleDiffusionCoeff(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_dt_[m])[particle_i] += diff_coeff_ij * phi_ij * surface_area_ij;
}
}
//=================================================================================================//
Expand Down Expand Up @@ -159,15 +159,15 @@ DiffusionRelaxationDirichlet<ParticlesType, ContactParticlesType, KernelGradient
//=================================================================================================//
template <class ParticlesType, class ContactParticlesType, class KernelGradientType>
void DiffusionRelaxationDirichlet<ParticlesType, ContactParticlesType, KernelGradientType>::
getDiffusionChangeRateDirichletContact(size_t particle_i, size_t particle_j, Vecd &e_ij,
getDiffusionChangeRateDirichlet(size_t particle_i, size_t particle_j, Vecd &e_ij,
Real surface_area_ij, const StdVec<StdLargeVec<Real> *> &gradient_species_k)
{
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);
Real diff_coeff_ij =
this->all_diffusions_[m]->getInterParticleDiffusionCoeff(particle_i, particle_i, e_ij);
Real phi_ij = (*this->gradient_species_[m])[particle_i] - (*gradient_species_k[m])[particle_j];
(*this->diffusion_dt_[m])[particle_i] += diff_coff_ij * phi_ij * surface_area_ij;
(*this->diffusion_dt_[m])[particle_i] += diff_coeff_ij * phi_ij * surface_area_ij;
}
}
//=================================================================================================//
Expand All @@ -189,7 +189,7 @@ void DiffusionRelaxationDirichlet<ParticlesType, ContactParticlesType, KernelGra

const Vecd &grad_ijV_j = this->contact_kernel_gradients_[k](index_i, index_j, dW_ijV_j_, e_ij);
Real area_ij = 2.0 * grad_ijV_j.dot(e_ij) / r_ij_;
getDiffusionChangeRateDirichletContact(index_i, index_j, e_ij, area_ij, gradient_species_k);
getDiffusionChangeRateDirichlet(index_i, index_j, e_ij, area_ij, gradient_species_k);
}
}
}
Expand All @@ -213,7 +213,7 @@ DiffusionRelaxationNeumann<ParticlesType, ContactParticlesType, KernelGradientTy
//=================================================================================================//
template <class ParticlesType, class ContactParticlesType, class KernelGradientType>
void DiffusionRelaxationNeumann<ParticlesType, ContactParticlesType, KernelGradientType>::
getDiffusionChangeRateNeumannContact(size_t particle_i, size_t particle_j,
getDiffusionChangeRateNeumann(size_t particle_i, size_t particle_j,
Real surface_area_ij_Neumann, StdLargeVec<Real> &heat_flux_k)
{
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
Expand Down Expand Up @@ -241,7 +241,7 @@ void DiffusionRelaxationNeumann<ParticlesType, ContactParticlesType, KernelGradi
const Vecd &grad_ijV_j = this->contact_kernel_gradients_[k](index_i, index_j, dW_ijV_j_, e_ij);
Vecd n_ij = n_[index_i] - n_k[index_j];
Real area_ij_Neumann = grad_ijV_j.dot(n_ij);
getDiffusionChangeRateNeumannContact(index_i, index_j, area_ij_Neumann, heat_flux_k);
getDiffusionChangeRateNeumann(index_i, index_j, area_ij_Neumann, heat_flux_k);
}
}
}
Expand All @@ -268,7 +268,7 @@ DiffusionRelaxationRobin<ParticlesType, ContactParticlesType, KernelGradientType
//=================================================================================================//
template <class ParticlesType, class ContactParticlesType, class KernelGradientType>
void DiffusionRelaxationRobin<ParticlesType, ContactParticlesType, KernelGradientType>::
getDiffusionChangeRateRobinContact(size_t particle_i, size_t particle_j,
getDiffusionChangeRateRobin(size_t particle_i, size_t particle_j,
Real surface_area_ij_Robin, StdLargeVec<Real> &convection_k, Real &T_infinity_k)
{
for (size_t m = 0; m < this->all_diffusions_.size(); ++m)
Expand Down Expand Up @@ -299,7 +299,7 @@ void DiffusionRelaxationRobin<ParticlesType, ContactParticlesType, KernelGradien
const Vecd &grad_ijV_j = this->contact_kernel_gradients_[k](index_i, index_j, dW_ijV_j_, e_ij);
Vecd n_ij = n_[index_i] - n_k[index_j];
Real area_ij_Robin = grad_ijV_j.dot(n_ij);
getDiffusionChangeRateRobinContact(index_i, index_j, area_ij_Robin, convection_k, T_infinity_k);
getDiffusionChangeRateRobin(index_i, index_j, area_ij_Robin, convection_k, T_infinity_k);
}
}
}
Expand Down
6 changes: 3 additions & 3 deletions tests/2d_examples/test_0d_regression_test/regression_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(L + BW, H + BW));
//----------------------------------------------------------------------
// Global parameters on material properties
//----------------------------------------------------------------------
Real diffusion_coff = 1.0e-3;
Real bias_coff = 0.0;
Real diffusion_coeff = 1.0e-3;
Real bias_coeff = 0.0;
Real alpha = Pi / 4.0;
Vec2d bias_direction(cos(alpha), sin(alpha));
Real initial_temperature = 0.0;
Expand Down Expand Up @@ -103,7 +103,7 @@ class DiffusionMaterial : public DiffusionReaction<Solid>
public:
DiffusionMaterial() : DiffusionReaction<Solid>({"Phi"}, SharedPtr<NoReaction>())
{
initializeAnDiffusion<DirectionalDiffusion>("Phi", "Phi", diffusion_coff, bias_coff, bias_direction);
initializeAnDiffusion<DirectionalDiffusion>("Phi", "Phi", diffusion_coeff, bias_coeff, bias_direction);
};
};
//----------------------------------------------------------------------
Expand Down
6 changes: 3 additions & 3 deletions tests/2d_examples/test_2d_depolarization/depolarization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@ StdVec<Vecd> observation_location = {Vecd(0.3, 0.7)};
//----------------------------------------------------------------------
// Basic parameters for material properties.
//----------------------------------------------------------------------
Real diffusion_coff = 1.0;
Real bias_coff = 0.0;
Real diffusion_coeff = 1.0;
Real bias_coeff = 0.0;
Vec2d fiber_direction(1.0, 0.0);
Real c_m = 1.0;
Real k = 8.0;
Expand Down Expand Up @@ -83,7 +83,7 @@ int main()
SolidBody muscle_body(system, makeShared<MuscleBlock>("MuscleBlock"));
SharedPtr<AlievPanfilowModel> muscle_reaction_model_ptr = makeShared<AlievPanfilowModel>(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
muscle_body.defineParticlesAndMaterial<ElectroPhysiologyParticles, MonoFieldElectroPhysiology>(
muscle_reaction_model_ptr, TypeIdentity<DirectionalDiffusion>(), diffusion_coff, bias_coff, fiber_direction);
muscle_reaction_model_ptr, TypeIdentity<DirectionalDiffusion>(), diffusion_coeff, bias_coeff, fiber_direction);
muscle_body.generateParticles<ParticleGeneratorLattice>();

ObserverBody voltage_observer(system, "VoltageObserver");
Expand Down
6 changes: 3 additions & 3 deletions tests/2d_examples/test_2d_diffusion/diffusion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ BoundingBox system_domain_bounds(Vec2d(0.0, 0.0), Vec2d(L, H));
//----------------------------------------------------------------------
// Basic parameters for material properties.
//----------------------------------------------------------------------
Real diffusion_coff = 1.0e-4;
Real bias_coff = 0.0;
Real diffusion_coeff = 1.0e-4;
Real bias_coeff = 0.0;
Real alpha = Pi / 6.0;
Vec2d bias_direction(cos(alpha), sin(alpha));
//----------------------------------------------------------------------
Expand Down Expand Up @@ -44,7 +44,7 @@ class DiffusionMaterial : public DiffusionReaction<Solid>
public:
DiffusionMaterial() : DiffusionReaction<Solid>({"Phi"}, SharedPtr<NoReaction>())
{
initializeAnDiffusion<DirectionalDiffusion>("Phi", "Phi", diffusion_coff, bias_coff, bias_direction);
initializeAnDiffusion<DirectionalDiffusion>("Phi", "Phi", diffusion_coeff, bias_coeff, bias_direction);
};
};
using DiffusionParticles = DiffusionReactionParticles<SolidParticles, DiffusionMaterial>;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ BoundingBox system_domain_bounds(Vec2d(-BW, -BW), Vec2d(L + BW, H + BW));
//----------------------------------------------------------------------
// Basic parameters for material properties.
//----------------------------------------------------------------------
Real diffusion_coff = 1;
Real diffusion_coeff = 1;
std::array<std::string, 1> species_name_list{"Phi"};
//----------------------------------------------------------------------
// Initial and boundary conditions.
Expand Down Expand Up @@ -95,7 +95,7 @@ class DiffusionMaterial : public DiffusionReaction<Solid>
public:
DiffusionMaterial() : DiffusionReaction<Solid>({"Phi"}, SharedPtr<NoReaction>())
{
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coff);
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coeff);
}
};
using DiffusionParticles = DiffusionReactionParticles<SolidParticles, DiffusionMaterial>;
Expand Down
4 changes: 2 additions & 2 deletions tests/2d_examples/test_2d_heat_transfer/heat_transfer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ StdVec<Vecd> observation_location = {Vecd(0.0, DH * 0.5)};
//----------------------------------------------------------------------
// Global parameters on the material properties
//----------------------------------------------------------------------
Real diffusion_coff = 1.0e-3;
Real diffusion_coeff = 1.0e-3;
Real rho0_f = 1.0; /**< Density. */
Real U_f = 1.0; /**< Characteristic velocity. */
Real c_f = 10.0 * U_f; /**< Speed of sound. */
Expand Down Expand Up @@ -101,7 +101,7 @@ class ThermofluidBodyMaterial : public DiffusionReaction<WeaklyCompressibleFluid
ThermofluidBodyMaterial()
: DiffusionReaction<WeaklyCompressibleFluid>({"Phi"}, SharedPtr<NoReaction>(), rho0_f, c_f, mu_f)
{
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coff);
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coeff);
};
};
using DiffusionBaseParticles = DiffusionReactionParticles<BaseParticles, ThermofluidBodyMaterial>;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,8 @@ Real poisson = 0.4995;
Real bulk_modulus = 2.0 * a0[0] * (1.0 + poisson) / (3.0 * (1.0 - 2.0 * poisson));
/** Electrophysiology parameters. */
std::array<std::string, 1> species_name_list{"Phi"};
Real diffusion_coff = 0.8;
Real bias_coff = 0.0;
Real diffusion_coeff = 0.8;
Real bias_coeff = 0.0;
/** Electrophysiology parameters. */
Real c_m = 1.0;
Real k = 8.0;
Expand Down Expand Up @@ -72,7 +72,7 @@ class FiberDirectionDiffusion : public DiffusionReaction<LocallyOrthotropicMuscl
{"Phi"}, SharedPtr<NoReaction>(),
rho0_s, bulk_modulus, fiber_direction, sheet_direction, a0, b0)
{
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coff);
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coeff);
};
};
using FiberDirectionDiffusionParticles = DiffusionReactionParticles<ElasticSolidParticles, FiberDirectionDiffusion>;
Expand Down Expand Up @@ -363,7 +363,7 @@ int main(int ac, char *av[])
SharedPtr<AlievPanfilowModel> muscle_reaction_model_ptr = makeShared<AlievPanfilowModel>(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
physiology_heart.defineParticlesAndMaterial<
ElectroPhysiologyParticles, MonoFieldElectroPhysiology>(
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coff, bias_coff, fiber_direction);
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coeff, bias_coeff, fiber_direction);
(!system.RunParticleRelaxation() && system.ReloadParticles())
? physiology_heart.generateParticles<ParticleGeneratorReload>(io_environment, "HeartModel")
: physiology_heart.generateParticles<ParticleGeneratorLattice>();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ int main(int ac, char *av[])
SharedPtr<AlievPanfilowModel> muscle_reaction_model_ptr = makeShared<AlievPanfilowModel>(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
physiology_heart.defineParticlesAndMaterial<
ElectroPhysiologyParticles, MonoFieldElectroPhysiology>(
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coff, bias_coff, fiber_direction);
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coeff, bias_coeff, fiber_direction);
(!system.RunParticleRelaxation() && system.ReloadParticles())
? physiology_heart.generateParticles<ParticleGeneratorReload>(io_environment, "HeartModel")
: physiology_heart.generateParticles<ParticleGeneratorLattice>();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ Real poisson = 0.4995;
Real bulk_modulus = 2.0 * a0[0] * (1.0 + poisson) / (3.0 * (1.0 - 2.0 * poisson));
/** Electrophysiology parameters. */
std::array<std::string, 1> species_name_list{"Phi"};
Real diffusion_coff = 0.8;
Real bias_coff = 0.0;
Real diffusion_coeff = 0.8;
Real bias_coeff = 0.0;
/** Electrophysiology parameters. */
Real c_m = 1.0;
Real k = 8.0;
Expand Down Expand Up @@ -74,7 +74,7 @@ class FiberDirectionDiffusion : public DiffusionReaction<LocallyOrthotropicMuscl
{"Phi"}, SharedPtr<NoReaction>(),
rho0_s, bulk_modulus, fiber_direction, sheet_direction, a0, b0)
{
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coff);
initializeAnDiffusion<IsotropicDiffusion>("Phi", "Phi", diffusion_coeff);
};
};
using FiberDirectionDiffusionParticles = DiffusionReactionParticles<ElasticSolidParticles, FiberDirectionDiffusion>;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ int main(int ac, char *av[])
SharedPtr<AlievPanfilowModel> muscle_reaction_model_ptr = makeShared<AlievPanfilowModel>(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
physiology_heart.defineParticlesAndMaterial<
ElectroPhysiologyParticles, MonoFieldElectroPhysiology>(
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coff, bias_coff, fiber_direction);
muscle_reaction_model_ptr, TypeIdentity<LocalDirectionalDiffusion>(), diffusion_coeff, bias_coeff, fiber_direction);
(!system.RunParticleRelaxation() && system.ReloadParticles())
? physiology_heart.generateParticles<ParticleGeneratorReload>(io_environment, "HeartModel")
: physiology_heart.generateParticles<ParticleGeneratorLattice>();
Expand All @@ -147,7 +147,7 @@ int main(int ac, char *av[])
SharedPtr<AlievPanfilowModel> pkj_reaction_model_ptr = makeShared<AlievPanfilowModel>(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
pkj_body.defineParticlesAndMaterial<
ElectroPhysiologyReducedParticles, MonoFieldElectroPhysiology>(
pkj_reaction_model_ptr, TypeIdentity<DirectionalDiffusion>(), diffusion_coff * acceleration_factor, bias_coff, fiber_direction);
pkj_reaction_model_ptr, TypeIdentity<DirectionalDiffusion>(), diffusion_coeff * acceleration_factor, bias_coeff, fiber_direction);
pkj_body.generateParticles<NetworkGeneratorWithExtraCheck>(starting_point, second_point, 50, 1.0);
TreeTerminates pkj_leaves(pkj_body);
//----------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 5f36441

Please sign in to comment.