Discontinuous Galerkin method for steady transport problem

下面讨论如何使用 Discontinuous Galerkin 求解恒定对流问题。

1.简介

恒定状态对流方程

\[\begin{equation} a\cdot \nabla \mathbf{u} = f \end{equation}\]

出现在多种问题中，如海洋模型中求解连续方程计算垂向速度，明渠恒定流动问题等。

2.数值离散

首先需要将方程写为离散格式，以一维问题为例

\[\begin{eqnarray} \begin{aligned} & \frac{\partial u}{\partial x} = sin(x) \quad x\in[0, 2\pi] \cr & u(0) = 1 \end{aligned} \end{eqnarray}\]

将方程乘以实验函数（test function）并在控制单元内积分，利用分部积分将方程转化为

\[\begin{equation} \oint_{\Omega}(l_i \cdot u) n_x dx - \int_{\Omega} \frac{\partial l_i}{\partial x} u dx = \int_{\Omega} l_i f_h dx \end{equation}\]

用基函数线性组合近似变量 \(u \approx u_h = \sum_{j=1}^P l_j u_j\)，源项 \(f_h\) 也采用基函数线性近似 \(f_h = \sum_{j=1}^P l_j f_j\) 结果为

\[\begin{equation} \left( \oint_{\Omega}(l_i \cdot l_j) dx \right) u_j \cdot n_x - \left( \int_{\Omega} \frac{\partial l_i}{\partial x} l_j dx \right) u_j = \left( \int_{\Omega} l_i l_j dx \right) f_j \end{equation}\]

写成矩阵形式为

\[\begin{equation} J_s M_e \hat{u}_n - J M Dr \cdot u = J M \cdot f \end{equation}\]

其中边界通量 \(\hat{u}_n\) 采用如下方式计算

\[\hat{u}_n = u_L \cdot n_x \]

即始终采用边界左侧节点值计算数值通量。这主要是因为边界条件定义在计算域左侧，在计算时也将由左向右逐个单元进行计算，也就是数据从左向右进行传递。

3.联立求解

与包含时间项的非恒定方程不同，恒定对流方程需要将系数矩阵联立求解。

与普通系数矩阵构造不同，这里首先设多元函数 \(L:\mathbb{R}^{Np} \to \mathbb{R}^{Np}\)，

\[L(\mathbf{u}) = J_s M_e \hat{u}_n - J M Dr \cdot \mathbf{u} = (L_1, L_2, \cdots, L_{Np})^T \]

最终目标是寻找变量 \(\mathbf{u}_0\)，使得等式 \(L(\mathbf{u}_0) = JM \cdot f\) 成立。

这里令 \(\mathbf{e}_i\) 表示第\(i\)个分量为单位1，其余分量为0的解。假设函数满足线性关系：若

\[\begin{equation} u_0 = \sum_{i=1}^{Np} u_i \mathbf{e}_i \end{equation}\]

那么对应的函数有

\[\begin{equation} L(\mathbf{u}_0) = \sum_{i=1}^{Np} u_i L(\mathbf{e}_i) \end{equation}\]

那么我们通过构造系数矩阵 \(A = \left( L(\mathbf{e}_1), L(\mathbf{e}_2), \cdots, L(\mathbf{e}_{Np}) \right)\)，联立 \(Au_i = JM \cdot f\)，便可得到最终未知解 \(u_0 = \left( u_1, u_2, \cdots, u_{Np} \right)\)。

4.边界条件

在恒定输运方程中，在给定边界 \(\Gamma_D\) 上为Dirichlet边界 \(u = u_D\)。参考有限元方法，可以采用置大数法，即修改对应系数矩阵 \(A\)，与源项 \(f\) 来耦合 \(\Gamma_D\) 上已知解。具体方法请参考有限元边界 Dirichlet 条件处理

5.代码

SteadyConvectionDriver 负责构造计算所需的网格及标准线单元系数矩阵（刚度矩阵，质量矩阵等），方程源项 \(f\) 也在此给定。

function SteadyConvectionDriver
% solving steady convection problem by DGM
% \nabla u = sin(x)
% 
x1 = 0; x2 = 2*pi; % domain

N = 1; nElement = 20;
[~, VX, ~, EToV] = Utilities.Mesh.MeshGen1D(x1, x2, nElement);
BC = [2,1];

% 
line = StdRegions.Line(N);
mesh = MultiRegions.RegionLineBC(line, EToV, VX, BC);

f = sin(mesh.x);

u = SteadyConvectionSolver(mesh, f);

plot(mesh.x(:), u(:), 'b', mesh.x(:), cos(mesh.x(:)), 'r');
end% func

SteadyConvectionSolver 负责求解方程组，根据所给边界条件，采用置大数法修改对应系数矩阵及右端项系数，

function u = SteadyConvectionSolver(mesh, f)
% set up and solve the equation system
% Input:
%   mesh - mesh object
%   f - source term
% Output:
%   u - unknown variable

f = mesh.J.*(mesh.Shape.M*f);

%% set up and solve global matrix coeffcient

% get system global matrix coefficient
A = SteadyConvectionCoeffMatrix(mesh);

% boundary condition
u0 = 1; M = 1e8;
A(1, 1) = M; f(1) = u0*M;

solvec = A\f(:);
u = reshape(solvec, size(mesh.x) );

end% func

SteadyConvectionCoeffMatrix负责构造系数矩阵，每次计算单位变量 \(\mathbf{e}_{i}\) 对应的多元函数值 \(L(\mathbf{e}_{i})\)

function A = SteadyConvectionCoeffMatrix(mesh)
% set up symmetric matrix

A = zeros(mesh.nNode, mesh.nNode);
g = zeros(size(mesh.x));

% Build matrix -- one column at a time
for i = 1:mesh.nNode
    g(i) = 1;
    
    Avec = SteadyConvectionRHS(mesh, g);
    A(:, i) = Avec(:);
    g(i) = 0;
end% for

end% func

SteadyConvectionRHS计算函数 \(L(\mathbf{u})\)，根据信息传递方向，在边界处数值通量采用迎风格式计算。

function rhs = SteadyConvectionRHS(mesh, u)
% right hands of equation
% 

us = zeros( size(u(mesh.vmapM)) );
us(1, :) = u(mesh.vmapP(1, :));
us(2, :) = u(mesh.vmapM(2, :));

us = us.*mesh.nx;

rhs = (mesh.Shape.Mef * us) - mesh.J.*( mesh.rx .*( mesh.Shape.Dr'*(mesh.Shape.M*u) ));
end% func

6.计算结果

已知方程精确解为 \(u(x) = -cos(x) + 2\)，分别采用不同阶（N=1，2，3）与不同个数（Ne=10，20，40，80）单元进行计算，统计对应的 \(L_1\)、\(L_2\) 误差及收敛速率。

6.1.N=1

Ne	L1	rate
10	0.029536	\
20	0.007977	1.888594
40	0.002040	1.967584
80	0.000513	1.991309

Ne	L2	rate
10	0.037214	\
20	0.009872	1.914490
40	0.002506	1.978170
80	0.000629	1.994526

6.2.N=2

Ne	L1	rate
10	0.001776	\
20	0.000172	3.368547
40	0.000019	3.170322
80	0.000002	3.081502

Ne	L2	rate
10	0.002184	\
20	0.000233	3.227791
40	0.000028	3.073642
80	0.000003	3.020009

6.3.N=3

Ne	L1	rate
10	0.000175	\
20	0.000011	3.938733
40	0.000001	3.975718
80	0.000000	3.855050

Ne	L2	rate
10	0.000189	\
20	0.000012	3.946121
40	0.000001	3.978587
80	0.000000	3.871224