1 #pragma once 2 3 /// \file Conductivity.h 4 /// \brief Artificial thermal conductivity 5 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz) 6 /// \date 2016-2021 7 8 #include "quantities/Storage.h" 9 #include "sph/equations/DerivativeHelpers.h" 10 #include "sph/equations/EquationTerm.h" 11 12 NAMESPACE_SPH_BEGIN 13 14 /// \brief Artificial thermal conductivity 15 /// 16 /// See Price (2008). 17 class ArtificialConductivity : public IEquationTerm { 18 public: 19 class Derivative : public DerivativeTemplate<Derivative> { 20 private: 21 Float alpha, beta; 22 ArrayView<const Vector> r, v; 23 ArrayView<const Float> m, rho, u, p, cs; 24 ArrayView<Float> du; 25 26 SignalSpeedEnum sig; 27 28 public: Derivative(const RunSettings & settings)29 explicit Derivative(const RunSettings& settings) 30 : DerivativeTemplate<Derivative>(settings) { 31 alpha = settings.get<Float>(RunSettingsId::SPH_AC_ALPHA); 32 beta = settings.get<Float>(RunSettingsId::SPH_AC_BETA); 33 sig = settings.get<SignalSpeedEnum>(RunSettingsId::SPH_AC_SIGNAL_SPEED); 34 } 35 additionalCreate(Accumulated & results)36 INLINE void additionalCreate(Accumulated& results) { 37 results.insert<Float>(QuantityId::ENERGY, OrderEnum::FIRST, BufferSource::SHARED); 38 } 39 additionalInitialize(const Storage & storage,Accumulated & results)40 INLINE void additionalInitialize(const Storage& storage, Accumulated& results) { 41 tie(m, rho, u, p, cs) = storage.getValues<Float>(QuantityId::MASS, 42 QuantityId::DENSITY, 43 QuantityId::ENERGY, 44 QuantityId::PRESSURE, 45 QuantityId::SOUND_SPEED); 46 r = storage.getValue<Vector>(QuantityId::POSITION); 47 v = storage.getDt<Vector>(QuantityId::POSITION); 48 du = results.getBuffer<Float>(QuantityId::ENERGY, OrderEnum::FIRST); 49 } 50 additionalEquals(const Derivative & other)51 INLINE bool additionalEquals(const Derivative& other) const { 52 return alpha == other.alpha && beta == other.beta && sig == other.sig; 53 } 54 55 template <bool Symmetrize> eval(const Size i,const Size j,const Vector & grad)56 INLINE void eval(const Size i, const Size j, const Vector& grad) { 57 const Float eps = 1.e-6_f; 58 const Vector e = (r[i] - r[j]) / (getLength(r[i] - r[j]) + eps); 59 const Float rho_bar = 0.5_f * (rho[i] + rho[j]); 60 Float vu_sig; 61 if (sig == SignalSpeedEnum::PRESSURE_DIFFERENCE) { 62 vu_sig = sgn((p[i] - p[j]) * (u[i] - u[j])) * sqrt(abs(p[i] - p[j]) / rho_bar); 63 } else { 64 SPH_ASSERT(sig == SignalSpeedEnum::VELOCITY_DIFFERENCE); 65 vu_sig = abs(dot(v[i] - v[j], e)); 66 } 67 const Float a = alpha * vu_sig * (u[i] - u[j]); 68 const Float heat = a * dot(e, grad) / rho_bar; 69 du[i] += m[j] * heat; 70 71 if (Symmetrize) { 72 du[j] -= m[j] * heat; 73 } 74 } 75 }; 76 ArtificialConductivity(const RunSettings & settings)77 ArtificialConductivity(const RunSettings& settings) { 78 const SignalSpeedEnum sig = settings.get<SignalSpeedEnum>(RunSettingsId::SPH_AC_SIGNAL_SPEED); 79 const Flags<ForceEnum> forces = settings.getFlags<ForceEnum>(RunSettingsId::SPH_SOLVER_FORCES); 80 if (sig == SignalSpeedEnum::PRESSURE_DIFFERENCE && forces != ForceEnum::PRESSURE) { 81 throw InvalidSetup( 82 "Artificial conductivity with pressure-based signal speed cannot be used with forces other " 83 "than pressure gradient. Consider using the velocity-based signal speed instead."); 84 } 85 } 86 setDerivatives(DerivativeHolder & derivatives,const RunSettings & settings)87 virtual void setDerivatives(DerivativeHolder& derivatives, const RunSettings& settings) override { 88 derivatives.require(makeAuto<Derivative>(settings)); 89 } 90 initialize(IScheduler & UNUSED (scheduler),Storage & UNUSED (storage),const Float UNUSED (t))91 virtual void initialize(IScheduler& UNUSED(scheduler), 92 Storage& UNUSED(storage), 93 const Float UNUSED(t)) override {} 94 finalize(IScheduler & UNUSED (scheduler),Storage & UNUSED (storage),const Float UNUSED (t))95 virtual void finalize(IScheduler& UNUSED(scheduler), 96 Storage& UNUSED(storage), 97 const Float UNUSED(t)) override {} 98 create(Storage & UNUSED (storage),IMaterial & UNUSED (material))99 virtual void create(Storage& UNUSED(storage), IMaterial& UNUSED(material)) const override {} 100 }; 101 102 103 NAMESPACE_SPH_END 104