4

我一直在努力使用以下代码。它是Forward-Euler 算法的 F# 实现，用于对在引力场中移动的恒星进行建模。

let force (b1:Body) (b2:Body) = 
    let r = (b2.Position - b1.Position)
    let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
    if (b1 = b2) then
        VectorFloat.Zero
    else
        r * (b1.Mass * b2.Mass) / (Math.Sqrt((float)rm) * (float)rm)    

member this.Integrate(dT, (bodies:Body[])) = 

    for i = 0 to bodies.Length - 1 do
        for j = (i + 1) to bodies.Length - 1 do
            let f = force bodies.[i] bodies.[j]
            bodies.[i].Acceleration <- bodies.[i].Acceleration + (f / bodies.[i].Mass)
            bodies.[j].Acceleration <- bodies.[j].Acceleration - (f / bodies.[j].Mass)
        bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
        bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT

虽然这可行，但它并不完全是“功能性的”。它的性能也很糟糕，比等效的 c# 代码慢 2.5 倍。body 是 Body 类型的结构数组。

我正在努力解决的问题是 force() 是一个昂贵的函数，所以通常你为每对计算一次并依赖 Fij = -Fji 的事实。但这确实会打乱任何展开的循环等。

建议已收到！不，这不是家庭作业...

谢谢，

阿德

更新：澄清 Body 和 VectorFloat 被定义为 C# 结构。这是因为程序在 F#/C# 和 C++/CLI 之间进行互操作。最终，我将在 BitBucket 上获取代码，但这是一项正在进行的工作，我有一些问题需要解决，然后才能将其发布。

[StructLayout(LayoutKind.Sequential)]
public struct Body
{
    public VectorFloat Position;
    public float Size;
    public uint Color;

    public VectorFloat Velocity;
    public VectorFloat Acceleration;
              '''
}   

[StructLayout(LayoutKind.Sequential)]
public partial struct VectorFloat
{
    public System.Single X { get; set; }
    public System.Single Y { get; set; }
    public System.Single Z { get; set; }
}

向量定义了标准向量类所期望的运算符类型。对于这种情况，您可能可以使用 .NET 框架中的 Vector3D 类（我实际上正在研究切换到它）。

更新 2：根据以下前两个回复改进了代码：

    for i = 0 to bodies.Length - 1 do
    for j = (i + 1) to bodies.Length - 1 do
        let r = ( bodies.[j].Position -  bodies.[i].Position)
        let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
        let f = r / (Math.Sqrt((float)rm) * (float)rm)    
        bodies.[i].Acceleration <- bodies.[i].Acceleration + (f * bodies.[j].Mass)
        bodies.[j].Acceleration <- bodies.[j].Acceleration - (f * bodies.[i].Mass)
    bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
    bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT

  • force 函数中覆盖 b1 == b2 情况的分支是最严重的违规者。如果 softeningLength 始终非零，即使它非常小（Epsilon），您也不需要它。此优化在 C# 代码中，但不在 F# 版本中（doh！）。

  • Math.Pow(x, -1.5) 似乎比 1/ (Math.Sqrt(x) * x) 慢很多。本质上，这个算法有点奇怪，因为它的性能取决于这一步的成本。

  • 将力计算内联并消除一些分歧也带来了一些改进，但性能确实被分支杀死并且由 Sqrt 的成本主导。

WRT 在结构上使用类：在某些情况下（此代码的 CUDA 和本机 C++ 实现以及 DX9 渲染器），我需要将主体数组放入非托管代码或 GPU 上。在这些情况下，能够 memcpy 一个连续的内存块似乎是要走的路。不是我从一组 Body 中得到的东西。

4

2 回答 2

3

我不确定以功能样式重写此代码是否明智。我已经看到一些尝试以函数方式编写对交互计算的尝试，并且它们中的每一个都比两个嵌套循环更难遵循。

在查看结构与类之前（我相信其他人对此有聪明的说法），也许您可​​以尝试优化计算本身？

您正在计算两个加速度增量，我们称它们为 dAi 和 dAj：

dAi = r*m1*m2/(rm*sqrt(rm)) / m1

dAj = r*m1*m2/(rm*sqrt(rm)) / m2

[注：m1 = body.[i].mass，m2=body.[j].mass]]

质量除法取消如下：

dAi = r m2 / (rm sqrt(rm))

dAj = r m1 / (rm sqrt(rm))

现在您只需计算每对 (i,j) 的 r/(rm sqrt(rm))。 这可以进一步优化，因为1/(rm sqrt(rm)) = 1/(rm^1.5) = rm^-1.5，所以如果让r' = r * (rm ** -1.5)，那么 编辑：不，它不能，那是过早的优化（见评论）。计算 r' = 1.0 / (r * sqrt r) 是最快的。

dAi = m2 * r'

dAj = m1 * r'

您的代码将变成类似

member this.Integrate(dT, (bodies:Body[])) = 
    for i = 0 to bodies.Length - 1 do
        for j = (i + 1) to bodies.Length - 1 do
            let r = (b2.Position - b1.Position)
            let rm = (float32)r.MagnitudeSquared + softeningLengthSquared            
            let r' = r * (rm ** -1.5)
            bodies.[i].Acceleration <- bodies.[i].Acceleration + r' * bodies.[j].Mass
            bodies.[j].Acceleration <- bodies.[j].Acceleration - r' * bodies.[i].Mass
        bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
        bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT

看，妈，别再分裂了！

警告：未经测试的代码。尝试自担风险。

于 2010-01-31T09:44:59.200 回答
1

我想使用您的代码，但这很困难，因为缺少 Body 和 FloatVector 的定义，而且您指向的原始博客文章似乎也缺少它们。

我猜测您可以使用 F# 的惰性计算以更实用的方式改进性能并重写：http: //msdn.microsoft.com/en-us/library/dd233247 (VS.100).aspx

这个想法非常简单，您将任何可以在惰性 ( ... ) 表达式中重复计算的昂贵计算包装起来，然后您可以根据需要强制计算多次，并且只会计算一次。

于 2010-01-31T07:55:27.663 回答