draft/en/note/elsa.org

   1 * <lisp>(setq Title "elsa开发手册")</lisp>
   2 <contents depth="4">
   3
   4 * 1 空间离散后为ODE,采用(伪)unsteady时间积分求解器,进入(伪)时间循环
   5 ** 1.1在这个循环内部
   6  - flux和source计算
   7  - bc计算
   8  - 附属参数计算,如p,黏性
   9  - 时间步长计算,收敛加速计算
  10
  11  稳态计算 || 时间循环直到收敛
  12  瞬态计算 || 循环到指定的时间步结束
  13
  14 * 2 实现
  15  - elsA使用UML建模
  16  - elsA kernel 采用C++的OOP设计,只有很费cpu时间的计算采用Fortran,并且在不影响OO设计的情况下
  17
  18 * 3 输入数据
  19  - 几何数据,mesh 坐标
  20  - 拓扑结构数据, block的连接方式
  21  - 物理数据,初场
  22  - BC
  23
  24 * 4 用户界面
  25  - PyGelsA,为GUI界面
  26  - Python脚本,复杂仿真的首选,用法在elsA用户手册上
  27   elsA可以被看为一个标准的Python模块,elsA.py . 将c++代码转换的Python 工具为swig.
  28
  29 * 5 并行
  30  - 采用MPI库,使用coarse-grained并行策略.
  31
  32 * 6 多目标耦合
  33  - 气体固体
  34
  35 * 7 优化
  36  - Opt模块实现了离散Adjoint方法.用来嵌入气动外形自动优化过程.
  37
  38 * 8 后处理
  39  -
  40
  41 * 9 Kernel 设计
  42 ** 9.1 职责分类
  43  - 流体物理属性
  44  - 构建数值仿真的目标空间
  45  - 构建方程组:计算空间离散项(flux,source terms)和边界条件
  46  - 控制时间递进
  47
  48 ** 9.2 layers
  49  elsA kernel包括400个class,分配到26个modules.这些模块又归于6个layer,上层只依赖下层,单向职责.
  50   - Base
  51   - Geometry
  52   - Physical model
  53   - Space Discretization
  54   - Solver
  55   - Factory(顶层)
  56
  57 *** 9.2.1 Base由以下模块组成
  58  - Fld: 数据存储类,封装了存储任何计算数据所需要的动态数据管理,提供了monitoring 方法,优化了array-like的语义(类似fortran 90语法),并提供于fortran的交互方法.Fld提供最有效率的cpu和内存.STL的容器没法与其效率相比.
  59  - Pcm: 并行实现,封装了消息传递界面(目前只有MPI)
  60
  61 *** 9.2.2 Geometry Layer
  62  - Blk: 定义了block概念,block为一个定义了网格的离散了的物理区域,包括边界条件和初始数据.考虑了grid motion,ALE,chimera和HMR.
  63  - Geo: 定义了抽象的计算网格,提供了所有的FV需要的geometry要素.
  64    - Metris: cell的体积,surface of cell interface
  65    - Topological relations: cells, interfaces和nodes分别的关系.ghost cells 减少了计算循环的indirections.提高了cpu效率
  66  - Dtw: 保存说有的distance和边界层积分厚度
  67  - Mask: 定义用在Chimera 技术里的概念
  68  - Join: 处理multi-block计算.可以为1-1连接,1-n连接或者mismatched.
  69
  70 *** 9.2.3 Physical model层
  71  - Eos: 计算p,t,层流年星系数等变量
  72  - Tur: 处理湍流模型和转淚预测
  73
  74 *** 9.2.4 空间离散 layer
  75  - Oper: 计算一个单独的CFD物理量
  76        - Fxc: convective fluxes
  77        - Fxd: diffusive fluxes
  78        - Sou: source terms
  79  - Bnd: 处理边界条件
  80
  81 *** 9.2.5 Solve layer, kernel最复杂的部分,
  82  - Rhs: 计算方程组的Right Hand side
  83  - Lhs: 汇总implicit方法: 每种方法都要从线性化方程组中构建和转折矩阵
  84  - Tmo: 时间积分,管理＝总时间循环
  85
  86 *** 9.2.6 Factory (elsA 顶层layer)
  87  - 负责动态创建说有的kernel对象,从用户的输入构建说有的相应对象
  88
  89 ** 10 部件详述
  90 *** 10.1 FLD 部件
  91 **** 10.1.1基本的数值容器
  92   Fields是elsA最基本的操作对象,作为real,integer,boolean的容器.用后缀表示版本,F代表Float,I,B.
  93  - FldArray: 保存数值,没有任何的位置信息
  94  - FldField: 保存在grid上的数值信息.有FldNode,FldCell,FldInt三种.可以用如下详细指定
  95 <source lang="c">
  96   typedef FldFieldF FldCellF
  97 </source>
  98
  99 **** 10.1.2 接口
 100   fileld定义为2D结构.
 101  - 一个为节点的个数 0 ~ _size-1
 102  - 一个为节点内保存的变量个数 1 ~ _nfld
 103
 104 *****    10.1.2.1. Construction of a field which stores the unknowns of the CFD problem (ro, rou, rov, row, roE):
 105
 106        E_Int nfld = 5;
 107        FldCellF wCons(ncell, nfld);
 108        Construction of a field which stores fluxes:
 109        FldIntF flux(3*ncell, nfld);
 110
 111 *****    10.1.2.2. Construction of a field which stores mesh coordinates:
 112        FldNodeF x(ncell, 3);
 113        FldNodeF y(ncell, 3);
 114        FldNodeF z(ncell, 3);
 115
 116 *****    10.1.2.3. FldArray or FldField can be used to store values without geometric links, such as:
 117        FldArrayF TurKO::getModConst() const
 118        {
 119           FldArrayF modConst (7);
 120           modConst[0] = _kappa;
 121           modConst[1] = _sigma1;
 122           modConst[2] = _sigmae1;
 123           modConst[3] = _beta1;
 124           modConst[4] = _wsig1;
 125           modConst[5] = _betae;
 126           modConst[6] = _Sr;
 127        }
 128
 129
 130 *****    10.1.2.4. To access individual elements, a syntax similar to Fortran is used:
 131        FldArrayF f(100,2);
 132        f(3,2)=3.14159; // assigns pi to the fourth element of component 2
 133        FldArrayF g(100);
 134        g[0] = 2.22;
 135
 136 *** 10.1.3 内存
 137   Fld 比c/c++数组的优势
 138   - 内存使用的check control
 139   - 完全的数据初始化控制
 140
 141 *** 10.1.4 传递field数据给Fortran
 142   传递数据给Fortran77函数.通过给FldField动态分配的内存地址
 143 **** 10.1.4.1 FldArray 内部结构
 144   FldArray在一个连续的内存空间存储数据,为动态分配的.我们可以认为它是c 指针管理内存的wrapper._data指向这块内存.这种1D的布局和传统的Fortran和c代码一致.
 145
 146 <source lang="c">
 147 //返回第fld物理变量的第l个节点数据,按_nfld存
 148 inline E_Float
 149 FldArrayF::operator()(E_Int l, E_Int fld) const
 150 {
 151   return (_data[l + (fld-1)*_size]);
 152 }
 153
 154 //We would have the transpose (or swapped) implementation:
 155 //另外一种的存储索引方式,按l存
 156 inline E_Float
 157 FldArrayF::operator()(E_Int l, E_Int fld) const
 158 {
 159   return (_data[fld-1 + l*_nfld]);
 160 }
 161
 162 </source>
 163
 164
 165 **** 10.1.4.2 例子
 166   在fortran中定义
 167 <source lang="c">
 168 extern "C"
 169 {
 170 void
 171 denconvec_(const E_Int& ncell, const E_Int& neqtot, const E_Int& neq,
 172            const E_Int& ro, const E_Int& rou, const E_Int& roe,
 173            const E_Int& rog, const E_Int& roug, const E_Int& roeg,
 174            const E_Float* consvar, const E_Float* press,
 175            E_Float* fdx, E_Float* fdy, E_Float* fdz);
 176 }
 177
 178 </source>
 179
 180   在C++中使用
 181 <source lang="c">
 182 denconvec_(ncell, neqTot, nbEqMoyComp,
 183            rho, mom, ene, rhoG, momG, eneG,
 184            wCons.begin(), press.begin(),
 185            fdx.begin(), fdy.begin(), fdz.begin() );
 186 </source>
 187
 188  - 我们使用第一个index先增长方式,即c方式,按行存储,先o行,再1行
 189 <source lang="shell">
 190 fdX(      0,1), fdX(0,2), fdX(0,3),..., fdX(      0,neq),
 191 ...,
 192 fdX(ncell-1,1), ...,                    fdX(ncell-1,neq)
 193 </source>
 194  - 而不是
 195 <source lang="shell">
 196 fdX(0, 1), fdX(1,1),fdX(2,1),..., fdX(ncell-1, 1),
 197 ...,
 198 fdX(0,neq),...,                   fdX(ncell-1,neq)
 199 </source>
 200
 201   - 最终的结果为
 202 <source lang="fortran">
 203       SUBROUTINE denconvec(ncell, neqtot, neq,
 204     &                          ro, rou, roe,
 205     &                          rog, roug, roeg,
 206     &                          w, p,
 207     &                          fdx, fdy,     fdz)
 208          IMPLICIT NONE
 209 C_IN
 210          INTEGER_E ncell, neqtot, neq
 211          REAL_E        w(0:ncell-1,neqtot)     ! Conservative Variables
 212          REAL_E        p(0:ncell-1)            ! Pressure
 213 C_OUT
 214          REAL_E        fdx(0:ncell-1,neq)      ! Convective Flux      X-Component
 215          REAL_E        fdy(0:ncell-1,neq)      ! Convective Flux      Y-Component
 216          REAL_E        fdz(0:ncell-1,neq)      ! Convective Flux      Z-Component
 217          [.......]
 218          DO icell = 0, ncell-1
 219             roi = ONE / w(icell,rog)
 220             fdx(icell,ro)      = w(icell,roug)
 221             fdy(icell,ro)      = w(icell,rovg)
 222             fdz(icell,ro)      = w(icell,rowg)
 223          [.......]
 224 </source>
 225
 226 *** 10.2 GEO 部件
 227   GeoGrid在kernel中广泛使用,因为很多的CFD类需要一个指向GeoGrid的指针.GeoGrid为GeoGridMetrics和GeoConnect的组合.他们负责
 228  - metrics， : GeoGridMetrics
 229  - index counting 和position(connectivity) : GeoConnect
 230
 231 **** 10.2.1 Ghost 几何元素
 232   通过如下控制Ghost
 233  - GHOST_I1 and GHOST_I2 in the "I" direction,
 234  - GHOST_J1 and GHOST_J2 in the "J" direction,
 235  - GHOST_K1 and GHOST_K2 in the "K" direction.
 236
 237 ***** 10.2.1.1 Ghost cell numbering
 238     • GHOST_I1 ghost cells in IMIN(代表i的最小index), GHOST_I2 ghost cells in IMAX(代表i的最大index),
 239     • GHOST_J1 ghost cells in JMIN, GHOST_J2 ghost cells in JMAX,
 240     • GHOST_K1 ghost cells in KMIN, GHOST_K2 ghost cells in KMAX.
 241
 242 ***** 10.2.1.2 Ghost interface numbering
 243     • GHOST_I1 ghost interfaces in IMIN, GHOST_I2-1 ghost interface in IMAX,
 244     • GHOST_J1 ghost interfaces in JMIN, GHOST_J2-1 ghost interface in JMAX,
 245     • GHOST_K1 ghost interfaces in KMIN, GHOST_K2-1 ghost interface in KMAX.
 246
 247 ***** 10.2.1.3 Ghost node (mesh points) numbering
 248     • GHOST_I1 ghost nodes in IMIN, GHOST_I2-1 ghost node in IMAX,
 249     • GHOST_J1 ghost nodes in JMIN, GHOST_J2-1 ghost node in JMAX,
 250     • GHOST_K1 ghost nodes in KMIN, GHOST_K2-1 ghost node in KMAX.
 251
 252 ***** 10.2.1.4 Ghost defaultvalues
 253        The default values are:
 254      GHOST_I1 = 2; GHOST_I2 = 2;
 255      GHOST_J1 = 2; GHOST_J2 = 2;
 256      GHOST_K1 = 2; GHOST_K2 = 2;
 257
 258 ***** 10.2.1.5 cell/interface/node的分辨
 259    上述定义的好处就是我们在cell,interface和node的index之间关系简化,使得我们在他们中循环很方便.我们在i,j,k方向上的cells,nodes和interface数目完全相同.
 260 <source lang="fortran">
 261 DO n=n0cell, nfcell                  ! loop on cells
 262   ni1 = n                            ! "i" interface (left)
 263   nj1 = n +   ncell                  ! "j" interface (down),ncell是所有的cell个数,因为先存i,后存j,再存k
 264   nk1 = n + 2*ncell                  ! "k" interface (back)
 265   ....
 266   ni2 = n +           inccell(1,0,0) ! "i" interface (right
 267   nj2 = n +   ncell + inccell(0,1,0) ! "j" interface (up)
 268   nk2 = n + 2*ncell + inccell(0,0,1) ! "k" interface (front)
 269
 270 </source>
 271
 272  - 所以所有的interfaces的数目为cells的3倍.
 273 <source lang="c">
 274 inline E_Int
 275 GeoConnect::getNbCell() const
 276 {
 277   return ((_im - 1 + GHOST_I1 + GHOST_I2)*
 278           (_jm - 1 + GHOST_J1 + GHOST_J2)*
 279           (_km - 1 + GHOST_K1 + GHOST_K2));
 280 }
 281 inline E_Int
 282 GeoConnect::getNbInti() const
 283 {
 284   return getNbCell();
 285 }
 286 </source>
 287
 288 *** 湍流模型
 289
 290 **** 抽象类
 291  - 1 计算湍流涡黏性系数
 292  - 2 计算总应力张量(viscosity tensor+Reynolds tensor)
 293  - 3 转淚处理
 294
 295  其中所有的方法2,3的计算方式都一样
 296
 297 **** 计算方式分类
 298  - 基于算术方程: 计算渦黏性系数只需要守恒变量和壁面距离
 299  - 基于输运方程: 需要积分一个方程组.要计算source term,turbulent diffusive flux coefficients
 300
 301 *** Oper 部件
 302   Operator的类主要为空间离散.每一个operator都必须计算一个定义在特殊几何体(cell或者interface)上的Fld-FieldF对象.
 303  - Oper: 抽象类
 304  - Fxc:  对流fluxes(中心差分,迎风格式,人工耗散)
 305  - Fxd:  扩散fluxes(层流,平均流,湍流)
 306  - Sou:  源项(湍流封闭关系,moving frame,dual time stepping method)
 307
 308 **** Oper 模块
 309 ***** OperBase抽象类
 310   OperBase为抽象类.比较重要的属性:
 311  - _geoEntity: 计算几何体的类型
 312  - _borderDepth: 边界区域的宽度.需要特殊计算的区域
 313  - _elemSysOper: 控制方程的id
 314  - EosSysEq* _sysEq: 描述目前的问题
 315  - EosIdealGas* _eos: 目前的状态方程
 316  - GeoGridBase* _grid: 目前的计算网格
 317
 318 <source lang="c">
 319  virtual void prepare(const EosSysEq&, EosElemSysType,EosIdealGas* eos, const GeoGridBase* grid);
 320 </source>
 321  prepare函数用于初始化这些指针.
 322
 323  继承自OperBase的类有两种
 324  - OperTerm 用于计算方程右边值
 325  - Utility  用于附属的计算,OperTerm调用这些类,只实现了两个,OperGrad和OperGradInt.
 326
 327 ***** OperGrad类
 328   继承自OperBase.
 329  - OperGrad    用于在cell center计算梯度
 330  - OperGradInt 用于在interface center计算梯度
 331   这些操作符类不拥有计算结果
 332  - 他们的属性和OperBase一样
 333  - 有两个重载的compute函数计算
 334 <source lang="c">
 335   // compute the gradient of the conservative variables
 336   void compute( FldCellF& fldOut,
 337                 const list<BndPhys*>& listBnd,
 338                 const list<JoinBase*>& listJoin,
 339                 AuxField identOfField = MISC);
 340
 341   // compute the gradient of non conservative variables (?)
 342   virtual void compute( FldCellF& fldOut,
 343                         const list<BndPhys*>& listBnd,
 344                         const list<JoinBase*>& listJoin,
 345                         const FldCellF& fldIn,
 346                         AuxField identOfField);
 347 </source>
 348    梯度和flux计算在所有的cell和interfaces上使用相同的数值算法计算.但是标准的计算公式会在border区域给出错误的结果.所以operator对象需要和boundary对象协作.
 349  - 策略1: 分两步计算,先在所有的cell或interfaces进行计算,然后再在border上修正,OperGrad使用这个策略
 350  - 策略2: 一步计算:通过boundary 将Ghost cell填充在border区域周围,然后用标准算法就会得到正确的结果.OperGradInt使用这个策略.
 351
 352 ***** OperTerm抽象类
 353  计算OperTerm分为两步
 354  1. 计算所有需要的附属场值(通过守恒变量)
 355  2. 计算fluxes和source term(使用第一步计算的数值)
 356
 357  OperTerm::_setOfIdent为list<AuxField>类型,用于识别什么样的附属计算需要进行.computeAllField调用computeFiled()函数计算所有的注册过的附属场.
 358 <source lang="c">
 359   virtual void computeAllField(const list<BndPhys*>* listBnd = E_NULLPTR,
 360                               const list<JoinBase*>* listJoin = E_NULLPTR);
 361
 362   virtual void computeField(AuxField identOfField,
 363                            const list<BndPhys*>* listBnd = E_NULLPTR,
 364                            const list<JoinBase*>* listJoin = E_NULLPTR);
 365 </source>
 366
 367   Flus和source在这一个阶段还没有计算,下一节介绍.
 368
 369 ***** OperFlux抽象类
 370   继承自OperTerm,使用上面介绍的策略1实现了主要的compute()方法.
 371
 372 <source lang="c">
 373   void OperFlux::compute( FldIntF& fldOut,
 374                           const list<BndPhys*>& listBnd,
 375                           const list<JoinBase*>* listJoin )
 376 {
 377    compFlux(...);
 378    computeAllBorders(...);
 379 }
 380
 381 </source>
 382
 383   定义了flux类最重要的抽象方法compFlx.知道了cell center上的守恒变量,这个方法提供了interfaces上的flux.
 384
 385 <source lang="c">
 386   virtual void compFlux(const FldCellF& fldIn,FldIntF& fldOut) = 0;
 387 </source>
 388
 389   OperFlux还定义了另外一个重要的方法:computeAllBorders(),这个方法用于计算边界区域的flux,并且考虑边界条件.所有的物理边界,使用下述的思想:
 390  1. 定义在边界interfaces上的vector wb1的计算代理给了BndPhys* 对象,它负责处理边界
 391  2. wb1 然后被用来修正边界interface的flux计算.具体的实现由特定的flux类决定.
 392
 393 ***** OperSou 抽象类
 394   OperSou是source term的抽象类.这个类继承自OperTerm类.
 395
 396 **** Fxc 模块
 397
 398 **** Fxd 模块
 399
 400 **** Sou 模块
 401
 402 ****