Pointwise

网格导出流程

• [ ] Select Solver : OpenFOAM 3D
• [x] Set BC types : make sure “all boundaries” are set or you may encounter --> FOAM FATAL ERROR: Continuity error cannot be removed by adjusting the outflow. Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
• [x] Select blocks
• [ ] Export CAE to some “$CASE/constant/polyMesh”
• [x] Rescale to meter : transformPoints -scale ’(0.001 0.001 0.001)’
• [ ] checkMesh [注：这里的skewness感觉是有量纲的(1/L^3)，从mm变到m之后，skewness从O(1e-9)变到O(1)]
• [ ] serial run to check

配置OpenFOAM周期条件

语境
Pointwise里面网格为domain，OpenFOAM里面网格为patch

曾经天真的以为按照pointwise里面BC设定时选用cyclic，在网格文件生成完成后在constant/polyMesh/boundary里加上neighbourPatch就好，至少我在18.2R1这个版本不行。按照Maddalenna的提示。行得通的流程如下：

1. pointwise中要生成完全一样的domain，复制平移都不够（不够是指精度，重要的是点对点的对应关系），要达到这样的标准：create + periodic 这样生成的domain就是原先的twin，在domain list里面也可以看到domain之间的对应关系。这里就有了neighbourPatch的基础，没有这一步，极有可能会遇到face 1 area does not match neighbour *** by ***% -- possible face ordering problem.

2. 按照pointwise基本流程导出网格文件，但在设置BC的时候设置成patch而不是cyclic

3. 设置createPatchDict来完成patch到cyclic的转换. 默认createPatch 将生成新的网格文件，-overwrite可以改写原先网格文件，副产物还有一系列*.obj文件. 正确的Dict配置如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

//  File system/createPatchDict
//  姑且在这个语境中用“面”代替patch

    {
        // Name of new patch
        name front_cyclic;     // 操作完成后会是cyclic类的面1的name

        // Type of new patch
        patchInfo
        {
            type            cyclic;      // 新面1，类型
            matchTolerance  0.001;
            neighbourPatch  back_cyclic; // 对应的新面2
        }

        // How to construct: either from patches or set
        constructFrom patches;

        // If constructFrom = patches : names of patches. Wildcards allowed.
        patches ( front );    // 旧面1
    }

    {
        name back_cyclic;      // 新生成的会是cyclic类的面2

        patchInfo
        {
            type            cyclic;
            matchTolerance  0.001;
            neighbourPatch  front_cyclic;
        }

        constructFrom patches;

        patches ( back );     // 旧面2
    }

这里需要解释一下：新面1,2是由旧1,2生成，生成的条件有“类型，matchTolerance,neighbourPatch”，当新2还没有生成的时候怎么将新1与新2对应呢？当然，这里如果把第一个大括号当成了新1 patch的constructor，它将找不到新2的定义。不过回想一下，其实旧12和新12都互为twin的关系，互相依赖，这样写并无道理，只是…程序具体怎么实现的就不知道了.
注1：为了保证无误，新旧1和新旧2的名字没有取一样；想改成新旧一样？走完这个流程，最后把constant/polyMesh/boundary里面改过来就好
注2：createPatch在这里是由原先pointwise的twin domains/patches由类型patch变成了类型cyclic，据utility的功能介绍还可以将patches变成一个patch，或者将faceSet变一个patch

OpenFOAM

Simulation 完整流程

0. 环境变量
   a) purge modules
   b) set Foam environment : alias of301_intel
1. 网格
   a) transformPoints -scale '(0.001 0.001 0.001)'
   b) checkMesh [need a time dir (empty dir is OK)]
   c) 推荐使用格式上mesh变成binary，无论计算时的读写还是paraview读入都会更快 : foamFormatConvert -constant/-noZero [OF-2.3.1和OF-3.0.1在binary上不兼容，得通过ascii转换]
2. if needed, compile new Boundary Condition (which is a "*.so") for this simulation
3. 初始条件+边界条件
   a) IC (mapFields/topoSet)[OF-2.3.1 的 mapFields 就是个bug，推荐OF-4.x]
   b) BC (changeDictionary) [推荐此处IC,BC均用ascii，便于查错和纠正]
4. system/controlDict :
   a) check startTime 对不对, 保证startTime和endTime不相同 [如果相同OpenFOAM不会报错，log的末尾仍旧是”Finalising parallel run”]
   b) 如果有自己编译的模块，加入libs ("*.so");
5. serial run check
6. if 5 is correct, prepare for parallel run :
   a) decomposeParDict 指定MPI并行进程数
   b) 清除startTime目录里面’uniform’
   c) decomposePar -time 'startTime' : double check流场是否正确地被decompose. 首先写入硬盘的是constant/polyMesh的划分，后field transfer是流场的划分
   d) BC : double check 一下BC有没有被正确写入processor*[value可能被改写，但原则上member一定要都在，在BC编写时write(Ostream&)写对了就没有问题]
   e) if any warning : 要警惕，问题应当就出在当前步
7. 如果测试节点足够，可考虑interactive parallel test；如不够，配置slurm file:
   a) #SBATCH --job-name : sqeue 只能显示前8个字符
   b) 队列, 节点数#SBATCH --nodes，每个节点task数#SBATCH --ntasks-per-node，总task数需要在srun或者mpirun那一行用到（可以小于节点数*单个节点task数）
   c) 预估计算长#SBATCH --time
   d) 允许含有 bash variable
   e) 默认定向输出：标准输出#SBATCH --output；error输出#SBATCH --error
   f) [仅occigen] slurm执行任务时的计算环境初始化配置 : 为避免module相关输出到error(竟然module purge也会输出到error我也是醉了，而我python monitor目前只扫描error文件里面有没有输出，认为没有输出才是正常运行状态)，这里slurm file 就不再加入任何的环境配置，通过第一个步骤即步骤0里面实现slurm任务提交时正确的环境配置[登陆之后第一个load的module为默认slurm提交环境，如果想换个环境，得logout然后重新登陆]
   g) 注意 --exclusive 是否必要
   h) 注意 --mem 是否足够
   i) 主要任务执行行：即srun或mpirun那一行. 注意切忌行末加&幻想后台运行[例如后面一行还有其他executable的情况，如果当下行能后台运行，其后的command会被继续执行也许会有便利…但这样做的结果是&之后就没有然后了，而且还可能error message都没有，得不偿失]
   j) 还在主要任务执行行 : 通常将标准输出改到 > logFile : 注意最好把所有的log都留下，也就是每一个任务换一个文件名 [改到：因为前面谈到#SBATCH --output，用这个选项会由slurm在机群上的任务提交顺序来命名，可读性不强，通常改写]
8. [optional] submit job chain via python :
   a) [仅occigen] 因为$SCRATCH文件数目限制 : check-list python laundary
   b) check-list python auto submit
9. paraview : 有条件的话（因为通常processor文件数目很多）优先读decomposed case (因为internalField没有影响，但经测试reconstruct可能boundaryField的value会被篡改)

注：IC 指初始条件；BC 指边界条件

preProcessing

mapFields

此为大坑，尤其是OF-2.x的版本，存在一些bug(不是fatal，但会让mapFields运行得无比慢，比fatal还可恶。按照帖子改了还是不行)，但用同样的mapFieldsDict试一试OpenFOAM/4.0-foss-2016b或者OpenFOAM-5.x不仅速度快而且不会有莫名其妙的报错

1
2
3
4
5

#!/bin/bash
source=/some/case/to/be/mapped
#taget 这里target就是"."

mapFields $sourcedir -case . -noFunctionObjects -fields '(U p T)' -sourceTime '8' -targetRegion region0 > log.mapFields_8 2>&1

0. 回避
   a) $source里empty的BC（有试过，会有报错）
   b) 如果$source里面constant/polyMesh/boundary里面有mappedPatch，且如果$target里面没有相应的BC配置，可能会在mapFields最后写入数据的时候报错后果是例如U文件的写入遇到错误而被跳过].在$source有mappedPatch的情况下，$target里面constant/polyMesh/boundary也得改成相应BC
1. 如果不是consistent，目标case里面要编辑好文件mapFieldsDict
2. 确认$source里面startTime，它会是mapFields完成后的时间目录，输出格式改为ascii
3. 检查在controlDict里面libs (...)是否在当下环境中配置
4. mapFields (等待时间可能很长)
5. 手动检查mapFields是否无误地完成：检查目标case里面是否有-nan
6. 把映射后的场的BC由calculated改成相应的物理BC，这样才可以续算，这个步骤可以通过changeDictionary来完成
7. 串行试运行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

# 一个长度为5D的圆管映射到一个长度为10D的圆管，inlet*完全对应
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

patchMap        ( inlet0 inlet0         // Source 和 Target 都有的 patch
                  inlet1 inlet1         // 且一一对应
                  inlet2 inlet2
                  inlet3 inlet3
                  inlet4 inlet4 );

cuttingPatches  ( wall                  // 剩下 wall patch 和 outlet*
                  outlet0               // 此处为空->(U p)在 wall 和 oulet* 会出现异常大数甚至nan!!
                  outlet1               // 只填 wall， wall上面插值没有异常大数出现
                  outlet2               // 填上 wall patch 和 outlet*, outlet*都乖乖地保持了“最简单状态” (虽然U在 
                  outlet3               // outlet*上面不应该是 uniform (0 0 0)，但毕竟是 type calculated
                  outlet4 );

// ************************************************************************* //

# 两个T型圆管，几何完全一样，但网格不同
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

patchMap        ( Port1 Port1        
                  Port2 Port2
				  Port3 Port3
				  wall  wall  );

cuttingPatches
();
// ************************************************************************* //

旋转data然后mapFields

想要rotate data，transformPoints声称可以rotate polyMesh里面的points(也就是网格)，也可以rotate vector field.试过了，确实可以，paraview上就看出来转了90度，但是将rotate过后的再mapFields就不成功了，经过反复测试始终还是rotate之前的data被map过去了的感觉

postProcessing

管理postProcessing

rsync.py将远程的postProcessing同步到本地目录，并通过json将相关数据存入一个文本文件database.txt:

1. sourceDir:远程目录
2. targetDir:本地目录
3. name2plot:画图的时候用到的legend
   还有alias，即一个对case简单描述的string

画图的时候，在配置文件parameters_*.py里面有：读database.txt，通过alias来索引database[alias]['targetDir']和database[alias]['name2plot']

一旦按照以上操作“1.获取数据”并“2.抓取数据画图”，只需要确认有一致的alias就可以确认“画的数据是自己想要的”，具体的数据信息也可以通过database.txt来最终索引，只要这个索引不出错，就不会错

uv : reynoldsStress

1
2
3
4
5
6
7
8
9

# toolchain for uv_mean @12
userAverageVector U -time '6:12' # U_mean
userReyTensor_plusAvg U 12 -noIntermediateWrite -time '6:12' # reyTensor_mean, 不输出中间的reyTensor速度还是比较快的
userRMSuvw_fromReyTensor_mean -time '12' # uv_mean
sample -dict $sampleDicts/sets/set4uv -time '12' > log.set4uv 2>&1  # 需要sampleDict : set4uv

# 需要注意的在这里
# 以上最后一步输出有128个sample的结果，但只用其中64个来算uv，因为上下壁面会相互抵消(uv的定义相对于壁面，v的定义在两个壁面是反的)
# 这就是[parameters_](https://github.com/snow-stone/python/blob/master/userSetUV/parameters_t_r-2a_1_gradP0p703125.py)里面`'raw_sample_size':64`的原因

userProbeByLabel_noMean

or userProbeByLabelVector_noMean its version for vectors. $sliceStore is where you put your slice+number file (class labelList). userProbeByLabel_noMean don’t work on slices but on labelGroup which is itself a labelList containing all cell id of probes of interest (its defaut location is $sliceStore.

0. check-list userFindClosestInLabelList
1. userProbeByLabel_noMean T $sliceStore -time ‘0:xxxx’

userFindClosestInLabelList

前提是在input argument里面sliceStore处已经有sliceNumberList和refVectors，这个程序找的是遍历在sliceNumberList里面所有的slices找到离refVectors（注：这里仅有一个vector，但仍旧用vectors，为得是日后可以拓展）最近的cell，并将每个slice上面最近的cellID记下来写入到一个文件labelGroup.这个程序应当在一个OpenFOAM case中执行，其实就是提取网格信息，因此对于相同网格下的算例仅仅需要运行一次，于是就有了最后一个步骤将labelGroup移动到对应同一个网格的公共的目录下共享，详细check-list如下:

0. 为了最终写出的labelList人类可读:system/controlDikt change format to ascii
1. userFindClosestInLabelList xx xx
2. cd constant
3. mv file to a shared place (for all cases with the same mesh structure)

reconstructPar

在reconstructPar -fields '(U p)'之后，reconstructPar -fields '(phi)'会将phi添加到对应的时间目录里面.

sample

基本流程

1
2
3
4

#!/bin/bash
sampleDict=system/sampleDict

sample -dict $sampleDict/sets/Dai_lines_5cutsBetween2and3_typeFace_cell -time '8:10' > log.sample-Dai_lines_5cutsBetween2and3_typeFace_cell 2>&1

1. 编辑sampleDict：可能需要用python脚本来写入一系列sets的描述，例如160条线就不能全部手写
   a) 检查fields，如果写得不对，OpenFOAM并不会报错
   b) 检查是否有header（没有header会有报错），object填sampleDict
   c) sets或者surfaces的输出文件名都可以个性化编辑
2. 重命名sampleDict，是sets放目录$sampleDict/sets，是surfaces放目录$sampleDict/surfaces
3. sample

注：检查不出sample是否有结果？rm -rf postProcessing ; sample 这样会比较明确

用python写sampleDict

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

# file writeSampleDict_2Diagonals.py

import numpy as np
import os

def writeStartEnd(whichLineGroup,f,z,nPoints):
    if whichLineGroup == 0:
        f.write("\t\t"+ "type\t"  + "uniform;" +"\r\n")
        f.write("\t\t"+ "axis\t"  + "distance;" +"\r\n")
        f.write("\t\t"+ "start\t" + "("+z+ " 0.004 0.004);" +"\r\n")
        f.write("\t\t"+ "end\t" + "("+z+ " -0.004 -0.004);" +"\r\n")
        f.write("\t\t"+ "nPoints\t" + str(nPoints) + ";" +"\r\n")
    elif whichLineGroup == 1:
        f.write("\t\t"+ "type\t"  + "uniform;" +"\r\n")
        f.write("\t\t"+ "axis\t"  + "distance;" +"\r\n")
        f.write("\t\t"+ "start\t" + "("+z+ " -0.004 0.004);" +"\r\n")
        f.write("\t\t"+ "end\t" + "("+z+ " 0.004 -0.004);" +"\r\n")
        f.write("\t\t"+ "nPoints\t" + str(nPoints) + ";" +"\r\n")
    else:
        print "big problem"

def writeLines(f,z_value_List,N,Nb_lineGroup,resolution,dictName):
    for i in range(Nb_lineGroup*N):
        whichLineGroup = i // N  # quotient  0, 1 ,2... Nb_lineGroup
        j = i % N                # remainder 0, 1, 2... Ns
        
        f.write("\t"+dictName+'-'+str(i)+"\r\n")
        f.write("\t{\r\n")
        writeStartEnd(whichLineGroup,f,z_value_List[j],resolution)
        f.write("\t}\r\n")
        f.write("\r\n")        
    
def write_LineSet(f,z_value_List,Nb_slices_along_z,Nb_lineGroup,resolution,dictName):
    f.write("sets\r\n")
    f.write("(\r\n")
    writeLines(f,z_value_List,Nb_slices_along_z,Nb_lineGroup,resolution,dictName)
    f.write(");\r\n")

def write_headerSetting_sampleLines(f,fieldName,interpolationScheme):
    f.write("setFormat raw;\r\n")
    f.write("\r\n")
    f.write("interpolationScheme "+interpolationScheme+";\r\n")
    f.write("\r\n")
    f.write("fields\r\n")
    f.write("(\r\n")
    f.write("\t"+fieldName+"\r\n")
    f.write(");\r\n")

def write_sampleLines(dictName,fieldName,interpolationScheme,z_value_List,Nb_lineDisplacement,Nb_lineGroup,resolution):
    f= open(dictName,"w")
    write_headerSetting_sampleLines(f,fieldName,interpolationScheme)
    write_LineSet(f,z_value_List,Nb_lineDisplacement,Nb_lineGroup,resolution,dictName)
    f.close()
    
def NLines2Dict(Nb_lineDisplacement):
    l = np.linspace(-0.08,-0.04,num=Nb_lineDisplacement)
    listStr = [str(x) for x in l]
    print "A list of all line displacements :\n"
    print listStr
    
    dictName="lineOn2Diagonals"
    fieldName="U_mean"
    interpolationScheme="cell"
    write_sampleLines(dictName,fieldName,interpolationScheme,listStr,Nb_lineDisplacement,2,resolution=200)
    
    os.system('dos2unix '+dictName)
    print "\n"
    print "sampleDict written. Think to add a file header!"

NLines2Dict(64)

1. configure python script for sampleDict : dictName, fieldName, interpolationScheme… Sampled data will be named as `postProcessing/sets/time/dictName-n_fieldName. “dictName” may also includes information like “interpolationScheme”.

2. run python script

3. add header to sampleDict

4. run sample utility

5. verify generated files in postProcessing

6. begin plotting using python

7. edit plot parameters : a spatial statistic will require “time”, “number of samples” etc. This file is the footprint of the statistics.

8. configure plot script : ajust “dataShape” or likewise. In other word, not all data are eligable because not all data are of the same shape. We can interpolate but I d rather not.

9. run plot script

10. ajust plot and save fig to figure

11. write plot data (x,y) to txt file in foler data for reporducing use or later ajustement

python

laundary

为了防止occigen $SCRATCH关于文件数量的quota溢出，用python在后台nohup长时间地监控，对userDefinedLog里面removedTimes与dataWritingHistory的补集的时间步做reconstructPar和rm -rf processor*/timeStep

0. 拷贝laundary.py和reconstruct_occigen-OF301.sh
1. 保证userDefinedLog/removedTimes存在
2. 设置python脚本的时间参数（预估simu时长，留出laundary余量）和processor个数的参数
3. 检查上一个laundary是否结束:ls -l log.clean_*看时间
4. 判断reconstructPar的进程:ls processor0 | xargs ls -l > log.ls；如果有的目录没有顺利reconstructPar会有内容(非std output)输出到终端，至于目录对应上了但场是否完全(比如U p nu phi) reconstructPar 只能看log.ls
   a) reconstructPar有很多时间步积压，手动完成laundary，然后进行b)
   b) reconstructPar运行顺利，没有太多时间步积压在processor*里，那么准备续算：rm userDefinedLog/dataWritingHistory并清空userDefinedLog/removedTimes里面的内容
5. nohup python laundary.py OF301 > log.clean_*
6. 记录是在哪个登陆节点，需在对应的节点用ps -eaf | grep $USER才能找到nohup的任务

cluster

occigen

算例提交

0. openfoam env : 一句话如果用python来做slurm多个连续任务提交，必须在提交前做好这一步(即check-list第0步)，具体来说如果忘记了提交slurm会在第一个job报错，然后你发现slurm找不到mpirun，这时候补上环境变量再重新提交，新的提交第一个job能通过，但第二个job仍然会报错找不到 mpirun；似乎slurm的默认环境是你登陆某login节点后第一次提交任务的环境，也就是说：login之后得首先做这个事情再提交任何算例(虽然直接使用sbatch **倒是不影响，但对通过python来提交的job chain会被迫终断）
1. nohup python submit.py > log.submit : 提交job，通过squeue的返回值离散地监控job状态，通过返回值来检验一个算例完成后(这里还有个bug)修改system/controlDict里面startTime和endTime用于自动提交下一个算例(限制：simuLog不停地被改写，应该把每个job的log都留下来才好)；log.submit取好名字以方便查看
2. nohup python watchDog.py > log.watchDog : 设置一个运行最大时长，在这个时间内用reconstructPar来保留计算输出数据，删除processor*里面已经reconstructPar完成地时间步大幅度削减文件个数，保证不超过scratch的限额；log取好名字以方便查看
   G3. 笔记本上面记录下来是哪个login，什么算例，连续提交多少个job，每个job预计计算的物理的间隔是多少(用于更新startTime和endTime)：因为ps -eaf | grep $USER只能找出相应登陆节点上面的job
3. 检查log.submit，查看jobs的情况

数据同步

千万千万要注意para0，在dir里面并没有用到，一定要double check !!!!

1
2
3
4
5
6
7
8

#!/bin/bash

para0=inlet_0p3
para1=a_0p08
para2=setT_St_1
log=sync_log.$para1
dir=/store/lmfa/fct/hluo/occigen/caseByGeometry/T/shape_square/2a_3_T/BirdCarreau/inlet_0p3/$para1rsync -av occigen:/scratch/cnt0028/mfa0464/hluo/caseByGeometry/T/shape_square/2a_3_T/BirdCarreau/$para0/$para1/$para2/* DATA --exclude processor* &&
echo "DATA sync BirdCarreau : $para0 $para1 $para2 ended with success" >> $dir/$log

newton

1. No #SBATCH --exclusive (if not necessary)
2. Naming of the file with exactly 8 chars with -2 indicates for example 2nd run. Ex : p10D_gP5-2
3. consider remove #SBATCH --mem-per-cpu=4000 because there are nodes which bigger memory available

error

1. error while loading shared libraries: libpsm_infinipath.so.1: cannot open shared object file: No such file or directory : check #SBATCH --mem, this is an error saying you are maybe asking for more memory than the machine can offer
2. slurmstepd: error: Detected 1 oom-kill event(s) in step... : check #SBATCH –mem`, this is an errory saying that the memory you are asking for is not enough for the program

old note

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

#!/bin/bash
# FILE : p10D_gP5-2

#BATCH --job-name=MeshConvergence                    # Is this even useful ? SABTCH no?
#SBATCH --output=job.%j.out
#SBATCH --error=job.%j.err 
#SBATCH --mail-user=haining.luo@doctorant.ec-lyon.fr
#SBATCH --mail-type=ALL
#
#SBATCH --partition=mononode
#SBATCH --mem-per-cpu=4000
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=16
#SBATCH --cpus-per-task=1
#SBATCH -t 1-00:00:00

# number of simulation
nos=2
logFile=logSimulation$nos

# 在newton上module purge load部分可以保留
module purge
# Opt/Debug mode on
module load OpenMPI/1.6.5-GCC-4.8.3 ; \
. /home/lmfa/hluo/LocalSoftware/OpenFOAM/OpenFOAM-2.3.x/etc/bashrc.Opt

echo "NodeName = " $SLURMD_NODENAME >> $logFile
echo "JobID = " $SLURM_JOB_ID >> $logFile

mpirun -np 16 pimpleFoam -parallel >> $logFile

3. run mapFields need a big memory. There’s now a template for big-memory job. Especially, don’t put & at the end of commands or mapFields will only last for 2s and there’s no error message returned so could be confusing.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

#!/bin/bash
#
##SBATCH --job-name=T-1b_mirrorMerge_mapped
#SBATCH --output=job.%j.out
#SBATCH --error=job.%j.err 
#SBATCH --mail-user=haining.luo@doctorant.ec-lyon.fr
#SBATCH --mail-type=ALL
#
#SBATCH --partition=test
#SBATCH --mem=32000                  # asking for 32000Mb ~ 30Gb of memory 
#SBATCH --nodes=1                    # 1 node
#SBATCH --ntasks-per-node=1          # one task for one node
##SBATCH --cpus-per-task=1           # one cpu per task
##SBATCH -t 1:00:00

# number of simulation
nos=1
logFile=logSimulation$nos

module purge
# debug mode on
module load OpenMPI/1.6.5-GCC-4.8.3 ; \
. /home/lmfa/hluo/LocalSoftware/OpenFOAM/OpenFOAM-2.3.x/etc/bashrc.Opt

# run solver 
#decomposePar  # Do this manually

echo "NodeName = " $SLURMD_NODENAME >> $logFile
echo "JobID = " $SLURM_JOB_ID >> $logFile

#mpirun -np 64 icoFoam -parallel >> $logFile

sourcedir=/home/lmfa/hluo/LocalSoftware/OpenFOAM/hluo-2.3.x/run/1b_mirrorMerge_mapped_NearestFace
mapFields $sourcedir -case . -noFunctionObjects -fields '(U p)' -sourceTime '6.01' -targetRegion region0 > logMapFields_new 2>&1  # No "&" here

Here maybe a reason that zaurak well perform likely 2 times better than cpus on newton when the constraint on memory is not there.

1
2
3
4
5
6
7
8
9
10
11

hluo@zaurak $ cat /proc/cpuinfo | grep model
#...
model name	: Intel(R) Xeon(R) CPU E5-1620 v4 @ 3.50GHz

hluo@compil-haswell $ cat /proc/cpuinfo | grep model
#...
model name	: Intel(R) Xeon(R) CPU E5-2660 v3 @ 2.60GHz

hluo@visu $ cat /proc/cpuinfo | grep model
#...
model name	: Intel(R) Xeon(R) CPU E5-2640 v4 @ 2.40GHz

Inkscape

导出

• [ ] draw sketch
• [x] Edit-> Resize page to selection
• [ ] Export PNG image -> Export as -> Export

在最后一步如果遇到”the chosen area to be exported is invalid”,可以在Export PNG Image标签页里面四个标签Page,Drawing,Selection,Custom（默认是最后一个），选成比如说Drawing再点Export.

clip

1. import image
2. F4 (select “rectangular tool”) draw a rectangular (marked by red)
3. F1 (select “select tool”), ctrl+A select all
4. Object -> clip -> set

gimp

这个软件可以得到图片里面的坐标

paraview

显示几何外形

如果想要显示一下外形，但又不想要网格视角:
Clean to grid -> Extract surface -> change opacity

scripting

load state

load state也可以被载入到script里面(no legend)

I. 用paraview生成*.psvm

1. paraview -> start trace -> open *.foam -> …
2. saveFig (very possible I will do this. otherwise its nonsense.)
3. save state
4. stop trace

II. 拿出上一步得到的script，把关键行(例如)

1
2
3
4

LoadState('/store/T_c/1j/D2-NN-1j_test_from0p3_forcingSinus_St3p2_A_eq_0p05/hluo15_T_c_vorticity_z_300.pvsm', LoadStateDataFileOptions='Use File Names From State',
    DataDirectory='/store/T_c/1j/D2-NN-1j_test_from0p3_forcingSinus_St3p2_A_eq_0p05',
    OnlyUseFilesInDataDirectory=0,
    D2NN1k_syn_forcingfoamFileName='/store/T_c/1j/D2-NN-1j_test_from0p3_forcingSinus_St3p2_A_eq_0p05/D2-NN-1k_syn_forcing.foam')

转换为

1
2
3
4
5
6
7

SaveScreenshot(dirName+'/'+'hluo15_T_c_vorticity_z_300.png', renderView1, ImageResolution=[2754, 1838],
    FontScaling='Scale fonts proportionally',
    OverrideColorPalette='',
    StereoMode='No change',
    TransparentBackground=0, 
    # PNG options
    CompressionLevel='5')

尾巴上再加入

1
2
3
4
5

import os
print "Finalizing "+ os.path.basename(__file__) +" " + "@ dir : " + dirName 
import datetime
print datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
print "==============================="

III. 对不同目录（不同case）批量操作:

1. 拷贝*.psvm到相应目录（如果pvpython找不到*.psvm也会出图，白板）
2. 修改*.psvm里面带有关键字.foam的行，对应当下目录
3. 来到run.sh，用绝对路径作为参数开始批量运行

出图程序里面修改变量的range

LUT.RGBPoints和PWF.Points里面上下界都得改，PWF.RescaleTransferFunction也要改，当然LUT.UseLogScale也要改：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# get color transfer function/color map for 'k_mean_nonD'
# ...
k_mean_nonDLUT.RGBPoints = [1e-6, 0.231373, 0.298039, 0.752941, 0.00014459355054441403, 0.865003, 0.865003, 0.865003, 1, 0.705882, 0.0156863, 0.14902]

# get opacity transfer function/opacity map for 'k_mean_nonD'
# ...
k_mean_nonDPWF.Points = [1e-06, 0.0, 0.5, 0.0, 1, 1.0, 0.5, 0.0]

# Rescale transfer function
k_mean_nonDPWF.RescaleTransferFunction(1e-06, 1)

# Properties modified on k_mean_nonDLUT
k_mean_nonDLUT.UseLogScale = 0

纯粹读数据，摸索数据结构

似乎是默认读最后一个时间步，为免出错，每个case就一个时间步

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()

# create a new 'OpenFOAMReader'
cavityfoam = OpenFOAMReader(FileName='/home/hluo/OpenFOAM/hluo-2.3.1/run/tutorials/incompressible/icoFoam/cavity_org/cavity.foam')

print type(cavityfoam)
print cavityfoam.CellArrays
print type(cavityfoam.CellArrays)
print "GetAvailable()"
print cavityfoam.CellArrays.GetAvailable()
print "GetData()"
print cavityfoam.CellArrays.GetData()

data=servermanager.Fetch(cavityfoam)
print type(data)
print "GetNumberOfBlocks() : " , data.GetNumberOfBlocks()
#print "getBlock() : " , data.GetBlock(1)
obj = data.GetBlock(0)
print obj.GetNumberOfCells()

cellData = obj.GetCellData()
N = obj.GetNumberOfCells()
#print cellData
obj_p = cellData.GetArray('p')
print type(obj_p)
print "obj_p.GetDataTypeValueMin() : ", obj_p.GetDataTypeValueMin() # 这里得出的最大最小有问题
print "obj_p.GetDataTypeValueMax() : ", obj_p.GetDataTypeValueMax()

print obj_p.GetValue(10)

print "GetRange() : ", obj_p.GetRange() # 这里貌似是正确的
for i in range(N):
    obj_p.SetValue(i, obj_p.GetValue(i)/10)

print obj_p.GetValue(10)

export screen shot

export 有什么好说的？嗯……如果要对一系列算例做同样的图一般会选用load state来加载*.pvsm，似乎注意改一下*.foam对应的path就好，但是如果要让export screen shot输出同样像素的图的话，记得一定要全屏之后再export！！

do not skip time 0

1. paraview data.foam
2. uncheck all “Cell Array” ; remain defaut “Skip Zero Time” then apply : geometry will be visualized but no “Cell Array” (field data) is there. Make sense.
3. uncheck “Skip Zero Time” -> apply : cell array will then appear ; select field of interest.

color map

• 如果是对成的数据[-a,a]，用红白蓝diverging挺好，能分辨出0对应白
• 如果是[0,a]，用黑白灰最好，但paraview好像默认可以从Edit Color Map选项卡中带桃心的小按钮Choose Preset里面有X-ray和GrayScale，选择后apply（默认就会变到RGB color Space）；如果想要恢复diverging，最下面有个恢复默认按钮；如果想要自定义，参见RBG自定义

用load state来复现camera视角

target 目标视角：想要复制的视角，对应的case叫目标case
working 工作视角：想要在工作case下复现目标视角

0. 目标视角的存储通过目标case(Visu)里save state来实现(默认读取了一个绝对路径的但其实为空白的target.foam文件)，得到target.pvsm
1. 复制目标case下的target.pvsm到工作case(Visu)里，编辑查找关键字target.foam并替换成工作路径/working.foam
2. 工作路径下创建working.foam
3. 工作路径下打开paraview
4. load state 选择编辑后的target.pvsm-> “Load State Data File Options” 选Use File Names From State
5. 等待复现

注意：

1. 还涉及一个working case里面时间步是不是和state里面一致的事情，我的测试刚好target和working case有相对应的同一时刻的data
2. 在paraview-5.4.1测试成功

根据x坐标截取一个cell slice的label

cell slice : 这里的网格是正交的六面体，所以严格意义上一个x坐标对应“一层网格”，通过在paraview里面find Data中选取xcoord between xx and xx就可以选定这一系列cell
label : 一系列的cell的label合起来就是个labelList，而且是全局的，因此可以取出U中例如U[labelList_slice0]子集进行操作

下为流程：

0. paraview 可视化 case
1. calculator : xcoord (defaut : point data)
2. Filters -> Point data to cell data 这里我们只想对xcoord操作 但注意 这一步会将U(如果有的话)也插值一遍，值会跟先前有略微差别
3. Edit -> fidn Data -> (cells) (Pointdatatocelldata) (criteria : xcoord between [xmin, xmax] -> run selection
4. close “find Data” (selection is still effective)
5. split horozontal (very smalll icon on the right above render view) -> spread sheet -> show only selected elements -> attribute : cells 注意一定要选cells要不然得到的spread sheet行数不是cell number -> check 行数是否等于画网格时定下的个数
6. Toggle cell visibility (eliminate other colons : left with only labelIDs) -> export spread sheet *.csv -> optional depending on paraview version (filter colons by visibility)
7. modify *.csv : add OpenFOAM header (shown below), add ( and ) for list
8. use OpenFOAM to read modified file via class IOList<label>

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

// example of OpenFOAM header for file named "slice16"
$ head -10 slice16
FoamFile
{
    version     2.0;
    format      ascii;
    class       labelList; // IMPORTANT
    object      slice0;    // not important
}

16900
(
/*
 ...
*/
)

scripting

Base

利用Tools->Start Trace选上all properties Fully Trace Supplemental Proxies Show incremental Trace

0. Open -> some *.foam file (no need for the real file when running script afterwards, but necessary here for GUI use, to generate the script to run)
1. clip… slice…
2. export as
3. Tools -> End trace
4. pvpython script.py

Note

1. 值得注意的是legend，利用以上python trace的结果不加改动的情况下，会默认变得巨大，所以得加上
2. 加上auto-rescale是个必须得有的好习惯
3. 以上1和2的改动放在renderView1之前，经测试是有效的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

"""
some operations : read some certain array
"""
# Properties modified on t_meanLUTColorBar    调整legend nu_mean** 需要根据实际object名字替换
nu_meanLUTColorBar.TitleFontSize = 4
nu_meanLUTColorBar.LabelFontSize = 4
nu_meanLUTColorBar.ScalarBarThickness = 5
nu_meanLUTColorBar.ScalarBarLength = 0.5

# rescale color and/or opacity maps used to exactly fit the current data range  等同于auto-rescale
slice1Display.RescaleTransferFunctionToDataRange(False, True)

# current camera placement for renderView1    
renderView1.CameraPosition = [0.013999999035149813, -0.08839503035519257, 0.0]
renderView1.CameraFocalPoint = [0.013999999035149813, -0.0020000000949949026, 0.0]
renderView1.CameraViewUp = [0.0, 0.0, 1.0]
renderView1.CameraParallelScale = 0.022360679233547745

# save screenshot

4. customize range, logscale, use another color mapping 它们之间按照以下顺序排列是有效的

1
2
3
4
5
6
7
8
9
10
11

# Rescale transfer function
nu_meanPWF.RescaleTransferFunction(2e-06, 0.0003)

# convert to log space
nu_meanLUT.MapControlPointsToLogSpace()

# Properties modified on nu_meanLUT
nu_meanLUT.UseLogScale = 1

# Apply a preset using its name. Note this may not work as expected when presets have duplicate names.
nu_meanLUT.ApplyPreset('Yellow - Gray - Blue', True)

5. 如果有多个不同位置的slice，那么一定要记得加上下面这一行(在save screenshot前就行)以保证图”有效部分”的大小不会改变

1
2
3
4

# reset view to fit data
renderView1.ResetCamera()

# save screenshot

possible bug

1
2

raise ValueError("%s is not a valid value for attribute %s." % (value, name))
ValueError:  is not a valid value for attribute ScaleArrayComponent.

可能是最初没有读入后面要用到的数组，也有可能后面取slice什么的超出了计算域，感觉报的是个找不到目标数据的错

blender

一句话：通过paraivew做初期可视化，包括color map，导出后可通过blender调整camera和lamp. Paraview里面也有这个功能，但相对鸡肋.

blender version : 2.8 beta compatibility issue : glibc too old in CentOS

1. paraview -> visu -> File -> save Data -> *.ply (check coloring; check alpha : even if you dont have any transparency in paraview setting. This is important for blender/ or shall we say bug)
2. blencer -> File -> import -> *.ply
3. (optional but often needed) s0.01（导入的data通常会显得很大，需要缩小）: scale by 0.01;s90x以x为轴旋转90度；blender里面朝上通常为z
4. Now we have a render vie for data. Care now for coloring. -> choose blender render or cycles render
5. Add -> input -> attribut (the color attribute) : (i) name = Col ; (ii) connect “color” to “base color” of BSDF
6. 如果import的数据是contour，可能会有尖角出现，光滑一下：select mesh->object->shade smooth (need to be in object mode : normally by defaut)
7. 接下来就是调整camera和lamp了，最终blender出图是camera的视角，lamp可选太阳光什么的

Note : Gratitude to Juan Ignacio Polanco

svn

branching

video reference createBranch,workWithBranches,resolvingConflicts,svn_resolve_tree_conflict_in_merge

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

# Note a good pratice : make sure trunk is updated to the lastest version
# create branch by copying trunk. 
$ svn copy https://subversion.renater.fr/jnnf/trunk https://subversion.renater.fr/jnnf/branches/haining -m 'create branch haining'

# confirm existance of branch/haining in svn server
$ svn list https://subversion.renater.fr/jnnf/branches/ 

# get branch data from svn server
[branches] $ svn update # better than : svn checkout https://subversion.renater.fr/jnnf/branches/haining. This will give problems when doing svn delete at the end of checkList

# make changes on branch/haining, running tests
# if new file added : svn add
[haining] $ svn commit -m 'information on modification'

# make sure branch itself is updated 
[haining] $ svn update 

# now we are going to push changes to trunk. BUT
# firstly we have to merge from trunk to this branch (necessary step):
# possibly we may get some conflicts for THIS merge.
[haining] $ svn merge ^/trunk

# verify merge by running tests
# commit merge
[haining] $ svn commit -m 'merge from trunk to branch/haining is sucessful. ready to reintegrate into trunk'

# go to trunk
[trunk] $ svn update

# now re-integrate from branch/haining to trunk then managing possible conflicts (usually choose postpone and vim)
[trunk] $ svn merge --reintegrate ^/branches/haining

# verify again changes by running tests
# If all is well, we commit
[trunk] $ svn commit -m 'merged trunk from branches/haining'

# branches/haining is now useless
# thus remove them from version control meaning on the svn server. 
[branches] $ svn rm haining 
[branches] $ svn commit 'removing branches/haining'

# wont be able to see 'haining' anymore
$ svn list https://subversion.renater.fr/jnnf/branches/ 

# checking status :
[branches] $ svn status
?       haining        # ? means no longer in version control

resolving conflict

!M : svn rm file --force (if it is rm * not by svn remove)

out of date : svn update; svn resolved fileName (one by one) ; svn commit -m 'fix out of date'

svn和git之间相互转换

一定要非常注意！

0. 如果是我的thesis，先make clean
1. cp -r拷贝整个文件夹（最好不要是根目录，因为包含.svn或者.git）到目的地文件夹
2. 如果是根目录继续往下看
3. 如果是svn拷贝到git，一定要删掉.svn ; 如果反之，一定记得删掉 .git (.gitignore倒是无所谓啦)
4. 然后才开始进行其他操作

ubuntu in win10

Terminal默认不能输入中文，vim打开也会乱码，需如此配置:

1. Terminal上端right-click->properties->font，改为NSimSun : 这样在终端就可以敲出中文
2. vim 打开含中文字符的文本，:set encoding?如果不是utf-8`se:set encoding=utf-8`: 这样vim里面的中文字符就会正确解码
3. 将以上设为默认right-click->defaut;.vimrc

IO

IO的基本特征是硬盘读写，耗时间。

后处理里面的IO

目前我见过的OF版本里面后处理的编程范式，仍旧是读一个时间步进行操作再读下一个时间步，这样一段IO之后一点操作操作完了输出一个时间步的结果，然后下一个时间步再一段IO……可能会效率很低：也许来得不如先集中读入所有数据（内存足够大），然后再处理后集中输出。不过要能够这样处理的前提除去内存，对OpenFOAM里field类读写操作的类需要在时间维度上变成动态数组。

我的认识：读入field的时候确实是动态的数组，但在时间维度不是。

Serial or parallel ?

早期OpenFOAM

进行并行计算时要decomposePar将计算域分为多个部分，网格文件实际上也分开了。在写入的时候也是写入到每个processor文件中，这是早期OpenFOAM版本的做法，弊端是产生的文件量特别大。有些机器上面，比如occigen就有文件个数的限制，只允许200000个文件。如果500个processor，一个processor里面最多只能400个文件，而一个时间步里面至少有U, p, phi，所以这个限制还挺严格。

如何解决？要么定期reconstructPar然后把那些processor里面对应的U, p, phi删掉，需要注意的有两点：
1.第一个操作完成，第二步删除再进行
2.不要对OpenFOAM正在读写的文件进行操作
第2点比较好规避，选择所有时间步，避开倒数2个，剩下的一般都写好了（因为一般输出数据间隔都不会很小）；而第1点，就需要保证第一步的返回值为0，经由判断后再进行相应删除

较新的OpenFOAM

openfoam-plus_1712和openfoam5提供了一种直接输出一个整体流场U, p, phi的方式。(吐槽一下：这应该是一个并行CFD的基本功能才对)

经尝试想要在早期版本中编译并行输出的模块很复杂，弃。

数据格式

binary or ascii
模拟的数据输出可以选择，各有优劣：

• binary数据存储体积更小，推荐在simu中间使用，而且在paraview读取数据的时候（20M算例，两个时间步）至少比ascii要快7倍
• 在初始场中推荐使用ascii，因为万一要改边界条件binary格式不是那么好下手(vimdiff或者meld都对ascii支持更好)
• 再者OpenFOAM支持一个simu里面两种格式文件共存，utility和solver在读数据的时候都是以文件头上format相应格式读取

binary或者ascii，如何在二者中转换？
复杂解且不完全：比如从binary转为ascii，可以用decomposePar然后reconstructPar，第二步的时候把输出格式改为ascii，internalField的读写应该没问题，但是BC呢还是得check一下
简单解且官方：foamFormatConvert，用system/controlDict里面的格式来写输出（不过这会overwrite原数据）,默认会将constant/polyMesh里面的数据也按照格式重写

foamFormatConvert也不是一直都行得通：ascii格式的文件是OF-2.3.x/OF-3.0.1兼容的，但binary就不是，遇到以下报错需考虑版本问题(另外一点窍门是转换的时候如果方便，把constant和system移动到与时间步分离的目录下运行foamFormatConvert -constant，这样可以避免时间步的数据被影响)

1
2
3
4
5
6
7
8
9

--> FOAM FATAL IO ERROR: 
Expected a ')' while reading binaryBlock, found on line 20 an error

file: synthetic_phasedStepFrom_test_from_0/From0p3_3_of3/constant/polyMesh/faces at line 20.

    From function Istream::readEnd(const char*)
    in file db/IOstreams/IOstreams/Istream.C at line 111.

FOAM exiting

相关的类

UListIO.C里有底层的operator<<，判断了是否输出的是ascii后判断是不是uniform的UList，如果是uniform就用{}，如果不是就用()，这就是为什么OpenFOAM里面只要输出List就一定带有大括号或者小括号.

典型读写范例

bug

据测试…读的时候文件可以用绝对路径来且可以较长，但写的时候得注意绝对路径不一定奏效（即使确定了路径存在），为啥要说这一点，因为“不奏效”指无warning且不报错（如果debug不开的话）但就是找不到输出文件，写的时候用constant作为输出路径比较保险（当下目录.也需质疑）

读写范例

通过IOobject读入写场，注意，OpenFOAM里面场是*Field而Field是一个List，所以其实是List读写(上面写到还跟UList有关)

读

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

// (始) 读一个名为velocityFieldName的向量场
volVectorField velocityField
(
    IOobject
    (
        velocityFieldName,
        runTime.timeName(),   // 所在目录
        mesh,
        IOobject::MUST_READ
    ),
    mesh
);

// (结) 后面引用velocityField即可，注意这里只读了runtime.timeName()里面的，文件名为velocityFieldName这个文件

// （始）读一个名为"labelGroup"的labelList 
IOList<label> labelGroup
(
    IOobject
    (
        "labelGroup",
        position,             // 所在目录
        runTime,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
);

// 附上"labelGroup"文件的样本格式
FoamFile
{
    version     2.0;
    format      ascii;
    class       labelList;
    location    "constant";
    object      labelGroup;
}

21
(
26524615
25273885
24395085
23786685
23330385
22975485
22671285
22417785
22198085
21995285
21826285
21657285
21657285
21031985
20355985
19865885
19494085
19189885
18919485
18699785
18496984
)

// (结) 读一个名为"labelGroup"的labelList

Note : topoSet好像可以做labelList的操作，它的功能还很多！包括{label/zone/face/point/box/rotatedBox/cylinder/sphere/...}toCell

写

此例用runTime而没用用mesh作为objectRegistry，但需要system/controlDict

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

// IOobjectWriter.C
// List<label>

#include "Time.H"
#include "argList.H"

#include "IOList.H"    
#include "OFstream.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main( int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    IOList<scalar> IOScalarList 
    (
        IOobject
        (
            "IOSList",
            "constant",
            runTime,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        )
    );

	List<scalar> sList(8);
    sList[0] = 7.0;
    sList[1] = 9.0;
    sList[2] = 1.0;
    sList[3] = 2.1;
    sList[4] = 4.0;
    sList[5] = 7.0;
    sList[6] = 4.0;
    sList[7] = 0.0;
	IOList<scalar> IOScalarList1
	(
		IOobject
	   	(
			"IOSList1",
		    "constant",
			runTime,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
		),
		sList
	);

    Info<< "Writing " << IOScalarList.name() << " to " << IOScalarList.objectPath() << endl;

    OFstream os(IOScalarList.objectPath());   // 如果要用writeHeader，需要OFstream
    OFstream os1(IOScalarList1.objectPath());

	//writeHeader only
	IOScalarList.writeHeader(os);
	IOScalarList1.writeHeader(os1);

	//output as list                          // 没有header，就是"( )" 样式的输出
	//os << IOScalarList;
	//os1 << IOScalarList1;

	//output and header                       // 含有header，可以通过OpenFOAM再次读入，如上例volVectorField读入的时候文件必须含有header
	//IOScalarList.write();
	//IOScalarList1.write();
	
    Info<< "\nEnd\n" << endl;
}

// ************************************************************************* //


// Make/files
IOobjectWriter.C

EXE = $(FOAM_USER_APPBIN)/IOobjectWriter

// Make/options
EXE_INC = \
    -I$(LIB_SRC)/OpenFOAM/lnInclude

EXE_LIBS = \
    -lOpenFOAM

这里没有用到fvCFD.H这个容量超级大的*.H文件，也没有用到class fvMesh，仅仅是List的读入，算是个minimum code.但如果要加入读volScalarField呢？下面是上面代码加入volScalarField的变种，但如果想要保持尽量少的头文件并不容易.下面这段代码编译错误.

1
2
3
4
5
6
7
8
9
10
11
12
13

$ wmake
SOURCE=userProbeByLabel.C ;  OMPI_CXX="g++" mpicxx -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O2 -march=native -fuse-ld=bfd  -DNoRepository -ftemplate-depth-100 -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/finiteVolume/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/meshTools/lnInclude -IlnInclude -I. -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/userProbeByLabel.o
userProbeByLabel.C: In function ‘int main(int, char**)’:
userProbeByLabel.C:59:7: error: variable ‘Foam::volScalarField mean’ has initializer but incomplete type
       scalarFieldName+"_mean",
       ^
userProbeByLabel.C:76:14: error: variable ‘Foam::volScalarField scalarField’ has initializer but incomplete type
              scalarFieldName,
              ^
userProbeByLabel.C:84:18: error: variable ‘Foam::volScalarField sPrime’ has initializer but incomplete type
   volScalarField sPrime = scalarField - mean;
                  ^
make: *** [Make/linux64GccDPOpt/userProbeByLabel.o] Error 1

说明volScalarField并没有正确地被初始化，查userProbeByLabel.dep发现没有对应的模板类GeometricField的身影，但加入#include "GeometricField.H"的最终结果也是莫名其妙一堆错.当然一个肯定可行的解法是头文件加上fvCFD.H，基于测试我发现#include "fvc.H"或者#include "fvm.H"均可通过编译，而且能够正常运行.代码如下，关键行#include "fvc.H"被注释.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

#include "Time.H"
#include "argList.H"
#include "fvMesh.H"
#include "timeSelector.H"

//#include "fvc.H"  关键行

#include <fstream>

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main( int argc, char *argv[])
{
    timeSelector::addOptions();
	argList::validArgs.append("scalarFieldName");
    argList::validArgs.append("timeOfAverageField");
    argList::validArgs.append("position");

    #include "setRootCase.H"
    #include "createTime.H"

	word scalarFieldName(args.additionalArgs()[0]);
	word timeOfAverageField(args.additionalArgs()[1]);
	word position(args.additionalArgs()[2]);

	instantList timeDirs = timeSelector::select0(runTime, args);

    IOList<label> labelGroup
    (
        IOobject
        (
            "labelGroup",
            position,
            runTime,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

	//output and header
	//labelGroup.write();
	Info<< labelGroup << endl;

	mkDir("userDefinedLog");
	std::ofstream fluctuationLog
	(
	   	fileName(string("userDefinedLog")/string("fluctuation_labelGroup_"+scalarFieldName)).c_str(),
		ios_base::app
	);

    #include "createMesh.H"

	volScalarField mean
	(
	    IOobject
	    (
		    scalarFieldName+"_mean",  // 难道是 IOobject naming ???
		    timeOfAverageField,
		    mesh,
		    IOobject::MUST_READ
	    ),
		mesh
	);

	forAll(timeDirs, timeI)
	{
		runTime.setTime(timeDirs[timeI], timeI); 
	    Info<< "Time = " << runTime.timeName() << endl;

		volScalarField scalarField
		(
        	IOobject
        	(
            	scalarFieldName,
           	    runTime.timeName(),
           		mesh,
            	IOobject::MUST_READ
        	),
			mesh
		);

		volScalarField sPrime = scalarField - mean;

        fluctuationLog << runTime.timeName() << " ";
        forAll(labelGroup, i)
        {
            //fluctuationLog << labelGroup[i] << " " ;
            fluctuationLog << sPrime.internalField()[labelGroup[i]] << " " ;
        }
        fluctuationLog << std::endl;
	}
	
    Info<< "\nEnd\n" << endl;
}

// ************************************************************************* //

// Make/files
userProbeByLabel.C

EXE = $(FOAM_USER_APPBIN)/userProbeByLabel

// Make/options
EXE_INC = \
    -I$(LIB_SRC)/OpenFOAM/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lOpenFOAM \
    -lfiniteVolume \
    -lmeshTools

计算和写（标准）

1
2
3
4
5
6
7
8
9
10
11
12
13

volScalarField strainRate
(
   IOobject
    (
        "strainRate_"+velocityFieldName,
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    Foam::sqrt(2.0)*mag(symm(fvc::grad(velocityField))) // 这里就继承了量纲
);

计算和写（避开量纲）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

scalarField userCalcNu(scalarField strainRate) // 这里想要偷懒，全部无量纲
{
    scalar nuInf_ = 2e-6;
    scalar nu0_   = 3e-4;
    scalar k_     = 1;
    scalar n_     = 0.326;
    scalar a_     = 2.0;  // defaut BirdCarreau

    return
        nuInf_
      + (nu0_ - nuInf_)
        *pow(scalar(1) + pow(k_*strainRate, a_), (n_ - 1.0)/a_);
}

// 先初始化一个变量nu
volScalarField nu
(
    IOobject
    (
        outputFieldName,
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
	//userCalcNu(strainRate) Not a volScalarField ! 如何由scalarField变成volScalarField??
    mesh,
    dimensionedScalar
    (
        outputFieldName,
        dimless,                               // 续前面偷懒的地方，算个无量纲
        scalar(0.) // 很奇怪如果是vector这里可以用vector::zero，但scalar::zero就不行
    )
);

// 然后再赋值

nu.internalField()=userCalcNu(strainRate); // 如果用strainRate.internalField()会有什么区别？编译倒是可以通过

// 输出

nu.write();

计算一个跟壁面相关的场（也就是仅在边界上才有值）并输出

可以参照wallShearStress，这样输出的场（yPlus也是类似）可以在paraview可视化，也许需要去掉internal mesh的勾选

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

//初始化变量，依旧是volVectorField而不是什么boundaryField

volVectorField wallShearStress
(
    IOobject
    (
        "wallShearStress",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedVector
    (
        "wallShearStress",
        sqr(dimLength)/sqr(dimTime),
        vector::zero
    )
);

//赋值过程对着wallShearStress.boundaryField()赋值

//输出过程

wallShearStress.write();

初始化一个单位张量场

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

volTensorField identity
(
   IOobject
    (
        "identity",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedTensor // dimensioned<tensor>
    (
        "identity",
        dimless,
        tensor::I
    )
);

编写边界条件

fixedValueFvPatchVectorField

Dirichlet边界(可以随时间变化)，通常继承自fixedValueFvPatchVectorField，编写边界条件的时候要注意：

1. data member 申明(*.H)和初始化(*.C)的顺序要一致
2. 一定要要让所有constructor里面data member都完成初始化(*.C)
3. *.C里面的virtual void write(Ostream&) const务必要包含所有必要的data member，因为在decomposePar时写入processor*里边界条件用的就是这个write
4. *.C里面的makePatchTypeField似乎是个宏函数，一定要修改成makePatchTypeField(fvPatchVectorField, SameAsClassName);
5. 可以写一些non-member function用来做简单的计算

窍门
所在边界的名字：patch().name()
当前时间步的时间参数：this->db().time().timeOutputValue().在这个类runTime不可见

常用Tips

object

timeSelector : 通常以 -time 起，后接 ':500,1200:1300,3000:'

args : argList类 —— argument list. 构造函数对应 #include "setRootCase.H"

runTime : 可以用write()方法，将所有注册的object都写出来（如果不是NO_WRITE），但一般的流场也可以自己用write()，例如U.write(). 构造函数对应 #include "createTime.H"

mesh : fvMesh类，可以用来遍历：forAll(mesh.C(), celli)，也可以通过mesh.V()来找网格体积

函数

书写文件头

OpenFOAM里面文件头怎么来的？ : writeHeader

max/min

场(GeometricField，一般是某volScalarField)的最大值，最小值:（虽然我不确认并行的时候对不对）

1
2
3
4
5
6
7
8
9
10

max(someField); //这个给出的是internalField()和boundaryField()中数值上最大的值
max(someField.internalField()); // 这是网格内的体积元volume中的最大值
max(someField.boundaryField()); // 这个你会阴沟里翻船，因为这里相当于max(allBoundaryGeometricFields)，意味着这首先是“边界_i/patch_i”的list，虽然我不知道findMax会用什么方法来比较这几个list的大小，但从输出的index可以判断对应的就是边界的指标

//因此要求得boundaryField()上的最大值应当循环
forAll(mesh.boundaryMesh(), patchI)
{
	max(someField.boundaryField()[patchI]) // 这里才是对某一个边界上的所有面积元face上的值构成的list求最大值
}
//并且在求得每个边界上最大值后再从从中求最大值

findMax/findMin

这个紧接前面，要注意的是findMax返回的index是所操作数组的index而非全局index
面积元的全局index在哪里呢？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

//  通过start()来找到全局faces中，所在patch第一个face的全局index
forAll(mesh.boundaryMesh().names(), patchName) 
{
    polyPatch cPatch = mesh.boundaryMesh()[patchName]; //虽然觉得这样有点啰嗦，为啥要用name嘛
    forAll(cPatch, faceI) 
    {
        label faceId = cPatch.start() + faceI;
	    Info << faceId << " -> "  << endl;
	    Info << "  cPatch[faceI]        = " << cPatch[faceI] << endl;   // 此处输出的就是face本身，face是什么？四个点的label就确定一个面
	    Info << "  mesh.faces()[faceId] = " << mesh.faces()[faceId] << endl; // 全局输出看看是不是一样，mesh.faces()就是存所有faces的地方
    }
}

//输出
/*
...
1633 -> 
  cPatch[faceI]        = 4(733 734 755 754)
  mesh.faces()[faceId] = 4(733 734 755 754)
1634 -> 
  cPatch[faceI]        = 4(754 755 776 775)
  mesh.faces()[faceId] = 4(754 755 776 775)
...
*/

findPatchID

通过字符串来找到一个patchList里面对应的patchID

1

label patchID = mesh.boundaryMesh().findPatchID("walls");

boundary value还是离网格最近的cell value?

U.boundaryField()[patchi]这个是boundary value.U.patchInternalField()这是离boundary最近的cell些的U的cell value

faceCells

cell到face是很自然的，因为：cell是什么？多个face闭合起来成为cell。那怎么反过来找呢？

1
2
3

const fvPatchList& patches = mesh.boundary();

patches[somePatch].faceCells[facei]; // 返回值就是patches[somePatch][facei]对应的唯一cell的cellID

类

VectorSpace, Vector, Tensor: mag

VectorSpace 显然是这里最具有一般性最底层的模板类了，Vector和Tensor在它眼里仅是3个元素和9个元素的差别，求mag都一视同仁。

mag在VectorSpace/Vector里面就是所有元素平方和magSqr然后开平方，在Tensor变成9个元素平方和然后开平方。因此一般的Tensor不同重写magSqr，SymmTensor因为其形式特别就重写了magSqr，如此替换掉最一般的magSqr,对SymmTensor仍然可以进行mag操作，为什么呢？因为SymmTensor是Tensor是VectorSpace里面的元素，因此一定就有mag，此时母类mag调用的是SymmTensor版本的magSqr。

外积

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

#include "vector.H"
#include "tensor.H"
#include "IOstreams.H"

using namespace Foam;

int main()
{
/*
    Info<< vector::zero << endl
        << vector::one << endl
        << vector::dim << endl
        << vector::rank << endl;
*/

    vector a(1, 2, 3);
    Info << "a = " << a << endl;
    Info << "a*a = " << a*a << endl;

    return 0;
}

如果将#include "tensor.H"注释掉会报一个比较长的错：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

$ wmake
Making dependency list for source file Test-vector.C
SOURCE=Test-vector.C ;  OMPI_CXX="g++" mpicxx -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O2 -march=native -fuse-ld=bfd  -DNoRepository -ftemplate-depth-100  -IlnInclude -I. -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/Test-vector.o
OMPI_CXX="g++" mpicxx -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O2 -march=native -fuse-ld=bfd  -DNoRepository -ftemplate-depth-100  -IlnInclude -I. -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OSspecific/POSIX/lnInclude   -fPIC -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPOpt/Test-vector.o -L/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib \
      -lOpenFOAM -ldl   -lm -o /home/hluo/OpenFOAM/hluo-2.3.1/platforms/linux64GccDPOpt/bin/Test-vector
[hluo@zaurak userVector]$ vim Test-vector.C 
[hluo@zaurak userVector]$ wmake
Making dependency list for source file Test-vector.C
SOURCE=Test-vector.C ;  OMPI_CXX="g++" mpicxx -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O2 -march=native -fuse-ld=bfd  -DNoRepository -ftemplate-depth-100  -IlnInclude -I. -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude -I/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/Test-vector.o
Test-vector.C: In function ‘int main()’:
Test-vector.C:18:23: error: no match for ‘operator*’ (operand types are ‘Foam::vector {aka Foam::Vector<double>}’ and ‘Foam::vector {aka Foam::Vector<double>}’)
  Info << "a*a = " << a*a << endl;
                       ^
Test-vector.C:18:23: note: candidates are:
In file included from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpace.H:168:0,
                 from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/Vector.H:44,
                 from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/vector.H:39,
                 from Test-vector.C:1:
/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:552:13: note: template<class Form, class Cmpt, int nCmpt> Form Foam::operator*(Foam::scalar, const Foam::VectorSpace<Form, Cmpt, nCmpt>&)
 inline Form operator*
             ^
/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:552:13: note:   template argument deduction/substitution failed:
Test-vector.C:18:24: note:   cannot convert ‘a’ (type ‘Foam::vector {aka Foam::Vector<double>}’) to type ‘Foam::scalar {aka double}’
  Info << "a*a = " << a*a << endl;
                        ^
In file included from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpace.H:168:0,
                 from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/Vector.H:44,
                 from /home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/vector.H:39,
                 from Test-vector.C:1:
/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:565:13: note: template<class Form, class Cmpt, int nCmpt> Form Foam::operator*(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, Foam::scalar)
 inline Form operator*
             ^
/home/hluo/.local/easybuild/software/OpenFOAM/2.3.1-foss-2016a/OpenFOAM-2.3.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:565:13: note:   template argument deduction/substitution failed:
Test-vector.C:18:24: note:   cannot convert ‘a’ (type ‘Foam::vector {aka Foam::Vector<double>}’) to type ‘Foam::scalar {aka double}’
  Info << "a*a = " << a*a << endl;
                        ^
make: *** [Make/linux64GccDPOpt/Test-vector.o] Error 1

报错的意思是我们找Foam::Vector<double>与Foam::Vector<double>之间的operator*找不到，只能有后面vector前面scalar的operator*（分别在VectorSpaceI.H:552和VectorSpaceI.H:565），但因为vector,tensor和VectorSpace是由一种复杂的方式耦合在一起的，所以并没有报一个找不到俩vector的外积。但如果加上”tensor.H”一切都对了，输出结果：

1
2
3

$ Test-vector 
a = (1 2 3)
a*a = (1 2 3 2 4 6 3 6 9)

List

List and UList are different. List a(10) then a.append(something) will append to the tail of the initialized a resulting a a[11]=something

相互引用的数据

1
2
3
4
5
6
7
8
9
10

// 通过U找mesh()
const fvMesh & mesh  =  U.mesh();

// 通过runTime(Foam::Time is an objectRegistry)找db()，通过db()找phi
const surfaceScalarField& phi = runTime.db().lookupObject<surfaceScalarField>("phi"); // 能通过编译，但runTime能不能从objectRegistry找到是另一回事

// 在边界条件里面(fvPatchField)找db()
this->db()
// 再通过db()找计算time的值
const scalar t = this->db().time().timeOutputValue()

几个坑

scalarExplicitSetValue : 源fvOptions/constraints/general/explicitSetValue，居然连field不设置都可以不报错..

fixedTemperatureConstraint : 是针对能量方程的constraint，即compressible solvers

一些继承关系疏理

singlePhaseTransportModel ： nonNewtonianIcoFoam里面通过singlePhaseTransportModel来初始化流变模型，与viscosityModel没有继承关系，因此singlePhaseTransportModel不能用Foam::tmp<Foam::volScalarField> Foam::viscosityModel::strainRate()这个方法。然而BirdCarreau是viscosityModel的子类，它是如何被singlePhaseTransportModel包装起来的呢？

Pstream : 继承自UPstream(内有方法：parRun(), master(), procNo()…)，不过通常用的时候这样用Pstream::master()，用于比如输出到文件而不想因为mpi多次重复输出

Pirating

timeDir

如果小心输出时间步全部刚好差了0.01怎么办，例如[5,7]变成了[5.01,7.01]，timeDir里面的U,p都有个冠冕堂皇的location的值对应时间。但其实可以把目录5.01改成5，7.01改成7，这时比如要续算，solver便认为对应的时间步其实是5和7。似乎sample这些utility也可以这样骗[待确认]。

一句话：外面的目录名是个“货真价实”的幌子，可以欺骗solver和utility

setFields

setFields用于在internalField里面写一个值。选择区域最常用的估计是boxToCell，但定义box有格式

1
2
3

box (x1 y1 z1) (x2 y2 z2);

// 这里必须有 x1 < x2, y1 < y2, z1 < z2

如果不符合上述规范，不会有报错，但不会在internalField里面有写任何的值

这里给出一个符合规范的setFieldsDict

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

defaultFieldValues
(
    volScalarFieldValue T 0
);

regions
(
    boxToCell
    {
        box (-0.004 -0.08 -0.004) (0.004 -0.004 0.004);
        fieldValues
        (
            volScalarFieldValue T 1
        );
    }
);

probes

OpenFOAM本身的probes如果刚好碰上网格边上会有warning
属于functionObjects，在simu的同时运行，如果想要postProcssing，用execFlow…，但是：需要在constrolDict里面加入,自己写的BC会有奇怪报错…
对策：用wyldckat写的，摘干净的最通用的ExecFunctionObjects，把所有的functionObject都放在system/system/controlDict.functions，这样就可以分开了

reconstructPar与decomposePar

这俩不互为反函数，如果要用这种方式完成binary和ascii的转换，一定double check BC，因为BC可能被改写

一句话：reconstructPar/decomposePar 一定要看好BC

reconstructed case or decomposed case

paraview中有这个选项，然而，经历reconstructPar或者decomposePar后BC可能被篡改，所以单纯看边界的话有可能会被吓倒在地，怎么可能！同一个时间步decompose情况下用paraview看BC上的值好好的，结果reconstructPar之后paraview一看……面目全非！这种差异在numberOfSubdomains更容易出现(例如120)，而设置成4就有可能没问题，可能可以算作bug。当然也许是我写出来的bug，毕竟有的边界value要生成随机数，我随机种子又没变。

一句话：当你看到reconstruct之后的BC不对的时候，即使processor*已经被清空了，此时也不用慌，再一次decomposePar帮你恢复BC（至于numberOfSubdomains设置为多少，最好跟并行算的时候一样吧）。即使又出了什么幺蛾子BC的值不能恢复，至少internalField的值不会被动到

上图，这个就是我并行120个核算出来之后想都没想就reconstructPar的效果，一脸懵

但把原数据也就是放在processor*里面的数据用paraview来看的时候就变成了下图，用numberOfSubdomains=120再重新decomposePar也是一样的效果

同一个类的不同变种编译为不同名字的lib

如上，如果想要将多个lib混用在一个case里面，在system/controlDict里面都加入lib是必须的，但其实这样不行，做不到混用。如何做到呢？需要修改类的名字，并与makePatchTypeField(fvPatchVectorField, pVFvPatchVectorField2Dpf_Port1)这个macro function对应，这样不仅lib名字不同其实里面的类也不同，这样就可以混用了。

cyclic

U,p设置周期性条件后，在生成的数据里type为cyclic你会发现U,p都没有value，仅有一个type (phi却有，也是cyclic类型).那为啥我用patchIntegrate U inlet又不是零呢？经过仔细地看源码和比对

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

// 读入U，找到cyclic对应的patchLabel

Info << U.boundaryField()[patchLabel] << endl;
/*
你猜输出什么？
type            cyclic;
仅有type！没有value？不过，读入的就是这样的数据(volVectorField.boundaryField()[cyclic patch])，输出这样的数据也在情理之中
*/

vectorField& Ub_cyclic = U.boundaryField()[patchLabel]
Info << Ub_cyclic << endl;
/*
看输出...就有了，不过vectorField有type，这个似乎也不是特别有道理，毕竟vectorField仅仅是List<vector>
type            cyclic;


16900
(
...
...
...
)
*/

// 参照patchIntegrate.C，里面没有输出整个List而仅仅输出sum；于是我在同样的patch上面求sum输出，计算结果是一样的
Info << sum(U.boundaryField()[patchLabel]) << endl;
Info << sum(Ub_cyclic) << endl;
/*
(11560.70065 80.71560542 3.732255509)
(11560.70065 80.71560542 3.732255509)
*/

// 所以可能是GeometricBoundaryField类对cyclic type输出的问题？？

这里得出的结论是：U和p在数据文件里面，boundaryField项的其中一个patch上有可能仅有type（例如cyclic），但没有value，这不意味着U.boundaryField()[cyclic patch]上面没有值。简单地来说：当你看数据文件里面边界上没有值，并不是真的没有值，边界上还可以算通量呢！
PS : phi同为cyclic类型，但边界上却有值，且value写在数据文件里

如何做出成体系好debug的后处理tool chain？

答：先用bash把我的user*系列utility连起来，通过文件夹好好组织；后续的python后处理一定要主程序条理清楚，加上合适的输出，去掉冗余输出，有exception control；一个类别的后处理，尽量用同一个python module去实现通用。往具体的说：先做好牛顿流动的case，在一个个case中成长，在痛苦中净化和进化。

如何组织DNS数据？

答：

1
2
3
4
5

|- raw data      # 原始数据（用于备份） + 文件夹良好组织的小型后处理文件
|- visu data     # 特征时间步的数据（可能比raws小不了多少），用于生成大型的后处理文件(场），paraview save state 等等
|- image         # 把图片分开放
|- pythonScript  # python 都放这里

海量数据的存放？

答：
zaurak :
+++ 固态硬盘IO速度快效率高，包括查看文件大小也是速度奇快
/store/
– 容积小[2T]
/home/
– 容积小[170G]

srv-data-1 :
没有空间了，基本弃用

visu/newton :
++ 容积大[4.8T / 66T]
++ 用于读写各种操作的内存很足，当zaurak内存不够用则考虑开始启用
– 硬盘IO太慢
/store/lmfa/fct/
++ 硬盘存储空间足
++ 跟zaurak之间传输速度（待估计）比较快
$HOME
++ 有备份，可放重要数据
– 空间太小[1T / 2T]

occigen :
– $STORE[1.5T]只支持*.tar.gz，只能暂时存储
– $SCRATCH[4T && fileNb 600000]
$HOME基本不用

推公式哪本书好？

答：
Turbulence en mecanique des fluides - P. Chassaing 这本里面的公式灰常好，还带有很多物理解释

远程可以进行的操作

1. 提交算例
2. zaurak, occigen, newton上都有paraview，可以用script来做一些简单的可视化

bash

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

# kill一系列进程，ps 可能搞出很多其他不相干的进程号，一定要注意，而且前三行一定不是，怎样通过awk滤过还不清楚
$ ps | awk '{print $1}' | xargs kill


# check for all jobs including "nohup" (at the same login noed)
$ ps -eaf | grep $USER


# check for graphic card status
# both Xorg and gnome-shell stores for cache (computer has already run for many days). After rebooting, there will be much less.
$ nvidia-smi


# screen
# name a screen session
$ screen -S someName
# detache from screen
Ctrl + A then D
# see running screen
$ screen -ls
# re-attach to a screen
$ screen -r ** 

# rsync

# --dry-run
$ rsync --dry-run -av --progress sourcefolder /destinationfolder --exclude thefoldertoexclude
# show progress
$ rsync -av --progress sourcefolder /destinationfolder --exclude thefoldertoexclude
# rsync multiple files at one line, passwd input only once
$ rsync -av remote:dirToCase/{0.165,0.18,0.195,0.21,0.225} /some/local/dir/


# mkdir / rm / regular expressions / using {} for multiple-file usage

# create multiple dirs
$ mkdir processor{0..11}
# rm directories under processor0, processor1, processor2, processor3 which is named 0.1
$ rm -r processor{0..3}/0.1
# rm directories under processor* named '2.8*' '2.9*' using "regular expression"
$ rm -r processor*/2.[89]*
# this is cool, delete files in any directory named as XEq*-cellPoint_p.xy or XEq*-cellPoint_U.xy
$ rm */XEq*-cellPoint_{p,U}.xy   


# sort

$ sort file.csv -n -t "," -k 22 -k 23 -o output 
# sort file *.csv by colon 22 then by colon 23


# find

# counts for the number of directories including it self "." This also counts the hidden directories prefixed by '.'. This is actually a count for "newlines" because wc -l counts for newlines and -print produces output separated by newlines. Changing -print to -print0 which will separates output by null, then there will be no newline and wc -l will always give 0 as result.
$ find . -maxdepth 1 -type d -print | wc -l 

# find directories in . and grep the str (**fast**)
$ find . -type d | grep "2" 
# (**slower**) but only will not produce results of sub directory of '2' as the previous one did.
$ find . -type d -name '2' -print  # dry run 
$ find . -type d -name '2' -delete # this wont work !! instead : xxx | xargs rm -r 是有效的

# find file matching pattern (dry run) then delete them.
$ find . -type f -name 'plotoverline.*.csv' -print 
$ find . -type f -name 'plotoverline.*.csv' -delete 
# find in somedir all subdirectories and chmod to defaut 775
$ find somedir -type d -exec chmod 775 {} \;  


# chmod

# this will change subdirectories and files ! Which is normally not GOOD. You dont want to have files all green/executable.
chmod -R 755 somedir  


# file

$ file *.so        shared object
$ file executable  dynamically/statically linked
$ file *.a         ar archive

vim

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# search 
:/some_string\c    # case insensitive
:/some_string\C    # case sensitive

# git rid of ^M
# 有时候用vim打开文件会发现"^M"，如何去掉？
$ dos2unix fileName

# 多行选中
ctrl + v (visual block) -> 
use cursor (h j k l) to move and select a multi-line-block -> 
s (insert) -> 
input "#" -> 
Esc -> 
all select block is then replaced by "#"

# 1到10行，有caption关键词的行首加上#
:1,10s/.*caption/#\0/

# encoding (other than utf8)
:set fileencoding=iso-8859-1

# paste without autoindent
:set paste

# add cursorline
:set cursorline

# reload file
:e
:edit

# 遇到难缠的(可能是硬的tab，也就是当下如果vim配置是以空格来做tab，也能通过/\t，搜索得到，且vim光标摁'j'移动不到行首，会卡住在这些tab的位置)
# IndentationError: unindent does not match any outer indentation level
:retab

git

1
2
3
4
5
6
7
8

# 密码缓存300秒

$ git config credential.helper 'cache --timeout=300'

# 不再输入密码

$ git config --global credential.helper store

latex

建议做个makefile方便操作（见下）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

new=actualVersion
old=olderVersion

all:
    pdflatex ${new}.tex
    bibtex ${new}
    pdflatex ${new}.tex
    pdflatex ${new}.tex

old:
    pdflatex ${old}.tex
    bibtex ${old}
    pdflatex ${old}.tex
    pdflatex ${old}.tex

bibs:
    bibtex ${new}
    bibtex ${old}

diff:
    latexdiff ${old}.tex ${new}.tex > diff.tex
    latexdiff ${old}.bbl ${new}.bbl > diff.bbl
    pdflatex diff.tex
    pdflatex diff.tex

clean:
    rm *.bbl *.aux *.blg *.log

完整编译${new}.tex：make
完整编译${old}.tex：make old
用latexdiff完成对整篇new与old版本的对校（可以包括reference也可以不包括）: make diff

用到reference对校功能的时候要注意：

1. 比较的是*.bbl而不是*.bib
2. 对校reference的时候，用make bibs ; make diff
3. 对校成功的话，old为红色，new为蓝色
4. 改new.bib的时候要注意尽量不要增加新的field：例如old.bib里面没有pages,volume，如果在new.bib中加入，有可能让latexdiff old.bbl new.bbl > diff.bbl prompt，可以通过enter来跳过这一系列warning但最终得到的diff.pdf在reference部分会在增加field的当前项开始后面全篇都会被蓝色override干扰校对。经测试：增减1~2个field似乎在允许范围内
5. 让latexdiff old.bbl new.bbl > diff.bbl prompt有可能是Stop和NoStop，这种warning可以忽略

beamer

overlays :

1
2
3
4
5
6
7
8
9

\only<1>{}

\onlyenv{}

\begin{overprint}
    \onslide<1>
\end{overprint}

%...

检查表制度

我开始认识到“检查表”(check-list)是源于傅盛的文章，他提到让工作变得可靠高效，检查表是个保障时效的好东西，看似机械但是可靠。随后我在《黑匣子思维》里面读到，进一步了解到纽约哈德孙河上迫降成功的背后（美当局还是对机组进行了详细地调查）起到积极作用的有无数次飞行的经验，过硬的心理素质，还有副机长分担责任在机长思考的同时协力沟通，副机长在短时间内把空客给的双发停车的check-list走完。检查表制度已经完全融入了这些航空工业优秀从业者的工作中，如作者所说：人的因素可能是正面的也可能是负面的，然而可靠性的保证，需要仰仗人训练自己“黑匣子”一般的思维。回到傅盛，自己定下的check-list就一定要严格执行，需要像”黑匣子“一样完全记录，要不然你将看不到当前check-list的缺陷，改进也就无从谈起，低级和已知的错误还是会犯。check-list的存在意义在于，当你头脑不是特别清醒，状态不好的时候，也能给你一个完整严谨的流程，最大程度规避已知的一系列“低级错误”，不单纯依靠记忆和头脑，复杂系统尤为需要如此作业。

为什么我要用在OpenFOAM上?

同fluent或其他CFD软件，OpenFOAM可能更称得上“复杂系统”，初学者容易犯错：不像fluent通过按钮和已经编程的GUI来辅助构建组织求解一个问题的逻辑链，OpenFOAM从边界条件、初始条件、网格文件、流体物性、离散格式、求解器（U求解器，p求解器）还有utility到用到的例如sampleDict，changeDictionaryDict，refineMeshDict，topoSetDict等，里面的逻辑关系通过一系列的文本文件来完成。信息的碎片化很可怕，想想你需要用vim去不同的目录打开这一系列的文本文件来检查，这样的过程费事、费力且艰辛，更重要的是极容易犯错。为什么呢？

1. 文本文件有特定格式
2. 打错字
3. OpenFOAM在tutorial里面给的范本参照并不是MVP（Minimum Viable Product），也就是说通常都有不必要的部分，初学者者并不清楚哪些信息必要，哪些重要，哪些完全没用

总之细节很多，细节一旦繁多而又难以从中抽离就难有大局观。而check-list就是提醒自己“从哪里出发”，“现在在哪里”，“还要去往哪里”。举例：像飞机起飞的check-list一样，一个check-list可以是对于某算例从无到有的按步骤一步步的配置。“从哪里出发”包括最基本的0，constant和system的文本文件结构，包括边界条件、初始条件，那么“项目的出发点”就是自己剔除上面提到的不必要部分后得到的MVP case，基于MVP case的基础可以做加密网格或者其他的变式。顾名思义，check-list是地图，“现在在哪里”有check-list指明当前步骤，作为机长，“要去往哪里，飞行姿态如何保持”才是自己心里应该操心的核心的事，换句话说check-list在这里担当了类似副机长的辅助职责。

记录“已知在某个步骤上会出现的问题”，然后据此设计check-list来规避已知问题，这是高于平时“工作的自己”层次的思维方式。工作的时候我要debug，要进入多层嵌套的代码，是细节的把控，而check-list让我跳出来，让我看到自己在地图的哪个地方，在某某步骤上的经验和教训是什么，接下来可能的问题是什么，问自己是主要矛盾吗，时间成本会是多少，应该如何决策，这是做决策的思维层次(《原则》)。因为没有一个流程或者说没有包装好的GUI会告诉你“现在到了哪一步，下一步需要做什么”。这里多说一句，GUI同check-list，也仅仅是辅助驾驶的功能，目前自动驾驶系统也没聪明到跟数千小时飞行经验的飞行员一样，能根据经验快速反应做出正确的决策。

总结一下：相较于fluent，OpenFOAM的劣势是文本文件之间的自然逻辑没有GUI来完成，优势是对一个功能如果充分了解，就可以有很大自由度进行操作，可以按照上面说到的check-list来完成步骤，实现步骤的“自动化”，因而劣势可以弥补; OpenFOAM最大的优势还是完全开源，只有不懂的地方才是黑箱，经手的所有步骤都在掌握，所有的变量都在掌握。

个性化检查表以配合工作需求

要知道，检查表可繁可简，为节约时间成本阿波罗13号出事之后检查表都被大规模简化，也就是“跳步”，而明白跳步的风险是简化的必须。所以检查表到底应该多繁多简，视情况而定。我根据自己的经验制作了一系列检查表，有的是配置算例流程，有的是关于某一个utility，还有NotToDo-list，总之对我怎么有用怎么来，还不止OpenFOAM，因为尤其跟CAE相关的软件都很大型，按钮特别多，记忆是不可靠的。加上如果在不同的计算机群上面计算，队列是不一样的，对用户的要求是不一样的（比如occigen上面对文件个数要求特别苛刻），用户的需求是不一样的（比如考虑到priority的问题，我又时候可能会选择提交一个计算时间较短的算例），情形很多，未必都需要检查表化，但按照自己的逻辑记录下要点，比起每次决策的时候再去网站上找零散的信息来得更高效。

环境

都说python是很棒的语言，入门学习曲线可以说是相当亲民，但各种环境的“崎岖”让人望而却步，我到底应该用pip呢还是anaconda这样的问题让人挠头，见图便知。

我的需求

easybuild依赖于python，简单介绍一下用easybuild：它属于Modules上一级的工具，可以（自动搜寻依赖关系）安装各种版本的python，安装相应的python库，配置相应的清晰干净的python环境，尤其在多用户高性能计算机群上广泛使用。我需要管理手头多个版本的OpenFOAM和与它们依赖的库（其实主要用到的也就mpi和scotch）。用户体验：当你的工作环境很单一不用Modules没问题，不过一旦复杂度上去了，用Modules工具可以做到事半功倍，”模块化的环境”不仅让工作环境清爽整洁、条理清楚，还可以不费劲地将自己熟悉的整个工作环境在另一台机器上复现出来。我手头与python依赖的环境有：

1. Linux CentOS 系统自带 /usr/bin/python
2. easybuild (基于1的安装)
3. /usr/bin/ipython
4. /usr/bin/pip
5. anaconda

我要使用easybuild安装的module的时候用这个alias：

1
2
3
4

alias Easybuild='EASYBUILD_PREFIX=$HOME/.local/easybuild; \
                 module use $EASYBUILD_PREFIX/modules/all; \
                 module load EasyBuild; \
                 EASYBUILD_MODULES_TOOL=EnvironmentModules'

根据我的设置，默认不加载easybuild的环境。

1
2
3
4
5

$ echo $PYTHONPATH

$ Easybuild
$ echo $PYTHONPATH
/home/hluo/.local/easybuild/software/EasyBuild/3.2.1/lib/python2.7/site-packages

看到上图里面右上角那个至少有三个箭头朝外的环境变量了吗？就是它。

而anaconda有它自己一套的环境变量:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

$ Anaconda                  # 加载anaconda环境
$ which python
~/bin/anaconda2/bin/python  # 此处有python
$ which ipython
~/bin/anaconda2/bin/ipython # 此处有ipython
$ which pip
~/bin/anaconda2/bin/pip     # 此处有pip，是不是开始复杂起来了。。

$ ipython
In [1]: import sphinx

In [2]: sphinx.__path__
Out[2]: ['/home/hluo/bin/anaconda2/lib/python2.7/site-packages/Sphinx-1.4.1-py2.7.egg/sphinx']     # 把库放在了$HOME/bin/anaconda2
In [3]: import numpy as np

In [4]: np.__path__
Out[4]: ['/home/hluo/.local/lib/python2.7/site-packages/numpy']

为啥两个库存放的地方不同？按理说这是pip安装的目录(见下文)，我肯定之前干了什么然后失忆了，促成了此“姻缘”。。。总之我在尝试安装h5py的时候遇到了各种问题，最终因为盘根错节根本厘不清！想搞清楚？重装吧。其实每次重新安装都有“做得更好”的可能，前提是厘清楚，然后养成好的习惯在使用的时候用哪一个就哪一个，切忌同时使用pip和anaconda。

一个程序就做一件事情，然后把它做好，身边的计算机牛人们说过的话，受用了。

pip

实验室哥们推荐我单用pip，他给我的使用建议：
安装scipy，pip install --user scipy，所安装的包就在$HOME/.local/lib/python2.7/site-packages
更新scipy，pip install scipy --user -U，所安装的包就在$HOME/.local/lib/python2.7/site-packages
移除所安装的包直接删除对应位置的文件即可，如果想用另一个版本直接将当前的包删除即可，这样干净整洁不容易出错，此例中rm -rf $HOME/.local/lib/python2.7/site-packages/scipy*
查看包的版本

1
2
3
4
5
6

$ pip install scipy==
Collecting scipy==
  Could not find a version that satisfies the requirement scipy== (from versions: 0.8.0, 0.9.0, 0.10.0, 0.10.1, 0.11.0, 0.12.0, 0.12.1, 0.13.0, 0.13.1, 0.13.2, 0.13.3, 0.14.0, 0.14.1, 0.15.0, 0.15.1, 0.16.0, 0.16.1, 0.17.0, 0.17.1, 0.18.0rc2, 0.18.0, 0.18.1, 0.19.0, 0.19.1, 1.0.0b1, 1.0.0rc1, 1.0.0rc2, 1.0.0, 1.0.1, 1.1.0rc1, 1.1.0)
No matching distribution found for scipy==
You are using pip version 9.0.3, however version 10.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

如果是python3的库，那么就改用pip3，库函数就会放在$HOME/.local/lib/python3.4/site-packages（这里是3.4）

安装特例

需要指出的是一般python包按照以上傻瓜安装就行，并不需要操心这么多，这里仅为一个特例。

h5py的安装需要用到devel包（for developper 或者叫 no-binary 包），因为基于h5py的特性需要基于已有的mpi和hdf5重新编译，比如安装的时候报错Python.H not found，那就说明python包里面该有*.H的地方没有源文件，也许安装的不是开发者安装包（仅仅是pre-compiled后的库函数，类似于apt-get install安装时从源上拷贝到本地的packages）；当然另一种可能的情况是给的头文件寻找路径不对。

小结一下，如果报错找不到mpi.h或者hdf5.h，那么就是环境中找不到mpi或hdf5包里的头文件，也对应前面两种可能：包里面不包含*.H或者搜索路径不对。这种时候，意味着h5py需要在本地重新编译。

这里说起来有点绕，这是因为h5py是python对hdf5库的一个wrapper，也就是基于用C或者C++的hdf5库的接口重新包装的结果，因而需要保证在当前环境中hdf5库已经配置完整（包括CPATH指定头文件的路径，LD_LIBRARY_PATH指定*.so库的搜索路径），别忘了hdf5依赖于mpi，那么mpi也应该完成配置。这样才能保证重新编译的基本要求，在python里面no-binary安装对应：

1
2

$ pip install --user h5py # 预编译
$ pip install --no-binary=h5py h5py --user # no-binary 版本

anaconda

anaconda也可以安装在自己$HOME下且包很全，有个conda命令来管理包的历史、更新、安装、卸载。我在初期就用anaconda，试图用conda安装h5py失败了（安装h5py的动机源自一篇关于OpenFOAM湍流进口条件的学士论文，可惜这个人写的库只适用于structured rectilinear grid，也就是进口得是方形，圆管就不行了。对于我来说，直接用用不上，但eddylicious库的编写里得更多的是面向过程，加上只有一个轴的方向y是边界层，蛮好读懂；另一方面他改写了timeVaryingMappedFixedValue，把输出格式改成了HDF5，按照要求编译：要注意Make/options里面通过HDF5_DIR而不是CPATH和LD_LIBRARY_PATH来找，且在环境中设置不够，需要写入到Make/options里面HDF5_DIR=$HOME/.local/easybuild/software/HDF5/1.8.17-foss-2016a，这里的HDF5用easybuild安装，跟使用的OpenFOAM/2.3.1-foss-2016a匹配），就h5py来说，pip网上的资料还是多些。

IDE

我在windows和linux下也都有anaconda，主要用在包里的spyder，可以实时查看数组变量（函数内部的变量在IDE里面查看不了，因为没有暴露在当前执行脚本的环境中），这在python学习初期很有用，在初期学习中可以有意识地不写函数，这样所有的变量都在可追溯的Variable explorer里面便于实时查看。想要查看某变量的类型，在里面集成的ipython console通过type(variable)输出即可。
我选spyder的主要原因：

1. 用matplotlib画图的时候就在ipython console即时输出图片，避免另弹出画图的新窗口+手动关闭窗口，看到图片就可以调整legend的位置还有字体大小什么的，很直观；如果要拷贝图片另存为png格式直接右键save image as就可以，所见即所得（为何这样说？如果legend在图片box外面，savefig不会将legend包括进去，也就说savefig最终不会包含box外的legend）。
2. 库函数自动补全
3. 注释和去注释方便，缩进和去缩进方便，高亮易辨识
4. 鼠标点击函数，在不同的库之间跳转

福利: spyder -> File -> Print to File (PDF) 有语法高亮

局限: 上面1提到的legend在外的情况可能会有点麻烦；在3d plot里面有些图（比如3D云图）是可以鼠标拖着旋转的，但在spyder里面就是个静止的图（可能原因在于就其本源可能还不如货真价实的ipython console）；因此animation当然也在spyder里面看不了，这个时候就都得回到无IDE下的python解释器环境

基本数据类型

Sequence Types — str, unicode, list, tuple, bytearray, buffer, xrange

1. str : string, written in single or double quotes
2. unicode : unicode string, much like string
3. list : constructed by square brackets, its element seperated by comma
4. tuple : constructed by the comma operator (not within square brackets), with or without enclosing parentheses, but an empty tuple must have the enclosing parentheses, such as a, b, c or (). A single item tuple must have a trailing comma, such as (d,).
5. bytearray : created with the built-in function bytearray()
6. buffer : Buffer objects are not directly supported by Python syntax, but can be created by calling the built-in function buffer()
7. xrange : similar to buffers in that there is no specific syntax to create them, but they are created using the xrange() function

Operators on list

1
2
3
4
5
6

>>> myList = ["The", "earth", "revolves", "around", "sun"] # initialize a list with double quotes
>>> myList
 ['The', 'earth', 'revolves', 'around', 'sun'] # output are single quoted
>>> myList = myList + ["sure"]  # need a list to perform operator+ on list objects
>>> myList
 ['The', 'earth', 'revolves', 'around', 'sun', 'sure'] # list + list = list

Sequence Types – dict
A mapping object maps hashable values to arbitrary objects. Mappings are mutable objects (in c++ there are also mutable objects, what are they?). There is currently only one standard mapping type, the dictionary.
Dictionaries can be created by placing a comma-separated list of key: value pairs within braces, for example: {‘jack’: 4098, ‘sjoerd’: 4127} or {4098: ‘jack’, 4127: ‘sjoerd’}, or by the dict constructor.

以上是抄写，这里是个人使用心得：

1. list里面可以放不同类型的变量，很有用；

2. tuple可以作为函数return多个值的时候的选项例如fig, ax = plt.subplots()，当然dictionary也可以；

3. 我在OpenFOAM后处理脚本里面用到dict，因为一个dictionary就可以包含算例的绝对路径、后处理时间区间，后处理数据shape…更重要的是在程序读取并处理数据时，主程序里暴露出从key到value的过程，比如如果我画平均值就取名为mean，如果画均方根就取名为rms，这样后处理程序可读性会大大提高；numpy里面array当然用得很多，list可以转化为array

库的使用

python强大的地方当然是有各种很棒的库的支持(numpy, scipy, matplotlib等)，使用库用关键字import就行。如何查询所用库的路径？进入ipython

1
2
3
4
5
6
7
8
9
10

In [1]: import numpy as np
In [2]: np.__                   # 自动补全
np.__
np.__NUMPY_SETUP__   np.__delattr__       np.__getattribute__  np.__new__           np.__repr__          np.__version__
np.__all__           np.__dict__          np.__git_revision__  np.__package__       np.__setattr__       
np.__builtins__      np.__doc__           np.__hash__          np.__path__          np.__sizeof__        
np.__class__         np.__file__          np.__init__          np.__reduce__        np.__str__           
np.__config__        np.__format__        np.__name__          np.__reduce_ex__     np.__subclasshook__
In [2]: np.__path__
Out[2]: ['/home/hluo/.local/lib/python2.7/site-packages/numpy'] # 可以看到输出的是一个list，里面路径变量是str

自定义库

自定义库也很简单，写了一个脚本叫Dai_thesis.py，里面有个函数叫Fig4p8a，那么在主程序（同一目录下）只需import Dai_thesis，然后就可以Dai_thesis.Fig4p8a()，在执行主程序的时候Dai_thesis.pyc将会自动生成，这是个二进制文件，表明Dai_thesis.py以库的形式被使用，这样就使用了Dai_theis这篇论文里面Fig4.8(a)里面的数据。当然库的使用也有嵌套，例如在上述基础上加一个层级:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

# 高层级的库 reference_database.py ： 我画图的时候作为对照的的数据来源汇总：杂志文章、博士论文、解析表达式

import Niewstadt_article_1995

import Dai_thesis

import analytic_functions

import Eggels_thesis

import Eggels_article_1994



# 低一个层级的 Dai_thesis.py : 功用是去特定的地方读特定的文件，返回numpy array组成的一个tuple

import numpy as np

def Fig4p8a(fluid):
    string='/home/hluo/Pictures/Dai_T/meanProfile_velocity/XEqM4_debitMin/profile_XEqM4_debitMin.csv'
    data=np.genfromtxt(string,skip_header=1,delimiter=',')
    
    switcher={
        'EAU':1,  # 选择流体类型
        'PAA':2,
        'XG':3
    }
    
    x = data[:,0]
    y = data[:,switcher[fluid]]
    
    return x, y

def Fig4p9a(fluid):
    string='/home/hluo/Pictures/Dai_T/meanProfile_velocity/XEqP12_debitMin/profile_XEqP12_debitMin.csv'
    data=np.genfromtxt(string,skip_header=1,delimiter=',')
    
    switcher={
        'EAU':1,
        'PAA':2,
        'XG':3
    }
    
    x = data[:,0]
    y = data[:,switcher[fluid]]
    
    return x, y
   


# parameters_T_RES1d_MethodMapped_subMethod_NearestFace.py ： 用dict存储一个OpenFOAM算例后处理里的复合数据类型

# physical parameters

physics={
        'R':0.004,
        'nu':1.0e-6,
        'uTau':0.0473
        }

# output of sample utility

sampling={
        'raw_sample_size':160,
#        'dataShape':(199,4)     # uniform 
#        'dataShape':(188,4)     # face 3 
        'dataShape':(195,4)     # face 2     
#        'dataShape':(195,4)     # face 1    

# 这里展示了这个后处理过程的复杂性：dataShape有多个值
# 因为OpenFOAM sample sets操作后的数组用不同mode输出的数据长度不同
# 例如: uniform 和 face 两个mode下输出的长度不同，在同一mode下不同的地方进行sample操作也可能得(对于
# 复杂的几何外形的算例，网格在空间不均匀，换个sample的位置，几乎一定会变)的不同长度的数据，例如这里的
# face 1,2,3

# 为何要确定长度？因为如果做空间统计（这个例子是时间，只要网格不变，sets不变那数据长度就不变，因为取
# 样的labelList不变），比如在周期条件的圆管流动里，里面沿着半径方向取一系列数学上有对称性的sets，然后
# 想求对于sets的平均（即空间平均），每个sets在这里对应一个数据文件。这时候即使是相当规则的网格（butterfly
# 分成5块），我的测试结果是:这一系列看起来数学上完全对称的sets有时候也会出现个别的特例（要么输出数据
# 长度为0，要么输出数据跟大部分的shape不同，OpenFOAM当然不会报错警告你，所以我自己建立了一个检查机制，
# 这样呢就需要一个目标shape和读取实际数据shape的对照的过程，这样每次后处理我就知道是否有exception。

# 于是在这里，通过这样一个dictionary完成了参数parameters['sampling']['dataShape']的传递）
        }

# data entry parameters

dataEntry={
        'startTime':5.5,  # KinecticEnergy stationary
        'endTime':7.2, 'chunkStep':50,
        'NbOfFiles':171,
        'path':"/store/caseByGeometry/T/new-mesh/pointwise/postProcessing/1d_mapped_NearestFace",
        }

parameters={
        'physics':physics,
        'sampling':sampling,
        'dataEntry':dataEntry,
        }