A primal-dual finite element method for first-order transport problems

TL;DR

This paper introduces a primal-dual weak Galerkin finite element method for first-order transport problems, providing a stable, mass-conserving, and accurate numerical solution under low regularity conditions.

Contribution

It develops a novel PD-WG finite element method that handles low regularity solutions and ensures local mass conservation for first-order transport equations.

Findings

Achieves optimal error estimates in various norms.

Provides numerical results demonstrating accuracy and stability.

Ensures local mass conservation on each element.

Abstract

This article devises a new numerical method for first-order transport problems by using the primal-dual weak Galerkin (PD-WG) finite element method recently developed in scientific computing. The PD-WG method is based on a variational formulation of the modeling equation for which the differential operator is applied to the test function so that low regularity for the exact solution of the original equation is sufficient for computation. The PD-WG finite element method indeed yields a symmetric system involving both the original equation for the primal variable and its dual for the dual variable (also known as Lagrangian multiplier). For the linear transport problem, it is shown that the PD-WG method offers numerical solutions that conserve mass locally on each element. Optimal order error estimates in various norms are derived for the numerical solutions arising from the PD-WG method with weak regularity assumptions on the modelling equations. A variety of numerical results are presented to demonstrate the accuracy and stability of the new method.

Tables

Table 1. Table 1: Numerical error for exact solution u = 1 𝑢 1 u=1 on the unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	$\Vert {\lambda }_0\Vert$	$\Vert {\lambda }_b\Vert$
1	0.00E+00	7.54E-17	4.71E-16
2	3.93E-17	2.58E-16	5.11E-16
4	8.33E-17	7.80E-17	1.39E-16
8	1.06E-16	3.85E-17	7.00E-17
16	1.60E-16	5.70E-17	8.36E-17
32	3.45E-16	1.99E-16	2.86E-16

Table 2. Table 2: Numerical error for exact solution u = 1 𝑢 1 u=1 on the unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	$\Vert {\lambda }_0\Vert$	$\Vert {\lambda }_b\Vert$
1	0.00E+00	5.34E-17	7.28E-17
2	1.04E-16	2.06E-16	4.94E-16
4	1.23E-16	4.72E-16	7.54E-16
8	1.70E-16	1.23E-16	1.92E-16
16	1.81E-16	1.35E-16	1.95E-16
32	3.32E-16	4.87E-17	7.04E-17

Table 3. Table 3: Numerical error for exact solution u = 1 𝑢 1 u=1 on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	$\Vert {\lambda }_0\Vert$	$\Vert {\lambda }_b\Vert$
1	0	7.92E-16	1.67E-15
2	3.64E-16	6.63E-16	1.19E-15
4	3.00E-16	2.76E-16	4.34E-16
8	4.73E-16	7.00E-16	1.02E-15
16	6.84E-16	8.89E-16	1.27E-15
32	1.16E-15	5.21E-16	7.42E-16

Table 4. Table 4: Numerical error for exact solution u = 1 𝑢 1 u=1 on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	$\Vert {\lambda }_0\Vert$	$\Vert {\lambda }_b\Vert$
1	3.14E-16	2.94E-15	5.26E-15
2	1.88E-16	1.37E-15	2.35E-15
4	2.92E-16	1.02E-15	1.53E-15
8	3.30E-16	9.58E-16	1.40E-15
16	8.18E-16	1.46E-15	2.09E-15
32	1.28E-15	1.30E-15	1.85E-15

Table 5. Table 5: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.06461		0.6836		0.7941
2	0.02966	1.123	0.1892	1.853	0.2238	1.827
4	0.01318	1.170	0.04600	2.040	0.05273	2.086
8	0.006401	1.042	1.11E-02	2.046	1.23E-02	2.098
16	0.003180	1.009	2.73E-03	2.029	2.94E-03	2.064
32	0.001589	1.001	6.75E-04	2.016	7.18E-04	2.036

Table 6. Table 6: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.06022		0.8152		0.9571
2	0.02685	1.165	0.2034	2.003	0.2454	1.964
4	0.01212	1.148	0.04709	2.111	0.05471	2.165
8	0.006096	0.9914	1.12E-02	2.066	1.26E-02	2.124
16	0.003092	0.9791	2.74E-03	2.036	2.98E-03	2.076
32	0.001561	0.9858	6.77E-04	2.019	7.23E-04	2.042

Table 7. Table 7: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.1701		1.114		1.769
2	0.08362	1.025	0.3369	1.725	0.5205	1.765
4	0.03665	1.190	0.1355	1.314	0.2020	1.365
8	0.01409	1.379	0.05085	1.414	0.07398	1.449
16	0.005947	1.245	0.01424	1.837	0.02045	1.855
32	0.002693	1.143	0.003720	1.936	0.005309	1.946

Table 8. Table 8: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , − 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,-1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.1590		1.364		2.210
2	0.06908	1.203	0.4459	1.614	0.7237	1.611
4	0.02846	1.279	0.1998	1.158	0.2995	1.273
8	0.01178	1.272	0.05795	1.786	0.08450	1.826
16	0.005410	1.123	0.01505	1.945	0.02164	1.965
32	0.002590	1.063	0.003816	1.979	0.005449	1.990

Table 9. Table 9: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = − 1 𝑐 1 c=-1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.3572		0.9412		1.115
2	0.2116	0.7556	0.9473	-0.009340	1.482	-0.4110
4	0.09173	1.206	0.3473	1.448	0.5146	1.526
8	0.03086	1.571	0.1101	1.657	0.1573	1.710
16	0.01171	1.398	0.02954	1.899	0.04127	1.930
32	0.005356	1.128	0.007508	1.976	0.01037	1.993

Table 10. Table 10: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = − 1 𝑐 1 c=-1 , τ = 0.001 𝜏 0.001 \tau=0.001 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.3582		0.9358		1.102
2	0.2116	0.7593	0.9461	-0.01592	1.480	-0.4250
4	0.09169	1.207	0.3470	1.447	0.5142	1.525
8	0.03084	1.572	0.1101	1.656	0.1572	1.709
16	0.01169	1.399	0.02953	1.898	0.04127	1.930
32	0.005350	1.128	0.007508	1.976	0.01037	1.993

Table 11. Table 11: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = − 1 𝑐 1 c=-1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.4448		0.6546		0.2871
2	0.2306	0.9474	0.7959	-0.2819	1.216	-2.082
4	0.09213	1.324	0.3166	1.330	0.4658	1.384
8	0.02928	1.654	0.1058	1.582	0.1509	1.627
16	0.01022	1.518	0.02920	1.857	0.04080	1.887
32	0.004612	1.149	0.007484	1.964	0.01034	1.980

Table 12. Table 12: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = − 1 𝑐 1 c=-1 , τ = 1000 𝜏 1000 \tau=1000 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.6863		0.1841		0.06168
2	0.3773	0.8632	0.2349	-0.3514	0.2919	-2.243
4	0.1802	1.066	0.1257	0.9026	0.1692	0.7866
8	0.08886	1.020	0.05015	1.325	0.06880	1.299
16	0.03389	1.391	0.01999	1.327	0.02760	1.318
32	0.008018	2.079	0.006678	1.581	0.009201	1.585

Table 13. Table 13: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.6665		1.9483		2.172
2	0.2840	1.231	2.267	-0.2188	3.394	-0.6439
4	0.1080	1.395	0.6346	1.837	0.9082	1.902
8	0.04060	1.412	0.1547	2.036	0.2158	2.074
16	0.01904	1.092	0.03706	2.062	0.05090	2.084
32	0.009381	1.021	0.009125	2.022	0.01244	2.033

Table 14. Table 14: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = 1 𝑐 1 c=1 ; τ = 0.001 𝜏 0.001 \tau=0.001 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.66831		1.9405		2.160
2	0.2840	1.235	2.265	-0.2230	3.390	-0.6501
4	0.1079	1.396	0.6342	1.836	0.9076	1.901
8	0.04055	1.412	0.1547	2.036	0.2157	2.073
16	0.01902	1.092	0.03706	2.061	0.05089	2.084
32	0.009374	1.021	0.009125	2.022	0.01243	2.033

Table 15. Table 15: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = 1 𝑐 1 c=1 ; τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.9362		1.076		0.5230
2	0.3222	1.539	1.774	-0.7216	2.616	-2.322
4	0.1133	1.507	0.5773	1.619	0.8230	1.668
8	0.03770	1.588	0.1515	1.930	0.2111	1.963
16	0.01745	1.111	0.03689	2.038	0.05067	2.059
32	0.008610	1.019	0.009107	2.018	0.01242	2.029

Table 16. Table 16: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ sin ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\sin(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ 1 , 1 ] 𝜷 1 1 {\boldsymbol{\beta}}=[1,1] , c = 1 𝑐 1 c=1 ; τ = 1000 𝜏 1000 \tau=1000 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	1.208		0.2749		0.1019
2	0.5886	1.037	0.3499	-0.3481	0.3809	-1.903
4	0.2668	1.142	0.1865	0.9081	0.2370	0.6846
8	0.1119	1.254	0.07380	1.337	0.1001	1.243
16	0.03311	1.756	0.02724	1.438	0.03748	1.417
32	0.009289	1.834	0.008355	1.705	0.01143	1.714

Table 17. Table 17: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ y − 1 , − x + 1 ] 𝜷 𝑦 1 𝑥 1 {\boldsymbol{\beta}}=[y-1,-x+1] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.2577		1.829		3.933
2	0.1345	0.9376	0.4039	2.179	0.7490	2.393
4	0.06931	0.9570	0.07443	2.440	0.1299	2.528
8	0.03677	0.9143	0.01586	2.231	0.02558	2.344
16	0.01970	0.9006	3.71E-03	2.096	0.005621	2.186
32	0.01050	0.9072	9.03E-04	2.038	0.001318	2.093

Table 18. Table 18: Numerical rates of convergence for exact solution u = sin ⁡ ( x ) ​ cos ⁡ ( y ) 𝑢 𝑥 𝑦 u=\sin(x)\cos(y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 = [ y − 1 , − x + 1 ] 𝜷 𝑦 1 𝑥 1 {\boldsymbol{\beta}}=[y-1,-x+1] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.1779		3.173		6.112
2	0.1147	0.6330	0.4859	2.707	0.8792	2.797
4	0.06512	0.8168	0.07896	2.622	0.1396	2.655
8	0.03560	0.8714	0.01633	2.273	0.02684	2.378
16	0.01929	0.8843	0.003771	2.115	0.005791	2.212
32	0.01033	0.9002	9.13E-04	2.047	0.001341	2.111

Table 19. Table 19: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on the cracked square domain Ω 3 subscript Ω 3 \Omega_{3} : 𝜷 = [ y , − x ] 𝜷 𝑦 𝑥 {\boldsymbol{\beta}}=[y,-x] , c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.5008		2.009		3.088
2	0.2900	0.7881	1.0270	0.9680	1.615	0.9350
4	0.1555	0.8994	0.3126	1.716	0.4735	1.770
8	0.07209	1.109	0.06884	2.183	0.09902	2.257
16	0.03540	1.026	0.01601	2.105	0.02202	2.169
32	0.01765	1.004	0.003888	2.041	0.005219	2.077

Table 20. Table 20: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on the cracked square domain Ω 3 subscript Ω 3 \Omega_{3} : 𝜷 = [ y , − x ] 𝜷 𝑦 𝑥 {\boldsymbol{\beta}}=[y,-x] , c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.3958		3.348		4.985
2	0.2562	0.6275	1.417	1.240	2.207	1.176
4	0.1477	0.7943	0.3559	1.994	0.5442	2.020
8	0.06988	1.080	0.07175	2.310	0.1042	2.385
16	0.03472	1.009	0.01628	2.140	0.02250	2.211
32	0.01740	0.9972	0.003922	2.053	0.005271	2.094

Table 21. Table 21: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 ​ ( x , y ) = [ 1 , − 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[1,-1] if y < 1 − x 𝑦 1 𝑥 y<1-x and 𝜷 ​ ( x , y ) = [ − 1 , 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[-1,1] otherwise, c = 1 𝑐 1 c=1 , τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.1411		0.6130		0.7236
2	0.1427	-0.01700	0.4904	0.3219	0.6096	0.2472
4	0.07178	0.9916	0.1188	2.046	0.1388	2.134
8	0.03058	1.231	3.33E-02	1.833	0.03937	1.818
16	0.01446	1.080	8.44E-03	1.981	9.86E-03	1.997
32	0.007131	1.020	2.11E-03	2.003	2.43E-03	2.018

Table 22. Table 22: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on an unit square domain Ω 1 subscript Ω 1 \Omega_{1} : 𝜷 ​ ( x , y ) = [ 1 , − 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[1,-1] if y < 1 − x 𝑦 1 𝑥 y<1-x and 𝜷 ​ ( x , y ) = [ − 1 , 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[-1,1] otherwise, c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.1953		0.7864		1.058
2	0.1136	0.7815	0.5297	0.5703	0.6675	0.6650
4	0.06042	0.9109	0.1236	2.099	0.1466	2.187
8	0.02824	1.097	0.03433	1.849	0.04123	1.830
16	0.01383	1.030	0.008656	1.988	0.01026	2.007
32	0.006899	1.003	0.002155	2.006	0.002524	2.023

Table 23. Table 23: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 ​ ( x , y ) = [ 1 , − 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[1,-1] if y < 1 − x 𝑦 1 𝑥 y<1-x and 𝜷 ​ ( x , y ) = [ − 1 , 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[-1,1] otherwise; c = 1 𝑐 1 c=1 ; τ = 1 𝜏 1 \tau=1 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.4095		2.047		2.881
2	0.2575	0.6693	1.484	0.4640	2.180	0.4020
4	0.1138	1.178	0.5139	1.530	0.7252	1.588
8	0.04564	1.319	0.1377	1.900	0.1912	1.923
16	0.02152	1.085	0.03337	2.045	0.04585	2.060
32	0.01060	1.022	0.008210	2.0231	0.01122	2.031

Table 24. Table 24: Numerical rates of convergence for exact solution u = sin ⁡ ( π ​ x ) ​ cos ⁡ ( π ​ y ) 𝑢 𝜋 𝑥 𝜋 𝑦 u=\sin(\pi x)\cos(\pi y) on the L-shaped domain Ω 2 subscript Ω 2 \Omega_{2} : 𝜷 ​ ( x , y ) = [ 1 , − 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[1,-1] if y < 1 − x 𝑦 1 𝑥 y<1-x and 𝜷 ​ ( x , y ) = [ − 1 , 1 ] 𝜷 𝑥 𝑦 1 1 {\boldsymbol{\beta}}(x,y)=[-1,1] otherwise, c = 1 𝑐 1 c=1 , τ = 0 𝜏 0 \tau=0 .

$1/h$	$\Vert e_h\Vert$	order	$\Vert {\lambda }_0\Vert$	order	$\Vert {\lambda }_b\Vert$	order
1	0.3716		3.493		5.315
2	0.2102	0.8219	1.907	0.8731	2.835	0.9066
4	0.1047	1.006	0.5779	1.723	0.8193	1.791
8	0.04487	1.222	0.1414	2.031	0.1967	2.059
16	0.02160	1.055	0.03373	2.067	0.04639	2.084
32	0.01070	1.013	0.008287	2.025	0.01133	2.034

Equations

$$\begin{array}{rl}\nabla \cdot (\mathbf{\beta }u)+cu=& f,\text{in }\Omega ,\\ u=& g,\text{on }{\Gamma }_{-},\end{array}$$

$${\Gamma }_{-}=\{\mathbf{x}\in \partial \Omega :\mathbf{\beta }(\mathbf{x})\cdot \mathbf{n}<0\},$$

$$(u,\mathbf{\beta }\cdot \nabla \sigma -c\sigma )=⟨g,\mathbf{\beta }\cdot \mathbf{n}\sigma {⟩}_{{\Gamma }_{-}}-(f,\sigma ),\forall \sigma \in H_{0,{\Gamma }_{+}}^1(\Omega ),$$

$$H_{0,{\Gamma }_{+}}^1(\Omega )=\{v\in H^1(\Omega ):v=0\text{on}{\Gamma }_{+}\}.$$

$$(u,\mathbf{\beta }\cdot {\nabla }_w\{\sigma \}-c\sigma )=⟨g,\mathbf{\beta }\cdot \mathbf{n}\sigma {⟩}_{{\Gamma }_{-}}-(f,\sigma ),\forall \sigma \in H_{0,{\Gamma }_{+}}^1(\Omega ),$$

$$(u_h,\mathbf{\beta }\cdot {\nabla }_{w,h}\sigma -c\sigma )=⟨g,\mathbf{\beta }\cdot \mathbf{n}\sigma {⟩}_{{\Gamma }_{-}}-(f,{\sigma }_0),\forall \sigma \in W_h^{0,{\Gamma }_{+}},$$

$$(v,\mathbf{\beta }\cdot {\nabla }_{w,h}{\lambda }_h-c{\lambda }_h)=0,\forall v\in M_h.$$

$$\{\begin{array}{rl}s({\lambda }_h,\sigma )+(u_h,\mathbf{\beta }\cdot {\nabla }_{w,h}\sigma -c\sigma )=& ⟨g,\mathbf{\beta }\cdot \mathbf{n}\sigma {⟩}_{{\Gamma }_{-}}-(f,{\sigma }_0),\forall \sigma \in W_h^{0,{\Gamma }_{+}},\\ (v,\mathbf{\beta }\cdot {\nabla }_{w,h}{\lambda }_h-c{\lambda }_h)=& 0,\forall v\in M_h,\end{array}$$

$${\int }_{\partial T}{\text{F}}_h\cdot \mathbf{n}ds+{\int }_{T}c{\tilde{u}}_{h}dT={\int }_{T}fdT,T\in {\mathcal{T}}_h.$$

$$\mathcal{W}(T)=\{v=\{v_0,v_b\}:v_0\in L^2(T),v_b\in L^2(\partial T)\}.$$

$$({\nabla }_{w}v,\mathbf{\psi }{)}_T=-(v_0,\nabla \cdot \mathbf{\psi }{)}_T+⟨v_b,\mathbf{\psi }\cdot \text{n}{⟩}_{\partial T},\forall \mathbf{\psi }\in [H^1(T){]}^d.$$

$$({\nabla }_{w,r,T}v,\mathbf{\psi }{)}_T=-(v_0,\nabla \cdot \mathbf{\psi }{)}_T+⟨v_b,\mathbf{\psi }\cdot \text{n}{⟩}_{\partial T},\forall \mathbf{\psi }\in [P_r(T){]}^d,$$

$$({\nabla }_{w,r,T}v,\mathbf{\psi }{)}_T=(\nabla v_0,\mathbf{\psi }{)}_T-⟨v_0-v_b,\mathbf{\psi }\cdot \text{n}{⟩}_{\partial T},\forall \mathbf{\psi }\in [P_r(T){]}^d,$$

$$W_j(T)=\{\{{\sigma }_0,{\sigma }_b\}:{\sigma }_0\in P_j(T),{\sigma }_b\in P_j(e),e\subset \partial T\}.$$

$$W_{j,h}^{0,{\Gamma }_{+}}=\{\{{\sigma }_0,{\sigma }_b\}\in W_{j,h}:{\sigma }_b{|}_e=0,e\subset {\Gamma }_{+}\}.$$

$$M_{k-1,h}=\{w:w{|}_T\in P_{k-1}(T),\forall T\in {\mathcal{T}}_h\}.$$

$$({\nabla }_w\sigma ){|}_T={\nabla }_{w,k-1,T}(\sigma {|}_T),\sigma \in W_{j,h}.$$

$$s(\rho ,\sigma )$$

$$b(v,\sigma )$$

$$s_T(\rho ,\sigma )$$

$$b_T(v,\sigma )$$

$$s({\lambda }_h,\sigma )+b(u_h,\sigma )$$

$$b(v,{\lambda }_h)$$

$$\mathbf{\beta }\cdot \nabla \Psi -c\Psi$$

$$\Psi$$

$$\Vert \Psi {\Vert }_{\gamma }\le C\Vert \theta \Vert ,$$

$$Q_{h}w=\{Q_{0}w,Q_{b}w\}$$

$${\nabla }_w(Q_{h}w)={\mathcal{Q}}_h(\nabla w),\forall w\in H^1(T).$$

$${\boldsymbol{\beta}}\cdot\nabla\lambda_{0}-c\lambda_{0}=0\quad\mbox{strongly in $\Omega$}.$$

$$b(u_h,\sigma )=0,\forall \sigma \in W_{j,h}^{0,{\Gamma }_{+}}.$$

$$\begin{array}{rl}0=& b(u_h,\sigma )\\ =& \sum _{T\in {\mathcal{T}}_h}(u_h,\mathbf{\beta }\cdot {\nabla }_w\sigma -c{\sigma }_0{)}_T\\ =& \sum _{T\in {\mathcal{T}}_h}-({\sigma }_0,\nabla \cdot (\mathbf{\beta }{u}_h){)}_T+⟨{\sigma }_b,\mathbf{\beta }{u}_h\cdot \mathbf{n}{⟩}_{\partial T}-(u_h,c{\sigma }_0{)}_T\\ =& -\sum _{T\in {\mathcal{T}}_h}({\sigma }_0,\nabla \cdot (\mathbf{\beta }{u}_h)+c{u}_h{)}_T+\sum _{e\subset {\mathcal{E}}_h\setminus {\Gamma }_{+}}⟨{\sigma }_b,[[\mathbf{\beta }{u}_h\cdot \mathbf{n}]]{⟩}_e,\end{array}$$

$$\begin{array}{r}0=\sum _{T\in {\mathcal{T}}_h}{h}_T^2\Vert \nabla \cdot (\mathbf{\beta }{u}_h)+c{u}_h{\Vert }_T^2+\sum _{e\subset {\mathcal{E}}_h\setminus {\Gamma }_{+}}{h}_T\Vert [[\mathbf{\beta }{u}_h\cdot \mathbf{n}]]{\Vert }_e^2,\end{array}$$

Full text

1 Introduction

In this paper we are concerned with the development of new numerical methods for first-order linear convection equations in divergence form. For simplicity, consider the model problem that seeks an unknown function $u$ satisfying

[TABLE]

where $\Omega$ is an open bounded and connected domain in $\mathbb{R}^{d}\ (d=2,3)$ with Lipschitz continuous boundary $\partial\Omega$, $\Gamma_{-}$ is the inflow portion of the boundary defined by

[TABLE]

where ${\mathbf{n}}$ is the unit outward normal vector to the boundary $\partial\Omega$ at a given point. Assume that the convection vector satisfies ${\boldsymbol{\beta}}=(\beta_{1},\cdots,\beta_{d})\in[L^{\infty}(\Omega)]^{d}$ and is locally as smooth as $W^{1,\infty}$, the reaction coefficient $c\in L^{\infty}(\Omega)$, the load function $f\in L^{2}(\Omega)$, and the inflow boundary data $g\in L^{2}(\Gamma_{-})$.

First-order linear partial differential equations (PDEs) of hyperbolic-type are called transport equations or linear convection equations. Problems of hyperbolic-type arise in many areas of science and engineering, such as fluid dynamics and neutron transport. In the past several decades, there have been increasing research activities devoted to the development of accurate and efficient numerical methods for hyperbolic problems. Due to the localized phenomena, such as propagating discontinuities and sharp transition layers, it has been a challenging task to develop effective numerical methods for hyperbolic problems. Due to largely the fact that linear hyperbolic PDEs admit discontinuous solutions for non-smooth boundary data, it is difficult to develop numerical methods that provide high-order accurate approximations in regions of smooth solution as well as sharp resolution of discontinuity while avoiding spurious oscillations at discontinuities [4]. Linear hyperbolic equations also form prototype equations for general hyperbolic equations, such as systems of nonlinear conservation laws [32] or transport equations in phase space [33]. It has been shown that successful numerical methods for linear convection equations can be used as building blocks for the numerical solution of complex hyperbolic PDEs [32]. In literature, a series of numerical methods have been developed for linear transport equations, including the streamline-upwind Petrov-Galerkin method [23], the residual distribution framework [29, 23, 1], the least-squares finite element methods [18, 5, 26, 36, 26, 3, 36], the stabilized finite element methods [13], and various discontinuous Galerkin finite element methods [35, 31, 19, 21, 27, 28, 24, 15, 10, 14, 16, 17].

Most of the existing study for the linear transport equation (1.1) typically assumes certain coercivity condition on the convection vector ${\boldsymbol{\beta}}$ and the reaction coefficient $c$ in the form of $c+\frac{1}{2}\nabla\cdot{\boldsymbol{\beta}}\geq\alpha_{0}$ or alike for some fixed positive number $\alpha_{0}$. In practical applications, these conditions are very restrictive and often rule out many important physics such as those with compressible flow fields and exothermic reactions [13]. One such exception is the stabilized finite element method developed by Burman in [13], in which the convection vector is assumed to satisfy ${\boldsymbol{\beta}}\in[W^{2,\infty}(\Omega)]^{d}$ and $c\in W^{1,\infty}(\Omega)$.

The objective of this paper is to develop a new numerical method for the linear convection problem (1.1) for piecewise smooth convection vector ${\boldsymbol{\beta}}$ and reaction coefficient $c$ without assuming any coercivity conditions on the equation. Our new numerical schemes will be devised by following the primal-dual weak Galerkin (PD-WG) framework introduced and studied in [39, 40, 41, 42, 38]. The PD-WG finite element method was originally developed for the second order elliptic PDEs in non-divergence form through a constraint optimization approach of the problem in which the constraint was given by a straightforward discretization of the PDE and the optimization was imposed to minimize the “discontinuity” of the approximating functions. The resulting Euler-Lagrange formulation reveals a symmetric problem involving both the primal (original) equation and the dual (adjoint) equation integrated through various stabilizers designed to provide certain “weak continuity or smoothness”. The framework of the primal-dual methods in numerical PDEs was also developed by Burman [12, 13] in other finite element contexts and was named stabilized finite element methods. The PD-WG finite element methods have shown great promises for PDE problems where no traditional variational formulations are readily available for analysis and discretization.

Let us briefly discuss the philosophy and the key ingredient in the PD-WG finite element method for the first-order hyperbolic problem (1.1). First of all, the solution of the model problem (1.1) can be characterized by seeking $u\in L^{2}(\Omega)$ such that

[TABLE]

where $\Gamma_{+}=\partial\Omega\setminus\Gamma_{-}$ is the outflow boundary satisfying ${\boldsymbol{\beta}}\cdot{\mathbf{n}}\geq 0$, and $H_{0,\Gamma_{+}}^{1}(\Omega)$ is the subspace of $H^{1}(\Omega)$ with vanishing boundary value on $\Gamma_{+}$; i.e.,

[TABLE]

Secondly, using the weak gradient operator $\nabla_{w}$ [43], we may reformulate (1.2) as follows:

[TABLE]

where $\{\sigma\}=\{\sigma|_{T},\sigma|_{\partial T}\}$ is understood as a weak function in the WG context. The weak function is discretized by piecewise polynomials in each element $T$ as well as on its boundary ${\partial T}$. The weak gradient operator $\nabla_{w}$ is then approximated by vector-valued polynomials, denoted as $\nabla_{w,h}$. The weak formulation (1.3) can then be approximated by seeking $u_{h}\in M_{h}$ (i.e., trial space) such that

[TABLE]

where $W_{h}^{0,\Gamma_{+}}$ is a test space for the weak functions with vanishing boundary value on $\Gamma_{+}$. However, the problem (1.4) is not well-posed unless the inf-sup condition of Babus̆ka [2] is satisfied. The primal-dual scheme overcomes this difficulty by coupling (1.4) with its dual equation which seeks $\lambda_{h}\in W_{h}^{0,\Gamma_{+}}$ satisfying

[TABLE]

A formal coupling between (1.4) and (1.5) makes an effective numerical scheme through the use of a stabilizer, denoted by $s(\lambda,\sigma)$. This stabilizer measures the level of “continuity” of $\sigma\in W_{h}$ in the sense that $\sigma\in W_{h}$ is of classical $C^{0}$-conforming if and only if $s(\sigma,\sigma)=0$. The resulting primal-dual weak Galerkin method for solving the hyperbolic model problem (1.1) seeks $u_{h}\in M_{h}$ and $\lambda_{h}\in W_{h}^{0,\Gamma_{+}}$, such that

[TABLE]

where $s(\cdot,\cdot)$ is known as the stabilizer or smoother which enforces a certain weak continuity for the numerical Lagrangian multiplier $\lambda_{h}$ in the weak finite element space $W_{h}^{0,\Gamma_{+}}$.

In this paper, we show that the PD-WG finite element method (1.6) has one and only one solution if the linear convection problem (1.1) admits at most one solution. The numerical scheme (1.6) will be shown to be conservative locally on each element in the sense that there exists a numerical solution $\tilde{u}_{h}$ and a numerical flux $\textbf{F}_{h}$ satisfying (see Theorem 4)

[TABLE]

Some optimal order error estimates will be derived for the numerical solution arising from the PD-WG scheme (1.6) under ultra-weak assumptions on the convection vector ${\boldsymbol{\beta}}$ and the reaction coefficient $c$.

The paper is organized as follows. In Section 2 we shall briefly review the weak gradient operator as well as its discretization. In Section 3, we give a detailed presentation on the primal-dual weak Galerkin algorithm for the linear hyperbolic problem (1.1). A discussion on the solvability (i.e., the solution existence and uniqueness) of the PD-WG scheme is given in Section 4. In Section 5 we show that the PD-WG method offers numerical solutions that conserve mass locally on each element. Next in Section 6, we shall derive an error equation for the PD-WG approximations. In Section 7, we establish some technical estimates/results useful in convergence analysis. Section 8 is devoted to the establishment of an optimal order error estimate for the PD-WG approximations in some discrete Sobolev norms. In Section 9, an error estimate in a weak $L^{2}$ topology is derived based on a local $H^{1}$-regularity assumption for the dual problem. In Section 10, a series of numerical results are reported to demonstrate the effectiveness and accuracy of the PD-WG method developed in the previous sections.

Throughout the paper, we follow the usual notations for Sobolev spaces and norms. For any open bounded domain $D\subset\mathbb{R}^{d}$ with Lipschitz continuous boundary, denote by $\|\cdot\|_{s,D}$ and $|\cdot|_{s,D}$ the norm and semi-norm in the Sobolev space $H^{s}(D)$ for any $s\geq 0$, respectively. The inner product in $H^{s}(D)$ is denoted by $(\cdot,\cdot)_{s,D}$. The space $H^{0}(D)$ coincides with $L^{2}(D)$, for which the norm and the inner product are denoted by $\|\cdot\|_{D}$ and $(\cdot,\cdot)_{D}$, respectively. When $D=\Omega$, or when the domain of integration is clear from the context, the subscript $D$ is dropped in the norm and the inner product notations.

2 Discrete Weak Gradient

The principle differential operator in the weak formulation (1.2) for the linear convection equation in divergence form (1.1) is given by the gradient operator. In this section we shall review the weak gradient operator as well as its discrete version introduced originally in [43].

Let $T$ be a polygonal or polyhedral domain with boundary $\partial T$. By a weak function on $T$ we mean a pair $v=\{v_{0},v_{b}\}$ such that $v_{0}\in L^{2}(T)$ and $v_{b}\in L^{2}(\partial T)$. The components $v_{0}$ and $v_{b}$ can be understood as the value of $v$ in the interior and on the boundary of $T$, respectively. Note that $v_{b}$ is not necessarily the trace of $v_{0}$ on $\partial T$, though taking $v_{b}$ as the trace of $v_{0}$ on $\partial T$ is a feasible option. Denote by ${\mathcal{W}}(T)$ the space of weak functions on $T$; i.e.,

[TABLE]

The weak gradient of $v\in{\mathcal{W}}(T)$, denoted by $\nabla_{w}v$, is defined as a continuous linear functional in the Sobolev space $[H^{1}(T)]^{d}$ satisfying

[TABLE]

Denote by $P_{r}(T)$ the space of polynomials on $T$ with degree $r$ and less. A discrete version of $\nabla_{w}v$ for $v\in{\mathcal{W}}(T)$, denoted by $\nabla_{w,r,T}v$, is defined as the unique polynomial-valued vector in $[P_{r}(T)]^{d}$ satisfying

[TABLE]

which, from the usual integration by parts, gives

[TABLE]

provided that $v_{0}\in H^{1}(T)$.

3 Primal-Dual Weak Galerkin Algorithm

Let ${\cal T}_{h}$ be a partition of the domain $\Omega\subset{\mathbb{R}}^{d}(d=2,3)$ into polygons in 2D or polyhedra in 3D which is shape regular described as in [43]. Denote by ${\mathcal{E}}_{h}$ the set of all edges or flat faces in ${\cal T}_{h}$ and ${\mathcal{E}}_{h}^{0}={\mathcal{E}}_{h}\setminus\partial\Omega$ the set of all interior edges or flat faces. Denote by $h_{T}$ the meshsize of $T\in{\cal T}_{h}$ and $h=\max_{T\in{\cal T}_{h}}h_{T}$ the meshsize of the partition ${\cal T}_{h}$.

For any integer $j\geq 0$, denote by $W_{j}(T)$ the local space of discrete weak functions; i.e.,

[TABLE]

Patching $W_{j}(T)$ over all the elements $T\in{\cal T}_{h}$ through a common value $v_{b}$ on the interior interface ${\mathcal{E}}_{h}^{0}$ yields a global weak finite element space $W_{j,h}$. Let $W_{j,h}^{0,\Gamma_{+}}$ be the subspace of $W_{j,h}$ with vanishing boundary values on $\Gamma_{+}$; i.e.,

[TABLE]

For any integer $k\geq 1$, let $M_{k-1,h}$ be the space of piecewise polynomials of degree $k-1$; i.e.,

[TABLE]

The discrete weak gradient $\nabla_{w,r,T}$ shall be taken in the polynomial subspace $P_{r}(T)$ with $r=k-1$. For simplicity of notation and without confusion, denote by $\nabla_{w}\sigma$ the discrete weak gradient $\nabla_{w,k-1,T}\sigma$ for any $\sigma\in W_{j,h}$ computed by (2.1) on each element $T$; i.e.,

[TABLE]

For any $\rho,\sigma\in W_{j,h}$ and $v\in M_{k-1,h}$, we introduce the following bilinear forms

[TABLE]

where

[TABLE]

with $\tau\geq 0$ being a parameter.

The numerical scheme for the linear convection equation (1.1) in divergence form based on the variational formulation (1.2) in the framework of primal-dual approach is given as follows:

Primal-Dual Weak Galerkin Algorithm 3.1**.**

Find $(u_{h};\lambda_{h})\in M_{k-1,h}\times W_{j,h}^{0,\Gamma_{+}}$ such that

[TABLE]

4 Solution Existence and Uniqueness

The following is the adjoint problem for the linear transport equation: For a given $\theta\in L^{2}(\Omega)$, find $\Psi$ such that

[TABLE]

The adjoint problem (4.1)-(4.2) is said to have the $H^{\gamma}$-regularity with some parameter $\gamma\in(0,1]$ if it has a solution $\Psi\in H^{\gamma}(\Omega)$ satisfying

[TABLE]

with a generic constant $C$.

For simplicity, we now introduce several $L^{2}$ projection operators. On each element $T$, denote by $Q_{0}$ the $L^{2}$ projection operator onto $P_{j}(T)$. For each edge or face $e\subset\partial T$, denote by $Q_{b}$ the $L^{2}$ projection operator onto $P_{j}(e)$. For any $w\in H^{1}(\Omega)$, we use $Q_{h}w$ to denote the $L^{2}$ projection of $w$ in the finite element space $W_{j,h}$ given by

[TABLE]

on each element $T$. ${\cal Q}_{h}$ is used to denote the $L^{2}$ projection operator onto the space $M_{k-1,h}$.

Lemma 1**.**

[43]** For $j\geq k-1$, the $L^{2}$ projection operators $Q_{h}$ and ${\cal Q}_{h}$ satisfy the following commutative property:

[TABLE]

For the convenience of analysis, in what follows of this paper, we assume that the convection vector ${\boldsymbol{\beta}}$ and the reaction coefficient $c$ are both piecewise constants with respect to the partition ${\mathcal{T}}_{h}$. However, all the analysis and results can be generalized and extended to piecewise smooth cases for the convection vector ${\boldsymbol{\beta}}$ and the reaction coefficient $c$.

Theorem 2**.**

Assume that the linear transport problem (1.1) and the dual problem (4.1)-(4.2) have unique solutions. If $\tau>0$, then the PD-WG algorithm (3.5)-(3.6) has one and only one solution for $j=k-1,k$.

Proof.

It suffices to show that the homogeneous problem of (3.5)-(3.6) has only the trivial solution. To this end, we assume $f=0$ and $g=0$. By letting $v=u_{h}$ and $\sigma=\lambda_{h}$ in (3.5) and (3.6), we obtain $s(\lambda_{h},\lambda_{h})=0$, which implies $\lambda_{0}=\lambda_{b}$ and ${\boldsymbol{\beta}}\cdot\nabla\lambda_{0}-c\lambda_{0}=0$ on each $\partial T$. It follows that

[TABLE]

Since, on the portion $\Gamma_{+}$ of the domain boundary, $\lambda_{0}=\lambda_{b}\equiv 0$ holds true, then from the solution uniqueness for the adjoint problem (4.1)-(4.2) we have $\lambda_{0}\equiv 0$ in $\Omega$. It follows from $\lambda_{0}=\lambda_{b}$ on each $\partial T$ that $\lambda_{b}\equiv 0$ and thus $\lambda_{h}\equiv 0$.

Next, substituting $\lambda_{h}\equiv 0$ into (3.5) gives

[TABLE]

From the usual integration by parts and (2.1) we obtain

[TABLE]

where we have used $\sigma_{b}=0$ on $\Gamma_{+}$ on the last line, and ${[\![}{\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}{]\!]}$ is the jump of ${\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}$ on $e\subset{\cal E}_{h}\setminus\Gamma_{+}$ in the sense that ${[\![}{\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}{]\!]}={\boldsymbol{\beta}}u_{h}|_{T_{1}}\cdot{\mathbf{n}}_{1}+{\boldsymbol{\beta}}u_{h}|_{T_{2}}\cdot{\mathbf{n}}_{2}$ for $e=\partial T_{1}\cap\partial T_{2}\subset{\cal E}_{h}^{0}$ with ${\mathbf{n}}_{1}$ and ${\mathbf{n}}_{2}$ being the unit outward normal directions to $\partial T_{1}$ and $\partial T_{2}$, respectively, and ${[\![}{\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}{]\!]}={\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}$ for $e\subset\Gamma_{-}$. By setting $\sigma_{0}=-h_{T}^{2}(\nabla\cdot({\boldsymbol{\beta}}u_{h})+cu_{h})$ on each $T\in{\cal T}_{h}$ and $\sigma_{b}=h_{T}{[\![}{\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}{]\!]}$ on each $e\subset{\cal E}_{h}\setminus\Gamma_{+}$, we may rewrite (4.6) as follows:

[TABLE]

which gives $\nabla\cdot({\boldsymbol{\beta}}u_{h})+cu_{h}=0$ on each $T\in{\cal T}_{h}$, ${[\![}{\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}{]\!]}=0$ on each $e\subset{\cal E}_{h}^{0}$, and ${\boldsymbol{\beta}}u_{h}\cdot{\mathbf{n}}=0$ on each $e\subset\Gamma_{-}$. This implies that $\nabla\cdot({\boldsymbol{\beta}}u_{h})+cu_{h}=0$ in $\Omega$ and $u_{h}=0$ on $\Gamma_{-}$. Thus, from the solution uniqueness assumption, we have $u_{h}\equiv 0$ in $\Omega$. This completes the proof. ∎

Theorem 3**.**

Assume that the linear transport problem (1.1) and the dual problem (4.1)-(4.2) have unique solutions. The following results hold true:

• •

if $j=k-1$, then the PD-WG algorithm (3.5)-(3.6) has one and only one solution for any non-negative value of the stabilizer parameter $\tau\geq 0$.

• •

if additionally the dual problem (4.1)-(4.2) has the $H^{\gamma}$- regularity (4.3) with some $0<\gamma\leq 1$, then the PD-WG algorithm (3.5)-(3.6) has a unique solution for any $\tau=0$ and $j=k-1,k$ provided that the meshsize $h<h_{0}$ holds true for a sufficiently small, but fixed $h_{0}>0$.

Proof.

It suffices to show that the homogeneous problem of (3.5)-(3.6) has only the trivial solution. To this end, we assume $f=0$ and $g=0$. By letting $v=u_{h}$ and $\sigma=\lambda_{h}$ in (3.5) and (3.6) we arrive at $s(\lambda_{h},\lambda_{h})=0$, which implies $\lambda_{0}=\lambda_{b}$ on each $\partial T$ for any $\tau\geq 0$. It follows from (3.6) and (2.2) that

[TABLE]

where we have used $\lambda_{0}=\lambda_{b}$ on each $\partial T$. By taking $v={\boldsymbol{\beta}}\cdot\nabla\lambda_{0}-c{\cal Q}_{h}\lambda_{0}$ we obtain

[TABLE]

on each element $T\in{\cal T}_{h}$. From $\lambda_{0}=\lambda_{b}$ on each $\partial T$, we have $\lambda_{0}\in H^{1}(\Omega)$ so that

[TABLE]

where $\theta=c({\cal Q}_{h}\lambda_{0}-\lambda_{0})$.

Case 1:  $j=k-1$. Since ${\cal Q}_{h}\lambda_{0}=\lambda_{0}\equiv 0$ in this case, then we have $\theta=0$. It follows from $\lambda_{0}|_{\Gamma_{+}}=0$ and the solution uniqueness for (4.7) that $\lambda_{0}\equiv 0$.

Case 2:  with the $H^{\gamma}$- regularity (4.3). In this case, we use the $H^{\gamma}$- regularity assumption (4.3) and the error estimate for the $L^{2}$ projection operator ${\cal Q}_{h}$ to obtain

[TABLE]

which gives

[TABLE]

This implies that $\lambda_{0}\equiv 0$ in $\Omega$ provided that the meshsize $h<h_{0}$ holds true for a sufficiently small but fixed $h_{0}>0$ such that $Ch^{\gamma}<1$.

For both cases, from $\lambda_{0}=\lambda_{b}$ on each $\partial T$, we obtain $\lambda_{b}\equiv 0$ so that $\lambda_{h}\equiv 0$ is verified.

The proof of $u_{h}\equiv 0$ in $\Omega$ can be easily carried out by using exactly the same argument for obtaining $u_{h}\equiv 0$ in Theorem 2. Details are thus omitted here. This completes the proof. ∎

5 Mass Conservation

The linear convection equation (1.1) can be rewritten in a conservative form as follows:

[TABLE]

On each element $T\in{\mathcal{T}}_{h}$, we may integrate (5.1) over $T$ to obtain the integral form of the mass conservation:

[TABLE]

We claim that the numerical solution arising from the primal-dual weak Galerkin scheme (3.5)-(3.6) for the linear convection problem (1.1) retains the local mass conservation property (5.3) with a numerical flux ${\textbf{F}}_{h}$. To this end, for any given $T\in{\cal T}_{h}$, by choosing a test function $\sigma=\{\sigma_{0},\sigma_{b}=0\}$ in (3.5) such that $\sigma_{0}=1$ on $T$ and $\sigma_{0}=0$ elsewhere, we obtain

[TABLE]

It follows from (2.1) and the usual integration by parts that

[TABLE]

where ${\mathbf{n}}$ is the outward normal direction to $\partial T$. The equation (5.4) implies that the primal-dual weak Galerkin algorithm (3.5)-(3.6) conserves mass with a numerical solution and a numerical flux given by

[TABLE]

It remains to show that the numerical flux ${\textbf{F}}_{h}\cdot{\mathbf{n}}$ is continuous across each interior edge or flat face. To this end, we choose a test function $\sigma=\{\sigma_{0}=0,\sigma_{b}\}$ in (3.5) such that $\sigma_{b}$ is arbitrary on one interior edge or flat face $e=\partial T_{1}\cap\partial T_{2}$, and $\sigma_{b}=0$ elsewhere, to obtain

[TABLE]

where we have used (2.1), ${\mathbf{n}}_{T_{1}}$ and ${\mathbf{n}}_{T_{2}}$ are the unit outward normal directions along $e=\partial T_{1}\cap\partial T_{2}$ pointing exterior to $T_{1}$ and $T_{2}$, respectively. This shows that

[TABLE]

and hence the continuity of the numerical flux along the normal direction on each interior edge or flat face.

The result can be summarized as follows.

Theorem 4**.**

Let $(u_{h};\lambda_{h})$ be the numerical solution of the linear convection model problem (1.1) arising from the primal-dual weak Galerkin finite element method (3.5)-(3.6). Define a new numerical approximation and a numerical flux function as follows:

[TABLE]

Then, the flux approximation $\textbf{F}_{h}$ is continuous across each interior edge or flat face in the normal direction, and the following conservation property is satisfied:

[TABLE]

6 Error Equations

Let $u$ and $(u_{h};\lambda_{h})\in M_{k-1,h}\times W_{j,h}^{0,\Gamma_{+}}$ be the exact solution of (1.1) and the numerical solution arising from the primal-dual weak Galerkin scheme (3.5)-(3.6), respectively. Note that the exact solution of the Lagrangian multiplier is $\lambda=0$. The error functions for the primal variable $u$ and the dual variable $\lambda$ are thus given by

[TABLE]

Lemma 5**.**

Let $u$ and $(u_{h};\lambda_{h})\in M_{k-1,h}\times W_{j,h}^{0,\Gamma_{+}}$ be the exact solution of (1.1) and the numerical solution arising from the primal-dual weak Galerkin scheme (3.5)-(3.6), respectively. Then, the error functions $e_{h}$ and $\epsilon_{h}$ satisfy the following equations:

[TABLE]

where

[TABLE]

Proof.

From (3.6) we have

[TABLE]

which gives rise to the equation (6.2).

Next, by subtracting $b({\cal Q}_{h}u,\sigma)$ from both sides of (3.5) we arrive at

[TABLE]

where we have used (2.2), the usual integration by parts, and the facts that $\nabla\cdot({\boldsymbol{\beta}}u)+cu=f$, $u=g$ on $\Gamma_{-}$, and $\sigma_{b}=0$ on $\Gamma_{+}$. This completes the proof of the lemma. ∎

7 Some Technical Estimates

Recall that ${\cal T}_{h}$ is a shape-regular finite element partition of the domain $\Omega$. For any $T\in{\cal T}_{h}$ and $\phi\in H^{1}(T)$, the following trace inequality holds true [43]:

[TABLE]

If $\phi$ is a polynomial on the element $T\in{\cal T}_{h}$, then from the inverse inequality we have [43],

[TABLE]

The following defines a semi-norm in the finite element space $M_{k-1,h}$:

[TABLE]

Lemma 6**.**

Assume that the solution of the linear transport problem (1.1) is unique. Then, the seminorm ${|||}\cdot{|||}_{M_{h}}$ given in (7.3) defines a norm in the linear space $M_{h}$.

Proof.

We shall only verify the positivity property for ${|||}\cdot{|||}_{M_{h}}$. To this end, we assume that ${|||}v{|||}_{M_{h}}=0$ for some $v\in M_{k-1,h}$. It follows from (7.3) that $\nabla\cdot({\boldsymbol{\beta}}v)+cv=0$ on each $T\in{\cal T}_{h}$, ${[\![}{\boldsymbol{\beta}}v\cdot{\mathbf{n}}{]\!]}_{e}=0$ on each $e\in{\cal E}_{h}^{0}$ and ${\boldsymbol{\beta}}v\cdot{\mathbf{n}}=0$ on each $e\subset\Gamma_{-}$. Therefore, we obtain $\nabla\cdot({\boldsymbol{\beta}}v)+cv=0$ in $\Omega$ and $v=0$ on $\Gamma_{-}$, which gives $v\equiv 0$ in $\Omega$ from the uniqueness assumption for the solution of the linear convection problem (1.1). This completes the proof of the lemma. ∎

We further introduce a semi-norm in the weak finite element space $W_{j,h}$ as follows. For any $\lambda=\{\lambda_{0},\lambda_{b}\}\in W_{j,h}$, define

[TABLE]

with $\tau\geq 0$ being the parameter in the stabilizer $s_{T}(\cdot,\cdot)$ given in (3). It is readily seen that

[TABLE]

Lemma 7** (inf-sup condition).**

For $j=k-1$ and $j=k$, there exists a constant $C>0$, such that for any $v\in M_{k-1,h}$ there exists $\tilde{\sigma}_{v}\in W_{j,h}^{0,\Gamma_{+}}$ satisfying

[TABLE]

Proof.

For any $v\in M_{k-1,h}$ and $\sigma\in W_{j,h}^{0,\Gamma_{+}}$, from the definition of the weak gradient (2.1) we have

[TABLE]

where we have used $\sigma_{b}=0$ on $\Gamma_{+}$ in the last line. By setting $\tilde{\sigma}_{v}=\{\tilde{\sigma}_{0};\tilde{\sigma}_{b}\}$ where $\tilde{\sigma}_{0}=-h_{T}^{2}(\nabla\cdot({\boldsymbol{\beta}}v)+cv)$ on each $T$ and $\tilde{\sigma}_{b}=h_{T}{[\![}{\boldsymbol{\beta}}v\cdot{\mathbf{n}}{]\!]}$ on each $e\subset{\cal E}_{h}\setminus\Gamma_{+}$, we have

[TABLE]

Now using the triangle inequality and the trace inequality (7.2) we obtain

[TABLE]

Moreover, we have from the inverse inequality that

[TABLE]

It follows from (7.8) and (7.9) that

[TABLE]

Thus, combining (7.7) and (7.10) gives the inf-sup condition (7.6). This completes the proof of the lemma. ∎

8 Error Estimates

The goal of this section is to derive some error estimates for the solution of the primal-dual weak Galerkin algorithm (3.5)-(3.6) by using the error equations (6.1)-(6.2).

Lemma 8**.**

[43]** Let ${\cal Q}_{h}$ be the $L^{2}$ projection operator onto the finite element space $M_{k-1,h}$ subordinated to the shape-regular finite element partition ${\cal T}_{h}$ for the domain $\Omega$. For any $0\leq s\leq 1$, $s\leq m\leq k$, $s-1\leq n\leq k$, there holds

[TABLE]

[TABLE]

Theorem 9**.**

Let $u$ and $(u_{h};\lambda_{h})\in M_{k-1,h}\times W_{j,h}^{0,\Gamma_{+}}$ be the exact solution of the linear convection problem (1.1) and its numerical approximation arising from the primal-dual weak Galerkin scheme (3.5)-(3.6). Assume that the solution $u$ of (1.1) is sufficiently regular such that $u\in\oplus_{i=1}^{J}H^{k}(\Omega_{i})$ where $\{\Omega_{i}\}_{i=1}^{J}$ is a non-overlapping partition of the domain $\Omega$. Then, the following error estimate holds true:

[TABLE]

Proof.

By letting $\sigma=\epsilon_{h}$ in (6.1) and $v=e_{h}$ in (6.2), we arrive at

[TABLE]

where

[TABLE]

From the relation (7.5) and the above two equations we have

[TABLE]

It follows from the Cauchy-Schwarz inequality, the trace inequality (7.1), (7.4), and (8.1) with $m=k$ that

[TABLE]

Analogously, we have from the standard Cauchy-Schwarz inequality and the estimate (8.1) with $m=k$ that

[TABLE]

Substituting (8.5)-(8.6) into (8.4) yields

[TABLE]

which asserts the estimate for $\epsilon_{h}$ in (8.3). Observe that the left-hand side of (8.6) is zero for $j=k-1$. Thus, the estimate (8.7) holds true without the factor $\tau^{-\frac{1}{2}}$.

To estimate $e_{h}$, from the inf-sup condition in Lemma 7 we have a function $\tilde{\sigma}\in W_{j,h}^{0,\Gamma_{+}}$ such that

[TABLE]

It follows from the error equation (6.1), the estimates (8.7) and (8.8) that

[TABLE]

which gives rise to the estimate for $e_{h}$ in (8.3). Combining (8.7) with (8.9) completes the proof of the theorem. ∎

9 Error Estimate in $L^{2}$

We shall use the usual duality argument to establish an error estimate in $L^{2}$. To this end, consider the problem of seeking $\phi$ such that

[TABLE]

where $\eta\in L^{2}(\Omega)$. The auxiliary problem (9.1)-(9.2) is said to satisfy a local $H^{1}$-regularity if there exists a non-overlapping partition of the domain $\Omega=\bigcup_{i=1}^{J}\Omega_{i}$ such that the solution $\phi$ exists, $\phi\in H^{1}(\Omega_{i})$ for all values of $i$, and

[TABLE]

where $C$ is a generic constant.

Denote by $X_{1}^{*}$ the set of all functions $\eta\in L^{2}(\Omega)$ so that the dual problem (9.1)-(9.2) has the local $H^{1}$-regularity satisfying the estimate (9.3).

Theorem 10**.**

Let $u_{h}\in M_{k-1,h}$ be the numerical solution of the linear convection problem (1.1) arising from the PD-WG algorithm (3.5)-(3.6), with $\lambda_{h}\in W_{j,h}^{0,\Gamma_{+}}$ being the numerical Lagrangian multiplier. Assume that the exact solution $u$ is sufficiently regular such that $u\in\oplus_{i=1}^{J}H^{k}(\Omega_{i})$. Under the local $H^{1}$-regularity assumption (9.3) for the dual problem (9.1)-(9.2), the following error estimate holds true:

[TABLE]

Proof.

For any $\sigma=\{\sigma_{0};\sigma_{b}\}\in W_{j,h}^{0,\Gamma_{+}}$, we use (3.5), (1.1), and the usual integration by parts to obtain

[TABLE]

where we have used $u=g$ on $\Gamma_{-}$ and $\sigma_{b}=0$ on $\Gamma_{+}$ in the last two lines. The last equation can be rewritten as follows

[TABLE]

For any $\eta\in X_{1}^{*}$, let $\phi\in\oplus_{i=1}^{J}H^{1}(\Omega_{i})$ be the solution of the dual problem (9.1)-(9.2). By letting $\sigma=Q_{h}\phi=\{Q_{0}\phi,Q_{b}\phi\}$ in (9.5) we have from (9.1), (4.4), and (9.5) that

[TABLE]

where $\{I_{j}\}_{j=1}^{3}$ represent the corresponding terms in the two lines above them. The rest of the proof is devoted to an analysis for the three terms $I_{1}$, $I_{2}$, and $I_{3}$.

For the term $I_{1}$, we use the Cauchy-Schwarz inequality, the triangle inequality, the trace inequality (7.1), the regularity assumption (9.3), and the estimate (8.2) with $n=0$, and (9.1) to obtain

[TABLE]

As to the second term $I_{2}$, from the Cauchy-Schwarz inequality, the trace inequality (7.1), the estimate (8.1) with $m=k$, the estimate (8.2) with $n=0$, and the regularity assumption (9.3), we obtain

[TABLE]

Finally for the term $I_{3}$, we have

[TABLE]

Substituting (9.7)-(9.9) into (9.6) yields

[TABLE]

which, together with (8.3), gives

[TABLE]

This completes the proof of the Theorem. ∎

10 Numerical Experiments

In this section, we shall numerically demonstrate the accuracy and order of convergence for the primal-dual weak Galerkin finite element method (3.5)-(3.6) through a series of numerical experiments. The finite element partition ${\mathcal{T}}_{h}$ is generated through a successive uniform refinement of a coarse triangulation of the domain by dividing each coarse triangular element into four congruent sub-triangles by connecting the mid-points of the three edges of the triangular element.

The numerical experiments are conducted on both convex and non-convex polygonal domains. The convex domain is given by the unit square $\Omega_{1}=(0,1)^{2}$. The non-convex domains include two cases: (i) the L-shaped domain $\Omega_{2}$ with vertices $A_{1}=(0,0)$, $A_{2}=(2,0)$, $A_{3}=(2,1)$, $A_{4}=(1,1)$, $A_{5}=(1,2)$, and $A_{6}=(0,2)$; (ii) the cracked square domain $\Omega_{3}=(-1,1)^{2}\setminus(0,1)\times{0}$ with a crack along the edge $(0,1)\times 0$. The inflow boundary $\Gamma_{-}$ is determined by using the condition of ${\boldsymbol{\beta}}\cdot{\mathbf{n}}<0$ with ${\mathbf{n}}$ being the unit outward normal direction to $\partial\Omega$. The right-hand side function $f$ and the inflow Dirichlet boundary data $g$ are chosen to match the exact solution $u$ (if the exact solution $u$ is known for the test problem).

The primal-dual weak Galerkin scheme (3.5)-(3.6) is implemented for the lowest order $k=1$ and $j=k$. The finite element spaces for the primal variable $u_{h}$ and the Lagrangian multiplier $\lambda_{h}$ are specified respectively as follows:

[TABLE]

[TABLE]

Denote by $u_{h}\in M_{0,h}$ and $\lambda_{h}=\{\lambda_{0},\lambda_{b}\}\in W_{1,h}$ the numerical solution arising from the PD-WG finite element scheme (3.5)-(3.6) with $j=k=1$. The approximation $u_{h}$ for the primal variable is compared with the exact solution $u$ on each element at the element center - known as the nodal point interpolation $I_{h}u$. Note that the numerical Lagrangian multiplier $\lambda_{h}$ approximates the exact solution $\lambda\equiv 0$. The error functions are denoted as

[TABLE]

The following $L^{2}$ norms are used to measure the error functions in the numerical tests:

[TABLE]

[TABLE]

Tables 1-4 illustrate the performance of the PD-WG finite element scheme (3.5)-(3.6) for the test problem (1.1) with exact solution $u=1$ on the convex domain $\Omega_{1}$ and the non-convex domain $\Omega_{2}$, respectively. The convection vector is given by ${\boldsymbol{\beta}}=[1,-1]$ and the reaction coefficient is $c=1$. The stabilization parameter $\tau$ assumes the values of $\tau=0$ and $\tau=1$, respectively. The numerical results in Tables 1-4 show that the errors are in machine accuracy which is consistent with the theory of $u_{h}=u$ and $\lambda_{h}=0$ for this particular test problem. The purpose of this test is to demonstrate the correctness of the code for the PD-WG algorithm (3.5)-(3.6) with $j=k=1$.

Tables 5 - 8 illustrate the numerical performance of the PW-WG method when $u=\sin(x)\cos(y)$ is employed as the exact solution in the numerical test. The stabilization parameter is taken to be $\tau=1$ and $\tau=0$. The convection vector and the reaction coefficient are respectively given by ${\boldsymbol{\beta}}=[1,-1]$ and $c=1$. Tables 5- 6 show that the convergence for $u_{h}$ in the discrete $L^{2}$ norm on the convex domain $\Omega_{1}$ is at the rate of the optimal order ${\cal O}(h)$ for both $\tau=1$ and $\tau=0$, which is consistent with the theory. Tables 7- 8 indicate that the PD-WG method appears to be convergent at a rate slightly higher than ${\cal O}(h)$ on the L-shaped domain $\Omega_{2}$.

Tables 9-12 show the numerical performance on the unit square domain $\Omega_{1}$ for various values of the stabilization parameter $\tau=0,0.001,1,1000$. The exact solution is given by $u=\sin(\pi x)\sin(\pi y)$, the convection vector is ${\boldsymbol{\beta}}=[1,1]$ and the reaction coefficient is $c=-1$. Tables 9-11 demonstrate that the convergence rates for $e_{h}$ in the discrete $L^{2}$ norm are a bit higher than the expected optimal order ${\cal O}(h)$ for the case of $\tau=0$, $\tau=0.001$, and $\tau=1$. Observe from Tables 9-10 that the absolute error of $e_{h}$ is almost the same as the case of $\tau=0$ and $\tau=0.001$ when compared with the results obtained from the same mesh. The numerical results in Table 12 indicate that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm is much higher than the optimal order ${\cal O}(h)$ for the stabilization parameter $\tau=1000$, while the absolute error of $e_{h}$ in this case is larger than the absolute errors of $e_{h}$ for the case of $\tau=0,0.001,1$ compared in the same mesh. For large values of the stabilization parameter $\tau$, the rate of convergence for the Lagrangian multiplier $\lambda_{h}$ seems to suffer. This test suggests a preference of the PD-WG method with moderate values of the parameter $\tau$.

In Tables 13-16, we report some numerical results for the PD-WG method on the L-shaped domain $\Omega_{2}$ with exact solution $u=\sin(\pi x)\sin(\pi y)$ and various values of $\tau=0,0.001,1,1000$ for the stabilization parameter. The convection vector is given by ${\boldsymbol{\beta}}=[1,1]$ and the reaction coefficient by $c=1$. Tables 13-15 show that the convergence rate for $e_{h}$ in the discrete $L^{2}$ is of ${\cal O}(h)$. Table 16 shows that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm is higher than ${\cal O}(h)$ when the stabilization parameter has the value $\tau=1000$. Note that the absolute error of $e_{h}$ in Table 16 with $\tau=1000$ is bigger than those in Tables 13-15 when compared with the result obtained using the same mesh.

Figure 1 illustrates the numerical performance of the PD-WG method for a test problem with the following configuration: the domain is the unit square $\Omega_{1}$, the exact solution is $u=\sin(x)\cos(y)$, the convection vector is ${\boldsymbol{\beta}}=[y-0.5,-x+0.5]$, the reaction coefficient $c=1$, and the stabilizer parameters $\tau=0,1,10000$. Figure 1 shows that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm is of order ${\cal O}(h^{0.9})$ which is a bit lower than the expected optimal order ${\cal O}(h)$ when $\tau=0$ (left figure) and $\tau=1$ (middle figure) are employed. We conjecture that the slight deterioration on the convergence rate is caused by the rotational nature of the flow. The right figure in Figure 1 indicates that the convergence rate for $e_{h}$ is of ${\cal O}(h^{1.3})$ when $\tau=10000$, which is better than the theoretical prediction.

Tables 17-18 show the numerical results on the L-shaped domain $\Omega_{2}$ for the stabilization parameters $\tau=1$ and $\tau=0$. The exact solution is given by $u=\sin(x)\cos(y)$. The convection vector is given by ${\boldsymbol{\beta}}=[y-1,-x+1]$ and the reaction coefficient is chosen as $c=1$. Tables 17-18 indicate that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm is of ${\cal O}(h^{0.9})$. This slight deterioration of the convergence rate might be caused by the rotational nature of the flow.

In Figure 2, we present some numerical results on the cracked square domain $\Omega_{3}$. The exact solution is chosen to be $u=\sin(x)\sin(y)$. The convection vector is ${\boldsymbol{\beta}}=[y,-x]$ and the reaction coefficient is $c=1$. Figure 2 indicates that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm arrives at an order of ${\cal O}(h^{1.1})$ when the stabilization parameters $\tau=1$ and $\tau=0$ are employed.

Tables 19-20 illustrate the numerical results on the cracked square domain $\Omega_{3}$ when the convection vector is given by ${\boldsymbol{\beta}}=[y,-x]$ and the reaction coefficient by $c=1$. The exact solution is given by $u=\sin(\pi x)\cos(\pi y)$. The numerical results in Tables 19-20 show that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm is of ${\cal O}(h)$ when $\tau=1$ and $\tau=0$.

Figure 3 illustrates the numerical performance of the PD-WG method on the unit square domain $\Omega_{1}$ for the exact solution $u=\sin(x)\cos(y)$. The convection vector ${\boldsymbol{\beta}}(x,y)$ is piece-wisely defined in the sense that ${\boldsymbol{\beta}}(x,y)=[y,-x]$ for $y<1-x$ and ${\boldsymbol{\beta}}(x,y)=[y-1,1-x]$ otherwise; and the reaction coefficient is $c=1$. The stabilization parameter assumes the values of $\tau=1$ and $\tau=0$. The numerical results in Figure 3 show the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm arrives at an optimal order of ${\cal O}(h)$, which is in good consistency with the theory.

Tables 21-24 illustrate the performance of the PD-WG method for the model problem (1.1) with piecewise constant coefficients. The exact solution is chosen as $u=\sin(\pi x)\cos(\pi y)$, and the domain is given by either the unit square $\Omega_{1}$ or the L-shaped domain $\Omega_{2}$. The convection vector ${\boldsymbol{\beta}}(x,y)$ is piece-wisely defined as follows: ${\boldsymbol{\beta}}(x,y)=[1,-1]$ if $y<1-x$ and ${\boldsymbol{\beta}}(x,y)=[-1,1]$ elsewhere. The reaction coefficient is given by $c=1$. The stabilization parameter has values $\tau=1$ and $\tau=0$. Tables 21-22 show that the convergence rate for $e_{h}$ in the discrete $L^{2}$ norm on the unit square domain $\Omega_{1}$ is of ${\cal O}(h)$, which is consistent with the theory. The numerical results in Tables 23-24 suggest that the convergence for $e_{h}$ in the discrete $L^{2}$ norm on the L-shaped domain $\Omega_{2}$ is of order ${\cal O}(h)$.

Fig. 4 illustrates the numerical performance of the PD-WG on the unit square domain $\Omega_{1}$. The convection vector is given by ${\boldsymbol{\beta}}=(\cos(\frac{\pi}{6}),\sin(\frac{\pi}{6}))=(b_{1},b_{2})$ and reaction coefficient is $c=0$. The exact solution is chosen as $u=\left((y-\frac{b_{2}}{b_{1}}x-\frac{1}{2})^{2}+\frac{1}{10}\right)^{-1}$. The stabilization parameter assumes the values of $\tau=1$ and $\tau=0$. The numerical results in Fig. 4 suggest a convergence of $e_{h}$ in the discrete $L^{2}$ norm at the rate of ${\cal O}(h^{1.3})$, which outperforms the expected optimal order of ${\cal O}(h)$.

Figure 5 shows the numerical results on the unit square domain $\Omega_{1}$. The convection vector is given by ${\boldsymbol{\beta}}=(\cos(\frac{\pi}{6}),\sin(\frac{\pi}{6}))=(b_{1},b_{2})$ and the reaction coefficient is given by $c=0$. The exact solution is

[TABLE]

The stabilization parameter assumes the values of $\tau=1$ and $\tau=0$. The numerical results in Fig. 5 suggest a convergence rate of order ${\cal O}(h^{1.3})$ for $e_{h}$ in the discrete $L^{2}$ norm, which is better than the theoretical result of ${\cal O}(h)$.

The rest of the numerical tests shall assume no explicit formulation on the exact solution for the linear transport problem (1.1). The plot of the numerical solution $u_{h}$ arising from the PD-WG scheme (3.5)-(3.6) will be shown for each numerical experiment. To produce a smooth plot, we occasionally apply a post-processing technique to generate a post-processed approximation for the plotting purpose. The post-processed approximation has values on the vertices and the midpoint of each edge for each triangular element $T\in{\cal T}_{h}$. The value of the post-processed approximation at each vertex point is calculated as the simple average of $u_{h}$ on all the elements sharing the same vertex. Similarly, the value of the post-processed approximation at the midpoint of each edge is computed as the simple average of $u_{h}$ on the elements sharing the same edge.

Figures 6-7 show the plots of the numerical solution $u_{h}$ arising from the PD-WG scheme (3.5)-(3.6) on the unit square domain $\Omega_{1}$ with and without using the post-processing technique. The configuration of the test problem is as follows: the convection vector is ${\boldsymbol{\beta}}=[1,-1]$; the reaction coefficient is $c=0$; the stabilization parameter is $\tau=0$; the load function is $f=0$; the inflow boundary data is given by $g=1$ on the inflow boundary edge $\{0\}\times(0,1)$ and by $g=-1$ on the inflow boundary edge $(0,1)\times\{1\}$. The exact solution is known to be $u=1$ for $x<1-y$ and $u=-1$ otherwise. The left ones in Figures 6-7 are the surface plots for the numerical solution and the right ones are the contour plots of the numerical solution. Figures 6-7 show that the numerical solution $u_{h}$ obtained from the PD-WG scheme (3.5)-(3.6) is consistent with the exact solution $u$ of the linear transport problem (1.1).

In what follows of this section, all the plots are based on the post-processed approximations for the PD-WG solution $u_{h}$. Fig. 8 illustrates the contour plots of the numerical solution $u_{h}$ on the unit square domain $\Omega_{1}$ with the following configuration: the convection vector ${\boldsymbol{\beta}}(x,y)$ is piece-wisely defined such that ${\boldsymbol{\beta}}(x,y)=[y+1,-x-1]$ if $y<1-x$ and ${\boldsymbol{\beta}}(x,y)=[y-2,2-x]$ otherwise, the reaction coefficient is $c=0$, and the inflow boundary data is given by $g=\cos(5y)$. The stabilization parameter is set as $\tau=0$. Fig. 8 presents the contour plots for the post-processed numerical solution $u_{h}$ with the load function $f=1$ and $f=0$, respectively.

Fig. 9 shows the contour plots of the numerical solution $u_{h}$ on the unit square domain $\Omega_{1}$ for a test problem with the following configuration: the convection vector ${\boldsymbol{\beta}}=[y-0.5,-x+0.5]$, the reaction coefficient $c=1$, the inflow boundary data $g=\cos(y)$, and the stabilization parameter $\tau=0$. The left one in Fig. 9 demonstrates the contour plot of the numerical solution corresponding to the load function $f=10000$ and the right one is the contour plot of the numerical solution with load function $f=0$.

Fig. 10 shows the contour plots of the numerical solution $u_{h}$ resulting from the PD-WG scheme on the cracked square domain $\Omega_{3}$ when the load function is given by $f=0$ and $f=10000$. The configuration of this test problem is as follows: (1) the convection vector is ${\boldsymbol{\beta}}=[y,-x]$, (2) the reaction coefficient is $c=0$, (3) the inflow boundary data is $g=\sin(x)$, and (4) the stabilization parameter is $\tau=0$.

Figure 11 shows the contour plots of the numerical solution $u_{h}$ of the PD-WG method on the L-shaped domain $\Omega_{2}$. The test problem has the following configuration: (1) the load function is given by $f=0$ and $f=10000$, respectively, (2) the convection vector is ${\boldsymbol{\beta}}(x,y)=[1,-1]$ for $y<1-x$ and ${\boldsymbol{\beta}}(x,y)=[-1,1]$ elsewhere, (3) the reaction coefficient is $c=1$, (4) the inflow boundary data is $g=\sin(x)\cos(y)$, and (5) the stabilization parameter is $\tau=0$. The left one in Fig. 11 is the contour plot of the numerical solution corresponding to the load function $f=10000$, and the right one is the contour plot of the numerical solution for the load function $f=0$.

In conclusion, the numerical performance of the PD-WG method (3.5)-(3.6) for solving the linear transport problem (1.1) is consistent with or better than what the theory predicts in earlier sections of this paper. The numerical results clearly reveal a convergence of the method at the optimal order. We thus claim that the PD-WG finite element method is a stable, convergent, and practically useful numerical scheme for linear transport problems.