Discontinuous Galerkin method for steady transport problem

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

1.简介

恒定状态对流方程

[egin{equation} acdot abla mathbf{u} = f end{equation}]

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

2.数值离散

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

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

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

[egin{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) 结果为

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

写成矩阵形式为

[egin{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} o 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的解。假设函数满足线性关系：若

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

那么对应的函数有

[egin{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}) ight))，联立 (Au_i = JM cdot f)，便可得到最终未知解 (u_0 = left( u_1, u_2, cdots, u_{Np} ight))。

4.边界条件

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

5.代码

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

function SteadyConvectionDriver
% solving steady convection problem by DGM
% 
abla 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 = Af(:);
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