请教：第三代涡识别方法Liutex的代码在linux中运行的问题

李东岳

@五好青年我把这个帖子发到刘老师那面有个群里面去了，但目前没啥反馈。你这面有啥结果，告诉我。搞不定的话，我给你调试。

五好青年

@李东岳
目前我将cavity算例改成三维，三个方向网格相同，试算了后进行Liutex后处理。记录如下：
30^3网格单核，Liutex处理时间：4s
35^3网格单核，Liutex处理时间：14s
40^3网格单核，Liutex处理时间：60s
45^3网格单核，单核，Liutex处理时间：131s
50^3网格 Liutex处理时间：327s
55^3网格单核，浮点溢出
58^3网格单核，浮点溢出
60^3网格单核，浮点溢出
70^3网格单核，浮点溢出
80^3网格 4核并行后处理，报错

300万网格以前的算例，扩大机器内存，单核处理3天了，还没出结果

李东岳

@五好青年

#include <chrono>

auto start = std::chrono::steady_clock::now();                                                   
// Functions here
auto end = std::chrono::steady_clock::now();                                                     
auto diff = end - start;
Info<< "Calculation using " 
    << std::chrono::duration <double, std::milli> (diff).count()                                 
    << " ms" << endl;

你可以用上段代码测试哪个函数运行时间比较长

yhdthu

@yfleung jhd上的代码有bug，且效率低，我自己重写可用。
https://cloud.tsinghua.edu.cn/d/05b0149e37904e7ca61c/

五好青年

@yhdthu
谢谢海东博士，请问您这个是在OpenFOAM的哪个版本实现的呢？

五好青年

已测试，OpenFOAM-2006以上版本
感谢

五好青年

@yhdthu
您好
我尝试往OF官网版本上移植，提示了一些报错，主要有：error: ‘complex’ was not declared in this scope
经过很久的折腾，发现您可能参考tensor.C文件来修改的，并且我对比了OpenFOAM低版本的tensor.C来逐个修改报错，但是始终没能解决。
可否帮忙帮忙编一个低版本的，或者给一些修改的思路
感谢

dyfluid@dyfluid:~/postVersion$ wclean
dyfluid@dyfluid:~/postVersion$ wmake
Making dependency list for source file Liutex.C
SOURCE=Liutex.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -I/home/dyfluid/OpenFOAM/OpenFOAM-2.4.x/src/finiteVolume/lnInclude -I/home/dyfluid/OpenFOAM/OpenFOAM-2.4.x/src/meshTools/lnInclude -IlnInclude -I. -I/home/dyfluid/OpenFOAM/OpenFOAM-2.4.x/src/OpenFOAM/lnInclude -I/home/dyfluid/OpenFOAM/OpenFOAM-2.4.x/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/Liutex.o
Liutex.C: In function ‘Foam::vector getZdEVec(const tensor&)’:
Liutex.C:44:18: error: ‘complex’ was not declared in this scope
     const Vector<complex> EVals(eigenValues(T));
                  ^
Liutex.C:44:25: error: template argument 1 is invalid
     const Vector<complex> EVals(eigenValues(T));
                         ^
Liutex.C:44:47: error: cannot convert ‘Foam::vector {aka Foam::Vector<double>}’ to ‘const int’ in initialization
     const Vector<complex> EVals(eigenValues(T));
                                               ^
Liutex.C:46:25: error: type/value mismatch at argument 1 in template parameter list for ‘template<class Cmpt> class Foam::Vector’
     const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
                         ^
Liutex.C:46:25: note:   expected a type, got ‘complex’
Liutex.C:46:42: error: ‘Zero’ was not declared in this scope
     const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
                                          ^
Liutex.C:46:63: error: type/value mismatch at argument 1 in template parameter list for ‘template<class PrimitiveType> class Foam::pTraits’
     const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
                                                               ^
Liutex.C:46:63: note:   expected a type, got ‘complex’
Liutex.C:46:75: error: expression list treated as compound expression in initializer [-fpermissive]
     const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
                                                                           ^
Liutex.C:47:25: error: type/value mismatch at argument 1 in template parameter list for ‘template<class Cmpt> class Foam::Vector’
     const Vector<complex> standardBasis2(Zero, Zero, pTraits<complex>::one);
                         ^
Liutex.C:47:25: note:   expected a type, got ‘complex’
Liutex.C:47:69: error: type/value mismatch at argument 1 in template parameter list for ‘template<class PrimitiveType> class Foam::pTraits’
     const Vector<complex> standardBasis2(Zero, Zero, pTraits<complex>::one);
                                                                     ^
Liutex.C:47:69: note:   expected a type, got ‘complex’
Liutex.C:47:75: error: expression list treated as compound expression in initializer [-fpermissive]
     const Vector<complex> standardBasis2(Zero, Zero, pTraits<complex>::one);
                                                                           ^
Liutex.C:48:25: error: type/value mismatch at argument 1 in template parameter list for ‘template<class Cmpt> class Foam::Vector’
     const Vector<complex> EVec
                         ^
Liutex.C:48:25: note:   expected a type, got ‘complex’
Liutex.C:50:30: error: request for member ‘x’ in ‘EVals’, which is of non-class type ‘const int’
         eigenVector(T, EVals.x(), standardBasis1, standardBasis2)
                              ^
Liutex.C:53:23: error: request for member ‘x’ in ‘EVec’, which is of non-class type ‘const int’
     vector EVecR(EVec.x().real(), EVec.y().real(), EVec.z().real());
                       ^
Liutex.C:53:40: error: request for member ‘y’ in ‘EVec’, which is of non-class type ‘const int’
     vector EVecR(EVec.x().real(), EVec.y().real(), EVec.z().real());
                                        ^
Liutex.C:53:57: error: request for member ‘z’ in ‘EVec’, which is of non-class type ‘const int’
     vector EVecR(EVec.x().real(), EVec.y().real(), EVec.z().real());
                                                         ^
Liutex.C: In function ‘Foam::tensor rotation(const vector&)’:
Liutex.C:62:14: error: ‘Zero’ was not declared in this scope
     tensor q(Zero); 
              ^
Liutex.C: In function ‘Foam::vector liutex(const tensor&, const vector&)’:
Liutex.C:100:19: error: ‘Zero’ was not declared in this scope
     vector liutex(Zero);
                   ^
Liutex.C: In function ‘int main(int, char**)’:
Liutex.C:156:48: error: ‘Zero’ was not declared in this scope
             dimensionedVector("zero", dimless, Zero)
                                                ^
Liutex.C:158:39: error: ‘Foam::volVectorField {aka class Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>}’ has no member named ‘primitiveFieldRef’
         vectorField& iLiutex = Liutex.primitiveFieldRef();
                                       ^
Liutex.C:174:43: error: ‘Foam::volTensorField {aka class Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh>}’ has no member named ‘primitiveField’
         const tensorField& iGradU = gradU.primitiveField();
                                           ^
Liutex.dep:674: recipe for target 'Make/linux64GccDPOpt/Liutex.o' failed
make: *** [Make/linux64GccDPOpt/Liutex.o] Error 1
dyfluid@dyfluid:~/postVersion$

李东岳

@yhdthu 你写的看起来比期刊上的那个neat很多

五好青年

测试了一下这个程序，给大家反馈下

OpenFOAM-v2006以上可用
OpenFOAM-8调试好了，编译没报错，但是处理算例报错
foam-extend4.1上，调试好了，但是不识别浸入边界

目前只能在高版本上用了，如果谁编译了低版本，记得回来分享下啊

五好青年

东岳老师，能不能分享一个空的Ubuntu虚拟机出来呢。我记得以前是有的呢，但是找不到啦。
这样我在空的ubuntu系统上，尝试装一下不同版本，安装完我再反馈提供给您。分享给CFD中文网的Fans

尚善若水

@五好青年请问老师，后来of8的可用了吗？

五好青年

@尚善若水可以啦，直接用刘老师的fortran程序

尚善若水

@五好青年好的，我去查下，谢谢。

李东岳

按道理来说这个C++跑起来应该没问题啊为什么后来大佬们都转fortran了呢

李东岳

我把 http://www.jhydrodynamics.com/en/download-of-liutex-code/ 这个链接里面的代码，小改了一下，放到了OpenFOAM-9里面可以编译成功。

下载：Rotex.zip

我目前来看这个代码是每个网格做了循环，每个循环里面做了大量的操作。因此如果网格越多，确实会越来越慢。这跟其他的变量算起来完全不一样，其他的变量比如Q之类都是在体场层面做操作，这个Rotex从代码上来看，要从cell做循环来操作。

哪个大佬可以提供个算例，我来测试下。

尚善若水

@五好青年再次来打扰下 seeeeeeeeeeer 老师，没 google到 fortran 代码，可以给我贴一个吗？感谢。

李东岳

@bestucan 确实，这个大佬竟然等了半个月。。。是不是家里服务器太多了，有一个算例跑起来忘了，半个月之后才想起来，发现还在算Rotex！

之前不是说有个土豪车太多了，扔某个车库一个玛莎拉蒂，几年之后才想起来，车上已经一层灰。

256³的网格算了半个月都没有算出一个时间步文件！

bestucan

@李东岳佛系计算，不给服务器压力，不push

coolhhh

@yhdthu 在请教：第三代涡识别方法Liutex的代码在linux中运行的问题中说：

@yfleung jhd上的代码有bug，且效率低，我自己重写可用。
https://cloud.tsinghua.edu.cn/d/05b0149e37904e7ca61c/

你好，这个下载地址已经失效了，能否再提供一下吗，非常感谢

李东岳

@bestucan 想当年2021年这个时候，我还在计算LES。3台服务器，每个64核心每个跑一个算例。600万网格也是算了半个月，每个算了1000秒。现在俺手握800核，LES跑的很血腥，3000万网格跑的跟单相流一样。

bestucan

@李东岳摩尔定律之东岳流体分律

coolhhh

@李东岳

1、李老师，我用OpenFOAM-8编译了这个程序，编译过程提示有一些未使用变量，但能正常编译。

Making dependency list for source file Rotex.C
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3  -DNoRepository -ftemplate-depth-100 -I/home/install/OpenFOAM/OpenFOAM-8/src/finiteVolume/lnInclude -I/home/install/OpenFOAM/OpenFOAM-8/src/meshTools/lnInclude -IlnInclude -I. -I/home/install/OpenFOAM/OpenFOAM-8/src/OpenFOAM/lnInclude -I/home/install/OpenFOAM/OpenFOAM-8/src/OSspecific/POSIX/lnInclude   -fPIC -c Rotex.C -o Make/linux64GccDPInt32Opt/Rotex.o
Rotex.C: In function \u2018int main(int, char**)\u2019:
Rotex.C:153:20: warning: unused variable \u2018eig1c_r\u2019 [-Wunused-variable]
  153 |             scalar eig1c_r=-0.5*(aaaa+bbbb)-aa/3.0;
      |                    ^~~~~~~
Rotex.C:154:20: warning: unused variable \u2018eig1c_i\u2019 [-Wunused-variable]
  154 |             scalar eig1c_i=0.5*Foam::sqrt(scalar(3.0))*(aaaa-bbbb);
      |                    ^~~~~~~
Rotex.C:155:20: warning: unused variable \u2018eig2c_r\u2019 [-Wunused-variable]
  155 |             scalar eig2c_r=-0.5*(aaaa+bbbb)-aa/3.0;
      |                    ^~~~~~~
Rotex.C:156:20: warning: unused variable \u2018eig2c_i\u2019 [-Wunused-variable]
  156 |             scalar eig2c_i=-0.5*Foam::sqrt(scalar(3.0))*(aaaa-bbbb);
      |                    ^~~~~~~
g++ -std=c++11 -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3  -DNoRepository -ftemplate-depth-100 -I/home/install/OpenFOAM/OpenFOAM-8/src/finiteVolume/lnInclude -I/home/install/OpenFOAM/OpenFOAM-8/src/meshTools/lnInclude -IlnInclude -I. -I/home/install/OpenFOAM/OpenFOAM-8/src/OpenFOAM/lnInclude -I/home/install/OpenFOAM/OpenFOAM-8/src/OSspecific/POSIX/lnInclude   -fPIC -fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/Rotex.o -L/home/install/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/lib \
    -lfiniteVolume -lfvOptions -lmeshTools -lOpenFOAM -ldl  \
     -lm -o /home/user1/OpenFOAM/user1-8/platforms/linux64GccDPInt32Opt/bin/Rotex

2、提供4个加载CDRFG脉动风的空风洞算例结果下载链接：
（1）3个算例用snappyHexMesh画的六面体网格，数量分别为9.6w，150w，364w。
（2）1个算例用fluentmeshing画的poly网格，网格数量为141w
（3）RotexCase下载链接：
链接：https://pan.baidu.com/s/1HQZkWVg2-ZKUF8SmqwkwOg
提取码：fgzl

3、测试结果
（1）9.6w SHM网格算例，Rotex用单核、并行都能正常计算，单核40 s，并行20核计算用0.36 s。但用paraview看后处理结果，虽然可以看到Rotex数据，但由于Rotex输出的结果为矢量，画contour图的时候，选项找不到Rotex，但能选择Lambda2和Q画等值面图。
问题是采用Liutex如何画出三维形状的涡？

这是选择Rotex画出来的纵向截面图：

这是选择Q准则画出的涡形状：

（2）150w SHM网格算例。Rotex并行20核计算用时127 s。

（3）364w SHM网格算例。Rotex并行20核计算用时1855 s。可以看出运行Rotex随着网格增加，运行时间是非线性增加的。

（4）141w poly网格算例。Rotex并行20核计算，报错Floating point exception，但运行Lambda2和Q画等值面图，能正常运行。

Time = 9
    Reading U
[4] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[4] #1  Foam::sigFpe::sigHandler(int) at ??:?
[4] #2  ? in "/lib64/libc.so.6"
[4] #3  ? at ??:?
[4] #4  __libc_start_main in "/lib64/libc.so.6"
[4] #5  ? at ??:?
[hpws14:21072] *** Process received signal ***
[hpws14:21072] Signal: Floating point exception (8)
[hpws14:21072] Signal code:  (-6)
[hpws14:21072] Failing at address: 0x3e800005250
[hpws14:21072] [ 0] /lib64/libc.so.6(+0x36400)[0x7fab592fb400]
[hpws14:21072] [ 1] /lib64/libc.so.6(gsignal+0x37)[0x7fab592fb387]
[hpws14:21072] [ 2] /lib64/libc.so.6(+0x36400)[0x7fab592fb400]
[hpws14:21072] [ 3] Rotex[0x4136e0]
[hpws14:21072] [ 4] /lib64/libc.so.6(__libc_start_main+0xf5)[0x7fab592e7555]
[hpws14:21072] [ 5] Rotex[0x4146f0]
[hpws14:21072] *** End of error message ***
[9] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[9] #1  Foam::sigFpe::sigHandler(int) at ??:?
[9] #2  ? in "/lib64/libc.so.6"
[9] #3  ? at ??:?
[9] #4  __libc_start_main in "/lib64/libc.so.6"
[9] #5  ? at ??:?
[hpws14:21077] *** Process received signal ***
[hpws14:21077] Signal: Floating point exception (8)
[hpws14:21077] Signal code:  (-6)
[hpws14:21077] Failing at address: 0x3e800005255
[hpws14:21077] [ 0] /lib64/libc.so.6(+0x36400)[0x7f02ac8ea400]
[hpws14:21077] [ 1] /lib64/libc.so.6(gsignal+0x37)[0x7f02ac8ea387]
[hpws14:21077] [ 2] /lib64/libc.so.6(+0x36400)[0x7f02ac8ea400]
[hpws14:21077] [ 3] Rotex[0x4136e0]
[hpws14:21077] [ 4] /lib64/libc.so.6(__libc_start_main+0xf5)[0x7f02ac8d6555]
[hpws14:21077] [ 5] Rotex[0x4146f0]
[hpws14:21077] *** End of error message ***
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
[6] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[6] #1  Foam::sigFpe::sigHandler(int) at ??:?
[6] #2  ? in "/lib64/libc.so.6"
[6] #3  ? at ??:?
[6] #4  __libc_start_main--------------------------------------------------------------------------
mpirun noticed that process rank 9 with PID 0 on node hpws14 exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------

这是poly网格选择Q准则画出的涡形状：

李东岳

@coolhhh 在请教：第三代涡识别方法Liutex的代码在linux中运行的问题中说：

由于Rotex输出的结果为矢量，画contour图的时候，选项找不到Rotex，但能选择Lambda2和Q画等值面图

你用calculator做一个计算器，计算Liutex的mag，然后做等值面图可以么。

这个可能要看一下liutex的sci的定义，看是不是通过这个方法来进行的。

我下载你150万网格的看下Liutex哪一块最耗费时间。

coolhhh

@李东岳在请教：第三代涡识别方法Liutex的代码在linux中运行的问题中说：

你用calculator做一个计算器，计算Liutex的mag，然后做等值面图可以么。

这个可能要看一下liutex的sci的定义，看是不是通过这个方法来进行的。

李老师，我看了刘老师的论文（刘超群, Liutex-涡定义和第三代涡识别方法. 空气动力学学报, 2020. 038(003): 第413-431,478页.），里边提到以下几点：

Liutex算出的是一个矢量，而Q是标量。

2.文章中提到Liutex为基础的几个涡识别方法，其中之一就是可以用计算mag来画等值面，但这个时候又需要指定阈值：

3.用paraview计算9.6w 网格的算例的mag(Rotex)，阈值分别取不同结果如下。Liutex采用幅值等值面因此会存在跟用Q准则一样指定阈值的问题，但根据论文意思是此时画出的等值面和Liutex向量线是大致平行的。
（1）mag(Rotex)=0.1结果

（2）mag(Rotex)=0.5结果

（3）mag(Rotex)=2结果

李东岳

我测试了一下，单核计算，对于这个150万网格的算例，每个网格的操作要占用30ms来计算，150万网格就是45000秒。

因为我对Liutex这个计算过程不熟悉，不清楚这个是否能做简化。