Hybridized Discontinuous Galerkin (HDG)

This part was treated at the end of the internship, and it was harder for me to understand. There may be some imprecision left on this page. For safer information about HDG method, you can report to the corresponding page in Feel++ mathematics [FeelppMath], or in the article [HDG2020].

1. HDG formulation

The HDG method uses discontinuous elements (contrary to usual, the elements are continuous in standard Galerkin methods).

It consists in rewrite a system of equations in a first-order system.

In the HDG method, we have an unknown, the potential field (for example the temperature) on which we add the associated flux, who becomes unknown too. It allows formulating the problem with a better cost, by introducing the trace.

The HDG method has some attractive features :

It gives an optimal approximation for the potential fiels and the flux
It requires less degrees of freedom than usual Discontinuous Galerkin methods for a comparable accuracy
It allow local post-processing to obtain better approximations and conservation properties.

See [HDG2020] for more complete details about this method.

2. Application to Bestest heat equation

Those are the equtation, tken from here :

\[\begin{align} \rho_m c_m \frac{\partial T_w}{\partial t}+\nabla \cdot\left(-k \nabla T_w\right) &= 0 \tag{$H_w$}\label{hw}\\ \rho_a c_a \frac{\partial T_r}{\partial t} - h_i \cdot A_w \cdot \left( T_w - T_r \right) - \eta \cdot \rho_a c_a \cdot V_a \cdot \left( T_e - T_r \right) &= H \tag{$H_r$}\label{hr} \end{align}\]

Where the unknowns of the equation are :

$T_w$ the temperature in the wall
$T_r$ the temperature in the room (supposed to be constant in space over the room)

We introduce the thermal flux :

\[\q_w=-k\nabla T_w\]

Replacing this in the equation $\ref{hw}$, we get the mixte problem, with first-order equations :

\[\begin{cases} \frac{1}{k}\q_w+\nabla T_w = 0\quad\text{on}\, K\\ \rho_mc_m\frac{\partial T_w}{\partial t} + \nabla \cdot \q_w = 0 \quad\text{on}\,K\\ T_w = \widehat{T}_w\quad\text{on}\,\partial K \end{cases}\]

This system stand on each element $K$.

The third equation introduce the trace $\widehat{T}_w$ of $T_w$.

Let’s partition the border $\Gamma$ of $\Omega$ in three disjoint subsets : $\Gamma = \Gamma_D \cup \Gamma_N \cup \Gamma_{ibc}$

The quantities $\q_w$ and $T_w$ are subject to those boundary conditions :

$T_w=0$ on $\Gamma_D$ (Dirichlet condition)
$\q_w\cdot\n = 0$ on $\Gamma_N$ (Neumann condition)
$T_w$ is a constant (unknown) on $\Gamma_{ibc}$ and $\displaystyle\int_{\Gamma_{ibc}}\q_w\cdot\n=I_\text{target}$ (Integral boundary condition, or IBC)

With such a condition, $T_w$ is constant but unknown. It makes sense because the temperature in the wall is supposed to be constant in space, but can vary in time. It it what tells the EDO $\ref{hr}$

The value of $\int_{\Gamma_{ibc}} \q_w \cdot \n$ is equal to the quantity of heat exchanged with the air, so :

\[\int_{\Gamma_{ibc}} \q_w\cdot \n = \int_{\Gamma_{ibc}} h_i(T_w-T_r)\]

We can deduce the value of $I_\text{target}$ which is therefore equal to $\displaystyle\int_{\Gamma_{ibc}} h(T_w-T_r)$, because $T_w$ is supposed to be constant over $\Gamma_{ibc}$.

The main idea of the method is to solve a small problem on $\widehat{T}_w$, then with a post processing we will deduce the value of $(\q_w,T_w)$ over each element.

Let $\bar{\varphi}\in H^{1/2}(\Gamma)$ such that $\bar{\varphi}|_{\Gamma_D}=0$ and $\bar{\varphi}|_{\Gamma_{ibc}}=1$. The variational problem of our equation system is :

Find $\q_w\in H(div,\Omega)$, $T_w\in L^2(\Omega)$ and $\widehat{T_w} \in \mathrm{span} \left<\bar{\varphi}\right> \oplus H^{1/2}_{00}(\Gamma_N)$ such as for all $\v\in H(div;\Omega)$, $w\in L^2(\Omega)$, $\mu\in\mathrm{span}\left<\bar{\varphi}\right>\oplus H^{1/2}_{00}(\Gamma_N)$ we have :

\[\begin{align} \int_\Omega k^{-1}\q_w \cdot \v - \int_\Omega p\nabla \cdot \v + \int_\Gamma \widehat{T_w} \v\cdot\n = 0\\ \int_\Omega\frac{\rho_m c_m}{\Delta t} T_w w + \int_\Omega \nabla\cdot\q_w w = \int_\Omega\frac{\rho_m c_m}{\Delta t} T_w^\text{prev} w\\ \int_\Gamma \q_w\cdot\n = \int_{\Gamma_N}g_N\mu + I_\text{target}|\Gamma_{ibc}|^{-1}\int_{\Gamma_{ibc}}\mu \tag{$E_\text{var}$}\label{eq:pb_var} \end{align}\]

We use those spaces :

$H^{1/2}(\Gamma) = \left\{\varphi\in L^2(\Gamma)\,\middle|\,\exists u\in H^1(\Omega), u|_\Gamma = \varphi\right\}$
$H^{1/2}_{00}(\Gamma) = \left\{\varphi\in H^{1/2}(\Gamma) \,\middle|\, \varphi=0 \,\text{on}\,\Gamma_D\cup_{\Gamma_{ibc}}\right\}$
$H(div;\Omega) = \left\{u\in L^2(\Omega)\,\middle|\, div(u)\in L^2(\Omega)\right\}$

We make the usual first order approximation $\dfrac{\partial T_w}{\partial t} = \dfrac{T_w-T_w^\text{prev}}{\Delta t}$, where $T_w^\text{prev}$ is the temperature from the previous step of the simulation.

We can show (cf. [HDG2020]) that if $g_N\in L^2(\Gamma_N)$, this problem admits a unique solution $(\q_w,T_w,\widehat{T}_w)$, satisfying the conditions given earlier.

The trace $\widehat{T}_w$ of $T_w$ allows reducing the cost of computation. We can show (cf. [Sala2019]) that the HDG method gives the best approximation error on potential and the flux.

3. Discretization

Let $\T_h$ be a triangulation of $\Omega$. We denote by $\F_h$ the set of all faces of $\T_h$, $\F_h^\Gamma$ the set of faces belonging to $\Gamma$ (and a similar notation for the subsets of $\Gamma$), and $\F_h^0$ the set of the faces belonging to the interior of $\Omega$.

If $\q\in H(div,K)$, then its normal trace on the boundary of $K$, $\q\cdot\n$ belongs to $H^{-1/2}(K)$. In the following, we note $\q^K := \q|_K$. Let’s suppose that $\q\in \left(L^2(\Omega)\right)^d$, and $\q^K\in H(div,K)$ for all $K\in\T_h$. Let $F := \partial K_1\cap\partial K_2$ the border between two elements $K_1,K_2\in \T_h$, this face belonging to $\F_h^0$. We define the jump of the normal trace of $\q$ across $F$ as :

\[[\![\q]\!] := \q^{K_1}\cdot\n_{\partial K_1}|_F + \q^{K_2}\cdot\n_{\partial K_2}|_F\]

We can prove (cf [HDG2020]) that

\[[\![\q]\!] = 0 \forall F\in\F_h^0 \Longleftrightarrow \q\in H(div,\Omega)\]

We introduce those three spaces :

$V_h = \displaystyle\prod_{K\in \T_h}\mathbf{P}_k(K)$
$W_h = \displaystyle\prod_{K\in \T_h} P_k(K)$
$M_h = \mathrm{span}\left<\varphi^*\right> \oplus \displaystyle\prod_{F\in\F_h^0 \cup\F_h^{\Gamma_N}} P_k(F)$

With $\varphi^*\in L^2(\F_h)$ such as $\forall F\in \F_h^{\Gamma_{ibc}}, \varphi^*|_{F}=1$ and $\forall F\in \F\setminus\F_h^{\Gamma_{ibc}}, \varphi^*|_{F}=0$. And $\mathbf{P}_k(K)=(\P_k(K))^d$.

We see that the spaces are richer because the functions belonging to them are in general discontinuous. The functions of $M_h$ play the role of connectors between adjacent elements, in order to reduce the time of calculation.

We define the numerical normal flux on $\partial K$ :

\[\widehat{\q}^{\partial K}_{w,h} \cdot \n_{\partial K} = \q^K_{w,h}|_{\partial K} + \tau_{\partial K} \left(T_{w,h}^K|_{\partial K}-\widehat{T}_{w,h}|_{\partial K}\right) \tag{$E_\text{flux}$}\label{eq:flux}\]

Where $\tau_K$ is a non-negative function that can vary on $\partial K$, and $\tau_K > 0$ on at least one face of $\partial K$. This parameter plays on the continuité of the potential, and to get a stabilized digital flow, we may have to scale this parameter or the flux. It is called stabilization parameter.

The discrete variational problem is therefore : find $\q_{w,h}\in V_h$, $T_{w,h}\in W_h$ and $\widehat{T}_{w,h}\in M_h$ such as $\forall \v_h\in V_h, \forall w_h\in W_h, \forall \mu_h\in M_h$

\[\begin{align} \sum_{K\in\T_h} \left[\int_K k^{-1}\q^K_{w,h} \cdot \v_h^K - \int_K T_{w,h}^K\nabla\cdot \v + \int_{\partial K} \widehat{T}_{w,h}\v_h^K\cdot \n_{\partial K}\right] &= 0\\ \sum_{K\in\T_h} \left[\int_K \frac{\rho_m c_m}{\Delta t}T_{w,h}^K-\int_K \q_{w,h}^K\cdot \nabla w_h^K + \int_{\partial K} \widehat{\q}_{w,h}^{\partial K}\cdot\n_{\partial K} w_h^K\right] &= \sum_{K\in\T_h} \int_K \frac{\rho_m c_m}{\Delta t}T_w^{\text{prev},K}w\\ \sum_{K\in\T_h} \int_{\partial K} \widehat{\q}_{w,h}^{\partial K} \cdot \n_{\partial K} \mu_h &= \int_{\Gamma_N} g_N\mu_h + I_\text{target} |\Gamma_{ibc}|^{-1}\int_{\Gamma_{ibc}} \mu_h \tag{$E_\text{disc}$}\label{eq:pb_var_disc} \end{align}\]

4. Static condensation

The discrete equation $\ref{eq:pb_var_disc}$ hlod in the interior of each $K\in \T_h$ and can b solved for each $K$ to eliminate $\q_w^K$ and $T_{w,h}^K$ in favor of the variable $\widehat{T}_{w,h}^K$. This method of elimination is called static condensation.

With the equation $\ref{eq:flux}$, it is possible to express the normal numerical flux as a function of $\widehat{T}_{w,h}^K$.

We are now going to recast the system $\ref{eq:flux},\ref{eq:pb_var_disc}$ to a global linear system where only the trace of the solution on the boundaries of the mesh appears. After solving the global system, the uknowns wan be recovered locally, in parallel.

First of all, we write the space for the numerical trace $\widehat{T}_{w,h}$ as a direct sum :

\[M_h = \widetilde{M}_h \oplus M_h^*\]

with :

$\widetilde{M}_h = \left\{\mu\in L^2(\F_h) : \forall F\in\F_h^0\cup F_h^{\Gamma_N} \mu|_F\in\P_k(F) \text{ and } \forall F\in \F_h\setminus (\F_h^0\cup\F_h^{\Gamma_N})\mu|_F = 0\right\}$. This is the sapce the the standard trace of the elements we find in usual HDG methods.
$M_h^* = \left\{\mu\in L^2(\F_h) : \mu|_{\Gamma_{ubc}}\in\RR \text{ and } \forall F\in \F_h\setminus \F_h^{\Gamma_{ibc}}\mu|_F = 0\right\}$. This space handles the condition $T_w$ is constant over $\Gamma_{ibc}$.

Let $\lambda_{1,h} = \widehat{T}_{w,h}|_{\widetilde{M}_h}$ and $\lambda_{2,h} = \widehat{T}_{w,h}|_{M_h^*}$. Making integrations by parts, the formulation $\ref{eq:flux},\ref{eq:pb_var_disc}$ can be rewritten as :

\[\begin{align} \sum_{K\in\T_h} \left[\int_K k^{-1}\q_{w,h}^K\cdot \v_h^K - \int_K T_{w,h}^K \nabla\cdot\v_h^K + \int_{\partial K} \lambda_{1,h}\v_h^K\cdot\n_{\partial K} + \int_{\partial K} \lambda_{2,h} \v_h^K\cdot \n_{\partial K}\right] &= 0\\ \sum_{K\in\T_h} \left[ \int_K\frac{\rho_mc_m}{\Delta t}T_{w,h}^K \int_K\nabla\cdot \q_{w,h}^K w_h^K + \int_{\partial K}\tau_{\partial K} T_{w,h}^K w_h^K - \int_{\partial K}\tau_{\partial K}\lambda_{1,h}w_h^K - \int_{\partial K}\tau_{\partial K}\lambda_{2,h}w_h^K\right] &= \sum_{K\in\T_h} \int_K \frac{\rho_mc_m}{\Delta t}T_{w}^{\text{prev},K}\\ \sum_{K\in\T_h} \left[\int_{\partial K} \q_{w,h}^K\cdot\n_{\partial K}\mu_{1,h} + \int_{\partial K} \tau_{\partial K} T_{w,h} \mu_{1,h} - \int_{\partial K} \tau_{\partial K}\lambda_{1,h}\mu_{2,h} \right] &= \int_{\Gamma_N}g_N \mu_{1,N}\\ \sum_{K\in\T_h} \left[\int_{\partial K}\q_{w,h}^K\cdot \n_{\partial K} + \int_{\partial K} \tau_{\partial K}T_{w,h} \mu_{2,h} - \int_{\partial K}\tau_{\partial K}\lambda_{2,h}\mu_{2,h}\right] &= I_\text{target} |\Gamma_{ibc}|^{-1} \int_{\Gamma_{ibc}}\mu_{2,h} \tag{$HDG_{ibc}$}\label{eq:hdg_ibc} \end{align}\]

for all $\v_h\in V_h, w_h\in W_h, \mu_{1,h}\in\widetilde{M}_h, \mu_{2,h}\in M_h^*$. We can notice that $M_h^*$ has a dimension 1, so $M_h^*\simeq\RR$. Because of it, the fourth equation of $\ref{eq:hdg_ibc}$ is a single scalar equation, enforcing the integral boundary condition.

5. Algorithm

Resolution of the system with HDG method

Input : $T_w^0, T_r^0$ and all the parameters for the simulation.

$t\gets0$

Note : in the following, $X^\text{prev}$ corresponds to the quantity $X$ at the previous step of the while loop.

while $t < t_\text{max}$ :

Solve the EDO $\ref{hr}$ with the first order approximation : $\quad\quad T_r = \dfrac{H + \eta\rho_ac_aV_aT_e + h_iA_wT_w + \frac{\rho_ac_a}{\Delta t}T_r^\text{prev}}{\frac{\rho_ac_a}{\Delta t} + h_iA_w + \eta\rho_ac_aV_a}$
$I_\text{target} \gets \int_{\Gamma_{ibc}} h_i(T_w^\text{prev} - T_r)$
Solve the EDP with the discretized problem $\ref{eq:hdg_ibc}$ using discotinuous elements. It gives the trace $\widehat{T}_{w,h}$
With the relation between $\widehat{T}_{w}, T_{w}$ and $\q_w$ $\ref{eq:flux}$, we calculate $T_w$ and $\q_w$ over each element
$t \gets t + \Delta t$

6. New section

From [FeelppMath].

We begin from the system given above, on each element :

\[\begin{align} \frac{1}{k}\q_w + \nabla T_w &= 0 & \text{in}\, K\\ \nabla\cdot \q_w + \rho_mc_m\frac{\partial T_w}{\partial t}&=0 & \text{in}\, K\\ T_w &= \widehat{T}_w & \text{on}\,\partial K \end{align}\]

Those equation lead to this variationnal formulation :

\[\begin{align} \int_K\frac{1}{k} \q_{w,h}\cdot \v - \int_K T_{w,h}\nabla\cdot\v + \int_{\partial K}\widehat{T}_{w,h} \v\cdot\n &= 0\\ -\int_K \q_{w,h} \cdot \nabla w + \int_{\partial K}\q_{w,h}\cdot \n w + \int_K \frac{\rho_mc_m}{\Delta t} T_{w,h}w &= \int_K \frac{\rho_mc_m}{\Delta t} T_{w,h}^\text{prev} \end{align}\]

for all $(\v,w)\in V_h\times W_h$.

The trace $\widehat{T}_w$ can be determined from these equations :

\[\begin{align} [\![\widehat{T}_w]\!]_F &= 0 & \text{if}\,F\in \F_h^0\\ [\![\widehat{T}_w]\!]_F &= h_N & \text{if}\,F\in \F_h^{\Gamma_N}\\ T_w|_F &= \widehat{T}_{w} & \text{if}\,F\in \F_h^{\Gamma_D}\\ \end{align}\]

The first of those three equation is a weak formulation of the continuity of the solution, and the two other corrspond to the boundary conditions.

The variationnal formulation of this set of equation is :

\[\begin{align} \int_{F}\mu[\![\q_{w,h}]\!] &= 0 & \text{if}\,F\in\F_h^0 \\ \int_{F}\mu[\![\q_{w,h}]\!] &= \int_F\mu h_N & \text{if}\,F\in\F_h^{\Gamma_N} \\ \int_{F}\mu\widehat{T}_{w,h} &= \int_F\mu h_D & \text{if}\,F\in\F_h^{\Gamma_D} \end{align}\]

Here we impose weak Dirichlet condition : $\widehat{T}_w$ is the projection $L^2$ of $h_D$ on $M_h$.

The whole system can be written as :

\[\begin{align} \int_{\Omega_h}\frac{1}{k}\q_w\cdot\v - \int_{\Omega_h} T_{w,h} \nabla\cdot\v + \int_{\partial\Omega_h} \widehat{T}_{w,h}\v\cdot\n &=& 0 & &\forall \v\in V_h\\ -\int_{\Omega_h} \q_{w,h}\cdot\nabla w + \int_{\partial\Omega_h}\q_{w,h}\cdot\n w + \int_{\Omega_h}\frac{\rho_mc_m}{\Delta t} T_{w,h}w &=& \int_{\Omega_h}\frac{\rho_mc_m}{\Delta t} T_{w,h}^\text{prev}w & & \forall w\in W_h\\ \int_{\partial\Omega_h\setminus\Gamma} \mu \q_{w,h}\cdot\n &=& 0 & & \forall\mu \in M_h^0\\ \int_{\Gamma_N} \mu \q_{w,h}\cdot\n &=& \int_{\Gamma_N}\mu h_N & &\forall \mu \in M_h^{\Gamma_N}\\ \int_{\Gamma_D}\mu T_{w,h} &=& \int_{\Gamma_D}\mu h_D & &\forall \mu \in M_h^{\Gamma_D} \end{align}\]

Using the formula of numerical flux given in $\ref{eq:flux}$, this set of equation gives this problem : we look for $(\q_{w,h}, T_{w,h}, \widehat{T}_{w,h} \in V_h\times W_h\times M_h$ satisfying the equations :

\[\begin{align} \int_{\Omega_h} \frac{1}{k} \q_{w,h}\cdot \v & & - &\int_{\Omega_h} T_{w,h} \nabla\cdot\v & & + &\int_{\partial\Omega_h} \widehat{T}_{w,h} \v\cdot \n & & = & 0\\ \int_{\Omega_h} \nabla\cdot\q_{w,h} w & & + &\int_{\partial\Omega_h} \tau T_{w,h}w + \int_{\Omega_h} \frac{\rho_mc_m}{\Delta t} T_{w,h}w & & - &\int_{\partial\Omega_h} \tau\widehat{T}_{w,h} w & &= & \int_{\Omega_h}\frac{\rho_mc_m}{\Delta t}T_{w,h}^\text{prev}\\ \int_{\partial\Omega_h\setminus\Gamma} \q_{w,h}\cdot \n \mu_1 & & + & \int_{\partial\Omega_h\setminus\Gamma} \tau T_{w,h}\mu_1 & & - &\int_{\partial\Omega_h\setminus\Gamma} & &= & 0\\ \int_{\Gamma_N} \q_{w,h}\cdot \n\mu_2 & & + & \int_{\Gamma_N}\tau T_{w,h}\mu_2 & & - &\int_{\Gamma_N}\tau\widehat{T}_{w,h}\mu_2 & &= & \int_{\Gamma_N} h_N\mu_2\\ & & & & & &\int_{\Gamma_D} \widehat{T}_{w,h}\mu_3 & &= &\int_{\Gamma_D}h_D\mu_3 \end{align}\]

for all $(\v,w,\mu_1,\mu_2,\mu_3) \in V_h\times W_h\times M_h^0\times M_h^{\Gamma_N}\times M_h^{\Gamma_D}$

The matrix representation of this system lead to this formula on each element :

\[\begin{bmatrix}\q_w^K\\T_w^K\end{bmatrix} = -(A^K)^{-1}B^K{T_w}|_{\partial K} + (A^K)^{-1}F^K\]

where :

\[A^K = \begin{bmatrix} \displaystyle\int_K\frac{1}{k}\q_{w,h}\cdot\v & -\displaystyle\int_K T_{w,h}\nabla\cdot\v\\ \displaystyle\int_K \nabla\cdot \q_{w,h} w & \displaystyle\int_{\partial\Omega_h} \tau T_{w,h}w + \int_{\Omega_h} \frac{\rho_mc_m}{\Delta t} T_{w,h}w \end{bmatrix} \quad B^K = \begin{bmatrix} \displaystyle\int_{\partial\Omega_h} \widehat{T}_{w,h} \v\cdot \n\\ -\displaystyle\int_{\partial\Omega_h} \tau\widehat{T}_{w,h} w \end{bmatrix} \quad F^K = \begin{bmatrix} \mathbf{0}\\ \displaystyle\int_{\Omega_h}\frac{\rho_mc_m}{\Delta t}T_{w,h}^\text{prev} \end{bmatrix}\]

With this formulation, we have a direct dependance of $(\q_{w,h},T_{w,h})$ on the trace $\widehat{T}_{w,h}$.

Then we can solve a small system to get $\widehat{T}_{w,h}$, which can be solved with the toolbox of Feel++.

6.1. Summary

From a big problem, we get a smaller problem to get $\widehat{T}_w$, then on each element we solve a small problem to get $(\q_w,T_w)$.

References

[CEN2007] EN 15026, Hygrothermal performance of building components and building elements - Assessment of moisture transfer by numerical simulation, CEN, 2007.
[FeelppMath] Feel++ Mathematics
[FppPicard] Feel, Non linear problems on http://docs.feelpp.org/math/fem/nonlinear/#_picard_strategy[Feel Mathematics]
[HDG2020] A HDG method for elliptic problems with integral boundary condition: Theory and Applications, Silvia Bertoluzza, Giovanna Guidoboni, Romain Hild, Daniele Prada, Christophe Prud’homme, R. Sacco, Lorenzo Sala, Marcela Szopos, In progess, 2020
[Sala2019] Lorenzo Sala. Mathematical modelling and simulation of ocular blood flows and their interactions.Numerical Analysis [math.NA]. Université de Strasbourg, 2019. English. NNT: 2019STRAD021 . tel-02284233v2
[Škerget2014] Škerget, L. Tadeu, A., BEM numerical simulation of coupled heat, air and moisture flow through a porous solid, Engineering Analysis with Boundary Elements, 2014, 40, p154-161