2020SummerWorkshopforScientiﬁcComputingSep.1,2020 Min-HungChen DiscontinuousGalerkinMethods

(1)

Discontinuous Galerkin Methods

Min-Hung Chen

Department of Mathematics, NCKU

2020 Summer Workshop for Scientific Computing Sep. 1, 2020

(2)

Outline

1 Introduction DG for ODE

2 Linear Hyperbolic System

3 Steady-State Problems Elliptic Problems Fenics and Colab

4 HDG

Elliptic Model Problem

Matlab tools for HDG - HDG3D

5 Summary

M.-H. Chen Discontinuous Galerkin Methods 2 / 47

(3)

Introduction

What is the Discontinuous Galerkin Method?:

a method between a finite element and a finite volume method local to the generating element

a practical framework for the development of high-order accurate methods using unstructured grids

(4)

Introduction

A short (and biased) historical overview of the DGM

¹

First DG method introduced in 1973 by Reed and Hill for linear transport.

First studied in 1974 by Lesaint and Raviart.

Extended to nonlinear hyperbolic conservation laws in the 90’s by B.

Cockburn and C.-W. Shu.

Extended to compressible flow in 1997 first by F. Bassi and S. Rebay.

New DG methods for diffusion appear and some old ones (the IP

methods of the late 70’s) are resuscitated. A unified analysis is proposed in 2002 by D. Arnold, F. Brezzi, B. Cockburn and D. Marini.

They clash with the well-established mixed and continuous Galerkin methods. In response, the HDG methods are introduced in 2009 by B.

Cockburn, J. Gopalakrishnan and R. Lazarov. The HDG methods are strongly related to the hybrid methods and to the hybridization

techniques of the mid 60’s introduced as implementation techniques for mixed methods.

1https://www.lacan.upc.edu/dg2017/

(5)

Notation

Jump Operator:

[[w]] =







w|_K+n⁺+ w|_K−n⁻ on ∂K⁺∩ ∂K⁻

wn on ∂Ω

[[v]] =







v|_K+· n⁺+ v|_K−· n⁻ on ∂K⁺∩ ∂K⁻

v · n on ∂Ω

Average Operator:

{w} =







1

2(w|_K++ w|_K−) on ∂K⁺∩ ∂K⁻

w on ∂Ω

{v} =



 1

2(v|_K++ v|_K−) on ∂K⁺∩ ∂K⁻

(6)

DG for ODE

ODE

²

Equation

d

dtu(t) = f (t), t ∈ (0, T ), u(0) = u0. Partition of (0, T ) = {tⁿ}^N_n=0, Iⁿ= (tⁿ, tⁿ⁺¹):

Local Basis: {φⁿ_j}^k_j=0;(for example, φⁿ_j(t) = P_j(^2t−t_t_n+1ⁿ⁺¹_−t^−t_nⁿ)) weak form:

Z

Iⁿ

d

dtu_h(s)v(s) ds = Z

Iⁿ

f (s)v(s) ds, or

− Z

Iⁿ

u_h(s)d

dtv(s) ds + u_hv|^t_tⁿ⁺¹n = Z

Iⁿ

f (s)v(s) ds, DG form:

− Z

Iⁿ

uh(s)d

dtv(s) ds +uchv|^t_tⁿ⁺¹n = Z

Iⁿ

f (s)v(s) ds,

2B. Cockburn, Discontinuous Galerkin methods, ZAMM Z. Angew. Math. Mech. 83

(7)

Numerical Flux

DG form: −R

Iⁿuh(s)_dt^dv(s) ds +cuhv|^t_tⁿ⁺¹n =R

Iⁿf (s)v(s) ds, Upwinding Numerical Flux u_ch(tⁿ) =

u0, if tⁿ= 0, uh(tⁿ⁻), otherwise.

Assuming Numerical Solution: uh|_In(t) =Pk

i=0aⁿ_iφⁿ_i(t) Numerical Scheme: On Iⁿ, for j = 0, ..., k,

− Z

Iⁿ k

X

i=0

aⁿ_iφⁿ_i(s)d

dtφⁿ_j(s) ds +

k

X

i=0

aⁿ_iφⁿ_iφⁿ_j(tⁿ⁺¹) −

k

X

i=0

aⁿ⁻¹_i φⁿ⁻¹_i φⁿ_j(tⁿ)

= Z

Iⁿ

f (s)φⁿ_j(s) ds.

Matrix Equation: M





 aⁿ₀

.. . aⁿ





= F

(8)

DG for ODE

Properties of the DG methods

Discontinuous from element to element.

Locally conservative

uch|^t_tⁿ⁺¹n =R

Iⁿf (s) ds,

Relation between the residuals and jumps

− Z

Iⁿ

uh(s)d

dtv(s) ds +uchv|^t_tⁿ⁺¹n = Z

Iⁿ

f (s)v(s) ds, Integration by parts,

Z

Iⁿ

d

dtu_h(s)v(s) ds + (uc_h− u_h)v|^t_tⁿ⁺¹n = Z

Iⁿ

f (s)v(s) ds, Or,

Z

Iⁿ

R(s)v(s) ds = (uh−uch)v|^t_tⁿ⁺¹n

Take v = 1,

Z

Iⁿ

R(s) ds = [[uh]] (tⁿ)

(9)

Example

d

dtu = e^tu(t), t ∈ (0, 2), u(0) = 1.

dg ode.zip

Rate of Convergence and super-convergence (P³) N L²-error order L^∞node order

2 0.50E+02 - 0.45E+02 -

4 0.66E+01 2.92 0.11E+01 5.33 8 0.63E-00 3.40 0.14E-01 6.34 16 0.44E-01 3.83 0.11E-03 6.93 32 0.28E-02 4.00 0.87E-06 7.01 64 0.17E-03 4.03 0.68E-09 7.02

(10)

DG for ODE

Numerical Solutions of P

³

(Left) and P

⁵

(Right)

(11)

Relation between the residuals and jumps

Ztⁿ⁺¹ tⁿ

R(s) ds = [[uh]] (tⁿ)

Interval 1 Interval 2

k Error [[u]] Ratio Error [[u]] Ratio

1 0.33E-01 0.19E+00 0.17 0.19E-00 0.11E+01 0.17 2 0.35E-02 0.26E-01 0.13 0.27E-01 0.20E-00 0.14 3 0.34E-03 0.31E-02 0.11 0.34E-02 0.30E-01 0.12 4 0.32E-04 0.33E-03 0.097 0.41E-03 0.40E-02 0.10 5 0.28E-05 0.32E-04 0.087 0.45E-04 0.50E-03 0.090 6 0.24E-06 0.30E-05 0.080 0.47E-05 0.57E-04 0.082

Interval 3 Interval 4

k Error [[u]] Ratio Error [[u]] Ratio

1 0.19E-00 0.11E+02 0.17 0.73E+02 0.96E+02 0.76 2 0.43E-01 0.30E+01 0.14 0.29E+02 0.17E+03 0.18 3 0.72E-02 0.60E-00 0.12 0.46E+01 0.37E+02 0.12 4 0.11E-03 0.11E-01 0.10 0.95E-00 0.87E+01 0.11 5 0.16E-04 0.17E-02 0.094 0.18E-00 0.18E+01 0.098

(12)

Linear scalar equation

Advection Equation

Equation (Ω ⊂ R^d)

ut+ ∇ · (au) = 0in Ω × (0, T ), u(t = 0) = u0,on Ω.

Triangulation: Th

Local Space: V (K) for K ∈ Th

weak form: For all v ∈ V (K)

Z

K

(uh)tv dx + Z

K

∇ · (au_h)v dA = 0,

or _Z

K

(uh)tv dx − Z

K

auh· ∇v dx + Z

∂K

auh· nv ds = 0,

DG form: _Z

K

(uh)tv dx − Z

K

auh· ∇v dx + Z

∂Kaudh· nv ds = 0,

(13)

Numerical Flux and Stability

u_t+ ∇ · (au) = 0,

Stability for the transport equation 1

2 Z

R^d

u²(x, T ) dx +1 2

Z T 0

Z

R^d

∇ · a(x)u²(x, t) dx dt = 1 2

Z

R^d

u²₀(x) dx Stable if ∇ · a ≥ 0

(14)

DG form

Z

K

(uh)tv dx − Z

K

auh· ∇v dx + Z

∂K

audh· nv ds = 0, Stability for the DG methods (Take v = u_h)

1 2 Z

R^d

u²_h(x, T ) dx +1 2

Z T 0

Z

R^d

∇ · a(x)u²_h(x, t) dx dt + Z T

0

Θh(t) dt =1 2 Z

R^d

u²_h,0(x) dx,

where

Θ_h(t) = X

K∈T_h

(1 2 Z

K

∇ · (a(x)u²_h)(x, t) dx) + Z

∂K

aud_h(x, t) · nu_h(x, t) ds).

= X

K∈T_h

Z

∂K

(audh· nu_h(x, t) −1

2au²_h· n) ds.

= X

e∈E_h

Z

e

(audh− a{u_h}) · [[u_h]] ds.

Numerical Flux: (Θh(t) ≥ 0)

I General: audh(tⁿ) = a{uh} + C [[u_h]]

I Up-winding: C =¹₂|a · n|Id ⇒audh(tⁿ) = a lim↓0uh(x − a)

I Lax-Friedrichs: C =¹₂|a|Id ⇒audh(tⁿ) = a{uh} +¹₂|a| [[u_h]]

(15)

DG Scheme

DG form _Z

K

(uh)tv dx − Z

K

auh· ∇v dx + Z

∂K

audh· nv ds = 0,

audh(tⁿ) = a{uh} + C [[uh]]

Assuming Numerical Solution: uh|K(t, x) = a_i(t)φ_i(x) Numerical Scheme: On K, for j = 0, ..., k,

Z

K

(ai)tφi(x)φj(x) dx − Z

K

aaiφi(x) · ∇φj(x) dx + Z

∂K

aud_h· nφj(x) ds = 0,

ODE

Md

dtU = N U

(16)

Example

ut+ ux= 0 in (0, T = 1) × (0, 1) u(t = 0) = sin 2πx on (0, 1) Rate of Convergence (P⁴with RK 5)

N L²-error order 10 0.17E-05 - 20 0.52E-07 5.01 40 0.16E-08 5.01 80 0.51E-10 5.00 160 0.16E-11 5.00

(17)

Procedures to implement the RK-DG methods

Conservation Form:

u_t+ ∇·f (u) = 0, DG space discretization

Z

K

(u_h)_tv_hdx − Z

K

f (u_h) · ∇v_hdx + Z

∂K

f (u[_h) · n_Kv_hds = 0.

Here, the proper definition of [f (uh)is essential for the stability and convergence of the method.

RK time discretizatin:

d

dtuh= L(uh)

(18)

Linear System

Wave Equation

Equation

utt− c²4u = 0in R^d× (0, T ).

First-order system

Ut+ ∇ · F (U ) = 0, in R^d× (0, T ),

where

U =





 q1

.. . qd

u







, F (U ) = −c







u . . . 0 ..

. . .. ... 0 . . . u q1 . . . qd





 .

Triangulation: Th

Local Space: U (K) = P^k(K) × . . . × P^k(K)for K ∈ Th

DG form: For all V ∈ U (K)

Z

K

(Uhi)tVi− Z

K

Fij(Uh)Vi,j+ Z

∂K

FcijnjVidx = 0,

(19)

Wave Equation

Equation

utt− c²4u = 0in R^d× (0, T ).

First-order system

Ut+ ∇ · F (U ) = 0, in R^d× (0, T ),

where

U =





 q1

.. . qd

u







, F (U ) = −c







u . . . 0 ..

. . .. ... 0 . . . u q1 . . . qd





 .

Triangulation: Th

Local Space: U (K) = P^k(K) × . . . × P^k(K)for K ∈ Th

DG form: For all V ∈ U (K)

Z

K

(Uhi)tVi− Z

K

Fij(Uh)Vi,j+ Z

∂K

FcijnjVidx = 0,

(20)

Linear System

Numerical Flux and Stability

DG form: for all V ∈ U (K)

Z

K

(Uhi)tVi− Z

K

Fij(Uh)Vi,j+ Z

∂K

FcijnjVidx = 0,

Stability for the DG methods (Take V = Uh)

1 2 Z

R^d

U_h²(x, T ) dx +1 2

Z T 0

Θh(t) dt =1 2 Z

R^d

U_h²(x, 0) dx,

where

Θh(t) = X

e∈E_h

Z

e

( cFij− {Fij}) · [[Uh, i]] dx.

Numerical Flux

Fcij= {Fij} + C_ijkl[[Uhk]]_l

(21)

Examples of the DG methods

Numerical Flux

Fcij= {Fij} + Cijkl[[Uhk]]_l

Up-winding:

Fcij = {Fij} + |c|

2 [[q]] δij+|c|

2 [[u]]_jδ_i(d+1) Lax-Friedrichs:

Fcij= {Fij} +|c|

2 [[qi]]_jδij+|c|

2 [[u]]_jδ_i(d+1) Generalization of the Up-winding flux :

Fcij= {Fij} + (C22[[q]] − C12· [[u]])δij+ (C12[[q]] + C11· [[u]])δij

(22)

Linear System

RK-DG Scheme

DG form _Z

K

(Uhi)tVi− Z

K

Fij(Uh)Vi,j+ Z

∂K

FcijnjVidx = 0,

Fcij= {Fij} + Cijkl[[Uhk]]_l

Numerical Solution: Uh|K(t, x) = a_i(t)φ_i(x) Numerical Scheme: On K, for j = 0, ..., k,

Z

K

(ai)tφi(x) · φj(x) dx − Z

K

F (aiφi(x)) · ∇φj(x) dx + Z

∂K

F · nφb j(x) ds = 0,

ODE

Md

dtU = N U

(23)

Numerical Experiments

1D Test: Wave equation with sound speedc(x) = ^√¹_µ

∂_tu − c(x)²∂_xv = 0, in (0, T ) × (−2, 2)

∂tv − ∂xu = 0, in (0, T ) × (−2, 2) c(x) =

2 if x ∈ (−1, 1) 1 otherwise Initial conditions:

u(x, t = 0) = φ(x), x ∈ (−2, 2) v(x, t = 0) = −φ(x), x ∈ (−2, 2) Boundary Condition: (Transparent BC)

u + v =0 at x = 2 u − v =2φ(x − t) at x = −2

(24)

Linear System

Profile

φ(x) = ψ(x + 1 + 0.5 0.5 ), ψ(y) =

(2y − 1)¹⁰(2y + 1)¹⁰ if |y| < .5,

0 otherwise.

φ(x − 1)

(25)

Rate of Convergence

h L^∞-error order L²-error order L¹-error order 0.1250 0.64E-01 0.00 0.13E-01 0.00 0.25E-01 0.00 0.0625 0.17E-02 5.25 0.34E-03 5.28 0.63E-03 5.32 0.0312 0.60E-04 4.81 0.52E-05 6.01 0.86E-05 6.19 0.0156 0.40E-05 3.90 0.26E-06 4.31 0.34E-06 4.68 0.0078 0.25E-06 3.97 0.16E-07 4.00 0.21E-07 4.02 0.0039 0.16E-07 3.99 0.10E-08 4.00 0.13E-08 4.00

P³-elements and SSP-RK4 scheme at T = 1.5

h L^∞-error order L²-error order L¹-error order 0.1250 0.92E-02 0.00 0.19E-02 0.00 0.38E-02 0.00 0.0625 0.12E-03 6.26 0.12E-04 7.40 0.23E-04 7.40 0.0312 0.37E-05 5.04 0.22E-06 5.70 0.26E-06 6.43 0.0156 0.12E-06 4.97 0.70E-08 5.00 0.76E-08 5.11 0.0078 0.37E-08 5.00 0.22E-09 5.00 0.24E-09 5.00 0.0039 0.12E-09 4.99 0.68E-11 5.00 0.74E-11 5.00

P⁴-elements, SSP-RK5 scheme at T = 1.5.

(26)

Elliptic Problems

Elliptic Equation

Consider a second-order elliptic model problem:

−∆u = f in Ω, u = uD on ∂Ω.

Introduce an auxiliary variable q and rewrite the equation as q = ∇u in Ω,

−∇·q = f in Ω, u = u_D on ∂Ω.

(27)

Definition

For each element K of the mesh T_h of the domain Ω, we define where n_K is the outward unit normal to K,

(u, v)K = Z

K

uv dx, (1)

hw, vi_∂K = Z

∂K

wv ds. (2)

(28)

Weak form and DG method

The approximate solution (q_h, uh)on the element K is taken in the space V (K) × W (K) and is defined as the solution, for all

(v, w) ∈ V (K) × W (K), of the equations

(q_h, v)_K+ (u_h, ∇·v)_K− hub_h, n_K· vi_∂K = 0, (q_h, ∇w)_K− hqb_h· n_K, wi_∂K = (f, w)_K, with Dirichlet boundary condition

buh= uD on ∂K ∩ ∂Ω

All the DG methods are generated by choosing the local spaces V (K) × W (K)and the numerical tracesqb_h· n_Kandub_h.

(29)

Numerical traces

The definition of the numerical tracesqb_h andub_hstrongly influences the properties of the corresponding DG method. In this context, we also require that the numerical traces be linear functions of the traces of q_h· n_K and u_h which are consistent and single valued.

Example (Cockburn and Shu 1998)

qb_h := {q_h} − C_qq[[q_h]] − C_qu[[u_h]] , ub_h := {u_h} − C_uu[[u_h]] − Cuq[[q_h]] .

Let us consider DG methods having a numerical tracebuh independent of q_h. This allows for the easy elimination of the variable q_h, which can now be expressed in terms of u_hin an elementwise manner, and results in the so-called primal formulation of the method.

(30)

Mixed Form of LDG Methods Find (q_h, u_h) ∈ M_h× V_hsuch that

a_h(q_h, v) + b_h(u_h, v) = G_h(v),

−b_h(w, qh) + ch(uh, w) = Fh(w),

The corresponding linear system has the form

A B

−B^t C

Q U

=

G F

Aand C are symmetric B is antisymmetric.

(31)

Numerical Result: Elliptic Interface Problem

Model Problem

−∇ · β∇u = f in Ω1∪ Ω2, u = gon ∂Ω, u|_Ω₁− u|Ω2 = aon ΓI, ((β∇u)|Ω₁− (β∇u)|Ω₂) · ne = bon ΓI,

where f , g, a, and b are functions of x and y, neis the outward unit normal vector to ∂Ω1, and β is a positive finite constant function on Ω1

and Ω2, separately.

(Chen, Wu, TJM2016 )

(32)

The Model Problem and LDG method

We assume that Ω = [−1, 1] × [0, 3] is the entire domain and Ω1is the open interior domain embedded in Ω with a complicated interface

ΓI(θ) =

0.6 cos θ − 0.3 cos 3θ 1.5 + 0.7 sin θ − 0.07 sin 3θ + 0.2 sin 7θ

for θ ∈ [0, 2π]. The exterior domain is Ω2= Ω \ Ω1.

Figure:Domain(left) and numerical solution (right) of test problem 2.

The media property is given by β =

1 in Ω1,

10 in Ω2. The discontinuous exact solution is u(x, y) =

e^x(y²+ x²sin y) in Ω1,

−(x²+ y²) in Ω2.

The corresponding jump conditions a and b, Dirichlet condition g, and source term f can be derived easily.

(33)

The Model Problem and LDG method

Table:p-convergence of u and q (25 elements) for test problem 2.

Basis function DoF L²-error of u L²-error of q Q¹ 300 6.615e − 01 2.641e + 01 Q² 675 1.599e − 01 1.129e + 01 Q³ 1200 4.450e − 02 4.346e + 00 Q⁴ 1875 6.711e − 03 9.273e − 01 Q⁵ 2700 2.652e − 03 4.010e − 01 Q⁶ 3675 3.512e − 04 7.475e − 02 Q⁷ 4800 1.095e − 04 2.204e − 02 Q⁸ 6075 2.148e − 05 5.815e − 03 Q⁹ 7500 3.297e − 06 7.967e − 04 Q¹⁰ 9075 9.766e − 07 3.258e − 04 Q¹¹ 10800 8.334e − 08 2.376e − 05 Q¹² 12675 3.172e − 08 1.316e − 05 Q¹³ 14700 3.406e − 09 1.621e − 06

(34)

p-convergence and Domain

Figure:p-convergence of u and q on a 25-element mesh (left). Corresponding curved-edge quadrilateral mesh with 25 elements (right).

(35)

Fenics Project

The FEniCS Project is a collection of free and open-source software components with the common goal to enable automated solution of differential equations.³

The components provide scientific computing tools for working with computational meshes, finite-element variational formulations of ordinary and partial differential equations, and numerical linear algebra.

FEniCS enables users to quickly translate scientific models into efficient finite element code. With the high-level Python and C++

interfaces to FEniCS, it is easy to get started. ⁴ FEniCS runs on a multitude of platforms.

3https://en.wikipedia.org/wiki/FEniCS Project

(36)

Fenics and Colab

Reference:

Solving PDEs in Python: The FEniCS Tutorial I (Free)

Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book (Free PDF)

Introduction to Automated Modeling with FEniCS by Ridgway Scott

FENICS Examples (John Burkardt)

(37)

Poisson problem (CG):R

Ω∇u · ∇v dx =R

Ωf v dx

%m a t p l o t l i b i n l i n e from f e n i c s import ∗

mesh = UnitSquareMesh ( 8 , 8 )

V = FunctionSpace ( mesh , ” Lagrange ” , 1 )

u0 = E x p r e s s i o n ( ” 1+x [ 0 ] ∗ x [ 0 ] + 2 ∗ x [ 1 ] ∗ x [ 1 ] ” , degree =2) bc = D i r i c h l e t B C ( V , u0 , ” on boundary ” )

f = Constant ( − 6 . 0 ) u = T r i a l F u n c t i o n ( V ) v = T e s t F u n c t i o n ( V )

a = i n n e r ( grad ( u ) , grad ( v ) ) ∗ dx L = f ∗ v ∗ dx

u = F u n c t i o n ( V ) s o l v e ( a == L , u , bc ) p l o t ( u )

(38)

Fenics and Colab

Google Colaboratory (Colab)

What is Colaboratory?

Colaboratory, or ”Colab” for short, allows you to write and execute Python in your browser, with

Zero configuration required Free access to GPUs Easy sharing

Some sample code:

MyPoissonCG link: CG method for Poisson problem.

MyPoissonDG link: DG method for Poisson problem.

Reference: https://github.com/leodenale/FenicsOnColab

(39)

Motivation for HDG methods

However, the DG methods (for second-order elliptic equations) have been criticized because:

For the same mesh and the same polynomial degree, the number of globally coupled degrees of freedom of the DG methods is much bigger than those of the CG method. Moreover, the orders of convergence of both the vector and scalar variables are also the same.

(40)

Elliptic Model Problem

Elliptic Equation

Consider a second-order elliptic model problem:

cq + ∇u = 0 in Ω,

∇·q = f in Ω, bu = u_D on ∂Ω.

Here c is a matrix-valued function which is symmetric and uniformly positive definite on Ω

(41)

The HDG methods

⁵

.

Rewrite the equations. If for each K ∈ Ω_h, we assume that we know the traceubon ∂K, we can obtain (q, u) inside K as the solution of Consider a second-order elliptic model problem:

cq + ∇u = 0 in K,

∇·q = f in K, u = ub on ∂K.

Thenubcan be determined as the solution, on each edge F ∈ E_h, of [[q · n]]b = 0 if F ∈ E_h^o, (interior edges)

ub = u_D if F ∈ E_h^∂, (exterior edges) Note: [[q · n]] :=b bq⁺· n⁺+bq⁻· n⁻.

5The Hybridizable Discontinuous Galerkin Methods, Proceedings of the

(42)

The local problems: A weak formulation on each element.

On each element K ∈ Ω_h, we define (q_h, u_h) ∈ V (K) × W (K)in terms of (ub_h, f )such that

(cq_h, v)K− (u_h, ∇·v)K+ hub_h, v · ni_∂K = 0,

−(q_h, ∇w)K+ hqbh· n, wi_∂K = (f, w)K, for all (v, w) ∈ V (K) × W (K), where

qb_h· n = q_h· n + τ (u_h−ub_h) on ∂K.

Note: (u, w)_K:=R

Kuw dx, hw, vi_∂K =R

∂Kwv ds

(43)

The global problem: The weak formulation for u b

_h

.

For each face F ∈ E_h^o, we takeub_h|_F in the space M (F ). We determine ubhby requiring that,

hν, [[qb_h]]i_F = 0, ∀ν ∈ M (F ) if F ∈ E_h^o, bu_h = u_D if F ∈ E_h^∂.

All the HDG methods are generated by choosing the local spaces V (K), W (K), M (F ) and the stabilization function τ .

(44)

Notice thatbu_h is the data of the local problems but is the unknown of the global problem. So, the only globally-coupled degrees of freedom are those ofbuh.

By solving the local problems, we express q_h, u_handqb_h in terms ofubhand f . With these expressions, we construct the matrix equation associated to the global problem.

After solving it, we can insert the actual values ofbu_hin the

expressions we had obtained for q_h, u_h andqb_h. Next, we describe this procedure more precisely.

(45)

The main features of the HDG methods.

The HDG methods are obtained by discretizing characterizations of the exact solution written in terms of many local problems, one for each element of the mesh Ω_h, with suitably chosen data, and in terms of a single global problem that actually determines them.

This permits an efficiently implementation since they inherit the above-mentioned structure of the exact solution. This is what renders them efficiently implementable, especially within the framework of hp-adaptive methods, as is typical of DG methods.

(46)

Matlab tools for HDG - HDG3D

HDG3D

Matlab implementation of the Hybridizable Discontinuous Galerkin method on general tetrahedrizations of polyhedra in three

dimensional space.

Developed by group Team Pancho at the Department of Mathematical Sciences at the University of Delaware.

Project website: https://team-pancho.github.io/HDG3D/

Remark:

Move “A simple example.m” up one level.

Move the directory “meshes” up one level.

(47)

Summary and Reference

Main features of the DG methods:

High-order accurate.

Locally Conservative.

Adaptivity

High parallelizability,

The HDG methods are obtained by constructing discrete versions (based on discontinuous Galerkin methods) of the above

characterization of the exact solution.

In this way, the globally coupled degrees of freedom will be those of the corresponding global formulations.

(48)

Reference

B. Cockburn, G. Karniadakis, C.-W. Shu, The development of

Discontinuous Galerkin methods, in Discontinuous Galerkin methods.

Theory, computation and applications, Lecture Notes in Computational Sicence and Engineering, Volume 11, Springer, 2000.

B. Cockburn and C.-W. Shu, Runge-Kutta Discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001), pp. 173-261.

D. Arnold, F. Brezzi, B. Cockburn and D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SINUM 39 (2002), pp. 1749-1779.

B. Cockburn, Discontinuous Galerkin methods, ZAMM Z. Angew. Math.

Mech. 83 (2003), pp. 731-754.

B. Cockburn, Discontinuous Galerkin methods for computational fluid dynamics, Encyclopedia of Computational Mechanics Second Edition (2018): 1-63.