如何在C++中使用Thrust和odeint求解简单的常微分方程

Question

我正在尝试创建一个简单的程序来熟悉Thrust的GPU计算能力和odeint的ODE求解能力。我希望能够在GPU上使用龙格-库塔方法解决简单的ODE（即dy/dx=3x^2y），希望以后能处理更复杂的问题

#include <boost/lambda/lambda.hpp>
#include <boost/numeric/odeint.hpp>
#include <iostream>
#include <iterator>
#include <algorithm>
using namespace boost::numeric::odeint;
using namespace std;
typedef std::vector<double> state_type;
void sys( state_type &y, state_type &dydx, double x){
    dydx[0] = 3*x*x*y[0];                           // dydx = 3*x^2*y
}
int main(){
    state_type y(3);
    runge_kutta4< state_type > rk4;
    y[0] = 2;                                       // y0 = 2
    double x = 1;                                   // x0 = 1
    double h = 0.1;                                 // h = 0.1
    for (int i = 0; i < 100; i++,x+=h){
        rk4.do_step(sys,y,x,h);
        cout << "(";
        cout << x+h;
        cout << ",";
        cout << y[0]; 
        cout << ")";
        cout << endl;
    }
}

然而，我很难理解推力是如何发挥作用的。我遇到的大多数在线资源都以洛伦兹参数研究为例，但我觉得这对我目前的水平来说太先进了。

我理解设备和主机矢量的概念，但我不明白如何使用GPU来解决我的问题。根据我自己的研究，我已经能够使用CUDA（而不是推力）求解简单的代数（非微分）方程。然而，事实证明，将我对odeint和thrust的知识结合起来比我预期的要困难。

特别是，我对以下内容感到困惑：

1） 龙格-库塔步进的改造

2） 调整系统函数本身（本例中dydx=3*x*x*y[0]）。

3） 将odeint和Thrust/boost目录都包含到程序中

如果这个问题太基本或要求太多，我深表歉意；我是StackOverflow的新手，还没有学会所有的"提问"协议，也没有学会我应该在多大程度上自己解决这个问题。

Answer 1

这个问题有点令人困惑。如果你想使用GPU，你通常有大型的微分方程系统。只有三个变量通常是不够的。GPU上的一条指令速度较慢，但它可以在一条指令中并行执行许多操作。Thrust设计用于处理大型数据结构，如GPU上有许多条目的矢量。

简而言之，要回答您的问题，您需要

1. 将thrust_algebra和thrust_operations添加到RK步进器的定义中
2. 利用推力实现系统功能，这是最困难的步骤，以及
3. 将#include <boost/numeric/odeint.hpp>和#include <boost/numeric/odeint/external/thrust.hpp>添加到源文件中。当然，您还需要链接到CUDA库，并使用nvcc编译所有内容。在odeint的examples目录中有一个makefile显示了它是如何工作的