我的天，它充满了星星：使用 Catch2 测试质点引力

欢迎回到 n 体问题系列。上次，我们设置了项目结构和 CMake。根据一些有用的 Reddit 评论，我稍微调整了我们的 CMake 设置 一下，但我仍然不确定目前的项目结构是否能支撑到项目结束。如果您错过了之前的任何一篇文章，可以在这里找到它们

1. 我的天，它充满了星星：n 体问题系列
2. 我的天，它充满了星星：然后有了 CMake

由于我们要练习 TDD，我们首先需要定义一些测试用例，然后才能实现我们的求解器库。一个简单的场景是两个质量相等（）、距离为（）且初始速度和加速度为零的质点。当然，我们可以用它来做数百个不同复杂的测试。但我希望从两个非常简单且具体的测试开始，即质点仅在 x 和 y 方向上加速和移动。这样，我们可以至少检查基础、假设、公式和实现是否正确。

我们还必须考虑我们用于计算加速度的引力方程的一个异常。当两个质量接近奇点（换句话说，发生碰撞）时会发生什么？如果两个质点处于同一空间位置，则两点之间的距离 r 为零。如果发生这种情况，加速度力 F 将趋于无穷大，因此两个质点的加速度也将趋于无穷大。两个质点将彼此分离，飞向虚空。

我们基本上有两种方法来解决这个问题。第一种解决方案是将两个质量合并为一个。任何黑洞爱好者可能都会同意这个解决方案，但有一个小问题。关于 合并恒星 的情况，我们完全不了解，因此很难对其进行建模。当然，我们的模型并不完美，有很多缺点，但模拟合并恒星会极大地增加其复杂性。第二个解决方案很简单，我们从不允许质量碰撞。这可以通过在引力中添加一个误差项 来实现。此外，我们必须添加 来获得力的正确符号。结果，我们得到了质点 j 对质点 i 施加的力 。

我们现在该如何开始？我们有几种选择来实现所谓的 显式和隐式 方法，例如 欧拉、欧拉-克罗默、蛙跳、Verlet 等等。我们将从最简单但相当不精确的显式欧拉方法开始。

让我们从 solverTest 可执行文件中一个名为“具有两个质点的显式欧拉算法”的简单测试用例开始。

TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) {

    Solver solver;

    ParticleBuilder particleBuilder;

    Particle particleA = particleBuilder
            .position()
            .mass()
            .build();

    Particle particleB = particleBuilder
            .position()
            .mass()
            .build();

    std::vector<Particle> particles = { particleA, particleB };

    const double EPSILON = 1e-3;

    SECTION( "Two still standing point mass are attracting each other" ) {
        std::vector<Particle> result = solver.solve(particles, EPSILON);

        REQUIRE(false);
    }
}

这个测试用例描述了我们迄今为止如何粗略地设计 solver。我们需要一个 Solver 来执行计算，一个 Particle 的向量，以及一个具有某种 流畅接口 的 ParticleBuilder，以便轻松为我们生成粒子。欧拉算法不仅需要计算粒子，还需要一个时间步长 。正如您在我第一篇 文章 中回忆的那样，我们将在源代码中称之为 EPSILON（不要与 混淆，我们需要它来解决奇点）。您可以下载 v0.2.0 来获取此初始状态。

我们的第一个测试将检查两个粒子在 x 方向测试中的以下结果：

在纯粹的 x 方向和 y 方向运动的质点的情况下，所有 3 个结果都是有效的。在我们实现了测试并使其通过后，编译 solverTest 文件如下所示：

double Inverse(double value) { return -value; }

Particle GenerateStandardParticle(double xPosition, double yPosition)
{
    ParticleBuilder particleBuilder;

    return particleBuilder
            .position({xPosition, yPosition})
            .mass(1e10)
            .build();
}

TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) 
{
    Solver solver;

    const std::vector<Particle> particlesX = { GenerateStandardParticle(0.5, 0),
                                               GenerateStandardParticle(-0.5, 0)};

    const std::vector<Particle> particlesY = { GenerateStandardParticle(0, 0.5),
                                               GenerateStandardParticle(0, -0.5)};

    const double EPSILON = 1e-3;

    //Solution
    const double acceleration = -0.6674079993; //m/s^2
    const double velocity = -0.06674079993; //m/s
    const double position = 0.48665184; //m

    SECTION( "Two still standing point mass are attracting each other in x-direction" ) {
        std::vector<Particle> result = solver.solve(particlesX, EPSILON);

        Particle &particle = result.front();
        REQUIRE(particle.acceleration.x == Approx(acceleration));
        REQUIRE(particle.velocity.x == Approx(velocity));
        REQUIRE(particle.position.x == Approx(position));

        particle = result.back();
        REQUIRE(particle.acceleration.x == Approx(Inverse(acceleration)));
        REQUIRE(particle.velocity.x == Approx(Inverse(velocity)));
        REQUIRE(particle.position.x == Approx(Inverse(position)));
    }

    SECTION( "Two still standing point mass are attracting each other in y-direction" ) {
        std::vector<Particle> result = solver.solve(particlesY, EPSILON);

        Particle &particle = result.front();
        REQUIRE(particle.acceleration.y == Approx(acceleration));
        REQUIRE(particle.velocity.y == Approx(velocity));
        REQUIRE(particle.position.y == Approx(position));

        particle = result.back();
        REQUIRE(particle.acceleration.y == Approx(Inverse(acceleration)));
        REQUIRE(particle.velocity.y == Approx(Inverse(velocity)));
        REQUIRE(particle.position.y == Approx(Inverse(position)));
    }
}

让我们先看看 SECTION。我们执行两个测试，一个测试纯 x 方向的加速度、速度和位置，另一个测试纯 y 方向。两个测试都使用 Inverse() 函数，该函数反转了预期结果，用于测试第二个粒子。 Catch2 的 Approx 包装类重载了比较运算符，除了涵盖大多数情况的标准配置外，还提供了 3 种不同的配置方法。

1. epsilon：设置允许的相对差值：100.5 == Approx(100).epsilon(0.01)
2. margin：设置允许的绝对差值：104 == Approx(100).margin(5)
3. scale：如果结果值的范围与预期值不同，例如，通过不同的单位系统或单位转换。

为了生成我们的粒子，我们使用 ParticleBuilder，它本身提供了一个 流畅接口，用于轻松配置和生成质点。对于如此简单的案例，使用带有流畅接口的生成器模式并不是必需的，但我们稍后将需要它来进行更高级的生成。

class ParticleBuilder {
public:
    ParticleBuilder() : mMass(0) {}

    ParticleBuilder & position(const Vector2D &position);
    ParticleBuilder & velocity(const Vector2D &velocity);
    ParticleBuilder & acceleration(const Vector2D &acceleration);
    ParticleBuilder & mass(double mass);
    Particle build() const;

private:
    double mMass;
    Vector2D mAcceleration;
    Vector2D mVelocity;
    Vector2D mPosition;
};

Vector2D 是迄今为止一个简单的数据结构，以后可能会扩展更多功能。我们的测试结果正如预期的那样失败了。

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
solverTest is a Catch v2.6.1 host application.
Run with -? for options

-------------------------------------------------------------------------------
Explicit euler algorithm with two point mass
  Two still standing point mass are attracting each other in x-direction
-------------------------------------------------------------------------------
/mnt/c/Develop/gravity/solverTest/src/solverTest.cpp:39
...............................................................................

/mnt/c/Develop/gravity/solverTest/src/solverTest.cpp:43: FAILED:
  REQUIRE( particle.acceleration.x == Approx(acceleration) )
with expansion:
  0.0 == Approx( -0.6674079993 )

-------------------------------------------------------------------------------
Explicit euler algorithm with two point mass
  Two still standing point mass are attracting each other in y-direction
-------------------------------------------------------------------------------
/mnt/c/Develop/gravity/solverTest/src/solverTest.cpp:53
...............................................................................

/mnt/c/Develop/gravity/solverTest/src/solverTest.cpp:57: FAILED:
  REQUIRE( particle.acceleration.y == Approx(acceleration) )
with expansion:
  0.0 == Approx( -0.6674079993 )

===============================================================================
test cases: 1 | 1 failed
assertions: 2 | 2 failed

下次，我们将开始实现相对简单的欧拉方法，并扩展我们的测试，加入一些更有趣的案例，例如奇点和性能测试。像往常一样，您可以通过 GitHub v0.3.0 获取源代码。

您喜欢这篇帖子吗？

您有什么想法吗？

欢迎评论和分享这篇帖子。