MOOSE多物理场耦合平台入门学习记录一(稳态热传导程序实例)


本人刚刚使用MOOSE进行热力耦合分析完成了毕业课题,想系统整理一下这段时间学习和使用MOOSE的一些过程,如果有错误的地方或者理解不到位的地方也希望能被指出。MOOSE是一个非常好的开源平台,现在使用的人越来越多,以后应该还会有更大的应用前景。

MOOSE的简介

MOOSE是美国爱达华国家实验室推出的开源多物理场耦合平台,从2014年开源至今已经有8年的开源历史,目前已经相对完善并有成熟的社区。MOOSE的突出特点是极强的扩展性与优秀的计算效率。所谓扩展性体现在:

1、同一物理扩散方程,可以直接应用在一维、二维以及三维网格上;

2、可以非常方便地引入其他模块,或者描述其他物理场的控制方程;

3、可以使用MultiApp功能实现多种并行模式,或者引入外部程序,实现数据传递与交互。

MOOSE底层支持openmp、mpich等多种并行库,并且提供指令可以指定计算资源的分配方式和并行模式,开发者不需要专门编写并行算法程序便可以获得极高的并行效率。

MOOSE的安装

Linux和Mac

MOOSE可以在Linux和Mac平台上原生运行,主要是基于miniconda安装一系列类库,提供虚拟化环境。

在国内,必须使用VPN,因为连接到INL的服务器非常慢,极有可能安装失败。

具体的安装教程,参照下面的网址:

https://mooseframework.inl.gov/getting_started/installation/conda.html

按照网址上的指令,依次执行。

注意的是,有可能在执行mamba init指令后,会提示报错,这可能是因为电脑中存在多个终端环境。一般如果之前的步骤没有问题,选定默认终端环境如bash,只需要重启终端就好。

如果在编译moose过程中报错,首先确定是否所有的类库都已经安装完成,因为MOOSE的下载速度特别慢,通常出现’cannot find ***'之类的错误就是类库缺失,只需要遵循之前的步骤重新安装所有的类库,然后再继续编译。

对于Mac而言,尤其是在Catalina及以前的系统中,因为默认Python版本为2.x会导致编译错误尤其是使用VS code中的终端环境,这种时候,一般使用系统自带的终端切换到moose环境下编译即可,如果还是不可以,建议安装zsh。

Windows

Windows安装MOOSE的方式非常多样,建议使用Windows自带的Ubuntu子系统安装,而不通过虚拟机或者Docker之类的容器技术,这样可以最大化减少安装复杂度和性能损失。

Windows子系统安装教程:

https://blog.csdn.net/user_san/article/details/120703568

建议安装Ubuntu 20.04版本,安装后,一定要先通过下述指令安装gcc和mpich

sudo apt-get install update
sudo apt-get install upgrade
sudo apt-get install gcc
sudo apt-get install mpich

上述指令中,可能需要输入Linux密码(用户自定义密码)
然后所有步骤,与Linux安装教程一致。

MOOSE程序的一般开发流程-以导热微分方程为例

简单问题的有限元处理

一个物理问题,最重要的是内部的物理控制方程、边界条件、初始条件。方程本身就需要封闭存在理论解。MOOSE采用有限元方法求解偏微分方程,这就需要对方程进行有限元离散处理,以导热微分方程为例:

ρ c ∂ T ∂ t = ∂ ∂ x ( λ ∂ T ∂ x ) + ∂ ∂ y ( λ ∂ T ∂ y ) + ∂ ∂ z ( λ ∂ T ∂ z ) + q ˙ \rho c\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {\lambda \frac{{\partial T}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\lambda \frac{{\partial T}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {\lambda \frac{{\partial T}}{{\partial z}}} \right) + \dot q ρctT=x(λxT)+y(λyT)+z(λzT)+q˙

为了便于推导,将上述方程改写成矢量形式,并转成稳态导热方程:

∇ ⋅ ( λ ∇ T ) + q ˙ = 0 \nabla \cdot \left( {\lambda \nabla T} \right) + \dot q = 0 (λT)+q˙=0

左右两侧同时乘以试函数 φ i {\varphi _i} φi,在全域上进行积分,改写成弱格式可以得到:

∫ ∇ ⋅ ( λ ∇ ϕ i ) d v + ∫ q ˙ ϕ i d v = 0 \int \nabla \cdot{(\lambda\nabla \phi_i)}d_v + \int \dot{q}\phi_i d_v = 0 (λϕi)dv+q˙ϕidv=0

假定导热系数 λ \lambda λ为常数,可以提出积分以外,并将热源 q ˙ \dot{q} q˙除以导热系数可以获得如下公式:

∫ ∇ ⋅ ∇ T d v + ∫ q λ ˙ ϕ i d v \int \nabla \cdot \nabla T d_v + \int \dot{\frac{q}{\lambda}}\phi_i d_v Tdv+λq˙ϕidv

令$ \dot{\frac{q}{\lambda}} = v$对以上公式进行简化,然后使用分部积分公式,对扩散项进行降阶:

尤其要注意的是,上式中是点乘符号,因为即便是标量,标量的梯度也是矢量,具有方向。对于分部积分公式及化简后的格式如下:

∫ ∇ ( f ⋅ g ) = ∫ ∇ f ⋅ g + ∫ f ⋅ ∇ g \int \nabla({f\cdot g}) = \int \nabla f\cdot g + \int f \cdot \nabla g (fg)=fg+fg

∫ ∇ ⋅ ( φ i ∇ T ) d v − ∫ ∇ φ i ⋅ ∇ T d v + ∫ v ˙ φ i d v = 0 \int \nabla \cdot (\varphi_i \nabla T)dv - \int \nabla \varphi_i \cdot \nabla T dv + \int \dot{v} \varphi_i dv = 0 (φiT)dvφiTdv+v˙φidv=0

对方程左边第一项采用散度定理:

∯ φ i ∇ T ⋅ n d s ⃗ − ∫ ∇ φ i ⋅ ∇ T d v + ∫ v ˙ φ i d v = 0 \oiint \varphi_i \nabla T\cdot n d\vec{s} - \int\nabla\varphi_i \cdot \nabla Tdv + \int \dot{v} \varphi_i dv = 0 φiTnds φiTdv+v˙φidv=0

第一项实际上就是第二类边界条件,对于外部绝热的物理问题,这一项为0,可以直接省略,将剩余的两项改写为内积形式。

− ( ∇ φ i , ∇ T ) + ( v ˙ , φ i ) = 0 -(\nabla \varphi_i,\nabla T) + (\dot{v},\varphi_i) = 0 (φi,T)+(v˙,φi)=0

这两项对应MOOSE中的两个Kernels,为什么是两个不同的Kernel呢?这是要看试函数的形态,前一项是试函数的梯度,后一项只是试函数与常数相乘。当推导到这一步,便可以作为MOOSE开发程序的依据了,并不需要继续下去根据试函数类型写出每个单元上的积分形式,MOOSE底层已经可以完成接下去的工作。这类问题比较简单,并不涉及到与其他物理量的耦合,接下来用户需要指定的仅是边界条件。

MOOSE程序的开发流程

MOOSE的开发使用的是C++语言,基于面向对象开发模式。所以需要用户本身对C++的类和对象的概念比较熟悉才能比较好上手MOOSE程序开发。

如果是按照官网上的安装教程,在projects文件下输入如下指令,可以创建一个初始的MOOSE程序(不介绍相关git操作)。

./moose/scripts/stork.sh ProjectName

ProjectName是程序名称,执行完上述指令后,会在projects文件夹下面生成一个全小写的程序名称的文件夹,使用cd命令进入文件夹,可以看到里面包含的MOOSE程序基本框架,非常齐全,包括Makefile文件,Makefile文件中不仅包含编译相关设置,同时也包含了是否启用MOOSE相关模块的选项,如果用到了MOOSE提供好的类(如结构力学、多孔介质传热、相场模拟、传热等领域),就必须在要文件中将选项由默认的"no"改为"yes"。


################################## MODULES ####################################
# To use certain physics included with MOOSE, set variables below to
# yes as needed.  Or set ALL_MODULES to yes to turn on everything (overrides
# other set variables).

ALL_MODULES                 := no

CHEMICAL_REACTIONS          := no
CONTACT                     := no
EXTERNAL_PETSC_SOLVER       := no
FLUID_PROPERTIES            := no
FSI                         := no
FUNCTIONAL_EXPANSION_TOOLS  := no
GEOCHEMISTRY                := no
HEAT_CONDUCTION             := no
LEVEL_SET                   := no
MISC                        := no
NAVIER_STOKES               := no
PERIDYNAMICS                := no
PHASE_FIELD                 := no
POROUS_FLOW                 := no
RAY_TRACING                 := no
REACTOR                     := no
RDG                         := no
RICHARDS                    := no
STOCHASTIC_TOOLS            := no
TENSOR_MECHANICS            := yes
XFEM                        := no

include $(MOOSE_DIR)/modules/modules.mk
###############################################################################

如上面代码所示。

最好对刚创建的程序进行编译,先进行一遍测试,看本机的编译环境和生成的程序是否正常。

cd projects/flowtest //这里项目名为flowtest
make -j8  //8核编译,如果编译成功执行下一步
./run_tests -j8

对于简单的二阶方程,只需要进行Kernels的开发,不涉及到边界条件、材料系数项、辅助项等等的开发,只需要创建kernels文件夹。

mkdir ./include/kernels //Kernels类的声明文件存放位置
mkdir ./src/kernels //Kernels类的源文件存放位置
mkdir problems //输入卡及结果文件存放位置

至此,一个简单MOOSE程序的创建准备工作已经完成,开始针对Kernels编写程序。

下面的代码显示了Kernel函数中的一些基本参数说明:

_u[_qp],_grad_u[_qp] //moose需要求解的变量值和梯度
_test[_i][_qp], _grad_test[_i][_qp] //试函数的值和试函数的梯度
_i //当前试函数的索引
_qp //当前积分点的索引,因为moose采用高斯正交积分法求解单元积分,所以积分点与节点并不一致,qp不是节点索引

MOOSE的Kernel基类

上表中,剪切自MOOSE的官方教程,类名带有AD开头,“AD”表示automatic differential 自动微分,极大地简化了MOOSE程序的开发流程。通常只需要重载表中override中的函数(都是依据弱格式方程写的残差求解函数)而不需要再额外编写Jocobian矩阵求解函数。
值得注意的是:

  • ADKernelValue的precomputeQpResidual()方法的数据类型是ADReal,在其中不需要再乘以_test[_i][_qp],因为已经默认乘上了。这个类一般用于源项,或者弱格式中只与 φ i \varphi_i φi相乘的项。
  • ADKernelGrad的precoputeQpResidual()方法的数据类型ADRealVectorValue,同样在其中不需要乘以_grad_test[_i][_qp],其中grad表示梯度,RealVectorValue是向量数据类型,_grad_test是试函数的梯度本身就是矢量,因此,这个类中,如果返回的是_u[_qp],就会报错,因为数据类型不匹配。这个类中需要返回的本身就是矢量类型,乘号代表着点乘运算。所以,如果返回的是_grad_u[_qp]*_test[_i][_qp]就不会报错,因为_grad_u是矢量类型。这个类适用于扩散项的开发,但是需要注意数据类型的匹配问题。

针对上面的简单的有限元问题,在ADKernelGrad和ADKernelValue类上进行继承和重载,完成程序开发。
首先在/include/kernels文件夹下面创建声明文件

#pragma once
#include "ADKernelGrad.h"
class FirstTerm : public ADKernelGrad //继承自ADKernelGrad类
{
public:
  static InputParameters validParams(); //静态成员函数,处理MOOSE输入卡中相关物理量

  FirstTerm(const InputParameters & parameters); //构造函数

protected:
  /// ADKernel objects must override precomputeQpResidual
  virtual ADRealVectorValue precomputeQpResidual() override; //残差求解函数,必须重载
};

然后创建源项的声明文件

#pragma once
#include "ADKernelValue.h"
class SecondTerm : public ADKernelValue //继承自ADKernelValue类
{
public:
  static InputParameters validParams(); //静态成员函数,处理MOOSE输入卡中相关物理量

  SecondTerm(const InputParameters & parameters); //构造函数

protected:
  /// ADKernel objects must override precomputeQpResidual
  virtual ADReal precomputeQpResidual() override; //残差求解函数,必须重载
  const Real & _value; //热源项的大小
};

注意源项与扩散项的precomputeQpResidual的返回数据类型差异,一个是ADReal(实数类型,相当于double类型),一个是ADRealVectorValue(向量类型)。
另外源项中,还定义了变量_value,这个地方是引用,并不是直接的定义变量。_value用来获取输入卡中的参数,改变源项的数值。
然后在src文件夹中创建源文件,FirstTerm.C和SecondTerm.C(命名简单粗暴),这里文件的扩展名为大写的.C后缀,而不是常见的.c文件,如果是.c文件,则在编译中默认使用c语言编译,会报错,而使用.C文件,则按照C++语言进行编译。
代码如下:

#include"FirstTerm.h"
registerMooseObject("flowtest3App", FirstTerm);
//注册类名,括号里前面为项目名称+App,后面为类名,通过这个将类注入程序使用
InputParameters
FirstTerm::validParams()
{
  InputParameters params = ADKernelGrad::validParams(); //初始化静态成员函数
  params.addClassDescription("The diffusion term of the equation"); //加入对象描述
  return params;
}

FirstTerm::FirstTerm(const InputParameters & parameters)
  : ADKernelGrad(parameters) //构造函数
    // Set the coefficients for the diffusion kernel
{
}

ADRealVectorValue
FirstTerm::precomputeQpResidual()
{
  return  -1.0 * _grad_u[_qp]; //实际上就是_grad_u[_qp] * _test[_i][_qp]
}

SecondTerm.C代码如下:

#include"SecondTerm.h"
registerMooseObject("flowtest3App", SecondTerm);
//注册类名,括号里前面为项目名称+App,后面为类名,通过这个将类注入程序使用
InputParameters
SecondTerm::validParams()
{
  InputParameters params = ADKernelValue::validParams(); //初始化静态成员函数
  params.addClassDescription("The diffusion term of the equation"); //加入对象描述
  params.addParam<Real>("value",1.0,"The source term of the equation");
  //定义输入卡中的参数value,默认数值为1.0,类型为实数类型
  return params;
}

SecondTerm::SecondTerm(const InputParameters & parameters)
  : ADKernelValue(parameters), //构造函数
  _value(getParam<Real>("value")) //将输入卡中的value值用来初始化_value
  //此处与扩散项的构造函数不同
    // Set the coefficients for the diffusion kernel
{
}

ADReal
SecondTerm::precomputeQpResidual()
{
  return _value; //实际上就是_value * _test[_i][_qp]
}

上面源码中,每一个kernels的残差求解函数都与前面所推导的弱格式相同,只是省去了直接_test[_i][_qp]和_grad_test[_i][_qp],因为在ADKernel类中已经乘过了。

最后是在problems文件夹中编写输入卡test.i文件

[Mesh]
  type = GeneratedMesh # Can generate simple lines, rectangles and rectangular prisms
  dim = 2              # 定义二维网格
  nx = 100             # x方向上的网格数目
  ny = 100              # y方向上的网格数目
  xmax = 1.0        # x方向上最长距离
  ymax = 1.0        # y方向上最长距离
[]

[Problem]
  type = FEProblem  # This is the "normal" type of Finite Element Problem in MOOSE
[]

[Variables]
  [T]
    #Adds a Linear Lagrange variable by default
    #这里定义形函数的阶数和类型,通过order和family参数来定义
    #如果不加定义,默认就是一阶拉格朗日方程
  []
[]

[Kernels]
  [diffusion]
    type = FirstTerm #
    variable = T  #这里物理量与variables中定义的一致
  []
  [source]
    type = SecondTerm #
    variable = T  #
    value = 100.0 #在程序中已经定义了value,所以value值可以直接被读入
  []
[]

[BCs]
  [inlet]
    type = ADDirichletBC # Simple u=value BC
    variable = T  # Variable to be set
    boundary = left      # Name of a sideset in the mesh
    value = 500         #定义左边界的值
  []
  [outlet]
    type = ADDirichletBC
    variable = T
    boundary = right
    value = 100            # 定义右边界的值
  []
[]

[Executioner]
  type = Steady       # Steady state problem
  solve_type = NEWTON # 这里可以选择求解方法,包括PJFNK、JFNK、NEWTON

  # Set PETSc parameters to optimize solver efficiency
  petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc option pairs with values below
  petsc_options_value = ' hypre    boomeramg'
[]

[Outputs]
  exodus = true # Output Exodus format
[]

以上输入卡,包括MOOSE中输入卡的基本组成部分。后面根据问题的需要,还可以加上Materials、AuxKernel等部分。
在程序文件夹中,编译程序,运行输入卡

make -j8
cd problems
mpiexec -n 8 ../flowtest-opt -i test.i //8核并行运行输入卡

当编译没有错误,执行完成后,在problems文件夹下面生成.e后缀的结果文件,使用paraview打开查看结果。

  • 9
    点赞
  • 27
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值