Hybrid Discontinue 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-processin 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} \q_w+k\nabla T_w = 0\\ \rho_mc_m\frac{\partial T_w}{\partial t} + \nabla \cdot \q_w = 0 \end{cases}\]

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}$.

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 \end{align} \tag{$E_\text{var}$}\label{eq:pb_var}\]

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 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 \end{align} \tag{$E_\text{disc}$}\label{eq:pb_var_disc}\]

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} \end{align} \tag{$HDG_{ibc}$}\label{eq:hdg_ibc}\]

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$

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