MATLAB Tools for HDG in Three Dimensions
===
---
Student : 戴嘉展, 張泳倉
Adivisor : 陳旻宏

---
## Introduction
Due to the DG methods were criticized for having too many degree of freedom and for not being easy to implement, the Hybridizable Discontinuous Galerkin (HDG) method was promoted to response these criticisms.

The HDG method was intorduced in the framework of diffusion problems. The HDG have the advantage of not using degrees of freedom to stabilize the discrete equations while keeping equal optimal order of vonvergence in all computed fields.

---
## reaction-diffusion problem
Let $\Omega$ ⊂ $R^3$ be a polyhedron with boundary $\Gamma_{D}$, divided into a Dirichlet and a Neumann part $\Gamma_{D}$ and $\Gamma_{N}$ such that each part is the union of faces of $\Gamma$. The unit outward-pointing normal vector field on $\Gamma$ is denoted $v$.
\begin{align*}
\kappa^{-1} q+\nabla u=0\\
\nabla \cdot q+cu=f\\
u=u_D\\
-q\cdot v=g_{N}\cdot v\\
\end{align*}

---
## rewrite the equation
Thw HDG method seeks to find approximation solution$(q_h,u_h,\hat u)$ which is determined as the only solution of the following weak statements:

$$(\kappa^{-1} q_{k},r)_K-(u_K,\nabla \cdot r)_K+(\hat{u}_h,r\cdot v_K)_{\partial K}=0$$
$$(\nabla \cdot q_K,w)_K+(cu_K,w)_K+<\tau_K(u_K-\hat{u_h}),w>_{\partial K}=(f,w)_K$$
$$<u,z>_{\partial K}=<u_D,x>_{\partial K}$$
$$-\mathbf{q}\cdot v=g_{N}\cdot v$$

---
## Galerkin formulation and the numerical traces
flux equilibrium on internal faces:
$<q_{K}\cdot v_{K}+\tau(u_{K}-\hat{u}|_{e}),\hat{v}>=<q_{K}\cdot v_{K}+\tau(u_{K}-\hat{u}|_{e}),v>$

normal numerical flux:
$\Phi :=q_{K}\cdot v_{K}+\tau(u_{\tilde{K}}-\hat{u}|_{e})$

The function $\tau_{K}$ appears in the local solver is called stabilization function which is a nonnegative piecewise constant function on the boundary of each tetrahedrizatoin.
$\tau_{K}|_{e}\in \rho_{0} \qquad \forall e\in \epsilon(K), \qquad \forall K\in \mathcal{T}_{h},\qquad and \qquad \tau_{K}\geq 0.$

---
## local solver

Locally, we will think of $(q_K,u_K)$, while $\hat{u}_K$ will be counted by global face numbering. The solution of PDE(1)-(2) is determined as the only solution $\hat{u_{h}}$ called local solver
The local solvers are related to the pair of first order PDEs.

$(\kappa^{-1} q_{k},r)_K$-$(u_K,\nabla \cdot r)_K$+$(\hat{u}_h,r\cdot v_K)_{\partial K}=0$
$(\nabla \cdot q_K,w)_K$+$(cu_K,w)_K$+$<\tau_K(u_K-\hat{u_h}),w>_{\partial K}$=$(f,w)_K$

---
## global solver
The global solver is used to compute $\hat{u_{h}}$ by boundary conditions. The other variables ($q_{h},u_{h}$) is done by solving local problems. By doing so, we can used $\hat{u_{h}}$ to solve the the global system.

---

# Solve Advection aquation with FR/CPR method

---

Student : 林伯鴻, 馬裕傑
Adivisor : 陳旻宏

---

[TOC]

---

## Introduction

---

Before starting solve traffic flow with FR/CPR method, we're going to practice solving advection with FR/CPR method.
In the next few parts,we will introduce advection equation and FR/CPR method and the choice of numerical flux.
After the introduction,we will compute error of the method with python code.
The figure and the form will be shown as the end.

---

## Advection

---

If the substance is simply carried along (advected) in a flow at some constant velocity $a$,then the flux function is $$f(u)=au.$$The local density $u(x,t)$ (in gams/meter, say) multiplied by the velocity (in meter/sec, say) gives the flux of material past the point $x$ (in grams/sec).

Since the total mass in $[x_1,x_2]$ changes only due to the flux at the endpoints,we have$$\dfrac{d}{dt}\int_{x_1}^{x_2}u(x,t)dx=f(u(x_1,t))-f(u(x_2,t)).$$The minus sign on the last term comes from the fact that $f$ is,by definition,the flux to the right.

If we assume that $u$ and $f$ are smooth functions,then this equation can be rewritten as $$\dfrac{d}{dt}\int_{x_1}^{x_2}u(x,t)dx=-\int_{x_1}^{x_2}\dfrac{\partial}{\partial x}f(u(x,t))dx.$$or,with some further modification,as$$\int_{x_1}^{x_2}[\dfrac{\partial}{\partial t}u(x,t)+\dfrac{\partial}{\partial x}f(u(x,t))]dx=0.$$Since this integral must be zero for all values of $x_1$ and $x_2$,it follows that integrand must be identically zero. This gives,finally,the differential equation$$\dfrac{\partial}{\partial t}u(x,t)+\dfrac{\partial}{\partial x}f(u(x,t))=0.$$This form of equation is called a $conservation\ law$.

For the case considered before, $f(u)=au$ with $a$ constant and this equation becomes **advection equation**.$$u_t+au_x=0.$$This equation requires initial conditions and possibly also boundary conditionsin order to determine a unique solution. The simplest case is the **Cauchy problem** on $-\infty\ <\ x\ <\ \infty$ (with no boundary),also called the pure initial value problem. Then we need only to specify initial data$$u(x,0)=\eta(x).$$This is the simplest example of a **hyperbolic** equation,and it is so simple that we an write down the exact solution,$$u(x,t)=\eta(x-at).$$One can verify directly that this is the solution. However,many of the issues that arise more generally in discretizing hyperbolic equations can be most easily seen with this equation.

---

## FR/CPR method

### Standard DG method

On each cell $E$ , denote$$(v,w) = \int_Ev(x)w(x)dx$$with any two function $v,w$.

The polynomial $u,f,\phi$ are in $\mathbf{P}_{K-1}$.

Let $\phi$, typically, is one of the basis functions, be the test function.
$(u_t,\phi) = (u,\phi)_t$
$$(u_t,\phi) + (f_x,\phi) = 0$$by IBP
$$(u,\phi)_t + (\hat{f}\phi)_{\partial E} - (f,\phi_x) = 0$$
which is called **weak form**.

Then by IBP again$$(u,\phi)_t + (\hat{f}\phi)_{\partial E} - (f\phi)_{\partial E} + (f_x,\phi) = 0$$ which is called **strong form**.

Here $\hat{f}$ is numerical flux which is the common term of two interface at the node.
We will introduce the numrical flux in next chapter.

### Flux Reconstruction

Switch to the local coordinate $\xi$ in $I=[-1,1]$, denote the **jumps** by$$f_L = \hat{f}(-1)-f(-1)\qquad f_R = \hat{f}(1) - f(1)$$

Noting that $dx = \frac{h}{2}d\xi$ and $f_xdx = f_\xi d\xi$ where $h$ is the width of the cell.
We have$$\frac{h}{2}(u_t,\phi)+f_R\phi(1)-f_L\phi(-1)+(f_\xi,\phi) = 0$$

Now we want to find a polynomial on $I$ which possesses the property that $$-\phi(-1) = (g_L^\prime,\phi)\qquad\qquad(*)$$for any $\phi\in\mathbf{P}_{K-1}$

Note that $g_L^\prime = (g_L)_\xi$ is called **correction function** for left boundary.

If $g_L$ is found, by reflection, $g_R(\xi) = g_L(-\xi)$ and $$\frac{h}{2}(u_t,\phi)+(f_Rg_R^\prime,\phi)+(f_Lg_L^\prime,\phi)+(f_\xi,\phi) = 0$$

IBP on $(*)$$$-\phi(-1)=g_L(1)\phi(1)-g_L(-1)\phi(-1)+(g_L,\phi^\prime)$$Hene $g_L$ must satisfy $$g_L(1)=0,\quad g_L(-1) = 1,\qquad(g_L,\phi^\prime)=0$$

Consider right **Radau polynomial** $R_{R,K}$ of degree $K$, in [1], it satisfy the condition above.

As a result, we choose $R_{R,K}$ to replace $g_L$.

Similarly, $g_R$ must satisfy$$g_R(1)=1,\quad g_R(-1) = 0,\qquad(g_R,\phi^\prime)=0$$and we choose $R_{L,K}$ to replace it.

$R_{R,K} = \frac{(-1)^K}{2}(\mathbf{P}_K-\mathbf{P}_{K-1})$
$R_{L,K} = \frac{1}{2}(\mathbf{P}_K+\mathbf{P}_{K-1})$

With the discuss above, we obtain$$\frac{h}{2}(u_t,\phi)+(f_Rg_R^\prime,\phi)+(f_Lg_L^\prime,\phi)+(f_\xi,\phi) = 0$$
and cancel the test function, we have$$u_t+\frac{2}{h}(f_Rg_R^\prime+f_Lg_L^\prime+f_\xi)=0$$

The recontruction flux $F=f_Rg_R+f_Lg_L+f$

Noticing that $F_x = \frac{2}{h}F_\xi = \frac{2}{h}F^\prime$

From now on, we can deal with $$u_t=-F_x$$
and march in time by Runge-Kutta method.

---

## Numerical flux

### Upwind flux

Define:
$u_h^\pm(x^n)=\lim_{\epsilon\to0}u_h^\pm(x^n\pm\epsilon)$
$\{u_h\}=\frac{1}{2}(u_h^++u_h^-)$ : Average
$[[u_h]]=u_h^--u_h^+$ : Jump

\begin{align*}
\widehat{u}_{h}(t^{n})
& = \Bigg\{ \begin{array}{l}
u_{0}, \quad \text{if} \quad x^{n} = 0, \\
\{u_{h}\} + \frac{1}{2} [[u_{h}]], \quad \text{if} \quad x^{n} \in(0,X),
\\
u_{h}^{-}(X), \quad \text{if} \quad x^{n} =X.
\end{array}
\end{align*}

### Roe flux

At the interface $x^n$
Let $\tilde{u}$ be the value satisfies MVT
Define:
\begin{align*}
a(\tilde{u})=
& \Bigg\{ \begin{array}{l}
\frac{f(u_h^+)-f(u_h^-)}{u_h^+-u_h^-}, \quad \text{if} \quad u_h^+\neq u_h^-, \\
a(u_h^+)=a(u^-_h), \quad \text{otherwise}
\end{array}
\end{align*}
and the Roe flux:$$\hat{f} = \frac{1}{2}[f(u^+_h)+f(u^-_h)]-\frac{1}{2}|a(\tilde{u})|(u^+_h-u^-_h)$$

---

## Numerical results

---

We introduce the way we calculate the error of the numercial solution before we do it. Since the equation depends on both **time** and **space**,we will get result in the form of vector if we just use numercial solution to minus exact solution. For obtaining the order successfully,we want to compare a set of disctrete values with a function.

### Errors

Finite difference methods do not produce a function $U(x)$ as an approximation to $u(x)$ . Instead they produce a set of values $U_i$ at grid points $x_i$.

We must first decide what the values $U_i$ are supposed to be approximating. Often the value $U_i$ is meant to be interpreted as an approximation to the pointwise value of the function at $x_i$,so $U_i\approx u(x_i)$. In this case it is natural to define a vector of errors $e=(e_1,e_2,...,e_N)$ by $$e_i=U_i-u(x_i).$$Once we have defined the vector of errors $(e_1,...,e_N)$,we can measure its magnitude using some norm. Since this is simply a vector with $N$ components,it would be tempting to simply use the vector norm e.g.,$$\|e\|_1=\sum_{t=1}^N\left|e_i\right|.$$However,this choice will give a very misleading ideal of the magnitude of the error. The quantity can be expected to be roughly $N$ times as large as the error at any single grid point and here $N$ is not the dimension of some physically relevant space,but rather the number of points on our grid. If we refine the grid and increase $N$,then the quantity might well **increase** even if the error at each grid point decreases,which is clearly not the correct behavior.

Instead we should define the norm of the error by discretizing the integral,which is motivated by considering the vector $(e_1,...,e_N)$ as a discretization of some error function $e(x)$. This suggests defining$$\|e\|_1=h\sum_{i=1}^N\left|e_i\right|$$with the factor of h corresponding to the $dx$ in the integral. Note that since $h\approx(b-a)/N$,this scales the sum by $1/N$ as the numer of grid points increases,so that $\|e\|_1$ is the average value of $e$ over the interval.The norm will be called a **grid function norm** and is distinct from the related vector norm. The set of values $(e_1,...,e_N)$ will sometimes be called a **grid function** to remind us that it is a special kind of vector that represents the discretization of a function.

Similarly,the q-norm should be scaled by $h^{1/q}$,so that the **q-norm** for grid functions is$$\|e\|_q=\Bigg(h\sum_{i=1}^N\left|e_i\right|^q\Bigg)^{1/q}$$Since $h^1/q\rightarrow1$ as $q\rightarrow\infty$,the **max-norm** remains unhanged,$$\|e\|_\infty=\max_{1\leq i\leq N}\left|e_i\right|$$.

### Analysis

Although we consider the process above,there is still a vector of the error if we calculated it every single time. In order to fix it,we just take the last vector in time to calculate our error. And here is the form
| | $dt$ | $\|e\|_1$ | $\|e\|_2$
--- | :---: | :---: | :---:
N=5 | $\frac{0.01}{5}$|9.314111e-02|8.255185e-02|
N=10 | $\frac{0.01}{10}$|8.486300e-04|1.179292e-03|
N=20|$\frac{0.01}{20}$|1.712560e-06|3.915992e-06|
N=40|$\frac{0.01}{40}$|2.045453e-08|6.399544e-08|
N=80|$\frac{0.01}{80}$|5.114468e-10|1.955048e-09|
N=160|$\frac{0.01}{160}$|1.251203e-10|7.132974e-10|
After we get the norm,we can define the **error ratio**$$R(h)=E(h)\ /\ E(h/2),$$we respect $$R(h)\approx 2^p,$$and hence$$p\approx log_2(R(h)).$$
Here refinement by a factor of 2 is used only as an example,since this choice is often made in practice. But more generally if $h_1$ and $h_2$ are any two grid spacings,then we can estimate $p$ basd on calculations on these two grids using$$p\approx \dfrac{log(E(h_1)/E(h_2))}{log(h_1/h_2)}.$$Then we have
|$p_{\|e_1\|}$|$N_{i-1}$|
:---:|:---:|
$N_1=5$||
$N_2=10$| 6.778139 |
$N_3=20$| 8.952837 |
$N_4=40$| 6.387590 |
$N_5=80$| 5.321692 |
$N_6=160$| 2.031268 |

|$p_{\|e_2\|}$|$N_{i-1}$|
:---:|:---:|
$N_1=5$||
$N_2=10$| 6.129301 |
$N_3=20$| 8.234329 |
$N_4=40$| 5.935265 |
$N_5=80$| 5.032693 |
$N_6=160$| 1.454628 |

### Figure

N = 5

N = 10

N = 20

N = 40

N = 80

N = 160

---

## Python

### Algorithm

1. Define the parameter, initial condition and exact solution.

2. Use Gauss-Legendre quadrature and the properties of Legendre polynomial to find the initial coefficients of initial condition.

3. Evaluate the numerical flux at each interface, then compute the jumps $f_L$ and $f_R$.

4. Obtain $F_x$ and use the corrected slope to calculate $\frac{du}{dt}$ and march in Runge-Kutta method.

### Code

import numpy as np
from scipy.special import legendre
import matplotlib.pyplot as plt
import math

xmin = 0 
xmax = 1 # range of x

a=1 # u_t+au_x=0

Nx = 20 # # of mesh
dx = ( xmax - xmin ) / Nx # width of mesh

deg = 4 # degree of LP
gau = 5

CFL=0.01
dt = dx * CFL

tmin = 0
tmax = 0.25

xc = 0.5 #center of x

def f(x):

    return np.exp( -200*( x - xc )**2 )

def u_exact(x,t):

    return np.exp( -200*( x - xc - a*t )**2 )


def gaulegf(a,b,n):
    x = list(range(n+1))
    w = list(range(n+1))
    eps = 3.0E-14
    m = (n+1)//2
    xm = 0.5*(b+a)
    xl = 0.5*(b-a)
    for i in list(range(1,m+1)):
        z = math.cos(math.pi*(i-0.25)/(n+0.5))
        while True:
            p1 = 1.0
            p2 = 0.0
            for j in list(range(1,n+1)):
                p3 = p2
                p2 = p1
                p1 = ((2.0*j-1)*z*p2-(j-1.0)*p3)/j

            pp = n*(z*p1-p2)/(z*z-1.0)
            z1 = z
            z = z1-p1/pp
            if abs(z-z1)<=eps:
                break

        x[i] = xm-xl*z
        x[n+1-i] = xm+xl*z
        w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
        w[n+1-i] = w[i]
    return x,w

t=np.array(gaulegf(-1,1,gau))
t.reshape(gau*2+2)

s=t[0,1:gau+1]
w=t[1,1:gau+1]

c = np.zeros((Nx,deg + 1))
cp = np.zeros((Nx,gau))

for i in range(Nx):
    xa=xmin+i*dx
    xb=xmin+(i+1)*dx
    
    def smm(k):
        for j in range(gau):
            cp[i,j]= w[j]* f( ((xb-xa)/2)*s[j] + (xa+xb)/2 ) * legendre(k)( s[j] )
        return sum(cp[i,0:gau])

    for m in range(deg+1):
        c[i,m]=((2*m+1)/2)*smm(m) 
        
u = c.transpose()

# 定義基底變換矩陣 T : P' to P
D = np.zeros((deg + 1, deg + 1))
for i in range(deg + 1):
    for j in range(deg + 1):
        if (i+j) % 2 == 1:
            if i<j:
                D[i,j] = 2*i+1


P_M = np.zeros(deg + 1)
P_m = np.zeros(deg + 1)

P_M[deg] = 1     #基底函數展開P_K
P_m[deg-1] = 1   #    ''     P_{K-1}


P_R = P_M + P_m  
P_L = P_M - P_m

g_L = (-1) ** deg / 2 * np.dot(D,P_L)
g_R = 1/2 * np.dot(D,P_R)

P_RB = np.ones(deg + 1)
P_LB = np.ones(deg + 1)

# P_RB = P(1) 、 P_LB = P(-1)

for i in range(deg + 1):
    P_LB[i] = P_RB[i] * ( (-1)**i )

    
dF = np.zeros((deg + 1 , Nx))
f_L = np.zeros(Nx)
f_R = f_L.copy()
flux = np.zeros((Nx+1))
A=np.zeros((Nx+1))
u_stage = np.zeros((deg + 1 , Nx))
dF_rk = np.zeros((deg + 1 , Nx))

while tmin < tmax:

    plt.clf()
 

    if (tmin + dt) > tmax:
        dt = tmax - tmin

    tmin = tmin + dt

    for i in range(1,Nx):
        if  abs(np.dot(u[:,i-1],P_RB)-np.dot(u[:,i],P_LB))<0.00000001:   

            A[i] = a*np.dot(np.dot(D,u[:,i-1]),P_RB)
    
        else:
        
            A[i] =a*(np.dot(u[:,i],P_LB)-np.dot(u[:,i-1],P_RB))/ (np.dot(u[:,i],P_LB)-np.dot(u[:,i-1],P_RB))
        
        flux[i]=(1/2)*(a*(np.dot(u[:,i],P_LB)+np.dot(u[:,i-1],P_RB))-abs(A[i])*(np.dot(u[:,i],P_LB)-np.dot(u[:,i-1],P_RB)))
    
    if  abs(np.dot(u[:,i-1],P_RB)-np.dot(u[:,i],P_LB))<0.00000001:   
        A[0] = a*np.dot(np.dot(D,u[:,Nx-1]),P_RB)
    else:
        A[0] =a*(np.dot(u[:,0],P_LB)-np.dot(u[:,Nx-1],P_RB))/ (np.dot(u[:,0],P_LB)-np.dot(u[:,Nx-1],P_RB))


    flux[0]=(1/2)*(a*(np.dot(u[:,0],P_LB)+np.dot(u[:,Nx-1],P_RB))-abs(A[0])*(np.dot(u[:,0],P_LB)-np.dot(u[:,Nx-1],P_RB)))
    flux[Nx]=flux[0]

    for i in range(Nx):
        f_L[i] = flux[i] -np.dot( a*u[:,i] , P_LB )

    for i in range(Nx):
        f_R[i] = flux[i+1] - np.dot( a*u[:,i] , P_RB )
            
    for p in range(Nx):
        dF[:,p] =( a*np.dot(D,u[:,p]) + f_L[p]*g_L + f_R[p]*g_R )
        
        u_stage[:,p] = u[:,p] + dt * (-2/dx*dF[:,p])
        dF_rk[:,p] = a*np.dot(D,u_stage[:,p]) + f_L[p]*g_L + f_R[p]*g_R
        u_stage[:,p] = (3/4) * u[:,p] + (1/4)*u_stage[:,p]+(1/4)*dt*(-2/dx)*dF_rk[:,p]
        dF_rk[:,p] = a*np.dot(D,u_stage[:,p]) + f_L[p]*g_L + f_R[p]*g_R

    u=(1/3)*u+(2/3)*u_stage+(2/3)*dt*(-2/dx)*dF_rk

te = tmin - int( tmin/1 )

for i in range(Nx):

    x_left = xmin + i*dx
    x_right = xmin + ( i + 1 )*dx

    x = np.linspace( x_left , x_right , 10 )


    def fh(x):

        ff=np.zeros(len([x]))
  
        for j in range(deg+1):
            ff = ff + u[j,i]*legendre(j) ( (2*x/( x_right - x_left ) )-(( x_right + x_left )/( x_right - x_left )) )

        return ff


    y_numerical = fh( x )
    y_exact = u_exact( x , te )

    plt.plot( x , y_exact , 'r' , label = "Exact solution" )
    plt.plot( x , y_numerical , 'bo-' , label = "Flux Reconstruction" )
    plt.axis( ( xmin - 0.12 , xmax + 0.12 , -0.2 , 1.4 ) )
    plt.grid( True )
    plt.xlabel( "Distance(x)" )
    plt.ylabel( "u" )
    plt.suptitle( "Time = %1.3f" % tmin )
  
plt.pause( 0.01 )


error = abs( y_numerical - y_exact )
e_1norm=dx*error.sum()
e_2norm = ( dx*(error**2).sum())**(1/2)

---

## Reference

[1] High-Order Methods Including Discontinuous Galerkin
by Reconstructions on Triangular Meshes , H. T. Huynh

[2] A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods, H. T. Huynh

[3] A Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin for Diffusion, H. T. Huynh

[4] https://hackmd.io/SVs6KZ8vR2ivQ6Ckq3XSAA?both

[5] https://hackmd.io/@2njaaOJLSAG_5EQNoHYykw/Sy32jbBUl?type=view

[6] https://www.mathworks.com/matlabcentral/fileexchange/44628-flux-reconstruction-fr-cpr?fbclid=IwAR2nnuJydMSeWR_i_iVVOLwJ4rxkBkidbGhNbkR7OG9Vq0NuaEC2oD5gwTQ

[7] R. J. LeVeque. $Finite\ Difference\ Methods\ for\ Ordinary\ and\ Partial\ Differential\ Equations.$ 2007.(cited on 201,202,210,251,252,255,256,315)