GCN 文件
概述
对于多物理场耦合问题,不同的物理场对应不同的偏微分方程,各个物理场之间存在着相互影响。对有限元语言来说,偏微分方程的表达式是由 PDE 文件和 FBC 文件给出的,为了继续描述不同场之间的相互作用关系,用户还需要编写两个文件:GCN 文件和 MDI 文件。这里给出 GCN 文件的相关说明。
GCN 文件名代表项目名,和 MDI 文件名保持一致。GCN 文件负责描述每个物理场使用的算法、不同物理场之间的耦合方式以及整个问题的求解流程。其编写内容由以下两段信息组成:
- 算法引入段,用于给出各个物理场使用的算法和耦合方式;
- 命令流段,用于给出具体的计算过程,组织和调用系统的组件命令。
下面将分别叙述这两段信息的语法规则。
算法引入段
该信息段起到指定不同物理场使用的算法模板的作用。
关键字 defi
语法规则
defi
场标识符 NFE文件名 场标识符 … 场标识符
场标识符 NFE文件名 场标识符 … 场标识符
场标识符 NFE文件名 场标识符 … 场标识符
… …
说明
- 第一行给出关键字 defi。从第二行开始,每一行表示一个场,有多少个场就写多少行。每一行的第一个字符串为场标识符,按从上到下的声明顺序,由英文字母 a、b、c … 来表示,即第一个场场标识符为 a,第二个场场标识符为 b,…,如此一一对应。第二个字符串为单场使用的算法模板,即 NFE 文件的名字(不含扩展名)。随后可以继续给出若干个场标识符,场标识符之间用空格隔开,以表明本场与其它场之间的耦合方式,也就是求解本场需要用到其它哪些场的计算结果。如果没有耦合关系,则 NFE 文件名之后留空。
- 需要注意的是,第一段信息只是用来初始化不同物理场的算法文件的。利用这些信息,FEtch 系统会调用指定的模板,生成项目专用的 NFE 文件,并命名为 "项目名+场标识符.nfe" 。初始化后,即当目录存在项目专用的 NFE 文件时,此信息段将会失效。用户后续可根据需要,继续修改 "项目名+场标识符.nfe" 文件。
- 对于 NFE 模板文件,系统会优先查找当前工作目录,之后是系统的 NFE 算法库。也就是说,如果工作目录有该 NFE 模板文件,系统就使用工作目录中的 NFE 模板文件;如果没有,则继续查找使用算法库中的 NFE 模板文件。
注:对于场的概念,这里需要解释一下。在有限元语言中,一个待解方程就是一个场。由已知量求梯度等也作为单独的一个场。每个场可以有多个自由度。求解之后,就可以得到场在网格节点上的值。
范例
假设 a 场计算温度,算法模板文件为 ell.nfe,计算结果为节点上的温度值;b 场计算热流,使用了 a 场的计算结果,算法模板文件为 str.nfe。那么,GCN 文件的算法引入段可以写为:
defi
a ell
b str a
将该 GCN 文件命名为 heat.gcn。运行 "脚本初始化" 后,FEtch 系统会生成项目专用的 NFE 文件。在当前工程目录,由算法库的 ell.nfe 生成 heata.nfe,由 NFE 库的 str.nfe 生成 heatb.nfe。如果当前目录已经存在了 heata.nfe 和 heatb.nfe 文件,则算法引入段自动失效。因此,后续只要根据工程需求修改 heata.nfe 和 heatb.nfe 文件即可。
机理解析
如果存在耦合关系,系统会基于本场调用的 NFE 模板文件,按以下规则生成本场的 NFE 文件。
- 如果模板文件中有 %0、%1、%2 … 符号,则依次替换为场标志符 0、场标志符 1、场标志符 2,等等。
- 如果模板文件中有 #0、#1、#2 … 符号,则依次替换为场标志符 0 对应的变量名、场标志符 1 对应的变量名、场标志符 2 对应的变量名,等等。
这些符号的应用情景主要包括与耦合场相关的 coef 语句、var 语句和 read 语句等,有兴趣的用户可以结合系统自带的 NFE 模板文件来学习和研究一下。
NFE 算法库
常见的线性偏微分方程主要有三大类:
椭圆型: $Lu=f$
抛物型: $M\frac{\partial u}{\partial t}+Lu=f$
双曲型: $M\frac{\partial^2 u}{\partial t^2}+C\frac{\partial u}{\partial t}+Lu=f$
对于非线性方程,对应以上格式,我们可以通过线性化方法将其转化成线性形式的表达式。
FEtch 系统已经把各种线性和非线性问题的常用算法做成了模板文件放入了算法库中,用户可以直接使用其中的算法,系统库中提供的算法如下。
算法文件 | 说明 |
---|---|
ell.nfe par.nfe parb.nfe wave.nfe newmark.nfe nell.nfe npar.nfe nparb.nfe nwave.nfe nnewmark.nfe str.nfe |
求解线性椭圆型方程 求解线性抛物型方程,时间离散采用 Crank-Nicolson 法 求解线性抛物型方程,时间离散采用向后差分法 求解线性波动方程,时间离散采用 Wave 速度法,只求解位移和速度 求解线性波动方程,时间离散采用 Newmark 法,求解位移、速度和加速度 求解非线性椭圆型方程 求解非线性抛物型方程,时间离散采用 Crank-Nicolson 法 求解非线性抛物型方程,时间离散采用向后差分法 求解非线性波动方程,时间离散采用 Wave 速度法 求解非线性波动方程,时间离散采用 Newmark 法 求解已知场的梯度场,采用最小二乘法 |
注:NFE 算法库中 NFE 文件的命名规则大致为: ell 代表椭圆,par 代表抛物,wave 和 newmark 代表波动,以上文件名前面加符号 n 表示非线性。
命令流段
该信息段由 Fortran 语言和 FEtch 系统组件命令组成,用于生成多场耦合计算的主程序。 组件命令是本系统的专有命令,较常用的组件命令如下:
求解器类型 | 初始化命令 | 求解命令 | 说明 |
---|---|---|---|
高斯消元法求解器 | startsin | solvsin | 对称矩阵 |
高斯消元法求解器 | startnin | solvnin | 非对称矩阵 |
非对称多波前求解器 | startumf | solvumf | 非对称矩阵 |
Intel MKL Pardiso 求解器 | startsinpds | solvsinpds | 对称矩阵 |
Intel MKL Pardiso 求解器 | startninpds | solvninpds | 非对称矩阵 |
梯度场求解器 | 无 | stress | 最小二乘和集中质量矩阵法 显式算法 |
对于系统的组件命令,以 # 开头作为标记,每个命令后空一格填写场标识符,以表示该组件命令所作用的场。随后,根据需要选择性地填写另一个场标识符,用于表明采用了哪个场的网格和材料数据。该功能常用于算子分裂法和求梯度场的情况。如果只采用本场的数据,则第二个场标识符可以缺省。
需要特别说明的是,梯度场求解器(stress)是一种特殊化了的求解器,专门用于求解已知场的梯度场。它所作用的场的方程组左端矩阵必须为集中质量矩阵。为满足这一要求,该场的 NFE 文件也需要按对应格式进行书写。为此,系统提供了 str.nfe 算法模板文件与之互相匹配。详见弹性静力学等具体算例。
标记符 #
语法规则
#组件命令 场标识符 场标识符
范例
#startsin a
#solvsin a
#stress b a
这表示对 a 场采用对称直接求解器模块进行先初始化,然后进行求解。在求解完 a 场的基础上,继续求解它的梯度场 b 场。
注意,对同一个场,组件命令字符串 start* 和 solv* 中的 * 必须相同,否则可能会出现数据传递错误。
stress 求解器不需要初始化,即没有对应的 start 语句。它需要使用与之耦合的场的网格和材料数据作为输入数据,因此,第二个场标识符是必不可少的。
命令流段最终会转化成项目的主程序,因此,它完全支持 Fortran 90 语句,以及使用全局的系统变量。例如
open(1,file="time.dat")
read(1,*) dt,tmax
close(1)
#startsin a
time = 0.0
do while (time <= tmax-1e-8)
time = time + dt
if (time>tmax) time=tmax
#solvsin a
end do
该命令流段实现了文本读取操作和时间步循环操作。这里的 dt 和 tmax 为全局变量,是系统保留字。详情请参考关于系统保留字的介绍。
注 :由于过于陈旧,原 FEPG 系统采用的 bft、post 等命令已被废除。如果参考了老版本的 GCN 文件,要注意针对性地进行修改。
综合案例
线弹性固体静力学问题
对于一个线弹性固体变形问题,需要解出固体的位移和应力,一个位移场和一个应力场共两个物理场。其中,求解应力时要用到位移场的值,即应力场计算与位移场相耦合。这里使用的是线弹性静力学的数学模型,所以计算位移时应选择隐式的椭圆型方程的算法,计算应力时应选择显式的最小二乘算法。求解顺序是先位移场,后应力场。按照上述分析编写成以下的 elas.gcn 文件。
elas.gcn 文件 | 说明 | |
---|---|---|
1 2 3 4 5 6 7 |
defi a ell b str a #startsin a #solvsin a #stress b a |
算法引入段 场 a,采用线性椭圆问题的算法 ell.nfe,不和其他场耦合 场 b,采用最小二乘算法 str.nfe,并采用 a 场的计算结果 空一行,开始命令流段 对 a 场进行初始化 求解 a 场 求解作为 a 场梯度场的 b 场,采用 a 场的网格和材料数据 |
弹塑性固体静力学问题
弹塑性固体静力学问题与弹性问题相比多了一个非线性迭代过程。容易编写成以下的 plas.gcn 文件。
plas.gcn 文件 | 说明 | |
---|---|---|
1 2 3 4 5 6 7 8 9 10 |
defi a nellp b str a #startsin a iend=0 do while (iend==0) #solvsin a enddo #stress b a |
算法引入段 场 a,采用非线性椭圆问题的算法 nell.nfe 场 b,采用最小二乘算法 str.nfe,耦合 a 场 空一行,开始命令流段 对 a 场进行初始化 对迭代标记 iend 赋初值,iend 是系统保留字,整数型全局变量 迭代循环开始 求解 a 场,计算过程中 a 场会更新 iend 的值 迭代循环结束 求解作为 a 场梯度场的 b 场,采用 a 场的网格和材料数据 |
瞬态热传导问题
如果想解决一个瞬态热传导问题,则需要引入时间步循环,并在每个时间步求解一次温度场的分布。按照上述分析编写成以下的 heat.gcn 文件。
heat.gcn 文件 | 说明 | |
---|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 |
defi a parb #startnin a tmax = 8 dt = 2 time = 0.0 do while (time <= tmax-1d-8) time = time + dt if (time > tmax) time = tmax #solvsin a enddo |
算法引入段 场 a,采用线性抛物问题的算法 parb.nfe 空一行,开始命令流段 初始化 a 场 结束时间 时间步长 初始值 循环开始 更新 time 越界修正 求解 a 场 循环结束 |
上面的文件将时间参数直接 "写死" 了,这么做通常仅仅是为了方便进行程序测试。为了更加灵活和实用,也可以通过读取磁盘文件来引入时间步参数,只需将上述 GCN 文件作简单修改即可。修改后的文件内容和说明如下。
heat.gcn 文件 | 说明 | |
---|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
defi a parb open(21,file='time.dat') read(21,*) dt,tmax close(21) #startsin a time = 0.0 do while (time <= tmax-1d-8) time = time + dt if (time > tmax) time = tmax #solvsin a enddo |
算法引入段 场 a,采用线性抛物问题的算法 parb.nfe 空一行,开始命令流段 从 time.dat 文件读取时间步参数 初始化 a 场 初始值 循环开始 更新 time 越界修正 求解 a 场 循环结束 |
可以看出,用户自定义了 time.dat 文件,来设置时间步长和终止时间。需要注意的是,在计算之前切记要提前准备好该文件,避免在计算时出现文件无法访问的情况。
非线性瞬态热传导问题
非线性瞬态热传导问题的计算流程较为复杂,主要包含两大循环过程,时间步循环和迭代步循环。同时,为了实现在不同时间段使用不同的时间步长,还可以额外增加一个文件读取循环。对应的 GCN 文件的内容和说明如下。
nheat.gcn | 说明 | |
---|---|---|
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 |
defi a nparb cc = 1.0 errb = 1.0e-8 itnmax = 30 #startsin a time = 0.0 it = 1 open(21,file='time.dat') do k=1,10 read(21,*,iostat=ios) dt,tmax if (ios \= 0) exit do while (time <= tmax-1e-8) time = time + dt if (time > tmax) then dt = dt-(time-tmax) time = tmax endif iend = 0 itn = 1 do while (iend==0) #solvsin a itn = itn+1 enddo it = it+1 if (itn < 5) dt = dt*1.2 enddo enddo close(21) |
算法引入段 场 a,采用线性抛物问题的算法 nparb.nfe 空行,开始命令流段 迭代变量赋值 误差阈值 最大迭代次数 初始化场 a 对时间 time 赋初始值 对时间步数 it 赋初始值 打开 time.dat 文件,建议I/O单元号>20 读取循环开始 读取时间步参数 读取错误时退出 时间步循环开始 更新 time 步长越界修正 迭代结束标志赋初值 迭代步数赋初值 迭代步循环开始 求解 a 场 更新 itn 迭代步循环结束 更新 it 变更时间步长 时间步循环结束 读取循环结束 关闭 time.dat 文件 |
注意,这里的 time.dat 文件在计算过程中一直是打开的。由于程序计算过程中还会读写其他的文件,为了避免标记文件的 I/O 单元号发生冲突,建议用户在插入的 Fortran 源代码中采用 20 以上的整数作为 I/O 单元号。