Blazek版非结构网格CFD求解器案例分析03：
前面铺垫了这么多，今天终于进入正题了。这部分内容比较长，是框架的整体介绍，对理解整个代码非常重要。程序启动的入口一般都是main函数，之后通过main函数调用不同的函数实现既定的功能。因此我们先从main函数的框架讲起。
1.完整的main函数概览
main函数出去注释外有200多行，完整代码如下所示。刚开始接触大家可能有一点头晕，其实笔者第一次看到这些也是一样的。此处先展示出来只是给大家一个直观的印象，不必一行一行读下去。本文后面部分会将main函数主要实现的几个功能尽量用简洁的语言描述清楚。随着后续对各个类的成员变量和方法的详细解读，大家一定会有更清晰的认识。
/// @file main.cpp
///
/// 非结构网格2D求解器
///
// Features:
// ~~~~~~~~~
// # unstructured finite-volume scheme of median-dual type
// # triangular elements only
// # ideal gas model (other models possible)
// # laminar flow (viscosity computed by Sutherland`s law)
// # Roe's flux-difference splitting scheme, Venkatakrishnan's limiter
// # explicit multistage time-stepping scheme (Runge-Kutta type)
// # preconditioning for low Mach numbers
// # global or local time steps
// # central implicit residual smoothing
// # characteristic boundary conditions for external and internal flows
// # special initial solution for compressor and turbine blades
//
//*****************************************************************************
//
// (c) J. Blazek, CFD Consulting & Analysis, www.cfd-ca.de
// Created February 15, 2014
// Last modification: September 22, 2014
//
# include <cstdio>
# include <cstdlib>
# include <ctime>
# include <stdexcept>
# include "error.h"
# include "output.h"
# include "solver.h"
# include "userInput.h"
using namespace std ;
void Info ( ) ;
/// Main program of the flow solver.
///
/// @param argc number of command line arguments
/// @param argv[] list of arguments
/// @return EXIT_SUCCESS or an error code
///
int main ( int argc , char * argv [ ] )
{
Output output ;
Solver solver ;
string str ;
clock_t startTime , endTime ;
cout << endl
<< " *************************************************" << endl
<< " * *" << endl
<< " * 2-D FLOW ON UNSTRUCTURED TRIANGULAR GRIDS *" << endl
<< " * *" << endl
<< " * (c) Jiri Blazek, CFD Consulting & Analysis *" << endl
<< " * www.cfd-ca.de *" << endl
<< " * *" << endl
<< " * Version 1.0 from 09/22/2014 *" << endl
<< " * *" << endl
<< " *************************************************" << endl
<< endl ;
if ( argc == 2 )
{
try
{
UserInput :: Read ( argv [ 1 ] , solver , output ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
}
else
Info ( ) ;
// initialize some constants
solver . InitConstants ( ) ;
output . StoreFlowType ( solver . bndConds . flowType ) ;
// print input parameters for checking
UserInput :: Print ( solver , output ) ;
// read grid dimensions, coordinates, and triangle nodes
cout << " Reading grid data ..." << endl
<< endl ;
try
{
solver . geometry . ReadGrid ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
// generate edge list
cout << " Generating edge list ..." << endl
<< endl ;
try
{
solver . geometry . GenerateEdgelist ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
// print some statistics
cout << " No. of interior nodes: " << solver . geometry . nndInt << endl ;
cout << " No. of dummy nodes : " << solver . geometry . nNodes - solver . geometry . nndInt << endl ;
cout << " No. of grid cells : " << solver . geometry . nTria << endl ;
cout << " No. of interior edges: " << solver . geometry . nedInt << endl ;
cout << " Total number of edges: " << solver . geometry . nEdges << endl ;
cout << " No. of boundary faces: " << solver . geometry . nBfaces << endl ;
cout << " No. of boundary nodes: " << solver . geometry . nBnodes << endl
<< endl ;
// compute face vectors & cell volumes and check them;
// set coordinates of dummy nodes = those of boundary nodes;
// correct face vectors at symmetry boundaries;
// compute projections of control volumes
cout << " Computing metrics ..." << endl
<< endl ;
try
{
solver . geometry . ComputeMetrics ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
// allocate remaining memory in SpaceDiscr and TimeDiscr
cout << " Allocating remaining memory ..." << endl
<< endl ;
try
{
solver . bndConds . AllocateMemory ( solver . geometry ) ;
solver . fluidProps . AllocateMemory ( solver . geometry . nNodes ) ;
solver . spaceDiscr . AllocateMemory ( solver . fluidProps . equsType , solver . geometry . nNodes ) ;
solver . timeDiscr . AllocateMemory ( solver . geometry . nNodes ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
// read / initialize flow field
if ( solver . restUse )
{
cout << " Reading initial solution ..." << endl
<< endl ;
try
{
solver . ReadSolution ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
}
else
{
cout << " Guessing initial solution ..." << endl
<< endl ;
solver . InitSolution ( ) ;
}
solver . fluidProps . DependentVarsAll ( solver . geometry . nNodes ) ;
if ( ! solver . restUse )
solver . bndConds . BoundaryConditions ( solver . geometry , solver . fluidProps , solver . precond ) ;
// compute limiter reference values
solver . spaceDiscr . LimiterRefvals ( solver . geometry , solver . fluidProps , solver . bndConds ) ;
// open file for convergence history
try
{
output . OpenConvergence ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
// ***************************************************************************
// iterate until steady state solution or max. number of iterations is reached
// ***************************************************************************
// 准备必要的输出信息
if ( solver . bndConds . flowType == FlowType :: External )
{
str . clear ( ) ;
str . insert ( 0 , 83 , '-' ) ;
cout << str << endl
<< " step resid resmax i-res cl cd cm"
<< endl
<< str << endl ;
}
else
{
str . clear ( ) ;
str . insert ( 0 , 69 , '-' ) ;
cout << str << endl
<< " step resid resmax i-res mass flow mass ratio"
<< endl
<< str << endl ;
}
if ( ! solver . restUse )
solver . iter = 0 ;
startTime = clock ( ) ;
do
{
solver . iter ++ ;
solver . timeDiscr . Solve ( solver . geometry , solver . fluidProps , solver . bndConds ,
solver . spaceDiscr , solver . precond ) ;
solver . Convergence ( output ) ;
if ( solver . iter % solver . outStep == 0 )
{
try
{
cout << endl
<< " Writing plot files ..." << endl ;
output . Surfaces ( solver . geometry , solver . fluidProps , solver . bndConds ,
solver . spaceDiscr , solver . iter ) ;
output . Flowfield ( solver . geometry , solver . fluidProps , solver . bndConds ,
solver . iter ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
}
} while ( ! solver . Converged ( ) ) ;
endTime = clock ( ) ;
// ***************************************************************************
str . clear ( ) ;
str . insert ( 0 , 85 , '-' ) ;
cout << endl
<< str << endl
<< endl ;
// close file with convergence history
output . CloseConvergence ( ) ;
// output the results
if ( solver . iter % solver . outStep != 0 )
{
try
{
cout << " Writing plot files ..." << endl
<< endl ;
output . Surfaces ( solver . geometry , solver . fluidProps , solver . bndConds ,
solver . spaceDiscr , solver . iter ) ;
output . Flowfield ( solver . geometry , solver . fluidProps , solver . bndConds ,
solver . iter ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
}
// store solution for restart
cout << " Writing solution file ..." << endl
<< endl ;
try
{
solver . WriteSolution ( ) ;
}
catch ( exception & e )
{
Error :: Message ( e . what ( ) ) ;
}
cout << " Finished in " << fixed << setprecision ( 1 )
<< ( endTime - startTime ) / ( ( double ) CLOCKS_PER_SEC )
<< " seconds." << endl
<< endl ;
return ( EXIT_SUCCESS ) ;
}
//*****************************************************************************
/// Displays program usage.
///
void Info ( )
{
printf ( "Usage:\n" ) ;
printf ( "unstruct2d <input_file>\n\n" ) ;
exit ( EXIT_FAILURE ) ;
}
2.main函数头文件部分
下面的代码为头文件部分。C++支持分离式编译机制，这样的好处是团队之间的协作非常方便，可以由不同的团队完成不同部分的编写之后再讲程序文件集合在一起进行编译；而且在部分程序发生变动时通过Make文件可以只针对变化影响部分进行编译，节省了编译时间。此处的头文件即是为了支持这一机制设计的。其中<>部分是系统内置库部分，“”部分是自定义的代码文件对应的头文件。下图中的头文件会分别实现一些标准输入输出、动态内存分配、数据类型转换、求解等，具体后面用到时会详细解释。
//系统头文件部分
# include <cstdio>
# include <cstdlib>
# include <ctime>
# include <stdexcept>
//自定义头文件部分
# include "error.h"
# include "output.h"
# include "solver.h"
# include "userInput.h"
3.辅助信息输出部分
第一行表示命名空间std内的关键词后续都不需要加作用域运算符（::）了（谷歌编程规范中对这种方式不推荐），第二行进行了函数Info的声明，顾名思义Info是information的缩写，主要用于输出一些必要的信息。这种函数声明前置的优点是结构比较清晰。在main函数的末尾为Info函数的定义部分，可以看到其实就是printf函数的包装。
using namespace std;
void Info(); //Info函数声明，后面main程序会用到，所以必须提前声明
//末尾部分Info函数的定义部分
void Info()
{
printf("Usage:\n");
printf("unstruct2d <input_file>\n\n");
exit(EXIT_FAILURE);
}
4.main函数体内：实例化一些关键类
接下来进入main函数的作用域部分，首先前四行分别创建了几个类的实例，Output和Solver是自定义的类，分别用于输出控制和求解；string和clock_t则是常见的字符串类型和typedef过得long类型，分别用于main函数执行过程的字符串操作以及程序执行时间记录。
int main(int argc, char *argv[])
{
Output output; //output用于输出
Solver solver; //网格信息存储、数据存储、求解器算法等
string str;
clock_t startTime, endTime; //计算程序运行时间
接下来的部分用于输出一些程序信息，如下图1所示，可以看到其输出部分与程序代码部分是一一对应的。
cout << endl
<< " *************************************************" << endl
<< " * *" << endl
<< " * 2-D FLOW ON UNSTRUCTURED TRIANGULAR GRIDS *" << endl
<< " * *" << endl
<< " * (c) Jiri Blazek, CFD Consulting & Analysis *" << endl
<< " * www.cfd-ca.de *" << endl
<< " * *" << endl
<< " * Version 1.0 from 09/22/2014 *" << endl
<< " * *" << endl
<< " *************************************************" << endl
<< endl;
图1.程序启动后输出信息（part 1）
5.main函数体内：读取字典文件channel_input
接下来执行的部分为读取控制文件部分，也可以称之为是字典文件。在前面头文件中包含了"userInput.h"，此处的UserInput就是自定义的用于处理文件读取的类。该类内有一个Read成员函数，用于读取**_input文件，此处指定的是channel_input文件主要信息如下图2所示。可以看到channel_input文件其实就是我们在商业软件中设置的各类求解器信息。
//读取字典文件channel_input
if (argc == 2)
{
try
{
UserInput::Read(argv[1], solver, output);
}
catch (exception &e)
{
Error::Message(e.what());
}
}
else
Info();
图2.channel_input文件
此外值得注意的是Blazek在代码中引入了异常处理部分，这样在程序执行过程中可以针对出现的问题进行通信并做出相应的处理。也就是说有时候有一些异常信息不必严重到必须停止程序运行的程度时可以采用这样的处理方式。这对CFD程序而言很有必要，因为有时候我们在提交程序到服务器运行后并不能一直待在电脑旁监测，如果程序没有设计合理的异常处理功能那么很有可能晚上提交计算第二天来看时发现程序早就因为种种原因停止运行了，这对按分钟计算的宝贵服务器资源是巨大的浪费。
图3.程序执行过程所需的各类文件
6.main函数体内：初始化计算部分
此部分是初始化部分，进行一些常用变量的初始化。粗略来讲就是前面读取的channel_input文件通过Read函数处理时已经将对象solver的一些成员变量初始化了，可以理解为是“基础常量”，如密度、温度等，但是有一些“衍生常量”还需要通过这些“基础常量”进行计算，如马赫数等。flowType枚举类型则定义了内流外流等，这些也是通过channel_input文件定义的，根据流动类型的不同这些会影响到后续的处理方法，此处不再赘述，后面讲到时会详细分析。
// 初始化计算部分，基于channel_input中给定的物理常量进行一些其他物理常量的计算
solver.InitConstants();
//存储流动类型信息
output.StoreFlowType(solver.bndConds.flowType);
读取完成之后会通过UserInput内的Print函数输出一些读取信息（图4），用于程序执行过程中工程师进行检查。
UserInput::Print(solver, output);
图4.程序启动后输出信息（part 2）
7.main函数体内：网格处理部分
此部分是网格处理部分，原始的网格文件为channel.ugr文件，该文件规定的信息如图5所示。可以看到ugr文件部分比较简单，主要是给定了点坐标、线组成信息（哪两个点连成线）、面组成信息（哪些线组合成面）、边界类型信息。显然这些信息是不足以支撑CFD计算的，比如面之间的相邻关系、线的中心信息等等，这些都是通过这部分信息完成的。处理完成之后会输出必要的网格信息，如图6所示。
图5. urg文件
// read grid dimensions, coordinates, and triangle nodes
cout << " Reading grid data ..." << endl
<< endl;
try
{
solver.geometry.ReadGrid();
}
catch (exception &e)
{
Error::Message(e.what());
}
// generate edge list
cout << " Generating edge list ..." << endl
<< endl;
try
{
solver.geometry.GenerateEdgelist();
}
catch (exception &e)
{
Error::Message(e.what());
}
// print some statistics
cout << " No. of interior nodes: " << solver.geometry.nndInt << endl;
cout << " No. of dummy nodes : " << solver.geometry.nNodes - solver.geometry.nndInt << endl;
cout << " No. of grid cells : " << solver.geometry.nTria << endl;
cout << " No. of interior edges: " << solver.geometry.nedInt << endl;
cout << " Total number of edges: " << solver.geometry.nEdges << endl;
cout << " No. of boundary faces: " << solver.geometry.nBfaces << endl;
cout << " No. of boundary nodes: " << solver.geometry.nBnodes << endl
<< endl;
图6. 程序启动后输出信息（part 3）
8.main函数体内：分配求解变量存储所需的动态内存
此处部分是在堆上用new的方法分配内存，这些内存用于存储求解过程中的因变量信息等，比如我们求解得到的压力、速度、温度等的分布就是通过这些变量存储。不过此处笔者发现Blazek的代码里用了大量的new方式，这样容易产生内存泄漏的问题。而在OpenFOAM中则通过tmp、autoPtr等方式实现了内存的自动管理，类似于C++ STL中提供的智能指针（shared_ptr和unique_ptr）。
cout << " Allocating remaining memory ..." << endl
<< endl;
try
{
solver.bndConds.AllocateMemory(solver.geometry);
solver.fluidProps.AllocateMemory(solver.geometry.nNodes);
solver.spaceDiscr.AllocateMemory(solver.fluidProps.equsType, solver.geometry.nNodes);
solver.timeDiscr.AllocateMemory(solver.geometry.nNodes);
}
catch (exception &e)
{
Error::Message(e.what());
}
9.main函数体内：流场初始化
该部分用于初始化流场，其中分为是restart方式还是直接读取的方式。相对较为直观，后面具体分析内容时再详细分析。
// read / initialize flow field
if (solver.restUse)
{
cout << " Reading initial solution ..." << endl
<< endl;
try
{
solver.ReadSolution();
}
catch (exception &e)
{
Error::Message(e.what());
}
}
else
{
cout << " Guessing initial solution ..." << endl
<< endl;
solver.InitSolution();
}
solver.fluidProps.DependentVarsAll(solver.geometry.nNodes);
if (!solver.restUse)
solver.bndConds.BoundaryConditions(solver.geometry, solver.fluidProps, solver.precond);
10.main函数体内：限制器参考值设置
此部分用于设定限制器的参考值，限制器的主要作用是针对可压缩流过程中产生间断的区域采用高分辨率捕捉格式，其原理的简要介绍可以参考本系列文章第二部分，后面讲解到具体代码时也会详细分析。
solver.spaceDiscr.LimiterRefvals(solver.geometry, solver.fluidProps, solver.bndConds);
11.main函数体内：求解部分
首先又是输出一些信息，输出的信息主要是求解过程中的检测量，对应执行exe文件后的提示界面如图7所示。
// 准备必要的输出信息
if (solver.bndConds.flowType == FlowType::External)
{
str.clear();
str.insert(0, 83, '-');
cout << str << endl
<< " step resid resmax i-res cl cd cm"
<< endl
<< str << endl;
}
else
{
str.clear();
str.insert(0, 69, '-');
cout << str << endl
<< " step resid resmax i-res mass flow mass ratio"
<< endl
<< str << endl;
}
此处进入了真正的求解部分，其核心就是solver.timeDiscre.Solve部分。通过该函数的实参可以比较直观的看到，前面定义和计算的一系列变量都被传递给Solve函数进行运算。这部分是求解的核心，后面需要仔细讲解。
if (!solver.restUse)
solver.iter = 0;
startTime = clock();
do
{
solver.iter++;
solver.timeDiscr.Solve(solver.geometry, solver.fluidProps, solver.bndConds,
solver.spaceDiscr, solver.precond);
solver.Convergence(output);
if (solver.iter % solver.outStep == 0)
{
try
{
cout << endl
<< " Writing plot files ..." << endl;
output.Surfaces(solver.geometry, solver.fluidProps, solver.bndConds,
solver.spaceDiscr, solver.iter);
output.Flowfield(solver.geometry, solver.fluidProps, solver.bndConds,
solver.iter);
}
catch (exception &e)
{
Error::Message(e.what());
}
}
} while (!solver.Converged());
endTime = clock();
12.main函数体内：收尾部分，计算结果处理及存储
此部分是程序的收尾部分，用于将计算信息从内存写到硬盘，后面再详细分析。
// ***************************************************************************
str.clear();
str.insert(0, 85, '-');
cout << endl
<< str << endl
<< endl;
// close file with convergence history
output.CloseConvergence();
// output the results
if (solver.iter % solver.outStep != 0)
{
try
{
cout << " Writing plot files ..." << endl
<< endl;
output.Surfaces(solver.geometry, solver.fluidProps, solver.bndConds,
solver.spaceDiscr, solver.iter);
output.Flowfield(solver.geometry, solver.fluidProps, solver.bndConds,
solver.iter);
}
catch (exception &e)
{
Error::Message(e.what());
}
}
// store solution for restart
cout << " Writing solution file ..." << endl
<< endl;
try
{
solver.WriteSolution();
}
catch (exception &e)
{
Error::Message(e.what());
}
cout << " Finished in " << fixed << setprecision(1)
<< (endTime - startTime) / ((double)CLOCKS_PER_SEC)
<< " seconds." << endl
<< endl;
return (EXIT_SUCCESS);
}
文章来源：清锡cae
默认 最新
当前暂无评论，小编等你评论哦！点赞 2 评论 收藏 1
关注
