OpenFOAM周期边界生成不成功，cyclic的问题

[李东岳]

转移原贴内容 http://www.cfd-china.com/topic/2188/25亿网格直接模拟

[TINGHAIK]

李老师，关于OF中网格check 的时候所输出来的各项指标的含义以及可能对求解产生什么影响，这方面有没有介绍或文献？

[random_ran]

不知道他们用的是 OF 吗？

我自己用 OF 4.1 瞎折腾过一次 2.5 亿的网格，生成周期边界的时候总是过不去。

$ checkMesh
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 4.1
Exec   : checkMesh
Date   : Jul 06 2018
Time   : 15:54:53
Host   : "cp0401"
PID    : 12872
Case   : /scratch/crazyuser/biggerMesh
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           257514000
    faces:            769512000
    internal faces:   766488000
    cells:            256000000
    faces per cell:   6
    boundary patches: 5
    point zones:      0
    face zones:       1
    cell zones:       1

Overall number of cells of each type:
    hexahedra:     256000000
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    FRONT               1000000  1002000  ok (non-closed singly connected)  
    INLET               256000   257257   ok (non-closed singly connected)  
    OUTLET              256000   257257   ok (non-closed singly connected)  
    CYLINDER            512000   514000   ok (non-closed singly connected)  
    BACK                1000000  1002000  ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-32 -32 0) (32 32 3.33332)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (-4.65606e-16 4.16736e-21 -4.97036e-16) OK.
    Max cell openness = 3.18122e-16 OK.
    Max aspect ratio = 29.4282 OK.
    Minimum face area = 1.3175e-06. Maximum face area = 0.0385213.  Face area magnitudes OK.
    Min volume = 1.71549e-08. Max volume = 0.000501578.  Total volume = 10720.6.  Cell volumes OK.
    Mesh non-orthogonality Max: 1.72876e-05 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.00266764 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

[youmengtian]

@random_ran 从check的结果来看，原始生成的网格没有什么问题。您说的周期性边界失败是指？可是尝试使用createPatch这个utility创建cyclic边界实现周期性边界。

[random_ran]

@youmengtian 在 25亿网格直接模拟 中说：

@random_ran 从check的结果来看，原始生成的网格没有什么问题。您说的周期性边界失败是指？可是尝试使用createPatch这个utility创建cyclic边界实现周期性边界。

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Reading createPatchDict

Adding new patch FRONT_CYC as patch 5 from 
{
    type            cyclic;
    neighbourPatch  BACK_CYC;
    matchTolerance  0.01;
}

Adding new patch BACK_CYC as patch 6 from 
{
    type            cyclic;
    neighbourPatch  FRONT_CYC;
    matchTolerance  0.01;
}


Moving faces from patch FRONT to patch 5
terminate called after throwing an instance of 'std::bad_array_new_length'
  what():  std::bad_array_new_length

cppreference 上说三种可能的原因，我觉得我的case是“too large”。

#include <iostream>
#include <new>
#include <climits>
 
int main()
{
    int negative = -1;
    int small = 1;
    int large = INT_MAX;
    try {
        new int[negative];           // negative size
        new int[small]{1,2,3};       // too many initializers
        new int[large][1000000];     // too large
    } catch(const std::bad_array_new_length &e) {
        std::cout << e.what() << '\n';
    }
}

我顺着 Moving faces from patch FRONT to patch 5 找到 createPatch.C 源代码:

====> createPatch.C:  Line 705~714
                forAll(pp, i)
                {
                    changePatchID
                    (
                        mesh,
                        pp.start() + i,
                        destPatchi,
                        meshMod
                    );
                }

似乎是死在了这个循环里:

changePatchID 的arguments：

(
    const polyMesh& mesh,
    const label faceID,
    const label patchID,
    polyTopoChange& meshMod
)

changePatchID 里嵌套的这个函数很可疑 meshMod.setAction， 特别是 polyModifyFace 这个函数

定义在

polyModifyFace.H 奇怪 polyModifyPoint.C 居然不存在.

看到这里，感觉无能为力了，想要找到类似这种 new int[large][1000000]; // too large
的声明看来是行不通了。

所以想请教大家的想法。

[youmengtian]

@random_ran polyModifyFace是一个类，而非函数。另外在polyModifyFace.H同级文件夹下有polyModifyPoint.H 。至于抛出的异常之前未见过，不过看样子应该是在分配内存时出现了异常，非常有可能是你说的问题。但是实际上可能和你演示的例子还是有差别，固定数组能够使用的栈本身大小有限很容易出现这种异常。但是实际中像createPatch之类需要对网格数据进行处理的往往是通过动态数组开辟内存空间的，其上线与计算机内存有直接关系。createPatch在运行中需要额外生成相关patch的数据信息，需要占用相应大小的内存，如果网格过大很有可能超出计算机内存上限抛出错误。

另，过于createPatch具体实现我也没有进行深入的分析，仅提供一个思路，正确与否请谨慎。

[youmengtian]

另，建议改小matchToTolerance。相对于你的网格尺寸，这个值明显偏大。

[random_ran]

嗨 youmengtian：

非常感谢阅读我的这个瞎折腾! 当时是想看看 OF 到底能跑多大的网格，所以才有了这些。

由于这是 cfd-china ， 默认大家都用中文，所以我有点担心会不会影响到别人。但是，这种边看代码，边写文档的东西，我目前只能在在 emacs 中高效地完成，所以就只能用英文来交流。 如果冒犯到你，请不要继续往下读。

I am using OF-dev version (don't know exactly the date I downloaded
it).

The main classes are found in the flowing path:

~/OpenFOAM-dev/src/dynamicMesh/polyTopoChange/polyTopoChange/modifyObject/

This for loop:

//Line 707
//~/OpenFOAM-dev/applications/utilities/mesh/manipulation/createPatch/createPatch.C

void changePatchID
(
  const polyMesh& mesh,
  const label faceID,
  const label patchID,
  polyTopoChange& meshMod
)

is the point that I feel something might wrong.

I don't know exactly what polyModifyFace is a function or a
class. I double checked the defination of polyModifyFace:

//Line 82
//~OpenFOAM-dev/src/dynamicMesh/polyTopoChange/polyTopoChange/modifyObject/polyModifyPoint.H
meshMod.setAction
  (
      polyModifyFace
      (
	  mesh.faces()[faceID],               // face
	  faceID,                             // face ID
	  mesh.faceOwner()[faceID],           // owner
	  -1,                                 // neighbour
	  false,                              // flip flux
	  patchID,                            // patch ID
	  false,                              // remove from zone
	  zoneID,                             // zone ID
	  zoneFlip                            // zone flip
      )
  );

I think polyModifyFace really looks like a function definition in
c++, but I might wrong.

Another thing that I feel strange is that the implementation of the
following classes are missing.

3.5k 12-14 23:52 polyModifyCell.H
7.5k 12-14 23:52 polyModifyFace.H
4.2k 12-14 23:52 polyModifyPoint.H

Is that becasue those classes are too easy?

I furthe study the class (~/OpenFOAM-dev/src/dynamicMesh/polyTopoChange/polyTopoChange/modifyObject/polyModifyFace.H):

//- Construct and return a clone
virtual autoPtr<topoAction> clone() const
{
    return autoPtr<topoAction>(new polyModifyFace(*this));
}

Is this clone() function mess up all the things? If so, what's the
limited of it? How big (size of the mesh) the OF can handle?

I agree with you. This mesh is converted from \*.msh file without any
issue. convertMesh seems have a much large memory requirement. I am
not sure exactly the order of it.

What I want to do, is to run some small test to investigate the
memory usage of convertMesh and creatPatch.

Finally, why make it (matchToTolerance) smaller? Do you have any
reference to this point? I might (can remember exactly) just use
some default number here.

[wwzhao]

@random_ran 我猜可能是 int 类型长度限制的原因。你看一下 echo $WM_LABEL_SIZE ，默认是 32，如果要用超大网格，需要在编译前将 WM_LABEL_SIZE 设为 64。

[youmengtian]

@random_ran 您说的太客气了，虽然咱们这儿叫CFD中文网但是也没说不让使用英文。你之前的阅读逻辑是没有什么错误的，但是我不知道你怎么读到polyModifyPoint.H这个头文件去，因为我觉得貌似跟这个头文件没有关系。但是我感觉你对错误的基本定位应该是准确的，应该就出现在polyModifyFace.H第190行开辟新的智能指针出错了

polyModifyFace.H      188
virtual autoPtr<topoAction> clone() const
{
    return autoPtr<topoAction>(new polyModifyFace(*this));
}

[random_ran]

因为我几次遇到过在网络上由于这种潜在规则的问题的不明了导致的一些不必要的麻烦，所以谨慎点应该的。看来之后的交流会更加便捷 ：）

@youmengtian

The reason why I checked polyModifyFace was that I wanted to study
the piece of code in:

//Line 82
//~OpenFOAM-dev/src/dynamicMesh/polyTopoChange/polyTopoChange/modifyObject/polyModifyPoint.H [wrong!]
///OpenFOAM-dev/applications/utilities/mesh/manipulation/createPatch/createPatch.C [correct]
meshMod.setAction
  (
      polyModifyFace
      (
	  mesh.faces()[faceID],               // face
	  faceID,                             // face ID
	  mesh.faceOwner()[faceID],           // owner
	  -1,                                 // neighbour
	  false,                              // flip flux
	  patchID,                            // patch ID
	  false,                              // remove from zone
	  zoneID,                             // zone ID
	  zoneFlip                            // zone flip
      )
  );

I don't know the reason why it crash, that's why I am asking.

return autoPtr<topoAction>(new polyModifyFace(*this));

Do you know how this function could trigger the following exception:

terminate called after throwing an instance of 'std::bad_array_new_length'
what():  std::bad_array_new_length

I don't know.

@wwzhao

thanks for helping. This test case was built Jul 06 2018,
so I am not sure the exactly envirmental varibles that I was playing with.

I've double checked the:

echo $WM_LABEL_SIZE
32

It looks like a right direction to further study this stuff. But, it
needs more time, since I need to rebuid OF. Which Version would you
like to suggest me to run? dev or v4.1?

Also, if the default

$WM_LABEL_SIZE</sub> is indeed 32

what is its theoretical upper bond, i.e. running pisoFoam (motobike tutorial?).

Can I change to any number (2^{5,6,7,8})?

[wwzhao]

@random_ran OpenFOAM中用label表示point、face和cell的size，默认情况下label为int32类型（WM_LABEL_SIZE=32）。抛出异常的原因很可能是由于对face size的计算超过了int32类型的最大限制，导致变为负数。若将label用int64类型表示（WM_LABEL_SIZE=64）则应该可以避免这个问题。不过所有label类型的变量占用内存是原来的两倍。

[wwzhao]

关于OpenFOAM版本，没必要追求最新，够自己用就行。

[wwzhao]

理论上限跟WM_LABEL_SIZE有关，若任何一个label数组长度超过2^(N-1)-1（对int32是2^31-1，对int64是2^63-1），则会抛出异常。

[random_ran]

好的，我再仔细研究研究 这个 LABEL_SIZE 的东西。有新的进展我再贴在这里。

(话说你的博客好赞啊！超喜欢这种风格。不知道是否接受 OF 热心读者的问题呢？直接留言在 Disqus 好还是在这里发帖，博主会看见呢？)

[wwzhao]

@random_ran

哈哈，多谢捧场。博客是用markdown写的，记录了一些自己对OpenFOAM的理解，初学者阅读可能比较吃力

有问题欢迎交流，Disqus跟论坛均可，我都能看到 :w: