# Non-Incremental Boundary Element Discretization of non-linear heat equation based on the use of the Proper Generalized Decompositions

Gaël Bonithon, Pierre Joyot, Francisco Chinesta, Pierre Villon

► To cite this version:

Gaël Bonithon, Pierre Joyot, Francisco Chinesta, Pierre Villon. Non-Incremental Boundary Element Discretization of non-linear heat equation based on the use of the Proper Generalized Decompositions. 11th International Conference on Advances in Boundary Element Techniques, Jul 2010, Berlin, Germany. hal-01580912

HAL Id: hal-01580912

https://hal.science/hal-01580912

Submitted on 3 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-Incremental Boundary Element Discretization of non-linear heat equation based on the use of the Proper Generalized Decompositions

G. Bonithon¹,⁴, P. Joyot¹, F. Chinesta² and P. Villon³

¹ ESTIA-Recherche, technopole izarbel, 64210 Bidart, France, pjoyot@estia.fr

² EADS Corporate Foundation International Chair, GEM CNRS-ECN, 1 rue de la Noë BP 92101, 44321 Nantes cedex 3, France, Francisco.Chinesta@ec-nantes.fr

³ UTC-Roberval UMR 6253, 60200 Compiègne, France, pierre.villon@utc.fr

⁴ EPSILON Ingénierie, 10 rue Jean Bart, BP 97431, France, 31674 Labège Cedex

**Keywords:** Boundary element method, Separated representations, Proper Generalized Decomposition

**Abstract.** In this work, we propose a new approach for solving the heat equation within the Boundary Elements method framework. This technique lies in the use of a separated representation of the unknown field that allows decoupling the space problem (that results steady state) from the temporal one (one dimensional that only involves the time coordinate).

## Introduction

The Boundary Elements Method (BEM) allows efficient solution of partial differential equations whose kernel functions are known. The heat equation is one of these candidates when the thermal parameters are assumed constant (linear model). When the model involves large physical domains and time simulation intervals the amount of information that must be stored increases significantly. We propose here an alternative strategy able to change the nature of the problem. Thus, the temperature field involved by the so called heat equation is approximated using a separated representation involving products of space and time functions. This kind of approximation is not new, in fact, proper orthogonal decomposition [1] allows such one decomposition, but in this case this decomposition must be applied a posteriori, i.e. on the transient solution of the considered model.

The technique that we propose in this paper allows to transform the transient model in a sequence of space problems (all of them steady state) and time problems (that only involve the time coordinate). This iteration procedure leads to a proper space-time generalized decomposition of the model solution. The efficiency of such one approach was proven in [2, 3, 4]. In our knowledge, this technique has never been coupled with a BEM for solving the resulting steady problem defined in the physical domain.

We start summarizing the main ideas of the Proper Generalized Decomposition and we will focus on the application of such technique in the context of the BEM. Finally, numerical example, with a non linear source term, will be presented and discussed.

## A Proper Generalized Decomposition Boundary Element Method

Let us consider the heat equation

$$ \frac{\partial u}{\partial t} - a\Delta u = f(u) \quad \text{in } \Omega \times (0, T_{\text{max}}] \qquad (1) $$

with homogeneous initial and boundary conditions, where a is the diffusion coefficient, $\Omega \subset \mathbb{R}^d, d \ge 1, T_{\text{max}} > 0$. The aim of the separated representation method is to compute N couples of functions
$\{(X_i, T_i)\}_{i=1,...,N}$ such that $\{X_i\}_{i=1,...,N}$ and $\{T_i\}_{i=1,...,N}$ are defined respectively in $\Omega$ and $(0, T_{max}]$ and the solution $u$ of this problem can be written in the separate form

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{N} T_i(t) \cdot X_i(\mathbf{x}) \qquad (2)$$

The weak formulation yields: Find $u(\mathbf{x}, t)$ such that

$$\int_{0}^{T_{\max}} \int_{\Omega} u^{*} \left( \frac{\partial u}{\partial t} - a \Delta u - f(\mathbf{x}, t) \right) d\mathbf{x}dt = 0 \qquad (3)$$

for all the functions $u^*(\mathbf{x}, t)$ in an appropriate functional space.

We compute now the functions involved in the sum (eq (2)). We suppose that the set of functional couples $\{(X_i, T_i)\}_{i=1,...,n}$ with $0 \le n < N$ are already known (they have been previously computed) and that at the present iteration we search the enrichment couple $(R(t), S(\mathbf{x}))$ by applying an alternating directions fixed point algorithm that after convergence will constitute the next functional couple $(X_{n+1}, T_{n+1})$. Hence, at the present iteration, $n$, we assume the separated representation

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{n} T_i(t) \cdot X_i(\mathbf{x}) + R(t) \cdot S(\mathbf{x}) \qquad (4)$$

The weighting function $u^*$ is then assumed as

$$u^* = S \cdot R^* + R \cdot S^* \qquad (5)$$

Introducing (eq (4)) and (eq (5)) into (eq (3)) it results

$$\int_{0}^{T_{\max}} \int_{\Omega} (S \cdot R^* + R \cdot S^*) \cdot \left( S \cdot \frac{\partial R}{\partial t} - a \Delta S \cdot R \right) d\mathbf{x}dt = \\ = \int_{0}^{T_{\max}} \int_{\Omega} (S \cdot R^* + R \cdot S^*) \cdot \left( f(\mathbf{x},t) - \sum_{i=1}^{n} X_i \cdot \frac{\partial T_i}{\partial t} + a \sum_{i=1}^{n} \Delta X_i \cdot T_i \right) d\mathbf{x}dt \qquad (6)$$

We apply an alternating directions fixed point algorithm to compute the couple of functions $(R, S)$:

* Computing the function $S(\mathbf{x})$.

First, we suppose that $R$ is known, implying that $R^*$ vanishes in (eq (5)). Thus, eq (6) writes

$$\int_{\Omega} S^* \cdot (\alpha_t S - a \beta_t \Delta S) d\mathbf{x} = \int_{\Omega} S^* \cdot \left( \gamma_t(\mathbf{x}) - \sum_{i=1}^{n} \alpha_t^i X_i + a \sum_{i=1}^{n} \beta_t^i \Delta X_i \right) d\mathbf{x} \qquad (7)$$
where

$$
\left\{
\begin{aligned}
\alpha_t &= \int_0^{T_{\max}} R(t) \cdot \frac{\partial R}{\partial t}(t) dt \\
\alpha_t^i &= \int_0^{T_{\max}} R(t) \cdot \frac{\partial T_i}{\partial t}(t) dt \\
\beta_t &= \int_0^{T_{\max}} R^2(t) dt \\
\beta_t^i &= \int_0^{T_{\max}} R(t) \cdot T_i(t) dt \\
\gamma_t(\mathbf{x}) &= \int_0^{T_{\max}} R(t) \cdot f(\mathbf{x}, t) dt; \quad \forall \mathbf{x} \in \Omega
\end{aligned}
\right.
\qquad (8)
$$

The weak formulation (eq (7)) is satisfied for all $S^*$, therefore we could come back to the associated strong formulation

$$ \alpha_t S - a \beta_t \Delta S = \gamma_t - \sum_{i=1}^{n} \alpha_t^i X_i + a \sum_{i=1}^{n} \beta_t^i \Delta X_i \quad (9) $$

that one could solve by using any appropriate discretization technique for computing the space function $S(\mathbf{x})$.

* Computing the function $R(t)$.

From the function $S(\mathbf{x})$ just computed, we search $R(t)$. In this case $S^*$ vanishes in (eq (5)) and eq (6) reduces to

$$
\begin{align}
& \int_0^{T_{\max}} \int_{\Omega} (S \cdot R^*) \cdot \left( S \cdot \frac{\partial R}{\partial t} - a \Delta S \cdot R \right) d\mathbf{x} dt = \nonumber \\
& = \int_0^{T_{\max}} \int_{\Omega} (S \cdot R^*) \cdot \left( f(x,t) - \sum_{i=1}^n X_i \cdot \frac{\partial T_i}{\partial t} + a \sum_{i=1}^n \Delta X_i \cdot T_i \right) d\mathbf{x} dt \tag{10}
\end{align}
$$

where all the spatial functions can be integrated in $\Omega$. Thus, by using the following notations

$$
\left\{
\begin{aligned}
\alpha_x &= \int_{\Omega} S(\mathbf{x}) \cdot \Delta S(\mathbf{x}) \, d\mathbf{x} \\
\alpha_x^i &= \int_{\Omega} S(\mathbf{x}) \cdot \Delta X_i(\mathbf{x}) \, d\mathbf{x} \\
\beta_x &= \int_{\Omega} S^2(\mathbf{x}) \, d\mathbf{x} \\
\beta_x^i &= \int_{\Omega} S(\mathbf{x}) \cdot X_i(\mathbf{x}) \, d\mathbf{x} \\
\gamma_x(t) &= \int_{\Omega} S(\mathbf{x}) \cdot f(\mathbf{x}, t) \, d\mathbf{x}; \quad \forall t
\end{aligned}
\right.
\qquad (11)
$$

equation (eq (10)) reads

$$ \int_{0}^{T_{\max}} R^{*} \cdot \left( \beta_{x} \frac{\partial R}{\partial t} - a \alpha_{x} R - \gamma_{x}(t) + \sum_{i=1}^{n} \beta_{x}^{i} \frac{\partial T_{i}}{\partial t} - \sum_{i=1}^{n} a \alpha_{x}^{i} \cdot T_{i} \right) dt = 0 \quad (12) $$

As eq (12) holds for all $S^*$, we could come back to the strong formulation
$$ \beta_x \frac{\partial R}{\partial t} = a \cdot \alpha_x \cdot R + \gamma_x(t) - \sum_{i=1}^{n} \beta_x^i \cdot \frac{\partial T_i}{\partial t} + \sum_{i=1}^{n} a \cdot \alpha_x^i \cdot T_i \quad (13) $$

which is a first order ordinary differential equation that can be solved easily (even for extremely small times steps) from its initial condition.

These two steps must be repeated until convergence, that is, until verifying that both functions reach a fixed point.

The BEM is used to solve eq (9). We can notice that this equation defines a steady state elliptic equation with constant coefficients.

## Numerical example

We considered a simple rectangular domain $\Omega = (0, 1) \times (0, 1)$ and a time interval $\mathcal{I} = (0, 1]$. The source term is set to $f(u) = u^2(1-u)$, the boundary conditions and the initial condition is set to an exact solution of this problem given by:

$$ u^{\text{ref}}(\mathbf{x}, t) = \frac{e^{\eta(\mathbf{x},t)}}{2 + e^{\eta(\mathbf{x},t)}} $$

$$ \text{with } \eta(\mathbf{x},t) = \frac{1}{\sqrt{2}} \left( x + \frac{t}{\sqrt{2}} \right) $$

The domain boundary $\Gamma$ consists of $n_{\Gamma} \times n_{\Gamma}$ segments $\Gamma_i$. The time interval $\mathcal{I}$ is discretized by using $n_{\tau}$ nodes, uniformly distributed.

First we are analyzing the convergence rates as a function of the space discretization (i.e. $n_{\Gamma}$). For all the space meshes the time discretization (i.e. $n_{\tau}$) is adapted in order to reach the maximum precision.

Figure 1 shows the evolution of the $L^2$ error in time and space as a function of the level of approximation, that is, as a function of the number of functional couples $X_i(\mathbf{x}) \cdot T_i(t)$ involved in the approximation of $u(\mathbf{x},t)$ for different meshes. This error is defined by:

$$ e_n = \frac{\left\| \sum_{i=1}^{n} X_i(\mathbf{x}) \cdot T_i(t) - u^{\text{ref}}(\mathbf{x}, t) \right\|_{L^2_{\Omega \times \mathcal{I}}}}{\| u^{\text{ref}}(\mathbf{x}, t) \|_{L^2_{\Omega \times \mathcal{I}}}} $$

We can notice that for a given number of functional couples the error $e_n$ decreases when $n_{\Gamma}$ increases, reaching an asymptotic value. For reducing the value of the error we must increase $n_{\Gamma}$ as well as the number of functional couples $X_i(\mathbf{x}) \cdot T_i(t)$ involved in the functional approximation. In the case of the example here addressed we must consider 4 functional couples for reaching a quadratic convergence rate for $4 \le n_{\Gamma} \le 16$.

Figure 3 depicts functions $\{X_1(\mathbf{x}), T_1(t)\}$, $\{X_2(\mathbf{x}), T_2(t)\}$, $\{X_3(\mathbf{x}), T_3(t)\}$ for $n_{\Gamma} = 16$ and $n_{\tau} = 256$. Finally, figure 2 depicts the unknown field $u(\mathbf{x},t)$.

## Conclusion

The proposed approach transforms an incremental BEM procedure into a decoupled one that needs the solution of some steady problems defined in space (Poisson equation in the case here addressed),
<figure><img src="image_2.png" /><figcaption>Fig. 1: Evolution of the error $e_n$ versus the space discretization for different levels of approximation $n$.</figcaption></figure>

and some ordinary differential equations that only involve the time coordinate. Significant reduction of CPU time is expected due to the non-incremental nature of the proposed technique, as well as a significant reduction of the amount of information to be stored. As shown by the numerical example, this technique seems specially well adapted for the treatment of non-linear transient BEM models.
<figure><img src="image_3.png" /><figcaption>Fig. 2: $u(\mathbf{x},t)$ and $u^{\text{ref}}(\mathbf{x},t)$ (dashed line) for $t = \{0.24s, 0.5s, 0.76s, 1s\}$, $n_{\Gamma} = 4$ and $n_{\tau} = 256$</figcaption></figure>
<figure><img src="image_4.png" /></figure>

$X_1(\mathbf{x})$

<figure><img src="image_5.png" /></figure>

$T_1(t)$

<figure><img src="image_6.png" /></figure>

$X_2(\mathbf{x})$

<figure><img src="image_7.png" /></figure>

$T_2(t)$

<figure><img src="image_8.png" /></figure>

$X_3(\mathbf{x})$

<figure><img src="image_9.png" /></figure>

$T_3(t)$

<caption>Fig. 3: Functional couples {$X_i(\mathbf{x})$, $T_i(t)$} for $n_{\Gamma} = 16$ and $n_{\tau} = 256$</caption>
References

[1] F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, and F. Boust. *Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization*. *Computer Methods in Applied Mechanics and Engineering*, 197(5):400–413, 2008.

[2] A. Ammar, F. Chinesta, and P. Joyot. *The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions*. *International Journal for Multiscale Computational Engineering*, 6(3):191–213, 2008.

[3] A. Ammar, B. Mokdad, F. Chinesta, and R. Keunings. *A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids*. *Journal of Non-Newtonian Fluid Mechanics*, 139(3):153–176, 2006.

[4] A. Ammar, B. Mokdad, F. Chinesta, and R. Keunings. *A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part ii: Transient simulation using space-time separated representations*. *Journal of Non-Newtonian Fluid Mechanics*, 144(2-3):98–121, 2007.