1 #pragma once 2 3 /// \file BruteForceGravity.h 4 /// \brief Simple gravity solver evaluating all particle pairs 5 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz) 6 /// \date 2016-2021 7 8 #include "gravity/IGravity.h" 9 #include "physics/Constants.h" 10 #include "quantities/Attractor.h" 11 #include "quantities/Storage.h" 12 #include "sph/kernel/GravityKernel.h" 13 #include "thread/ThreadLocal.h" 14 15 NAMESPACE_SPH_BEGIN 16 17 /// \brief Computes gravitational acceleration by summing up forces from all particle pairs. 18 /// 19 /// This implementation is not intended for high-performance code because of the O(N) complexity. Useful for 20 /// testing and debugging purposes. 21 class BruteForceGravity : public IGravity { 22 private: 23 ArrayView<const Vector> r; 24 ArrayView<const Float> m; 25 26 GravityLutKernel kernel; 27 Float G = Constants::gravity; 28 29 public: 30 /// \brief Default-construced gravity, assuming point-like particles 31 BruteForceGravity(const Float gravityContant = Constants::gravity) G(gravityContant)32 : G(gravityContant) { 33 SPH_ASSERT(kernel.radius() == 0._f); 34 } 35 36 /// \brief Constructs gravity using smoothing kernel 37 BruteForceGravity(GravityLutKernel&& kernel, const Float gravityContant = Constants::gravity) kernel(std::move (kernel))38 : kernel(std::move(kernel)) 39 , G(gravityContant) {} 40 build(IScheduler & UNUSED (scheduler),const Storage & storage)41 virtual void build(IScheduler& UNUSED(scheduler), const Storage& storage) override { 42 r = storage.getValue<Vector>(QuantityId::POSITION); 43 m = storage.getValue<Float>(QuantityId::MASS); 44 } 45 evalSelfGravity(IScheduler & scheduler,ArrayView<Vector> dv,Statistics & UNUSED (stats))46 virtual void evalSelfGravity(IScheduler& scheduler, 47 ArrayView<Vector> dv, 48 Statistics& UNUSED(stats)) const override { 49 SPH_ASSERT(r.size() == dv.size()); 50 SymmetrizeSmoothingLengths<const GravityLutKernel&> symmetricKernel(kernel); 51 parallelFor(scheduler, 0, r.size(), [&dv, &symmetricKernel, this](const Size i) { 52 dv[i] += this->evalImpl(symmetricKernel, r[i], i); 53 }); 54 } 55 evalAttractors(IScheduler & scheduler,ArrayView<Attractor> attractors,ArrayView<Vector> dv)56 virtual void evalAttractors(IScheduler& scheduler, 57 ArrayView<Attractor> attractors, 58 ArrayView<Vector> dv) const override { 59 SymmetrizeSmoothingLengths<const GravityLutKernel&> symmetricKernel(kernel); 60 // attractor-particle interactions 61 for (Attractor& a : attractors) { 62 parallelFor(scheduler, 0, r.size(), [&dv, &a, this, &symmetricKernel](const Size i) { 63 const Vector f = G * symmetricKernel.grad(r[i], setH(a.position, a.radius)); 64 dv[i] -= a.mass * f; 65 a.acceleration += m[i] * f; 66 }); 67 } 68 // attractor-attractor interactions 69 for (Size i = 0; i < attractors.size(); ++i) { 70 for (Size j = i + 1; j < attractors.size(); ++j) { 71 Attractor& a1 = attractors[i]; 72 Attractor& a2 = attractors[j]; 73 const Vector f = 74 G * symmetricKernel.grad(setH(a1.position, a1.radius), setH(a2.position, a2.radius)); 75 a1.acceleration -= attractors[j].mass * f; 76 a2.acceleration += attractors[i].mass * f; 77 } 78 } 79 } 80 evalAcceleration(const Vector & r0)81 virtual Vector evalAcceleration(const Vector& r0) const override { 82 struct NoSymmetrization { 83 const GravityLutKernel& kernel; 84 85 INLINE Vector grad(const Vector& r1, const Vector& r2) const { 86 return kernel.grad(r1 - r2, r1[H]); 87 } 88 }; 89 return this->evalImpl(NoSymmetrization{ kernel }, r0, Size(-1)); 90 } 91 evalEnergy(IScheduler & scheduler,Statistics & UNUSED (stats))92 virtual Float evalEnergy(IScheduler& scheduler, Statistics& UNUSED(stats)) const override { 93 SymmetrizeSmoothingLengths<const GravityLutKernel&> actKernel(kernel); 94 ThreadLocal<Float> energy(scheduler, 0._f); 95 parallelFor(scheduler, energy, 0, r.size(), [&actKernel, this](const Size i, Float& e) { 96 for (Size j = 0; j < m.size(); ++j) { 97 if (i != j) { 98 e += m[i] * m[j] * actKernel.value(r[i], r[j]); 99 } 100 } 101 }); 102 return 0.5_f * G * energy.accumulate(); 103 } 104 getFinder()105 virtual RawPtr<const IBasicFinder> getFinder() const override { 106 return nullptr; 107 } 108 109 private: 110 template <typename TKernel> evalImpl(const TKernel & actKernel,const Vector & r0,const Size idx)111 INLINE Vector evalImpl(const TKernel& actKernel, const Vector& r0, const Size idx) const { 112 SPH_ASSERT(r && m); 113 Vector a(0._f); 114 115 if (idx != Size(-1)) { 116 // do 2 for loops to avoid the if 117 for (Size i = 0; i < idx; ++i) { 118 a += m[i] * actKernel.grad(r[i], r0); 119 } 120 for (Size i = idx + 1; i < r.size(); ++i) { 121 a += m[i] * actKernel.grad(r[i], r0); 122 } 123 } else { 124 for (Size i = 0; i < r.size(); ++i) { 125 a += m[i] * actKernel.grad(r[i], r0); 126 } 127 } 128 return G * a; 129 } 130 }; 131 132 NAMESPACE_SPH_END 133