CN111105341A - 一种低功耗高运算性能求解计算流体动力学的框架方法 - Google Patents

一种低功耗高运算性能求解计算流体动力学的框架方法 Download PDF

Info

Publication number
CN111105341A
CN111105341A CN201911296165.7A CN201911296165A CN111105341A CN 111105341 A CN111105341 A CN 111105341A CN 201911296165 A CN201911296165 A CN 201911296165A CN 111105341 A CN111105341 A CN 111105341A
Authority
CN
China
Prior art keywords
data
substep
fpga
kernel
grid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911296165.7A
Other languages
English (en)
Other versions
CN111105341B (zh
Inventor
严伟安
丁雪海
童维勤
支小莉
程金凤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Transpacific Technology Development Ltd
Original Assignee
Beijing Transpacific Technology Development Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Transpacific Technology Development Ltd filed Critical Beijing Transpacific Technology Development Ltd
Priority to CN201911296165.7A priority Critical patent/CN111105341B/zh
Publication of CN111105341A publication Critical patent/CN111105341A/zh
Application granted granted Critical
Publication of CN111105341B publication Critical patent/CN111105341B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation

Abstract

本发明公开了一种低功耗高运算性能求解计算流体动力学的框架方法，属于计算流体动力学(CFD)领域，在FPGA端优化运行LBM计算流体力学模拟，可以达到低功耗高效能运算，其特征在于，该方法包括6个步骤：1.初始化网格数据和OpenCL环境；2.读取数据并计算密度及速度等宏观量；3.边界检测；4.格点碰撞计算；5.数据传播；6.根据迭代停止条件判断是否继续。其中1、6两个步骤在CPU端执行，其余步骤位于FPGA端执行，2、3两个步骤完全并行，其余各步骤内部也是并行计算。本专利使用的粒子密度分布数据采用乒乓缓存的方法进行存储，以减少片外内存的访问。

Description

一种低功耗高运算性能求解计算流体动力学的框架方法
技术领域
本发明涉及计算流体动力学领域，特别是通过格子玻尔兹曼方法来架构模拟及精确求解计算流体动力学问题及工程应用问题。
背景技术
为精确求解流体力学问题，通常要对流场进行网格划分，达到介观模拟尺度，运算量非常大、非常费时，如湍流的直接数值模拟等。现在通常采用Navier-Stokes(N-S)方程的隐式方法求解往往需要大规模的代数方程组求解，虽然N-S方程的显式解法可避免大规模代数方程组求解，具有较好的并行性，但其数值稳定性较差，收敛速度较慢；而隐式解法数值稳定性好，收敛速度快，但是其并行可扩展性较差。
在实际模拟求解流体力学问题或相互作用中，广泛使用的方法是使用基于多核CPU和GPU的大规模计算集群对LBM应用进行加速。在现有条件下，使用GPU对LBM进行加速可以达到目前的最高运算性能，但GPU的能耗非常高；而使用多核CPU集群时，加速效果远不如使用GPU加速效果，而且系统功耗也没有显著降低。
针对具体流体力学求解设计专用的集成电路(Application SpecificIntegrated Circuit,ASIC)在性能和能耗表现上来说都是最合适的，但ASIC内部电路的固化特性使得只能解决事先设定的问题，同时固化的集成电路与新应用算法不匹配，每次开发都需要较长周期、成本高。
开发FPGA需要特定的硬件描述语言，如Verilog等，还需具备精湛的数字电路设计经验，精通硬件描述语言的设计开发流程，还需要花费大量时间进行模拟仿真，验证设计的正确性和性能。因此，现在的方法使用硬件描述语言进行FPGA开发费时费力，对开发者要求太高。
发明内容
本发明提供一种低功耗高运算性能模拟及精确求解计算流体动力学的框架方法，采用的计算架构是不断迭代的模式。整个计算架构分为CPU端和FPGA端两部分，通过OpenCL框架进行组织调度，CPU端负责数据预处理、任务调度及数据的传输存储，FPGA端负责数据计算，通过FPGA板卡的板载内存进行数据交互。
包括以下6个步骤：
步骤1：初始化网格数据和OpenCL环境；
步骤2：读取数据并计算密度及速度等宏观量；
步骤3：边界检测；
步骤4：格点碰撞计算；
步骤5：数据传播；
步骤6：根据迭代停止条件判断是否继续；
其中步骤1、步骤6两个步骤在CPU端执行，其余步骤位于FPGA端执行，步骤2、步骤3两个步骤完全并行，其余各步骤内部的并行性在具体实施方式中阐述。
本发明基于FPGA计算部件，采用格子Boltzmann方法对计算流体动力学的应用进行高效能计算方法，方法框架如图1所示，包括一些步骤：首先，CPU端执行OpenCL环境初始化并根据应用的参数进行网格数据初始化，再将网格数据传输到FPGA端的外部内存中。其次，FPGA端中，数据读取核(Read Kernel)从板载内存中读取粒子密度分布和边界信息数据，并计算各格子对应的宏观量(密度，速度等),将处理后的信息通过FPGA中的数据通道内部格点碰撞核(Inside Cell Collision Kernel)和数据写出核(Write U Kernel)，然后对各个格点的数据进行边界检测，将边界格点的数据送入边界格点碰撞核(Edge CellCollision Kernel)。数据写出核用于将计算好的速度(U)数据写入板载内存中，供CPU端进行读取。内部格点碰撞核被设计为自动运行内核(Auto Run Kernel)，即当设计的二进制文件加载到FPGA芯片中时就开始运行，不需CPU端进行调度，并且不与外部数据接口有任何联系，只通过FPGA内部数据通道(channel)与其他内核进行数据交换，该内核读取来自数据读取核的数据，并进行计算，然后将计算结果传到数据传播核(Stream Kernel)。边界格点碰撞核也采用自动运行设计，读取来自边界检测核的数据并进行计算，然后将计算结果传到数据传播核。数据传播核从两个格点碰撞核的数据通道中读取数据，进行数据传播，并将新的粒子密度分布存储到FPGA板载内存中，供下一次迭代使用。至此，一次迭代执行结束，随后CPU端读取速度数据，并判断迭代结束条件，若未结束，则继续调度内核进行计算，需要注意的是粒子密度分布数据采用乒乓缓存的方法进行存储，以减少片外内存的访问。
本发明专利目的主要是通过下述技术方案得以实现的：
业界通常采用采用高级综合工具(High-Level Synthesis,HLS)进行FPGA开发，将高级语言(C/C++)自动的编译为低级硬件描述语言，在本专利实施过程中主要使用OpenCL语言进行FPGA开发。
OpenCL是面向异构系统通用目的并行编程的开放式标准，也具有统一的编程环境，是适用于多核心处理器(multi-core CPU)、图形处理器(GPU)以及数字信号处理器(DSP)等并行处理器组成的异构系统。
在本发明中还做出如下工作，在步骤1完成初始化网格数据和OpenCL环境工作，主要在CPU端实施，可分解为3个子步骤来实现，分别如下：
子步骤11：FPGA-OpenCL环境初始化；
子步骤12：网格数据初始化；
子步骤13：程序其他数据的初始化。
现场可编程门阵列(Field-Programmable Gate Array,FPGA)是一种新兴的计算部件，它内部电路的可编程，理论上它可以直接改变硬件电路结构对算法进行理想加速，是一个非常有应用前景的计算部件。在本发明子步骤11中要完成环境初始化，需要分解为5个子步骤，分别如下：
子步骤111：根据平台名称获取OpenCL平台对象(cl_platform_id)；
子步骤112：根据子步骤111获得的平台对象，获取设备对象(cl_device_id)；
子步骤113：根据子步骤112获得的设备对象，创建上下文对象(context)和程序对象(cl_program)；
子步骤114：根据上述步骤获得上下文对象和设备对象，创建命令队列(cl_queue)；
子步骤115：根据核函数名称和程序对象，创建核函数对象(cl_kernel)。
在本发明子步骤12中要网格数据初始化，需要分解为3个子步骤，分别如下：
子步骤121：根据程序输入的参数，生成网格计算域；
子步骤122：对子步骤121生成的计算域中的每个网格，根据程序输入参数，对密度和速度进行初始化；
子步骤123：根据子步骤122的密度和速度数据，使用平衡态分布函数对每个格点的粒子密度分布进行初始化，将其中的一个关键方程构造为：
Figure BDA0002320615260000021
迭代执行结束后CPU端读取速度数据，并判断迭代结束条件，若未结束，则继续调度内核进行计算，直到完成步骤1所需的工作。在本专利中粒子密度分布数据采用乒乓缓存的方法进行存储，以减少片外内存的访问。
时间、空间和相互作用的全离散格子玻尔兹曼方法(LBM)是与并行计算相适应的一种较新的求解算法，是一种适合于介观计算流体动力学求解方法。采用离散演化动力系统的计算模型，这些模型构建离散的空间格点和离散的时间步，在微观或细观层次上，基于局部动态平衡原理，描述个体自适应的演化运动，从而模拟复杂的流体动力学系统或相互作用。
在本发明专利中，FPGA端负责对格子数据进行计算，包括碰撞和传播两个过程。步骤2、步骤3并将LBM模型的完整迭代方程和拆分为碰撞、传播两个过程后的表达式进行实施。
在本发明步骤2中要完成读取数据并计算密度及速度等宏观量，主要在FPGA端实施，可分解为4个子步骤，分别如下：
子步骤21：根据步骤1中初始化的网格计算域，创建缓存对象(cl_mem)；
子步骤22：将数据传输到子步骤21创建的缓存中，设置kernel参数，设置迭代条件及读取数据的步长；
子步骤23：启动FPGA中的Kernel，等待执行完毕；
子步骤24：子步骤23所述的过程执行完毕后，判断迭代停止条件，判断读取速度的步长。
在本发明子步骤24中要判断迭代停止条件，判断读取速度的步长工作，可分解为3个子步骤，分别如下：
子步骤241：若未到达迭代停止条件，则继续调度内核；
子步骤242：若到达读取速度的步长，则从FPGA内存中读取速度，并保存到文件中；
子步骤243：若达到停止条件，则释放资源，包括cl_kernel对象、cl_program对象、cl_queue对象和cl_mem对象。
在本发明步骤3中要完成边界检测工作，主要在FPGA端实施，可分解为4个子步骤，分别如下：
子步骤31：数据读取核(Read Kernel，RK)以SWI模式运行，读取每个网格的粒子密度分布以及网格的边界信息；
子步骤32：根据子步骤31读取的各网格的粒子密度分布，计算网格的密度和速度，前后数据的对应关系为：
Figure BDA0002320615260000032
和
Figure BDA0002320615260000033
将这三个数据和该网格的位置数据打包为一个结构体，记为Cell，并通过Channel传输到内部格点碰撞核(Inside Lattice Collision Kernel,ILCK)；
子步骤33：如果迭代次数达到读取速度的步长，则将格点速度信息传输至速度写出核(WU)；
子步骤34：边界检测，根据边界信息与各网格的位置进行匹配，如果该网格属于边界网格，则将该Cell与邻近的Cell组成打包为一个结构体，记为Cell-Edge，并传输到边界格点碰撞核(Edge Lattice Collision Kernel,ELCK)。
在本发明步骤4中完成格点碰撞计算工作，主要在FPGA端实施，可分解为3个子步骤，分别如下：
子步骤41：ILCK和ELCK以SWI模式运行，以SWI模式运行，并且属于自启动(AutoRun Kernel)，所谓自启动内核，当设计的二进制文件加载到FPGA芯片中时就开始运行，不需CPU端进行调度，并且不与外部数据接口有任何联系，只通过FPGA内部数据通道(Channel)与其他内核进行数据交换；
子步骤42：ILCK中依据如下的方程进行碰撞计算，
Figure BDA0002320615260000031
格点中每个维度的数据计算完毕后，不需等待其他数据一起打包传输，直接将结果传输给传播核(Stream Kernel,SK)；
子步骤43：ELCK中使用的方程与ILCK不同，即边界格点的碰撞方程，定义如下：
Figure BDA0002320615260000041
平衡态函数函数的计算公式为：
Figure BDA0002320615260000042
非平衡态分布函数的计算公式为：
Figure BDA0002320615260000043
当边界的某些宏观量未知时可以使用邻近的内部格点的数据代替。边界格点的每个维度的数据碰撞完毕后，直接将结果传输至SK。
本发明专利步骤4的3个子步骤是并行工作的。
在本发明步骤5完成数据传播工作，主要在FPGA端实施，可分解为4个子步骤：
子步骤51：SK以SWI模式运行，是FPGA端计算架构的最后一部分，也是LBM的最后一个流程，SK读取来自ILCK和ELCK的新的粒子密度分布；
子步骤52：对于从ILCK传输来的数据，需要进行检测，将边界部分的数据过滤掉，只处理内部格点的数据，计算格点中各个维度分布的新位置，并直接存储到乒乓缓存的对应位置，传播公式为；
Figure BDA0002320615260000044
子步骤53：从ELCK传输来的数据不需要计算新位置，直接使用老位置，将新数据存储到乒乓缓存的对应位置即可。
子步骤54：WUK以SWI的模式运行，读取网格的尺寸，设置循环条件，进入循环后，监测Channel，读取Cell-U数据，并将相应格点的速度写入FPGA设备板载内存中，供CPU端读取。
步骤5的4个子步骤是并行执行的。
在本发明步骤6根据迭代停止条件判断是否继续迭代运算，主要在CPU端实施。
本专利发明提供了一种低功耗高运算性能模拟及精确求解计算流体动力学的框架方法，就可以避免在模拟或精确求解计算流体动力学时采用Navier-Stokes(N-S)方程的“两难”困境，并且该发明计算框架也具有更广泛的适应性。即兼顾了CPU和GPU大规模计算集群对LBM应用时的高运算性能，同时能够有效地避免CPU和GPU计算集群的系统功耗，并且在FPGA端的功耗显著低于CPU和GPU计算集群的系统功耗。采用OpenCL语言开发FPGA计算部件，并完成初始化工作，构建BLM模拟或精确求解计算流体动力问题工作环境，降低了开发FPGA的难度，省时省力。使用本发明计算框架，通过LBM在OpenCL计算框架优化后的FPGA计算部件运行，适合应用于处理以下工程问题，如解决多孔介质内的流动与传质问题，模拟计算气-固和流-固耦合问题，同时也能高效模拟各种复杂的非线性宏观现象，如高铁、飞机等高速运动物体和空气的作用界面，并且在这些领域研究应用中，达到低功耗高效能运算目的。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定，在附图中：
图1为本发明所述的一种低功耗高运算性能模拟及精确求解计算流体动力学的框架方法示意图
图2为本发明所属的优化FGPA计算框架及和CPU的相互作用示意图
图3基于图2本专利CPU与FPGA之间的通信通过OpenCL环境初始化示意图
图4基于图1本专利FPGA端进行LBM模型的计算架构及流程示意图。
具体实施方式
下面结合附图对本发明实施例通过优化FPGA设计实现用BLM方法高效能计算流体动力学的框架进行详细描述。显然，所描述的实施例是本发明一部分实施例。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互结合。
首先，CPU端执行OpenCL环境初始化并根据应用的参数进行网格数据初始化，再将网格数据传输到FPGA端的外部内存中。
其次，FPGA端中，数据读取核(Read Kernel)从板载内存中读取粒子密度分布和边界信息数据，并计算各格子对应的宏观量(密度，速度等),将处理后的信息通过FPGA中的数据通道内部格点碰撞核(Inside Lattice Collision Kernel)和数据写出核(Write UKernel)，然后对各个格点的数据进行边界检测，将边界格点的数据送入边界格点碰撞核(Edge Lattice Collision Kernel)。数据写出核用于将计算好的速度(U)数据写入板载内存中，供CPU端进行读取。
内部格点碰撞核被设计为自动运行内核(Auto Run Kernel)，即当设计的二进制文件加载到FPGA芯片中时就开始运行，不需CPU端进行调度，并且不与外部数据接口有任何联系，只通过FPGA内部数据通道(channel)与其他内核进行数据交换，该内核读取来自数据读取核的数据，并进行计算，然后将计算结果传到数据传播核(Stream Kernel)。
边界格点碰撞核也采用自动运行设计，读取来自边界检测核的数据并进行计算，然后将计算结果传到数据传播核。数据传播核从两个格点碰撞核的数据通道中读取数据，进行数据传播，并将新的粒子密度分布存储到FPGA板载内存中，供下一次迭代使用。
至此，一次迭代执行结束，随后CPU端读取速度数据，并判断迭代结束条件，若未结束，则继续调度内核进行计算，如图2所描述。需要注意的是粒子密度分布数据采用乒乓缓存的方法进行存储，以减少片外内存的访问。
如图1所示，本发明中的基于FPGA的LBM计算方法架构详细介绍如下：
本发明的计算架构是不断迭代的模式，下面对一次完整的迭代及迭代之间的联系进行详细阐述。整个计算架构分为CPU端和FPGA端两部分，通过OpenCL框架进行组织调度，CPU端负责数据预处理、任务调度及数据的传输存储，FPGA端负责数据计算，通过FPGA板卡的板载内存进行数据交互。
CPU端主要包含两个任务，分别是它与FPGA之间的通信和计算前后的数据处理。首先，CPU与FPGA之间的通信通过OpenCL框架进行，CPU首先需要进行OpenCL环境初始化。
具体来说，如图3所描述，第一步，根据平台名称获取平台对象(cl_platform_id)，如果平台对象获取失败，即如果平台对象为null，则程序初始化失败，退出程序；第二步，根据平台对象获取设备对象(cl_device_id)的数组；第三步，根据设备对象数组，创建上下文对象(cl_context),创建上下文对象时，研究人员需要根据应用是否具有可扩展性和资源调用策略，考虑应用需要同时使用多少块FPGA板卡；第四步，根据第一块设备(假设所有设备属于同一个类型)和FPGA二进制文件名获取二进制文件的字符串(std::string)，根据这个二进制文件的字符串、上下文对象和设备数组获取程序对象(cl_program),编译程序对象，此时，根据使用的设备数量，创建核函数对象(cl_kernel)和命令队列对象(cl_queue)数组，数组中的元素对应于每块设备，并将数据进行相应分割；第五步，根据上述数组的长度，分别创建其中的元素对象，具体来说，根据上下文对象和设备对象创建命令队列对象，根据程序对象和核函数名称创建核函数对象。完成上述步骤，并且步骤中没有发生错误，则OpenCL环境初始化成功，可以进行下一步工作，若发生错误，则程序退出。
接下来进行数据预处理，即根据程序的输入参数，进行格子数据初始化，具体来说，格子数据包含微观量和宏观量，微观量即每个格子的粒子密度分布，宏观量包括密度和速度，宏观量与微观量之间具有映射关系，即所有宏观量都可以用微观量通过公式换算出来。为了具体说明上述过程，下面将详细介绍初始化的过程以及所涉及到的公式。首先需要根据程序的输入参数生成网格计算域；其次对于计算域中的每个网格进行数据的初始化，密度和速度这两个宏观量的初始化也是根据程序的输入参数确定的；格子的粒子密度分布使用平衡态分布函数进行初始化，函数定义为：
Figure BDA0002320615260000061
其中ρ表示密度，u表示速度，这两个宏观量在上一步已经初始化完成，可以使用；c表示格子传播速度，由用户自定义；wi表示权重，ei表示离散速度向量，这两个参数由具体的格子Boltzmann模型决定，常用的有D2Q9、D3Q15和D3Q19三种模型，其中D2Q9模型具有9个速度方向(i＝0,…,8)用于模拟二维流体问题，D3Q15模型具有15个速度方向(i＝0,…,14)，D3Q19模型具有19个速度方向(i＝0,…,18)，这两个模型用来模拟三维流体问题。权重wi和离散速度向量ei的具体取值分别由表1和表2给出。
表1不同格子Boltzmann模型对应的离散速度向量ei的取值
Figure BDA0002320615260000062
表2不同格子Boltzmann模型对应的权重wi的取值
Figure BDA0002320615260000063
至此，格子数据的初始化工作完成。接下来，根据初始化后的格子大小，创建粒子密度分布和边界数据的缓存对象(cl_mem),由于宏观量可以通过微观量通过计算公式得到，且计算的时间远小于访问内存的时间，所以只传输微观量，宏观量通过计算获取。需要注意的是，由于每次迭代都需要使用上一次迭代后的数据进行计算，且每次迭代存在一个全局同步的过程，故需要使用两个缓存对象组成乒乓缓存来存储粒子密度分布数据，由于使用了乒乓缓存，减少了1/3的内存访问操作。然后将初始化数据传输到FPGA卡板载内存中，根据程序输入的参数设置迭代停止条件和速度数据的读取步长(即每隔多少次迭代读取一次速度)，设置kernel参数。进入迭代部分，启动所有Kernel(由于未使用NDRange模型，使用clEnqueueTask(…)对Kernel进行调度)，并等待Kernel执行结束(具体使用clWaitForEvents(…)和clFinish(…)这两个函数)，Kernel执行完毕后，判断迭代停止条件和读取步长，切换乒乓缓存开始下一次迭代。没当迭代步数到达读取步长时，CPU端通过clEnqueueReadBuffer(…)函数从FPGA设备板载内存中获取速度数据，并存储到文件中，已被后续作图分析。当到达迭代停止条件后，释放相关资源，包括cl_kernel对象、cl_program对象、cl_queue对象和cl_mem对象，分别使用clReleaseKernel(…)、clReleaseProgram(…)、clReleaseCommandQueue(…)和clReleaseMemObject(…)函数对上述资源进行释放。
至此，CPU端程序架构及执行流程全部描述完毕。
FPGA端负责对格子数据进行计算，本专利FPGA端进行LBM模型的计算架构及流程示意图如图4所示，包括碰撞和传播两个过程。为了后续对FPGA端的计算架构及流程进行更清晰的讲解，此处对LBM模型的完整迭代方程和拆分为碰撞、传播两个过程后的表达式进行实施。
整个模型的迭代方程为：
Figure BDA0002320615260000071
其中，fi(x,t)表示位于x处的网格点在t时刻的粒子密度分布函数，
Figure BDA0002320615260000072
是由系统当前的宏观量构造出来的局部平衡分布函数，τ表示松弛时间，与粘度系数(流体的粘性系数、热传导系数、质量扩散系数)有关，Δt表示时间步，ei为离散速度向量，上文中已进行详细介绍。其次，将上述方程拆分为碰撞和传播两个过程，并进行介绍。碰撞过程：
Figure BDA0002320615260000073
Figure BDA0002320615260000074
传播过程：
Figure BDA0002320615260000075
可以将
Figure BDA0002320615260000076
看做用来临时保存结果的中间变量，本发明中由于采用了FPGA内部的Channel数据通道，并对模型的碰撞和传播过程进行了解耦，即拆分到两个kernel中执行，所以无需专门在全局内存开辟额外的缓存区域，减少了芯片与片外存储的I/O通信，提高了性能。
首先，FPGA端的数据读取核(Read Kernel,RK)以单工作项(Single Work Item,SWI)的模式运行，所谓单工作项，即仅使用一个计算单元，并采用流水线并行模式运行，这种并行模式相对于NDRange模式来说，能够以很小的资源占用，达到很好的并行效果。RK首先读取网格计算域的尺寸，并依据此信息设置循环条件，每次循环读取一个格子的数据，每个格子的数据根据模型的不同，拥有固定维数(9、15或19)的粒子密度分布。读取时，采用unroll参数对内存访问进行优化，所谓unroll，是对固定次数的循环进行并行化展开，比如9维粒子密度分布需要循环9次进行数据读取，而使用unroll后，则将这9次循环展开成9个并行的数据访问操作，需要注意的是，如果这9个数在内存中存在的位置是邻近的，编译器会将这9次访问转化为1次大带宽访问，这无疑极大的减少了内存访问的成本。读取某个格子的每个维度粒子密度分布的同时，对该格子的宏观量进行计算，宏观量与微观量的映射关系如下：密度，ρ(x,t)＝∑ifi(x,t)；速度，ρ(x,t)·u(x,t)＝∑i eifi(x,t)，其中ρ(x,t)为位于x的网格，在t时刻的密度，u(x,t)为位于x的网格，在t时刻的速度。每个格点的粒子密度分布读取完，密度和速度计算完，则将这三个数据以及网格位置信息打包为一个结构体，记为Cell，并将Cell通过Channel(FPGA内部的一种通信通道，用于Kernel之间的数据交互)传输到内部格点碰撞核(Inside Lattice Collision Kernel,ILCK)；接下来，因为计算域的边界和内部格点的迭代规则不一样，所以需要对Cell进行边界检测，边界检测具体来说，就是根据从板载内存中读取到的边界数据与Cell中的位置信息进行匹配，如果是边界，则将Cell数据以及该格点邻近的格点数据打包为一个结构体，记为Edge-Cell，并将其通过Channel传输到边界格点碰撞核(Edge Lattice Collision Kernel,ELCK)；如果迭代步数达到读取步长的要求，则需要将格点的速度数据和位置信息打包为一个结构体，记为Cell-U，并将其通过Channel传入速度写出核(Write U Kernel,WUK)。
WUK以SWI的模式运行，读取网格的尺寸，设置循环条件，进入循环后，监测Channel，读取Cell-U数据，并将相应格点的速度写入FPGA设备板载内存中，供CPU端读取，CPU端的步骤上文已详细阐述。
ILCK和ELCK都以SWI模式运行，并且属于自启动(Auto Run Kernel)，所谓自启动内核，当设计的二进制文件加载到FPGA芯片中时就开始运行，不需CPU端进行调度，并且不与外部数据接口有任何联系，只通过FPGA内部数据通道(Channel)与其他内核进行数据交换。ILCK、ELCK和WUK三者之间是完全并行的，互不影响，没有数据依赖关系。ILCK和ELCK读取来自RK的Cell和Cell-Edge数据，并且以流水线并行的方式进行计算，结合unroll将固定循环并行化，有着很大的吞吐量。ILCK中依据如下的方程进行碰撞计算，
Figure BDA0002320615260000081
格点中每个维度的数据计算完毕后，不需等待其他数据一起打包传输，直接将结果传输给传播核(Stream Kernel,SK)。需要注意的是，ELCK中使用的方程前文没有介绍，即边界格点的碰撞方程，定义如下：
Figure BDA0002320615260000082
Figure BDA0002320615260000083
其中，fi(boundary,t)表示位于边界boundary的格点，t时刻下的粒子密度分布，由平衡态分布函数和非平衡态分布函数两部分组成，分别为
Figure BDA0002320615260000084
和
Figure BDA0002320615260000085
上文中已经详细阐述了平衡态分布函数的计算方程，非平衡态分布函数的计算公式为：
Figure BDA0002320615260000086
当边界的某些宏观量未知时可以使用邻近的内部格点的数据代替。边界格点的每个维度的数据碰撞完毕后，直接将结果传输至SK。
SK以SWI模式运行，是FPGA端计算架构的最后一部分，也是LBM的最后一个流程。SK读取来自ILCK和ELCK的新的粒子密度分布，对于ILCK来的数据，需要进行检测，将边界部分的数据过滤掉，只处理内部格点的数据，计算格点中各个维度分布的新位置，并直接存储到乒乓缓存的对应位置。ELCK不需要计算新位置，直接使用老位置，将新数据存储到乒乓缓存的对应位置即可。
至此FPGA端的计算架构、计算流程及并行优化思想全部讲解完毕。FPGA端执行完毕后，CPU进行迭代结束条件判断，数据读取，内核调度，资源释放等操作，上文已详述。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (12)

1.一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，采用的计算架构是不断迭代的模式，整个计算架构分为CPU端和FPGA端两部分，通过OpenCL框架进行组织调度，CPU端负责数据预处理、任务调度及数据的传输存储，FPGA端负责数据计算，通过FPGA板卡的板载内存进行数据交互，包括以下6个步骤：
步骤1：初始化网格数据和OpenCL环境；
步骤2：读取数据并计算密度及速度等宏观量；
步骤3：边界检测；
步骤4：格点碰撞计算；
步骤5：数据传播；
步骤6：根据迭代停止条件判断是否继续；
其中步骤1、步骤6两个步骤在CPU端执行，其余步骤位于FPGA端执行，步骤2、步骤3两个步骤完全并行，其余各步骤内部的也是并行性的。
2.如权利要求1所述的低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本专利实施过程中主要使用OpenCL语言进行FPGA开发。
3.如权利要求1和权利要求2所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在步骤1完成初始化网格数据和OpenCL环境工作，主要在CPU端实施，可分解为3个子步骤来实现，分别如下：
子步骤11：FPGA-OpenCL环境初始化；
子步骤12：网格数据初始化；
子步骤13：程序其他数据的初始化。
4.如权利要求3所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明子步骤11中要完成FPGA-OpenCL环境初始化，需要分解为5个子步骤，分别如下：
子步骤111：根据平台名称获取OpenCL平台对象(cl_platform_id)；
子步骤112：根据子步骤111获得的平台对象，获取设备对象(cl_device_id)；
子步骤113：根据子步骤112获得的设备对象，创建上下文对象(context)和程序对象(cl_program)；
子步骤114：根据上述步骤获得上下文对象和设备对象，创建命令队列(cl_queue)；
子步骤115：根据核函数名称和程序对象，创建核函数对象(cl_kernel)。
5.如权利要求3所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明子步骤12中要网格数据初始化，需要分解为3个子步骤，分别如下：
子步骤121：根据程序输入的参数，生成网格计算域；
子步骤122：对子步骤121生成的计算域中的每个网格，根据程序输入参数，对密度和速度进行初始化；
子步骤123：根据子步骤122的密度和速度数据，使用平衡态分布函数对每个格点的粒子密度分布进行初始化，将其中的一个关键方程构造为：
Figure FDA0002320615250000011
迭代执行结束后CPU端读取速度数据，并判断迭代结束条件，若未结束，则继续调度内核进行计算，直到完成步骤1所需的工作，并且在本专利中粒子密度分布数据采用乒乓缓存的方法进行存储，以减少片外内存的访问。
6.如权利要求1所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明专利中，FPGA端负责对格子数据进行计算，包括碰撞和传播两个过程；步骤2、步骤3并将LBM模型的完整迭代方程和拆分为碰撞、传播两个过程后的表达式进行实施。
7.如权利要求1和权利要求6所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明步骤2中要完成读取数据并计算密度及速度等宏观量，主要在FPGA端实施，可分解为4个子步骤，分别如下：
子步骤21：根据步骤1中初始化的网格计算域，创建缓存对象(cl_mem)；
子步骤22：将数据传输到子步骤21创建的缓存中，设置kernel参数，设置迭代条件及读取数据的步长；
子步骤23：启动FPGA中的Kernel，等待执行完毕；
子步骤24：子步骤23所述的过程执行完毕后，判断迭代停止条件，判断读取速度的步长。
8.如权利要求7所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明子步骤24中要判断迭代停止条件，判断读取速度的步长工作，可分解为3个子步骤，分别如下：
子步骤241：若未到达迭代停止条件，则继续调度内核；
子步骤242：若到达读取速度的步长，则从FPGA内存中读取速度，并保存到文件中；
子步骤243：若达到停止条件，则释放资源，包括cl_kernel对象、cl_program对象、cl_queue对象和cl_mem对象。
9.如权利要求1和权利要求6所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明步骤3中要完成边界检测工作，主要在FPGA端实施，可分解为4个子步骤，分别如下：
子步骤31：数据读取核(Read Kernel，RK)以SWI模式运行，读取每个网格的粒子密度分布以及网格的边界信息；
子步骤32：根据子步骤31读取的各网格的粒子密度分布，计算网格的密度和速度，前后数据的对应关系为：
Figure FDA0002320615250000021
和
Figure FDA0002320615250000022
将这三个数据和该网格的位置数据打包为一个结构体，记为Cell，并通过Channel传输到内部格点碰撞核(Inside Lattice Collision Kernel,ILCK)；
子步骤33：如果迭代次数达到读取速度的步长，则将格点速度信息传输至速度写出核(WU)；
子步骤34：边界检测，根据边界信息与各网格的位置进行匹配，如果该网格属于边界网格，则将该Cell与邻近的Cell组成打包为一个结构体，记为Cell-Edge，并传输到边界格点碰撞核(Edge Lattice Collision Kernel,ELCK)。
10.如权利要求1和权利要求6所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明步骤4中完成格点碰撞计算工作，主要在FPGA端实施，可分解为3个子步骤，分别如下：
子步骤41：ILCK和ELCK以SWI模式运行，以SWI模式运行，并且属于自启动(Auto RunKernel)，所谓自启动内核，当设计的二进制文件加载到FPGA芯片中时就开始运行，不需CPU端进行调度，并且不与外部数据接口有任何联系，只通过FPGA内部数据通道(Channel)与其他内核进行数据交换；
子步骤42：ILCK中依据如下的方程进行碰撞计算，
Figure FDA0002320615250000031
格点中每个维度的数据计算完毕后，不需等待其他数据一起打包传输，直接将结果传输给传播核(Stream Kernel,SK)；
子步骤43：ELCK中使用的方程与ILCK不同，即边界格点的碰撞方程，定义如下：
Figure FDA0002320615250000032
平衡态函数函数的计算公式为：
Figure FDA0002320615250000033
非平衡态分布函数的计算公式为：
Figure FDA0002320615250000034
当边界的某些宏观量未知时可以使用邻近的内部格点的数据代替，边界格点的每个维度的数据碰撞完毕后，直接将结果传输至SK；
本发明专利步骤4的3个子步骤是并行工作的。
11.如权利要求1和权利要求6所述的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明步骤5完成数据传播工作，主要在FPGA端实施，可分解为4个子步骤：
子步骤51：SK以SWI模式运行，是FPGA端计算架构的最后一部分，也是LBM的最后一个流程，SK读取来自ILCK和ELCK的新的粒子密度分布；
子步骤52：对于从ILCK传输来的数据，需要进行检测，将边界部分的数据过滤掉，只处理内部格点的数据，计算格点中各个维度分布的新位置，并直接存储到乒乓缓存的对应位置，传播公式为；
Figure FDA0002320615250000035
子步骤53：从ELCK传输来的数据不需要计算新位置，直接使用老位置，将新数据存储到乒乓缓存的对应位置即可。
子步骤54：WUK以SWI的模式运行，读取网格的尺寸，设置循环条件，进入循环后，监测Channel，读取Cell-U数据，并将相应格点的速度写入FPGA设备板载内存中，供CPU端读取；
步骤5的4个子步骤是并行执行的。
12.如权利要求1的一种低功耗高运算性能求解计算流体动力学的框架方法，其特征在于，在本发明步骤6根据迭代停止条件判断是否继续迭代运算，主要在CPU端实施。
CN201911296165.7A 2019-12-16 2019-12-16 一种低功耗高运算性能求解计算流体动力学的框架方法 Active CN111105341B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911296165.7A CN111105341B (zh) 2019-12-16 2019-12-16 一种低功耗高运算性能求解计算流体动力学的框架方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911296165.7A CN111105341B (zh) 2019-12-16 2019-12-16 一种低功耗高运算性能求解计算流体动力学的框架方法

Publications (2)

Publication Number Publication Date
CN111105341A true CN111105341A (zh) 2020-05-05
CN111105341B CN111105341B (zh) 2022-04-19

Family

ID=70423459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911296165.7A Active CN111105341B (zh) 2019-12-16 2019-12-16 一种低功耗高运算性能求解计算流体动力学的框架方法

Country Status (1)

Country Link
CN (1) CN111105341B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112100099A (zh) * 2020-09-28 2020-12-18 湖南长城银河科技有限公司 一种面向多核向量处理器的格子玻尔兹曼优化方法
CN112906887A (zh) * 2021-02-20 2021-06-04 上海大学 稀疏gru神经网络加速的实现方法和装置

Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101727512A (zh) * 2008-10-17 2010-06-09 中国科学院过程工程研究所 一种基于变分多尺度方法的通用算法及并行计算系统
CN102053945A (zh) * 2009-11-09 2011-05-11 中国科学院过程工程研究所 一种面向多尺度离散模拟的并行计算系统
CN102681972A (zh) * 2012-04-28 2012-09-19 浪潮电子信息产业股份有限公司 一种利用GPU加速格子-Boltzmann的方法
CN102779207A (zh) * 2012-06-19 2012-11-14 北京航空航天大学 基于OpenCL的并行差分进化算法的翼型优化设计方法
CN102945295A (zh) * 2012-10-15 2013-02-27 浪潮(北京)电子信息产业有限公司 一种格子玻尔兹曼方法的并行加速方法及系统
CN103064819A (zh) * 2012-10-25 2013-04-24 浪潮电子信息产业股份有限公司 一种利用MIC快速实现格子Boltzmann并行加速的方法
US20130198426A1 (en) * 2011-12-22 2013-08-01 Airbus Operations S.L. Heterogeneous parallel systems for accelerating simulations based on discrete grid numerical methods
CN103324531A (zh) * 2013-06-09 2013-09-25 浪潮电子信息产业股份有限公司 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法
CN103345580A (zh) * 2013-07-02 2013-10-09 上海大学 基于格子Boltzmann方法的并行CFD方法
CN103366045A (zh) * 2013-06-20 2013-10-23 华北水利水电大学 基于格子Boltzmann的流体可视化仿真方法
CN103698554A (zh) * 2013-12-17 2014-04-02 华中科技大学 一种流场实时精确测量系统及方法
CN104112291A (zh) * 2008-03-21 2014-10-22 柯斯提克绘图公司 用于光线追踪渲染的并行相交测试及着色的架构
CN104142845A (zh) * 2014-07-21 2014-11-12 中国人民解放军信息工程大学 基于OpenCL-To-FPGA的CT图像重建反投影加速方法
CN105278346A (zh) * 2015-11-06 2016-01-27 北京航空航天大学 一种基于离散格子Boltzmann双分布模型的热流体仿真方法
CN107122243A (zh) * 2017-04-12 2017-09-01 杭州远算云计算有限公司 用于cfd仿真计算的异构集群系统及cfd计算方法
CN107980118A (zh) * 2015-06-10 2018-05-01 无比视视觉技术有限公司 使用多线程处理的多核处理器设备
CN108597012A (zh) * 2018-04-16 2018-09-28 北京工业大学 一种基于cuda的医学图像的三维重建方法
US20190332869A1 (en) * 2017-04-17 2019-10-31 Intel Corporation Person tracking and privacy and acceleration of data using autonomous machines

Patent Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104112291A (zh) * 2008-03-21 2014-10-22 柯斯提克绘图公司 用于光线追踪渲染的并行相交测试及着色的架构
CN101727512A (zh) * 2008-10-17 2010-06-09 中国科学院过程工程研究所 一种基于变分多尺度方法的通用算法及并行计算系统
CN102053945A (zh) * 2009-11-09 2011-05-11 中国科学院过程工程研究所 一种面向多尺度离散模拟的并行计算系统
US20130198426A1 (en) * 2011-12-22 2013-08-01 Airbus Operations S.L. Heterogeneous parallel systems for accelerating simulations based on discrete grid numerical methods
CN102681972A (zh) * 2012-04-28 2012-09-19 浪潮电子信息产业股份有限公司 一种利用GPU加速格子-Boltzmann的方法
CN102779207A (zh) * 2012-06-19 2012-11-14 北京航空航天大学 基于OpenCL的并行差分进化算法的翼型优化设计方法
CN102945295A (zh) * 2012-10-15 2013-02-27 浪潮(北京)电子信息产业有限公司 一种格子玻尔兹曼方法的并行加速方法及系统
CN103064819A (zh) * 2012-10-25 2013-04-24 浪潮电子信息产业股份有限公司 一种利用MIC快速实现格子Boltzmann并行加速的方法
CN103324531A (zh) * 2013-06-09 2013-09-25 浪潮电子信息产业股份有限公司 一种基于格子Boltzmann理论CPU/MIC协同计算的大涡模拟方法
CN103366045A (zh) * 2013-06-20 2013-10-23 华北水利水电大学 基于格子Boltzmann的流体可视化仿真方法
CN103345580A (zh) * 2013-07-02 2013-10-09 上海大学 基于格子Boltzmann方法的并行CFD方法
CN103698554A (zh) * 2013-12-17 2014-04-02 华中科技大学 一种流场实时精确测量系统及方法
CN104142845A (zh) * 2014-07-21 2014-11-12 中国人民解放军信息工程大学 基于OpenCL-To-FPGA的CT图像重建反投影加速方法
CN107980118A (zh) * 2015-06-10 2018-05-01 无比视视觉技术有限公司 使用多线程处理的多核处理器设备
CN105278346A (zh) * 2015-11-06 2016-01-27 北京航空航天大学 一种基于离散格子Boltzmann双分布模型的热流体仿真方法
CN107122243A (zh) * 2017-04-12 2017-09-01 杭州远算云计算有限公司 用于cfd仿真计算的异构集群系统及cfd计算方法
US20190332869A1 (en) * 2017-04-17 2019-10-31 Intel Corporation Person tracking and privacy and acceleration of data using autonomous machines
CN108597012A (zh) * 2018-04-16 2018-09-28 北京工业大学 一种基于cuda的医学图像的三维重建方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
JAKSIC, ZORAN 等: "A highly parameterizable framework for Conditional Restricted Boltzmann Machine based workloads accelerated with FPGAs and OpenCL", 《FUTURE GENERATION COMPUTER SYSTEMS》 *
姜典坤: "基于异构处理器的深度卷积神经网络加速系统设计与实现", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *
崔冠哲 等: "基于格子Boltzmann方法土体CT扫描切片细观渗流场的数值模拟", 《岩土力学》 *
张云 等: "多松弛时间格子Boltzmann方法在GPU上的实现", 《计算机与应用化学》 *
张峰: "面向集成异构平台的负载分析与优化关键技术研究", 《中国博士学位论文全文数据库 信息科技辑》 *
张纲 等: "格子Boltzmann方法多GPU并行性能的研究", 《计算机与应用化学》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112100099A (zh) * 2020-09-28 2020-12-18 湖南长城银河科技有限公司 一种面向多核向量处理器的格子玻尔兹曼优化方法
CN112906887A (zh) * 2021-02-20 2021-06-04 上海大学 稀疏gru神经网络加速的实现方法和装置

Also Published As

Publication number Publication date
CN111105341B (zh) 2022-04-19

Similar Documents

Publication Publication Date Title
Jouppi et al. In-datacenter performance analysis of a tensor processing unit
KR20200069346A (ko) 대규모 병렬 소프트웨어로 정의된 하드웨어 시스템에서의 정적 블록 스케줄링
CN112433819A (zh) 异构集群调度的模拟方法、装置、计算机设备及存储介质
CN111105341B (zh) 一种低功耗高运算性能求解计算流体动力学的框架方法
CN105930201A (zh) 一种可重构专用处理器核的功能模拟器
CN102207904B (zh) 用于对可重构处理器进行仿真的设备和方法
Meng et al. Preliminary experiences with the uintah framework on intel xeon phi and stampede
CN109213587B (zh) GPU平台下的多Stream并行DAG图任务映射策略
Bosse et al. Structural health and load monitoring with material-embedded sensor networks and self-organizing multi-agent systems
Bernaschi et al. A factored sparse approximate inverse preconditioned conjugate gradient solver on graphics processing units
Fell et al. The marenostrum experimental exascale platform (MEEP)
Christou et al. Earth system modelling on system-level heterogeneous architectures: EMAC (version 2.42) on the Dynamical Exascale Entry Platform (DEEP)
Du et al. Feature-aware task scheduling on CPU-FPGA heterogeneous platforms
Luo et al. Design of FPGA-based accelerator for convolutional neural network under heterogeneous computing framework with OpenCL
Iványi CUDA accelerated implementation of parallel dynamic relaxation
Ducroux et al. Fast and accurate power annotated simulation: Application to a many-core architecture
CN115309502A (zh) 一种容器调度方法及装置
Amador et al. A CUDA-based implementation of stable fluids in 3D with internal and moving boundaries
CN114218874A (zh) 一种翼伞流固耦合并行计算方法
Liu et al. Parallel implementation and optimization of regional ocean modeling system (ROMS) based on sunway SW26010 many-core processor
US9098917B2 (en) Method and system for accelerating collision resolution on a reconfigurable processor
Cui et al. Design and Implementation of OpenCL-Based FPGA Accelerator for YOLOv2
Gan et al. Million-core-scalable simulation of the elastic migration algorithm on Sunway TaihuLight supercomputer
Chung et al. A compilation flow for the generation of CNN inference accelerators on FPGAs
Banerjee et al. Multi-stage parallel processing of design element access tasks in FPGA-based logic emulation systems

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Ding Xuehai

Inventor after: Yan Weian

Inventor after: Tong Weiqin

Inventor after: Zhi Xiaoli

Inventor after: Cheng Jinfeng

Inventor before: Yan Weian

Inventor before: Ding Xuehai

Inventor before: Tong Weiqin

Inventor before: Zhi Xiaoli

Inventor before: Cheng Jinfeng

GR01 Patent grant
GR01 Patent grant