1 #include "sph/kernel/GravityKernel.h" 2 #include "catch.hpp" 3 #include "math/Functional.h" 4 #include "tests/Approx.h" 5 #include "utils/SequenceTest.h" 6 7 using namespace Sph; 8 9 static_assert(IsGravityKernel<GravityKernel<CubicSpline<3>>>::value, "GravityKernel test failed"); 10 static_assert(IsGravityKernel<GravityKernel<ThomasCouchmanKernel<3>>>::value, "GravityKernel test failed"); 11 static_assert(IsGravityKernel<SolidSphereKernel>::value, "GravityKernel test failed"); 12 static_assert(!IsGravityKernel<GravityLutKernel>::value, "GravityKernel test failed"); 13 static_assert(!IsGravityKernel<CubicSpline<3>>::value, "GravityKernel test failed"); 14 static_assert(!IsGravityKernel<LutKernel<3>>::value, "GravityKernel test failed"); 15 16 using KernelImpl = decltype(getAssociatedGravityKernel(std::declval<CubicSpline<3>>(), std::declval<Size>())); 17 static_assert(IsGravityKernel<KernelImpl>::value, "GravityKernel test failed"); 18 19 TEST_CASE("Default kernel", "[kernel]") { 20 GravityLutKernel kernel; 21 REQUIRE(kernel.radius() == 0._f); 22 REQUIRE(kernel.value(Vector(2._f, 0._f, 0._f), 1._f) == -0.5_f); 23 REQUIRE(kernel.grad(Vector(2._f, 0._f, 0._f), 1._f) == Vector(0.25_f, 0._f, 0._f)); 24 } 25 26 TEST_CASE("M4 gravity kernel", "[kernel]") { 27 GravityLutKernel kernel = GravityKernel<CubicSpline<3>>(); 28 REQUIRE(kernel.value(Vector(1._f, 0._f, 0._f), 0.1_f) < 0._f); 29 REQUIRE(kernel.value(Vector(3._f, 0._f, 0._f), 5._f) < 0._f); 30 REQUIRE(kernel.value(Vector(0._f, 5._f, 0._f), 1._f) == -0.2_f); 31 REQUIRE(kernel.value(Vector(0._f, 5._f, 0._f), EPS) == -0.2_f); 32 33 REQUIRE_NOTHROW(kernel.grad(Vector(1._f, 0._f, 0._f), 1._f)); 34 REQUIRE(kernel.grad(Vector(0._f), 1._f) == Vector(0._f)); 35 REQUIRE(kernel.grad(Vector(0._f, 0._f, 5._f), 1._f) == approx(Vector(0._f, 0._f, 0.04_f))); 36 REQUIRE(kernel.grad(Vector(0._f, 0._f, 5._f), EPS) == approx(Vector(0._f, 0._f, 0.04_f))); 37 38 // check that value = int grad dx 39 __anon758436a00102(const Float x1, const Float x2, const Float h) 40 auto test = [&](const Float x1, const Float x2, const Float h) { 41 const Float lhs = integrate(Interval(x1, x2), [&](const Float r) { // 42 return kernel.grad(Vector(0._f, r, 0._f), h)[Y]; 43 }); 44 const Float rhs = kernel.value(Vector(0._f, x2, 0._f), h) - kernel.value(Vector(0._f, x1, 0._f), h); 45 REQUIRE(lhs == approx(rhs, 1.e-6_f)); 46 }; 47 test(0._f, 3._f, 1._f); 48 test(0.2_f, 0.25_f, 0.1_f); 49 test(0.2_f, 5._f, 0.5_f); 50 test(1._f, 6._f, 2._f); 51 } 52 53 TEST_CASE("M4 gravity kernel consistency", "[kernel]") { 54 GravityLutKernel kernel = GravityKernel<CubicSpline<3>>(); 55 CubicSpline<3> m4; 56 // check consistency of kernels; gravity kernel g must satisfy: 57 // W = 1/(4 pi r^2)*d/dr(r^2 dg/dr) 58 const Float x1 = 0.3_f; 59 const Float x2 = 2.5_f; 60 for (Float h : { 0.25_f, 0.5_f, 1._f, 2.3_f }) { __anon758436a00302(const Float r) 61 const Float lhs = integrate(Interval(x1, x2), [&](const Float r) { // 62 return 4._f * PI * sqr(r) * m4.value(Vector(r, 0._f, 0._f), h); 63 }); 64 const Float rhs = sqr(x2) * kernel.grad(Vector(x2, 0._f, 0._f), h)[X] - 65 sqr(x1) * kernel.grad(Vector(x1, 0._f, 0._f), h)[X]; 66 REQUIRE(lhs == approx(rhs, 1.e-6_f)); 67 } 68 } 69 70 TEST_CASE("Check getAssociatedGravityKernel", "[kernel]") { 71 GravityLutKernel kernel1 = GravityKernel<CubicSpline<3>>(); 72 GravityLutKernel kernel2 = getAssociatedGravityKernel(CubicSpline<3>()); 73 74 REQUIRE(kernel1.radius() == kernel2.radius()); 75 for (Float h : { 0.25_f, 0.5_f, 1._f, 2.3_f }) { 76 bool valueMatches = true; 77 bool gradientMatches = true; 78 for (Float x = 0; x < 3._f; x += 0.01_f) { 79 const Float w1 = kernel1.value(Vector(x, 0, 0), h); 80 const Float w2 = kernel2.value(Vector(x, 0, 0), h); 81 valueMatches &= w1 == approx(w2, 1.e-5_f); 82 83 const Float dw1 = kernel1.grad(Vector(x, 0, 0), h)[0]; 84 const Float dw2 = kernel2.grad(Vector(x, 0, 0), h)[0]; 85 gradientMatches &= dw1 == approx(dw2, 1.e-4_f); 86 } 87 REQUIRE(valueMatches); 88 REQUIRE(gradientMatches); 89 } 90 } 91