求解偏微分方程开源有限元软件deal.II学习--Step 37
求解偏微分方程开源有限元软件deal.II学习--Step 37
2016-12-26更新:
将方程求解过程更加细致地表述。
引子
本例展示了一种无矩阵(matrix-free)方法的使用,即不明确存储矩阵元素来求解二阶变系数Possion方程。同时求解时使用了多重网格算法。
应用无矩阵方法的原因是:当前科学计算的一个瓶颈是对于内存而不是高速缓存中的数据的读取:比如对于一个矩阵和一个向量的乘法运算,现在的CPU能很快地计算浮点数的乘法和加法,但通常需要很长时间等待数据从内存中传入。因此,如果我们不再从内存中寻找矩阵元素,而是重新计算它们,那么可能整体时间就能减少。
这个无矩阵方法中还涉及向量化/矢量化编程:把for循环的操作,用矩阵操作的形式代替。在向量化编程中,程序设计以向量为基本操作单位,采用向量运算代替循环操作以提高运行效率,这里的“向量”不同于一般数学中的概念, 它指的是数组或矩阵。该类具体进行向量化时,将一些单元合并成一个宏单元,这样用一条指令就可以同时对多个单元进行操作。
算例
这里要求解的是变系数Possion方程:
计算域是 Ω=[0,1]3,变系数是 a(x)=10.05+2|x|2
。系数虽然是以原点对称,但计算域不是,所以得到的结果也是不对称的。
方程的弱形式为:
所以,各项对应到程序中就是:
左端项:
phi.submit_gradient (coefficient(cell,q) * phi.get_gradient(q), q); |
右端项:
rhs_val += (fe_values.shape_value(i,q) * 1.0 * fe_values.JxW(q)); |
注意,上述左端项中仅有试探函数的梯度,而没有试探函数的值,所以,程序中是:
phi.evaluate (false,true,false); phi.integrate (false,true); |
如果仅用到试探函数的值,那么就得变成:
phi.evaluate (true,false,false); phi.submit_value(...); phi.integrate (true,false); |
矩阵与向量的乘法
求解一个方程,形如:
时,如果按照传统思路,需要求一个稀疏矩阵A与一个向量u的乘积。
先来看看有限元矩阵A是怎样组装的:
上式中,长方矩阵 Pcell,loc−glob定义当前单元从局部自由度到全局自由度的指标映射。
如果想要对上述矩阵A乘以一个向量u,即:
编程思路是这样的:先建立单元矩阵,即 Acell,然后用局部自由度和全局自由度之间的映射将真实空间中的u转换成单元上的量 ucell。再用vmult函数实现单元矩阵与单元向量的乘法,即:
,得到目标向量后,再用局部到全局的映射将该向量转换到真实空间中。
这个思路很正确,但是很慢。对于每个单元,都要用三个嵌套循环来构建单元矩阵,然后再用两个嵌套循环作乘法。一个改进的方法是意识到单元矩阵可在概念上视为三个矩阵的乘积:
这个形式跟力学中的单元刚度矩阵就是完全相同的了,矩阵B是形函数梯度矩阵,D是弹性矩阵。
当这个矩阵跟一个向量相乘时:
这样,就是从右往左做三次矩阵-向量乘法。这避免了在构建单元矩阵时的三次嵌套循环。
上述代码中一个瓶颈在于对每个单元都做FEValues::reinit,这个操作所用的时间可能跟其他所有操作加起来相同(至少对非规则网格是这样的,规则网格上的梯度通常不变)。这明显不理想,所以最好是优化一下。reinit所做的工作是:通过雅各比矩阵将参考单元上的梯度进行变换,从而计算真实空间上的梯度。这在每个单元上的每个积分点上对每个形函数都要操作。通常雅各比矩阵不依赖于形函数,但它在不同的积分点上不同。在之前的算例中,矩阵只构建一次,那么就没有必要对这个reinit函数做什么,因此在计算局部矩阵元素时必须得做这个雅各比变换。
然而,在应用无矩阵运算时,我们对只使用矩阵一次没有兴趣,而是想要对矩阵多次使用。所以,就想能不能缓存什么东西用于加速计算,但也不能缓存太多数据,否则就会陷入之前获取内存中数据缓慢的瓶颈。
这里用的方法就是识别出那个雅各比变换,然后仅在一个参考单元上应用一次。
即:
使用一个无矩阵、基于单元的有限元算子,需要使用跟之前代码不同的设计思路,deal.II中存储这种数据结构的类是MatrixFree类,然后FEEvaluation类来具体使用它。
无矩阵对象的使用
本例中,除了问题类的定义,还有一个类,LaplaceOperator,用来表示差分算子。为了通用性,将它设计成一个矩阵,即你可以获得它的大小,也可以将它用在一个向量上。跟真实的矩阵不同点在于:它不存储矩阵的元素,仅仅知道当它用于一个向量时应该怎么做。在这个类中用来存储数据的类是MatrixFree,它包含了局部自由度与整体自由度之间的映射关系,即雅各比矩阵,也能在并行计算时遍历所有单元,它能保证只对不共享自由度的单元进行操作,这就保证了线程安全性,它比之前提到的WorkStream类更加先进。
然后就是这个类中对MatrixFree类型数据的使用,包括以下功能:返回数据的维度、多种形式的矩阵-向量乘法、初始化数据等。这个类需要三个参数:维度(所以能处理不同维度的问题)、有限元的度(这是为了后续FEEvaluation的快速计算)、潜在的标量类型(我们想要对于最终矩阵使用double类型,而对于多重网格上的矩阵使用float类型)。
这个类的数据成员包括:真正使用的MatrixFree对象、在所有积分点上计算的变系数(这样在矩阵-向量乘法中就不用重复计算)、存储矩阵对角元素的容器(用于多重网格平滑)。
初始化数据时,最重要的是创建了MatrixFree的对象实例,同时计算了积分点上的系数:
additional_data.mapping_update_flags = (update_gradients | update_JxW_values | update_quadrature_points); data.reinit (dof_handler, constraints, QGauss<1>(fe_degree+1), additional_data); evaluate_coefficient(Coefficient<dim>()); |
注意,上面在MatrixFree的AdditionalData中需要设定要更新的flag。
然后就是这个类的主要功能所在:计算矩阵-向量乘积。注意:该类中的单元范围通常不等于triangulation的单元数。事实上,这里的“单元”可能是个错误的概念,因为它实际上是多个单元上的积分点的集合,MatrixFree对象将多个单元的积分点分组,变成一个块,形成一个新的向量化高度。这些“单元”的个数可以通过MatrixFree类的MatrixFree::n_macro_cells()获得。与之前的单元迭代器相比,MatrixFree类的所有单元都在同一层级的数组中,这样就可以直接通过整数指标来索引它们。
LaplaceOperator算子的使用步骤如下:
先创建一个FEEvaluation对象,用来进行后续对MatrixFree对象的计算。这个对象接收五个参数,分别是:维度、有限单元的degree、每个方向上积分点的个数(默认是fe.degree+1)、分量的个数(可以是向量,但这里是标量)、数据类型(因为想对多重网格预条件子只设置float类型)。然后就对给定的“单元”范围进行循环,在每个“macro cell”上具体所做的事情如下:
phi.reinit (cell); phi.read_dof_values(src); phi.evaluate (false,true,false); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_gradient (coefficient(cell,q) * phi.get_gradient(q), q); phi.integrate (false,true); phi.distribute_local_to_global (dst); |
(1)告诉那个FEEvaluation对象它要作用的单元,
(2)读入源向量的值,即上面分析中的ucell
,
(3)计算参考单元的梯度。因为FEEvaluation既能计算函数值,也能计算梯度,所以它提供了一个统一的界面来计算从0阶到2阶的梯度(0阶梯度即值本身)。因为这里只需计算梯度,所以就在第二个参数位置设置为true,在第一个和第三个参数位置设置为false。这一步的复杂度比传统的用FEValues来计算的复杂度要降低很多。
(4)然后就是雅各比矩阵变换、与变系数和积分权重的相乘。FEEvaluation有个函数getGradient能应用雅各比变换,同时获得真实空间中的梯度,那么就用它乘以变系数。再用submitGradient施加第二个雅各比矩阵和积分权重。注意:这里不要对积分点重复,否则会造成在此积分点上多次运算。
(5)然后就是积分。使用函数integrate,它的两个参数分别是说明是否对值、梯度进行积分。因为这里仅需对梯度积分,所以将第一个设置为false,第二个设置为true。
(6)最后就是将单元贡献叠加到整体上去。
程序解析
头文件
头文件需要增加两个:
#include <deal.II/matrix_free/matrix_free.h> #include <deal.II/matrix_free/fe_evaluation.h> |
驱动程序
using namespace Step37; LaplaceProblem<dimension> laplace_problem; laplace_problem.run (); |
创建问题类,注意这个类中不再使用稀疏矩阵,而是使用无矩阵方式:
typedef LaplaceOperator<dim,degree_finite_element,double> SystemMatrixType; typedef LaplaceOperator<dim,degree_finite_element,float> LevelMatrixType; SystemMatrixType system_matrix; MGLevelObject<LevelMatrixType> mg_matrices; |
其中无矩阵对象是LaplaceOperator类中的一个数据成员:
MatrixFree<dim,number> data; |
问题类中还提供一个用于输出每段程序进行多长时间的输出流:
ConditionalOStream time_details; |
在程序中其默认是关闭的,需要在构造函数初始化时更改一下让其起作用:
template <int dim> LaplaceProblem<dim>::LaplaceProblem () : fe (degree_finite_element), dof_handler (triangulation), time_details (std::cout, true) {} |
然后执行run函数。
运行函数
run函数控制程序运行的整个流程,包括创建并加密网格、创建系统、组装系统、组装多重网格、求解、输出结果:
if (cycle == 0) {GridGenerator::hyper_cube (triangulation, 0., 1.); triangulation.refine_global (3-dim); } triangulation.refine_global (1); setup_system (); assemble_system (); assemble_multigrid (); solve (); output_results (cycle); std::cout << std::endl; |
创建系统
注意是做一些初始化工作,初始化无矩阵对象、解向量、右端项等。
比如,初始化无矩阵对象:
system_matrix.reinit (dof_handler, constraints); |
此处只有一句话,实际做的东西很多:
template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>::reinit (const DoFHandler<dim> &dof_handler, const ConstraintMatrix &constraints, const unsigned int level) {typename MatrixFree<dim,number>::AdditionalData additional_data; additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::partition_color; additional_data.level_mg_handler = level; additional_data.mapping_update_flags = (update_gradients | update_JxW_values | update_quadrature_points); data.reinit (dof_handler, constraints, QGauss<1>(fe_degree+1), additional_data); evaluate_coefficient(Coefficient<dim>()); } |
从上可以看出,主要做的就是对无矩阵对象的data初始化,同时最后一句话还将算子类中的系数进行了赋值,具体来说就是先用一个FEEvaluation对象将data中的信息提取出来,用其取得积分点位置,这样就能通过系数函数确定积分点上的系数值,然后赋值:
template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: evaluate_coefficient (const Coefficient<dim> &coefficient_function) {const unsigned int n_cells = data.n_macro_cells(); FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data); coefficient.reinit (n_cells, phi.n_q_points); for (unsigned int cell=0; cell<n_cells; ++cell) {phi.reinit (cell); for (unsigned int q=0; q<phi.n_q_points; ++q) coefficient(cell,q) = coefficient_function.value(phi.quadrature_point(q)); } } |
组装系统
这一步不需要组装矩阵,只需要组装右端项,做法跟传统的组装右端项相同:
template <int dim> void LaplaceProblem<dim>::assemble_system () {Timer time; QGauss<dim> quadrature_formula(fe.degree+1); FEValues<dim> fe_values (fe, quadrature_formula, update_values | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) {cell->get_dof_indices (local_dof_indices); fe_values.reinit (cell); for (unsigned int i=0; i<dofs_per_cell; ++i) {double rhs_val = 0; for (unsigned int q=0; q<n_q_points; ++q) rhs_val += (fe_values.shape_value(i,q) * 1.0 * fe_values.JxW(q)); system_rhs(local_dof_indices[i]) += rhs_val; } } constraints.condense(system_rhs); setup_time += time.wall_time(); time_details << "Assemble right hand side (CPU/wall) " << time() << "s/" << time.wall_time() << "s" << std::endl; } |
组装多重网格
求解
用共轭梯度法求解:
cg.solve (system_matrix, solution, system_rhs, preconditioner); |
可以看出,求解形式跟之前的有矩阵的形式是相同的,这个地方是个坑,看似相同,但奥秘隐藏在算子类的vmult函数中。实际,共轭梯度法在计算时会调用这个函数。
template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>::vmult (Vector<double> &dst, const Vector<double> &src) const {dst = 0; vmult_add (dst, src); } |
然后,
template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>::vmult_add (Vector<double> &dst, const Vector<double> &src) const {data.cell_loop (&LaplaceOperator::local_apply, this, dst, src); const std::vector<unsigned int> & constrained_dofs = data.get_constrained_dofs(); for (unsigned int i=0; i<constrained_dofs.size(); ++i) dst(constrained_dofs[i]) += src(constrained_dofs[i]); } |
然后,
template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: local_apply (const MatrixFree<dim,number> &data, Vector<double> &dst, const Vector<double> &src, const std::pair<unsigned int,unsigned int> &cell_range) const {FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data); for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell) {phi.reinit (cell); phi.read_dof_values(src); phi.evaluate (false,true,false); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_gradient (coefficient(cell,q) * phi.get_gradient(q), q); phi.integrate (false,true); phi.distribute_local_to_global (dst); } } |
所以,这个地方是个连环调用。
扩展
要深刻理解方程中各项在程序中是怎样实现的,这样才能扩展。
本例的求解过程实际是prisms-pf中椭圆型方程的求解思路,仔细领会之。
求解偏微分方程开源有限元软件deal.II学习--Step 37相关推荐
- 求解二阶偏微分方程c语言,科学网—求解偏微分方程开源有限元软件deal.II学习--Step 3 - 亓欣波的博文...
完整版见:qixinbo.info. deal.II的程序结构 deal.II采用面向对象编程,其中包含了很多的Modules,各自实现不同的功能,并有机地结合起来.如上图所示.具体为:Triangu ...
- 流固耦合开源软件precice安装笔记(包括开源CFD软件OpenFOAM、插件swak4Foam,开源有限元软件CalculiX、deal.II)
安装环境:Ubuntu 20.04 LTS 1. 安装Python发行版Anaconda 可在Anaconda官网下载安装包,下载完成后在下载目录中鼠标右键打开终端,键入: bash Anaconda ...
- 开源的python有限元软件_想从Abaqus转用开源有限元软件有什么好的推荐吗?
推荐FreeCAD,Code-Aster,Calculix. (1)FreeCAD据说华为在尝试于FreeCAD基础上进行CAE开发,FreeCAD本是一款开源的CAD软件,具备一些简单的有限元分析功 ...
- 编译开源有限元软件calculix
calculix的编译过程 第一步:下载安装cygwin,选择gcc-core gcc-g++ gfortran perl注意版本一致,且为最低版本 第二步:修改makefile文件,因为最后没有生成 ...
- 开源的python有限元软件_python有限元
广告关闭 腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 节点坐标和单元的拓扑信息暂由有限元前处理软件导出后由python读入#nodes ...
- 【机器人学】机器人开源项目KDL源码学习:(5)KDL如何求解几何雅克比矩阵
这篇文章试图说清楚两件事:1. 几何雅克比矩阵的本质:2. KDL如何求解机械臂的几何雅克比矩阵. 一.几何雅克比矩阵的本质 机械臂的关节空间的速度可以映射到执行器末端在操作空间的速度,这种映射可以通 ...
- 免费开源Scada软件 RapidScada学习记录
免费开源Scada软件 RapidScada学习记录 由于工作需要,需要采集污水处理厂的数据,比如溶解氧,ph等数据,并且要求远程控制功能,设备数据异常手机推送功能,能通过网页查看摄像头当前设备及工作 ...
- 深度学习(Deep Ritz,Galerkin,PINN)求解偏微分方程(PDE)解读
深度学习在求解偏微分方程(PDE)中有三种主要方法: Deep Ritz, Galerkin 和 PINN. Deep Ritz 方法是通过使用深度神经网络来逼近解析解来求解 PDE 的. 它使用了 ...
- 深度学习(Deep Ritz,Galerkin,PINN)求解偏微分方程(PDE)实现代码地址
深度学习求解偏微分方程的实现代码可以在GitHub上查找.下面是几个可能有用的项目地址: Deep Ritz: https://github.com/yseop/DeepRitz Galerkin: ...
- 深度学习求解偏微分方程
深度学习求解偏微分方程 1. 稀疏回归解偏微分方程 2. 离散连续方程解偏微分方程 3. 物理神经网络解偏微分方程(PINN:物理激发的神经网络) 1. 稀疏回归解偏微分方程 论文:<Data- ...
最新文章
- SYNCHRO 4D可视化调度学习教程 SYNCHRO 4D: Visual Scheduling
- 什么是网络推广浅析如何提高搜索引擎的抓取频次?
- 火狐浏览器 xml 解析错误:文档元素后存有无效内容_五分钟了解浏览器工作原理...
- Java-开源工具类
- Java的原始字符串文字
- Linux C 内存管理
- centos7:塔建pure_ftpd虚拟用户
- 信息学奥赛一本通(1137:加密的病历单)
- 表的连接方式:NESTED LOOP、HASH JOIN、SORT MERGE JOIN【转】
- CentOS 7 最小化安装简单配置
- Spring MVC笔记 使用JdbcTemplate
- 上偏续关系哈斯图_偏序集的哈斯图G(A)跟A上的偏序关系≤的关系图G(≤)是一 一对应的,相互确定。...
- SHA256算法详解及python实现
- Discuz 手动添加 markdown 代码支持教程!
- 淘票票经典Python爬虫案例
- 橡胶支座抗压弹性模量计算公式_橡胶支座计算
- elementUI的 tree搜索过滤,可识别拼音,且不区分大小写
- java fn replace_JSTL fn:replace()函数替换 换行符
- python笑脸猫图案_酷叮猫编程课堂:python生成字符画
- 网飞文化-真正价值之沟通篇