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