求解偏微分方程开源有限元软件deal.II学习--Step 48

Posted on 2016-11-29   |   In computational material science   |   暂无评论

引子

本例提供了一个框架来应用MatrixFree类,既包括求解非线性偏微分方程过程,同时演示MatrixFree类怎样处理“constraints”以及如何在分布式节点上并行。这个算例显示基于单元的运算在六面体单元的二阶或更高阶插值上要比稀疏矩阵-向量乘法快得多,能达到后者10倍的浮点运算速率。
使用MatrixFree类,可以不用组装一个大型的稀疏矩阵,其运算都是基于单元。这里的并行也充分利用了现代超算机器的架构,分为三个层级:

  • 不同节点之间使用MPI并行
  • 单个节点内使用“动态任务规划”进行线程并行
  • 单个核心内使用处理器的向量单元进行显式向量化并行

很多通用有限元包,比如deal.II、libMesh等都将有限元计算和线性代数分开,然后依赖迭代型求解器的稀疏矩阵-向量乘法的内核,比如直接应用特定线性代数包PETSc或Trilinos等。这个算例想要挑战这种传统的将线性代数和有限元组装分开的思路。

具体来说,就是不组装一个全局的稀疏矩阵,仅仅存储参考单元形函数信息、自由度的枚举、从参考单元到真实单元的变换。通常情况下这个方法能够显著降低存储需求,在数值结算操作上比系数矩阵稍有增大,但有时甚至会更少。降低内存需求会提高浮点运算速率,因为稀疏矩阵运算通常受限于内存带宽,而不是浮点运算。尽管可以将稀疏矩阵优化成比如只运算非零元素的模式,但此时的浮点运算速率很少能超过峰值运算的2%-20%。
除了稀疏矩阵-向量乘法较差的计算性能,对于一个d维的(p-1)次的有限单元来说,矩阵中每排的非零元素的数目与p的d次方成正比,这样对于很高次的方法的开销就很大。如果我们将有限元算子分割成用参考单元的形函数值及其导数表示的函数计算和积分步,那么就可以通过一个张量乘法一次就对一个维度上的形函数信息同时计算,这使得计算复杂度变为每个自由度计算d^2p此,这在文献中称为sum-factorization方法。
之前的有限元方法应用中通常存储单元矩阵及其指标,这些单元矩阵再组成总体刚度矩阵,这是一种“一个元素一个元素”的存储方法。然后这种方法增大了内存消耗量。最近,不明确存储矩阵、基于单元的方法已经用在GPU运算和一些应用导向的软件包如SPEC-FEM 3D中。

算法

这里还是选择一个通用问题:算子A对一个向量u作用。
如果用了MPI并行,那么算法就是:
(1)从其他MPI进程中导入向量值,用于当前MPI进程在自己所拥有的单元上的计算。

update_ghost_values

(2)对当前拥有的单元进行循环(在当前MPI进程中进行线程并行)
(2.1)取得当前单元上向量的局部值
(2.2)通过积分计算当前单元上矩阵与向量的值
(2.3)将当前单元的贡献叠加到全局目标向量中
(3)与另一个MPI进程所拥有的单元交换所计算的信息

compress

第(2.2)步可以通过显式地构建单元矩阵A来实现,这就是通常的存储矩阵的所有元素的方法。为了避免这种显式地存储,就可以用无矩阵方式(还有另外一种方法,是FEniCs软件包中用的方式,但其局限于线性算子和简单几何上,即雅各比变换是常量,所以不如下面更通用):
在单元上通过积分计算“算子作用在向量上”:计算有限元函数u在所有积分点上的值和/或它的导数,然后用所有与此单元相关的试函数来测试。
下面是对方法的具体分析。

局部积分方法

为了方便说明,这里特化算子A为变系数的拉普拉斯算子:

−∇⋅K(x)∇

其中,K是d乘d的对称矩阵。
相应的有限元弱形式为:

(∇ϕj,K∇uh),j=1,…,n

其中, uh(x)=∑ni=1ϕi(x)u(i)是全局有限元函数u的离散插值, ϕj,j=1,…,n是试探函数。n是总自由度个数。
再定义参考单元上的基函数 ϕ^j(x^),j=1,…,pd。对应于全局的积分点 xq,其在参考单元上的对应的积分点就是 x^q。即在参考单元上的对应量都加上个冒,其中 (p−1)是有限单元的“度”,即p是每个方向上自由度的个数,d是维度。
,在某个单元k上,在积分点 xq上,有限元函数u的插值就是:

uhk(xq)=∑i=1pdϕ^i(x^q)u(i)k

在新方法中,显式地计算该局部有限元插值的梯度:

∇uhk(xq)=∑i=1pdJ−Tk(x^q)∇^ϕ^i(x^q)u(i)k=J−Tk(x^q)∑i=1pd∇^ϕ^i(x^q)u(i)k

其中, ∇代表真实坐标系中的梯度, ∇^代表参考单元上的梯度, J−Tk(x^q)是从参考单元到真实单元的逆转置的雅各比矩阵。注意:上面先求梯度,即先求和,再施加几何变换。
目标向量 vk=Akuk的每一个分量i都是一个积分,通过数值积分来获得:

v(i)k=∑q(∇ϕi(xq)TK(xq)∇uh(xq))wq|detJK(x^q)|=∑q(∇^ϕ^i(x^q)T(J−Tk(x^q))TK(xq)∇uh(xq))wq|detJk(x^q)|

因此,对于拉普拉斯算子的局部算法可以表述为:
(1)计算参考单元上所有积分点上的梯度值,即第一个公式的求和部分
(2)在每个积分点上:

  • 施加雅各比变换,从而根据第一个公式得到真实空间上的离散函数u的梯度
  • 乘以变系数矩阵
  • 再乘以雅各比行列式和积分权重
  • 再次施加雅各比矩阵变换

(3)乘以参考单元上形函数的梯度,然后遍历所有积分点求和

这种算法相比于传统的矩阵组装算法有如下特点:

  • 将梯度的计算分成两部分:在参考单元上计算梯度,然后再施加雅各比变换。而不是直接在真实空间中计算形函数梯度。
  • 参考单元上的计算对所有单元都是相同的,所以很多单元可以组合在一块同时计算
  • 对于张量积单元,可以应用下面的“sum-factorization”算法,其相比于传统的直接运算方法(对所有积分点和局部基函数嵌套循环)有更小的复杂度

FEEvaluation类

上面给出了在每个单元上的局部积分算法,这些基于单元的操作全部在FEEvaluation类中实现:

  • 读取源向量,然后再组装回全局目标向量中
  • 使用“sum-factorization”算法计算参考单元上的函数值、梯度及二阶导数,并对其积分
  • 应用雅各比变换、与变系数向量相乘等。

问题描述

这里要求解Sine-Gordon方程:

uttn⋅∇uu(x,t0)=Δu−sin(u)for(x,t)∈Ω×(t0,tf],=0for(x,t)∈∂Ω×(t0,tf],=u0(x).

这里采用显式的时间步进离散方法:

(v,un+1)=(v,2un−un−1−(Δt)2sin(un))−(∇v,(Δt)2∇un)

如果使用Gauss-Lobatto单元,将会产生一个对角质量矩阵M,那么更新方程就变成:

Un+1=M−1L(Un,Un−1)

其中,非线性算子 L(Un,Un−1)

就是上式右端项。

注意点

返回类型

FEEvaluation取得某个自由度和某个积分点上的函数值的返回类型都是value_type类型,其定义为:

typedef FEEvaluationAccess<dim,n_components_,Number> BaseClass;
typedef typename BaseClass::value_type value_type;

而此处用到的FEEvaluationAccess类是:

class FEEvaluationAccess<dim,1,Number>
typedef VectorizedArray<Number> value_type;

所以,此处得到返回值的类型都是VectorizeArray类型的。

同理,对于返回的梯度类型:

typedef FEEvaluationAccess<dim,n_components_,Number> BaseClass;
typedef typename BaseClass::gradient_type gradient_type;

typedef Tensor<1,dim,VectorizedArray<Number> > gradient_type;

即,返回的梯度是一个一阶张量,即矢量,它有dim个分量。

梯度计算

真实单元中的函数值和梯度与参考单元中的不同,两者要进行矩阵变换,同时梯度之间的变换还要再多一道,这在FEEvaluation的submit_gradient函数中能够明确地体现出来:

this->gradients_quad[0][d][q_point] = (grad_in[d] *
  this->cartesian_data[0][d] *
  JxW);

而在计算函数值时则不需要加中间的数值:

this->values_quad[0][q_point] = val_in * JxW;
#deal.II

求解偏微分方程开源有限元软件deal.II学习--Step 48相关推荐

  1. 求解二阶偏微分方程c语言,科学网—求解偏微分方程开源有限元软件deal.II学习--Step 3 - 亓欣波的博文...

    完整版见:qixinbo.info. deal.II的程序结构 deal.II采用面向对象编程,其中包含了很多的Modules,各自实现不同的功能,并有机地结合起来.如上图所示.具体为:Triangu ...

  2. 流固耦合开源软件precice安装笔记(包括开源CFD软件OpenFOAM、插件swak4Foam,开源有限元软件CalculiX、deal.II)

    安装环境:Ubuntu 20.04 LTS 1. 安装Python发行版Anaconda 可在Anaconda官网下载安装包,下载完成后在下载目录中鼠标右键打开终端,键入: bash Anaconda ...

  3. 开源的python有限元软件_想从Abaqus转用开源有限元软件有什么好的推荐吗?

    推荐FreeCAD,Code-Aster,Calculix. (1)FreeCAD据说华为在尝试于FreeCAD基础上进行CAE开发,FreeCAD本是一款开源的CAD软件,具备一些简单的有限元分析功 ...

  4. 编译开源有限元软件calculix

    calculix的编译过程 第一步:下载安装cygwin,选择gcc-core gcc-g++ gfortran perl注意版本一致,且为最低版本 第二步:修改makefile文件,因为最后没有生成 ...

  5. 开源的python有限元软件_python有限元

    广告关闭 腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 节点坐标和单元的拓扑信息暂由有限元前处理软件导出后由python读入#nodes ...

  6. 【机器人学】机器人开源项目KDL源码学习:(5)KDL如何求解几何雅克比矩阵

    这篇文章试图说清楚两件事:1. 几何雅克比矩阵的本质:2. KDL如何求解机械臂的几何雅克比矩阵. 一.几何雅克比矩阵的本质 机械臂的关节空间的速度可以映射到执行器末端在操作空间的速度,这种映射可以通 ...

  7. 免费开源Scada软件 RapidScada学习记录

    免费开源Scada软件 RapidScada学习记录 由于工作需要,需要采集污水处理厂的数据,比如溶解氧,ph等数据,并且要求远程控制功能,设备数据异常手机推送功能,能通过网页查看摄像头当前设备及工作 ...

  8. 深度学习(Deep Ritz,Galerkin,PINN)求解偏微分方程(PDE)解读

    深度学习在求解偏微分方程(PDE)中有三种主要方法: Deep Ritz, Galerkin 和 PINN. Deep Ritz 方法是通过使用深度神经网络来逼近解析解来求解 PDE 的. 它使用了 ...

  9. 深度学习(Deep Ritz,Galerkin,PINN)求解偏微分方程(PDE)实现代码地址

    深度学习求解偏微分方程的实现代码可以在GitHub上查找.下面是几个可能有用的项目地址: Deep Ritz: https://github.com/yseop/DeepRitz Galerkin: ...

  10. 深度学习求解偏微分方程

    深度学习求解偏微分方程 1. 稀疏回归解偏微分方程 2. 离散连续方程解偏微分方程 3. 物理神经网络解偏微分方程(PINN:物理激发的神经网络) 1. 稀疏回归解偏微分方程 论文:<Data- ...

最新文章

  1. MPB:南农成艳芬组-瘤胃微生物体外发酵过程与注意事项
  2. leetcode 两数之和 整数反转 回文数 罗马数字转整数
  3. Java里的容器 Collection 简介
  4. 【project】十次方-01
  5. 文本对抗攻击入坑宝典
  6. imgageJ开发【Java】
  7. 10个性鼠标指针主题包_游戏鼠标推荐
  8. PHP7通过yum源安装及性能测试
  9. 音频,视频合并软件——ffmpeg下载及使用
  10. 关于linux系统安装zabbix报错的解决方案
  11. 深蓝词库转换1.9发布——支持英库拼音、搜狗bin格式、FIT、中州韵等
  12. 如何使用高德地图 API 做一个路线规划应用,展示自定义路线
  13. IPv4地址--公网地址可以有多少
  14. AX4.0 SP2本地化的问题---汇兑损益报表打印
  15. RoboMaster无人机设计
  16. java 交流群 14187321 欢迎java爱好者参与
  17. 有冗余电源的台式计算机,台式机500w电源是否足够
  18. python内建函数是什么意思_python、什么是内建函数?
  19. 个人开公司的流程,以后用得着(经典)
  20. 中医副高级职称计算机考试试题,副高级职称计算机考试模拟题及答案.doc

热门文章

  1. 原码、反码、补码解析及其简单转化
  2. oracle plsql 字符串长度,PLSQL SQL
  3. 游戏开发之使用类封装双链表数据结构及双链表迭代器初版(C++基础)
  4. 华为交换机安全端口实验
  5. OSPFv3中LSA详解(八)——Type5类LSA详解
  6. IOCP中多次投递WSASend
  7. 学习 (2012.01)
  8. 权限管理----用户与角色关系
  9. Python 操作 protobuf 常见用法
  10. 服务器 --- 开发框架