开源风雷CFD软件多物理场耦合接口开发路线分享！！！-编程知识

本文将基于开发过程中积累的经验，介绍风雷如何基于preCICE开发适配器。

preCICE是一个开源的多物理场数值模拟耦合库，可以用于多个求解器联合求解一个复杂的多场问题，支持在大规模并行系统上应用，具有良好的并行效率。并且可以对各种支持二次开发的求解器进行快速适配，在开发方面具有高度的灵活性。

风雷适配器可以在风雷求解CFD问题的基础上继续拓展功能，最终目的是实现风雷和其他学科求解器的多物理场耦合计算。基于此目的开发了一个简单的风雷适配器以验证耦合可行性。

由于适配器的开发仅作为验证作用，功能并不完善，文章会在相应位置注明后续可拓展的方向。

下面会对风雷整体的开发思路进行介绍：

求解器功能分析

为了开发风雷适配器，在掌握preCICE相关API功能以外，还需要对风雷求解器的使用和代码架构进行学习。根据开发经验，主要需要学习以下信息：

🔹风雷的使用方法

🔹代码整体结构

🔹求解主流程

🔹网格数据表示方法

🔹物理场数据表示方法

风雷的使用方法

在风雷的算例中，通常包含三个文件夹，其分别为参数文件bin、网格文件grid和结果文件results。下图为风雷的算例文件结构：

图1 风雷算例文件结构

在完成具体配置信息的设置后，直接启动风雷即可。可以参照风雷提供的使用手册来学习。

代码整体结构

代码的整体结构如下图所示：

图2 风雷源代码结构

由图可知，代码主要分为Main、CFD、API、3rdparty、Mesh。其中有两个部分是与开发适配器密切相关的，一个是Main，另一个是API。后面会对这两部分进一步介绍。

求解主流程

求解主流程在Main文件夹的代码中，具体的流程如下图：

图3 风雷CFD求解主流程图

主流程主要分为三部分：

各种数据的初始化，包括各种物理场和边界条件的初始化。

完成初始化后，开始进行流场的迭代计算。

收集并输出流场数据和残差数据。

网格数据表示方法

由于风雷软件支持大规模并行计算，所以内部支持多块网格分区计算，风雷网格相关结构主要分为3个层级：

Region层级：
这一层级代表了一个进程，Region能够控制一个进程上的全部分区网格的计算，也就是说一个Region上可以包含多个网格块。

Zone层级：
Zone层级代表了一个Region上的一块分区网格，除此之外它还包含了求解器网格流程的求解器。

Grid层级：
Grid层级就是Zone包含的网格，在风雷中Grid是存放网格信息的数据结构。

物理场数据表示方法

在Grid内部存放了流场数据结构的成员变量：Data_Field *gField，在不可压缩求解器和可压缩求解器中，物理场的数据都可以通过统一的API接口调用。

目前耦合适配器是基于不可压缩求解器开发的，可压缩和不可压缩求解器获取流场信息时调用API相同，但传递的参数不同，如果对可压缩求解器进行拓展需要注意此问题。

开发思路

在完成风雷运行和代码架构的相关学习后，可以制定一个初步的方案来进行代码接口的设计，可以结合上文介绍内容来设计相应的实现方案。下图描绘了总体开发思路，以及开发思路的详细介绍：

图4 风雷适配器主要功能开发思路

1、耦合配置信息

基于风雷的算例文件结构，可以分析出风雷求解器配置文件的读取在bin文件夹，那么就可以在风雷的参数配置文件中直接添加耦合所需的配置信息。

这样做的目的是可以复用求解器内部已有的配置文件解析功能，不用自己再额外开发一套新的读取配置文件的代码或引入其他库来解析配置文件。

注：目前已经完成耦合配置信息的数据结构，以及使用的接口，暂未开发配置文件的读写功能，可以基于以上思路对配置信息进行拓展。

2、求解主流程

适配器代码主要分为三部分插入到风雷主流程：

第一个是对各种信息的初始化（包括网格、求解器并行信息等），此部分代码在求解器初始化完毕后、开始迭代计算前插入；

第二就是具体的耦合计算（耦合数据的读取和写入，以及推进preCICE进行耦合计算等），此部分代码需要添加在求解循环内部；

第三就是适配器资源的释放，这部分可以添加在风雷释放资源的函数代码处。

3、网格信息

基于风雷求解器的并行求解，需要对网格进行分区处理，由此可知在程序中需要分别处理每个Region上所有Zone的Grid信息。风雷求解器的数据是基于面心存取的，所以需要提供面心（faceCenter）的网格信息；

为了支持流固耦合作用，流固耦合会导致流场网格变形，需要更新节点（faceNode）坐标，所以需要提供网格节点信息。因为耦合功能需要向preCICE注册耦合网格信息，需要在程序内部知道耦合网格的信息，在此可以使用网格的边界名称来识别。

因为preCICE耦合原理是弱耦合，对于求解器而言，通常只需要边界处的网格信息。对于并行的情况，给出的处理办法是将每个进程内参与耦合的边界信息按边界名字区分，最后把同名边界信息合并到一起，提供给preCICE。

注：由于目前风雷求解器本身的限制——风雷暂不支持提供与流场计算相耦合的动网格功能，所以目前暂时无法支持节点坐标的更新，只能耦合流场数据。且暂时未提供结构网格信息提取的功能，后续可以对此功能进行拓展。

4、物理场信息

在完成网格信息提取后，风雷中的物理场信息的提取，就变得十分轻松。风雷中通过GetDataPtr接口，传入数据名称即可获取到相应的物理场。

在物理场中获取耦合网格的物理场数据指针，只需在获取耦合网格信息时额外储存一个网格的索引号，通过每个网格的索引号就可以快速在物理场数据指针获得相应网格块的数据。

在完成网格块数据的提取后，就可以通过preCICE进行相应数据的传输并执行preCICE内置的映射插值算法。

注：preCICE支持显式与隐式计算，目前仅开发了显式计算相关功能，隐式计算还需要对流场数据进行保存和回溯操作，在计算复杂问题时，需要开发隐式耦合功能提高数值稳定性。

在完成这些工作后，就可以直接开始适配器的开发工作了。

开发时使用的风雷API

下表为开发风雷适配器时，主要使用的风雷求解器的接口。

表1 开发使用的函数

适配器代码使用介绍

下面的Simulation::SolveSimulation() 是风雷实现主求解循环的函数，在此插入开发的适配器函数，即可实现在风雷求解过程中的耦合计算功能。

在SolveSimulation添加的耦合适配器代码主要分为三部分：

1. 适配器初始化阶段（代码24-29行）

2. 读取/发送耦合数据，推进耦合整体进行（代码37-39行）

3. 释放耦合适配器资源，完成耦合计算（代码49行）

void Simulation::SolveSimulation()  {      WriteLogFile("Start iterating ...");      bool continueSolve = true;      int maxSimuStep = PHSPACE::GlobalDataBase::GetIntParaFromDB("maxSimuStep");      int isOverset = PHSPACE::GlobalDataBase::GetIntParaFromDB("codeOfOversetGrid");    7:int sysGridType = PHSPACE::GlobalDataBase::GetIntParaFromDB("sys_gridtype");          //! Wait for all processors here!      if (!isOverset || sysGridType == STRUCTGRID)      {          PrintToWindow("Wait for all processors before iteration ... \n");          PH_Barrier();          PrintToWindow("Wait is over,start iteration ... \n");      }        int outIterStep = PHSPACE::GlobalDataBase::GetIntParaFromDB("outnstep");      if (maxSimuStep == 0)      {          continueSolve = false;      }        //创建适配器并初始化      Adapter* adapter = new Adapter(this->region);      adapter->GetCouplingBC();      adapter->GetCouplingMeshVertices();      adapter->CreatPreciceInterface();      adapter->SetMesh();      double dt = adapter->PreciceInit();        while (continueSolve)      {          UpdateOuterIterationStep(outIterStep);          SolveOneOuterStep();            //发送数据并推进时间前进          adapter->WriteData(dt);          dt = adapter->Advance(dt);          adapter->ReadData(dt);              if (outIterStep >= maxSimuStep && maxSimuStep > 0)          {              continueSolve = false;          }      }        //释放适配器对象      adapter->Finaliaze();  }