1 #include "quantities/Attractor.h"
2 #include "catch.hpp"
3 #include "gravity/BarnesHut.h"
4 #include "gravity/BruteForceGravity.h"
5 #include "gravity/CachedGravity.h"
6 #include "gravity/Moments.h"
7 #include "quantities/Quantity.h"
8 #include "quantities/Storage.h"
9 #include "sph/initial/Distribution.h"
10 #include "tests/Approx.h"
11 #include "thread/Pool.h"
12 #include "utils/SequenceTest.h"
13
14 using namespace Sph;
15
16 template <typename T>
17 AutoPtr<IGravity> createGravity();
18
19 template <>
createGravity()20 AutoPtr<IGravity> createGravity<BruteForceGravity>() {
21 return makeAuto<BruteForceGravity>(1._f);
22 }
23
24 template <>
createGravity()25 AutoPtr<IGravity> createGravity<BarnesHut>() {
26 return makeAuto<BarnesHut>(0.4_f, MultipoleOrder::OCTUPOLE, 25, 50, 1._f);
27 }
28
29 template <>
createGravity()30 AutoPtr<IGravity> createGravity<CachedGravity>() {
31 return makeAuto<CachedGravity>(0.5_f, createGravity<BruteForceGravity>());
32 }
33
34 template <typename T>
35 const Float gravityEps = EPS;
36 template <>
37 const Float gravityEps<BarnesHut> = 2.e-4_f;
38
39 TEMPLATE_TEST_CASE("Gravity with attractors", "[gravity]", BruteForceGravity, BarnesHut, CachedGravity) {
40 ThreadPool& pool = *ThreadPool::getGlobalInstance();
41
42 RandomDistribution distr(1234);
43 SphericalDomain domain1(Vector(0._f), 1.e6_f);
44 SphericalDomain domain2(Vector(1.e6_f, 0._f, 0._f), 5.e6_f);
45 Array<Vector> points1 = distr.generate(pool, 100, domain1);
46 Array<Vector> points2 = distr.generate(pool, 20, domain2);
47 const Float m1 = 3.e10_f;
48 const Float m2 = 1.5e10_f;
49
50 Storage storage1;
51 // combine particles and attractors in storage1
52 storage1.insert<Vector>(QuantityId::POSITION, OrderEnum::SECOND, points1.clone());
53 storage1.insert<Float>(QuantityId::MASS, OrderEnum::ZERO, m1);
54 for (const Vector& p : points2) {
55 storage1.addAttractor(Attractor(p, Vector(0._f), p[H], m2));
56 }
57
58 Storage storage2;
59 // put only particles in storage2
60 Array<Vector> allPoints;
61 allPoints.pushAll(points1);
62 allPoints.pushAll(points2);
63 storage2.insert<Vector>(QuantityId::POSITION, OrderEnum::SECOND, std::move(allPoints));
64 Array<Float> m(storage2.getParticleCnt());
65 for (Size i = 0; i < m.size(); ++i) {
66 m[i] = i < points1.size() ? m1 : m2;
67 }
68 storage2.insert<Float>(QuantityId::MASS, OrderEnum::ZERO, std::move(m));
69
70 Statistics stats;
71 stats.set(StatisticsId::RUN_TIME, 0._f);
72 ArrayView<Vector> dv1 = storage1.getD2t<Vector>(QuantityId::POSITION);
73 AutoPtr<IGravity> gravity = createGravity<TestType>();
74 gravity->build(pool, storage1);
75 gravity->evalSelfGravity(pool, dv1, stats);
76 gravity->evalAttractors(pool, storage1.getAttractors(), dv1);
77
78 ArrayView<Vector> dv2 = storage2.getD2t<Vector>(QuantityId::POSITION);
79 gravity->build(pool, storage2);
80 gravity->evalSelfGravity(pool, dv2, stats);
81
__anondcdac9d90102(const Size i) 82 auto test1 = [&](const Size i) -> Outcome {
83 if (dv2[i] != approx(dv1[i], gravityEps<TestType>)) {
84 return makeFailed("Incorrect acceleration: {} == {}", dv2[i], dv1[i]);
85 }
86 return SUCCESS;
87 };
88 REQUIRE_SEQUENCE(test1, 0, dv1.size());
89
90 ArrayView<const Attractor> attractors = storage1.getAttractors();
__anondcdac9d90202(const Size i) 91 auto test2 = [&](const Size i) -> Outcome {
92 const Vector acc1 = dv2[i + dv1.size()];
93 const Vector acc2 = attractors[i].acceleration;
94 if (acc1 != approx(acc2, gravityEps<TestType>)) {
95 return makeFailed("Incorrect acceleration: {} == {}", acc1, acc2);
96 }
97 return SUCCESS;
98 };
99 REQUIRE_SEQUENCE(test2, 0, attractors.size());
100 }
101