# Reflective Conditions for Radiative Transfer in Integral Form with H-Matrices

Olivier Pironneau, Pierre-Henri Tournier

► To cite this version:

Olivier Pironneau, Pierre-Henri Tournier. Reflective Conditions for Radiative Transfer in Integral Form with H-Matrices. 2023. hal-04123799

HAL Id: hal-04123799

https://hal.science/hal-04123799

Preprint submitted on 9 Jun 2023

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
# Reflective Conditions for Radiative Transfer in Integral Form with H-Matrices

Olivier Pironneau* and Pierre-Henri Tournier†

June 9, 2023

## Abstract

In a recent article the authors showed that the radiative Transfer equations with multiple frequencies and scattering can be formulated as a non-linear integral system. In the present article, the formulation is extended to handle reflective boundary conditions. The fixed point method to solve the system is shown to be monotone. The discretization is done with a $P^1$ Finite Element Method. The convolution integrals are precomputed at every vertices of the mesh and stored in compressed hierarchical matrices, using Partially Pivoted Adaptive Cross-Approximation. Then the fixed point iterations involve only matrix vector products. The method is $O(N^{3/\sqrt{N}} \ln N)$, with respect to the number of vertices, when everything is smooth. A numerical implementation is proposed and tested on two examples. As there are some analogies with ray tracing the programming is complex.

**Keywords** MSC classification 85A25, 37N30, 31A10, 35Q30, 68P30, 74S05, Radiative Transfer, Reflective boundaries, Integral equation, H-Matrix, Finite Element Methods.

## Introduction

The Radiative Transport Equations (RTE) describe the behavior of electromagnetic radiation in a domain $\Omega$ as it interacts with matter [14]. It is used to model a wide range of physical phenomena, including the propagation of light through plasma, tomography [18], atmospheric media [13], etc.

The RTE is derived from the basic principles of quantum and statistical mechanics; it is a partial differential equation (PDE) that describes the distribution of radiation intensity in space, time and frequencies, coupled with a budget balance equation (BBE) for the electronic temperature. The PDE takes into

<footnote>* LJLL, Sorbonne Université and Académie des Sciences</footnote>

<footnote>† LJLL, Sorbonne Université</footnote>
account both absorption and scattering of radiation by matter, as well as emission of radiation by sources, which, in the present case, will be restricted to the boundaries of the emitting material.

In [6],[7] the authors have shown that the PDE can be converted into an integral equation for the total radiation at each point in the domain and that the coupling with the BBE can be handled by fixed point iterations. The method leads also to a general proof of existence, uniqueness and regularity of the solution. The difference with earlier studies such as [4] is in the coupling with the equation for the temperature, the BBE, or even the PDE for the temperature when diffusion is important.

In [5] the authors have presented an implementation of the method using H-Matrix compression, a crucial ingredient which makes the evaluation of the integrals $O(N^{3/\sqrt{N}} \ln N)$ with respect to the number of vertices $N$ in the 3D mesh which discretizes the domain $\Omega$; $N \ln N$ is the complexity of the H-Matrix approximations but each element of the matrix requires an integral along a line in the domain. Compared with a brute force solution of the equations as in [10], the integral method keeps a manageable computing time for problems with frequency dependent parameters. However, it did not handle reflective boundary conditions [17].

H-Matrix compression [8], [1],[3], is a mathematical technique used to efficiently represent and manipulate large matrices that arise in a variety of applications. The technique uses a hierarchical structured representation of the matrices allowing fast and accurate numerical computations when the integrals have a convolution type integrand which decays with the distance.

H-Matrix compression works by first defining a hierarchical geometric partitioning of the matrix into smaller and smaller submatrices. This so-called hierarchical block tree is then traversed recursively and far-field interaction blocks which verify a geometric admissibility condition[1] are compressed by using a low rank approximation. The resulting H-Matrix allows for efficient matrix-vector multiplications, among other operations of linear algebra. The technique is particularly important and popular for computational electromagnetics in integral form such as boundary element methods.

With the Partially Pivoted Adaptive Cross-Approximation (ACA) [3] only the needed coefficients of the matrices are computed ($r$ rows and $r$ columns, where $r$ is the rank of the approximation). However the theory requires geometrical smoothness [2].

We have extended the implementation done in [5] using FreeFEM [9] and htool<fnref n="1" />; htool is a parallel C++ toolbox implementing H-Matrices, used in particular for the boundary element method in electromagnetism. FreeFEM is a popular open-source software package for solving PDE systems by the finite element method (FEM).

FreeFEM provides a wide range of pre-built FEM, as well as tools for mesh generation. It has a dedicated high level programming language that allows

<footnote n="1">¹https://github.com/htool-ddm/htool</footnote>
users to meet their specific needs. FreeFEM also supports parallel computing with mpi.

One of the main advantages of FreeFEM for the present study is its ability to handle complex geometries and boundary conditions, especially thanks to its powerful automatic interpolation from volume to surface meshes.

Adding reflective conditions (RC) to the FreeFEM code presented in [5] turned out to require solving the following difficulties:

* Integrate the RC into the integral formulation of the problem

* Show that the fixed point iterations are still monotone.

* Find a formulation compatible with the use of H-Matrices

* Implement the method in the FreeFEM language.

This paper presents the solutions found to overcome these four difficulties. It ends with a numerical test proposed in [11].

# 1 The Radiative Transfer Equations

The problem is formulated in a domain $\Omega \subset \mathbb{R}^3$ with boundary $\Gamma$. The unit sphere in $\mathbb{R}^3$ is called $\mathbb{S}^2$. One must find the radiation (called light from now on) intensity $I_\nu(\mathbf{x}, \omega)$ at all points $\mathbf{x} \in \Omega$, for all directions all $\omega \in \mathbb{S}^2$ and all frequencies $\nu \in \mathbb{R}_+$, satisfying:

$$ \omega \cdot \nabla I_{\nu} + \kappa_{\nu} I_{\nu} = \kappa_{\nu}(1 - a_{\nu})B_{\nu}(T) + \kappa_{\nu}a_{\nu}J_{\nu}, \quad J_{\nu} := \frac{1}{4\pi} \int_{\mathbb{S}^2} I_{\nu} d\omega, \quad (1) $$

$$ \int_{0}^{\infty} \kappa_{\nu}(1 - a_{\nu})(J_{\nu} - B_{\nu}(T))d\nu = 0, \quad (2) $$

$$ I_{\nu}(\mathbf{x}, \omega) = R_{\nu}(\mathbf{x}, \omega)I_{\nu}(\mathbf{x}, \omega - 2(\mathbf{n} \cdot \omega)\mathbf{n}) + Q_{\nu}(\mathbf{x}, \omega), \\ \text{on } \Sigma := \{(\mathbf{x}, \omega) \in \Gamma \times \mathbb{S}^2 : \omega \cdot \mathbf{n}(\mathbf{x}) < 0\}, \quad (3) $$

where $B_\nu(T) = \frac{\nu^3}{e^{\frac{\nu}{T}} - 1}$ is the (rescaled) Planck function. In the RC (3), $R_\nu$ is the portion of light which is reflected and $Q_\nu$ is the light source; $\mathbf{n}(\mathbf{x})$ is the outer normal of $\Gamma$ at $\mathbf{x}$. $\kappa_\nu > 0$ and $a_\nu \in [0, 1]$ are the absorption and scattering coefficients; in general they depend on $\nu$ and $\mathbf{x}$.

**Example 1** If an object $\mathcal{O}$ inside a box $\mathcal{B}$ radiates because it is at temperature $T_0$, then, $\Omega = \mathcal{B} \setminus \mathcal{O}$, $Q_\nu = Q^0[\omega \cdot \mathbf{n}] - B_\nu(T_0)$ on $\mathcal{O}$ and zero elsewhere and $\Sigma \subset \partial \mathcal{B} \times \mathbb{S}^2$. As usual $[f]_- = -\min(f, 0)$.

## 1.1 An Integral Formulation

For clarity we drop the subscript $\nu$ on $\kappa$, $a$ and $I$. Assume that $\Omega$ is bounded and convex (see remark 1). Let

$$ S_{\nu}(\mathbf{x}) = \kappa(1 - a)B_{\nu}(T) + \kappa a J_{\nu}, \quad (4) $$
For a given $\mathbf{x}$ and $\omega$, let $\tau_{\mathbf{x},\omega}$ be such that $(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega) := \mathbf{x} - \tau_{\mathbf{x},\omega}\omega, \omega) \in \Sigma$; the method of characteristics tells us that

$$I(\mathbf{x}, \omega) = I(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega), \omega)e^{-\int_{0}^{\tau_{\mathbf{x},\omega}} \kappa(\mathbf{x}-\omega s)ds} + \int_{0}^{\tau_{\mathbf{x},\omega}} e^{-\int_{0}^{s} \kappa(\mathbf{x}-\omega s')ds'} S_{\nu}(\mathbf{x}-\omega s)ds. \quad (5)$$

Notice that $\tau_{\mathbf{x},\omega} = |\mathbf{x}_{\Sigma} - \mathbf{x}|$ (see Figure 1), therefore, let

$$\begin{align}
J_{\nu}(\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} I(\mathbf{x}, \omega)d\omega = S_{\nu}^{E}(\mathbf{x}) + \mathcal{J}[S_{\nu}](\mathbf{x}) \quad \text{with} \nonumber \\
S_{\nu}^{E}(\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} I(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega), \omega)e^{-\int_{0}^{\tau_{\mathbf{x},\omega}} \kappa(\mathbf{x}-\omega s)ds}d\omega, \nonumber \\
\mathcal{J}[S](\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} \int_{0}^{\tau_{\mathbf{x},\omega}} e^{-\int_{0}^{s} \kappa(\mathbf{x}-\omega s')ds'} S(\mathbf{x}-\omega s)dsd\omega \tag{6} \\
&= \frac{1}{4\pi} \int_{\Omega} S(\mathbf{y}) \frac{e^{-\int_{[\mathbf{x},\mathbf{y}]} \kappa}}{|\mathbf{y}-\mathbf{x}|^2} d\mathbf{y}, \nonumber
\end{align}$$

$$\text{where } \omega'(\omega) := \omega - 2(\mathbf{n} \cdot \omega)\mathbf{n} \text{ and } \int_{[\mathbf{x},\mathbf{y}]} f := |\mathbf{y} - \mathbf{x}| \int_{0}^{1} f(s\mathbf{y} + (1-s)\mathbf{x})ds.$$

To justify the last formula we refer to the following lemma with $\Psi(\mathbf{x}, \mathbf{y}) = S(\mathbf{y})e^{-\int_{[\mathbf{x},\mathbf{y}]} \kappa}$. Again, for clarity, we drop the first argument $\mathbf{x}$.

**Lemma 1** Let $\Omega$ be a convex bounded open set of $\mathbb{R}^3$; let $\Gamma$ be its boundary. Let $\Psi: \Omega \to \mathbb{R}$ be continuous. Let $\tau_{\mathbf{x},\omega} \ge 0$ be such that $\mathbf{x} - \tau_{\mathbf{x},\omega}\omega \in \Gamma$, $\mathbf{x} \in \Omega$. Then

$$\int_{\mathbb{S}^2} \int_0^{\tau_{\mathbf{x},\omega}} \Psi(\mathbf{x} - \omega s) ds d\omega = \int_{\Omega} \frac{\Psi(\mathbf{y})}{|\mathbf{y} - \mathbf{x}|^2} d\mathbf{y}.$$

*Proof:* Denote $\tilde{\Psi}$ the extension of $\Psi$ by zero outside $\Omega$. Let $\omega = (\cos \theta \sin \varphi, \sin \theta \sin \varphi, \cos \varphi)^T$, $\theta \in (0, 2\pi)$, $\varphi \in (-\frac{\pi}{2}, \frac{\pi}{2})$. Consider a partition of the semi infinite line starting at $\mathbf{x}$ in direction $-\omega$ into segments of size $\delta s$ and denote $\mathbf{x}_n = \mathbf{x} - n\delta s\omega$. Then

$$\begin{align}
\int_{\mathbb{S}^2} \int_0^{\tau_{\mathbf{x},\omega}} \tilde{\Psi}(\mathbf{x} - \omega s) ds d\omega &= \lim_{\delta s \to 0} \sum_{n>0} \delta s \int_0^{2\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \tilde{\Psi}(\mathbf{x}_n) \cos \varphi d\varphi d\theta \tag{7} \\
&= \lim_{\delta s \to 0} \sum_{n>0} \int_0^{2\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{\tilde{\Psi}(\mathbf{x}_n)}{|\mathbf{x} - \mathbf{x}_n|^2} |\mathbf{x} - \mathbf{x}_n|^2 |\mathbf{x}_{n+1} - \mathbf{x}_n| \cos \varphi d\varphi d\theta. \nonumber
\end{align}$$

We note that $|\mathbf{x} - \mathbf{x}_n|^2|\mathbf{x}_{n+1} - \mathbf{x}_n| \cos \varphi d\theta d\varphi$ is the elementary volume in the sector $d\theta d\varphi$ between the spheres centered at $\mathbf{x}$ and of radii $|\mathbf{x} - \mathbf{x}_n|$ and $|\mathbf{x} - \mathbf{x}_{n+1}|$. Therefore the right-hand side is an integral in $\mathbf{y} \in \mathbb{R}^3$ of $\frac{\tilde{\Psi}(\mathbf{y})}{|\mathbf{x} - \mathbf{y}|^2} |\mathbf{x} - \mathbf{x}_n|^2$. $\square$

**Remark 1** When $\Omega$ is not convex, one may apply the lemma to its convex closure $\bar{\Omega}$ with $\kappa$ extended to $+\infty$ in $\bar{\Omega} \setminus \Omega$.

**Remark 2** When $R_{\nu} \equiv 0$, $S^{E}$ is given by (6) with $Q_{\nu}$ in place of $I$. As (2) defines a map $\mathcal{T}: J \to T$,

$$T(\mathbf{x}) = \mathcal{T}[J_{\nu}](\mathbf{x}), \quad \forall \mathbf{x} \in \Omega,$$
then, (4), (6) is a nonlinear integral formulation for J:

$$J_{\nu}(\mathbf{x}) = S_{\nu}^{E}(\mathbf{x}) + \mathcal{J}[\kappa(1-a)B_{\nu}(\mathcal{T}[J_{\nu}]) + \kappa a J_{\nu}^{k}](\mathbf{x}), \quad \forall \mathbf{x} \in \Omega. \quad (8)$$

The following fixed point method was shown in [6] to be monotone and convergent:

$$J_{\nu}^{k+1}(\mathbf{x}) = S_{\nu}^{E}(\mathbf{x}) + \mathcal{J}[\kappa(1-a)B_{\nu}(\mathcal{T}[J_{\nu}^{k}](\mathbf{x})) + \kappa a J_{\nu}^{k}](\mathbf{x}), \quad k = 0, 1, \dots \quad (9)$$

Let us extend these properties to the RTE with RC. For clarity let $\mathbf{x}_{\Sigma}$ be short for $\mathbf{x}_{\Sigma}(\mathbf{x}, \omega)$ and let

$$\omega'(\omega) := \omega - 2\omega \cdot \mathbf{n}(\mathbf{x}') \mathbf{n}(\mathbf{x}'), \quad \mathbf{x}_{\Sigma}' := \mathbf{x}_{\Sigma}(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega), \omega') \quad \text{with} \quad \omega := \frac{\mathbf{x} - \mathbf{x}'}{|\mathbf{x} - \mathbf{x}'|}.$$

Let us insert (5) and (3) in (6). Then,

$$\begin{align*}
S_{\nu}^{E}(\mathbf{x}) &= S_{\nu,1}^{E} + S_{\nu,2}^{E} + S_{\nu,3}^{E} \quad \text{with} \\
S_{\nu,1}^{E}(\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} Q_{\nu}(\mathbf{x}_{\Sigma}, \omega) e^{-\int_{0}^{\tau_{\mathbf{x},\omega}} \kappa(\mathbf{x}-\omega s)ds} d\omega, \\
S_{\nu,2}^{E}(\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} R_{\nu}(\mathbf{x}_{\Sigma}, \omega) Q_{\nu}(\mathbf{x}_{\Sigma}', \omega') \left[ e^{-\int_{0}^{\tau_{\mathbf{x}_{\Sigma},\omega'}} \kappa(\mathbf{x}_{\Sigma}-\omega's)ds} e^{-\int_{0}^{\tau_{\mathbf{x},\omega}} \kappa(\mathbf{x}-\omega s)ds} \right] d\omega, \\
S_{\nu,3}^{E}(\mathbf{x}) &:= \frac{1}{4\pi} \int_{\mathbb{S}^2} \left[ R_{\nu}(\mathbf{x}_{\Sigma}, \omega) e^{-\int_{0}^{\tau_{\mathbf{x},\omega}} \kappa(\mathbf{x}-\omega s')ds'} \int_{0}^{\tau_{\mathbf{x}_{\Sigma},\omega'}} e^{-\int_{0}^{s} \kappa(\mathbf{x}_{\Sigma}-\omega's')ds'} S_{\nu}(\mathbf{x}_{\Sigma}-\omega's)ds \right] d\omega.
\end{align*}$$

**Hypothesis 1** Let us rule out multiple reflections and focal points:

1. If $R_{\nu}(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega), \omega) > 0$, then $R_{\nu}(\mathbf{x}_{\Sigma}(\mathbf{x}_{\Sigma}(\mathbf{x}, \omega), \omega'), \omega) = 0$.

2. Given $\mathbf{x}$ and $\mathbf{y}$, there is only a finite number $M$ of $\mathbf{x}'_n \in \Gamma$ such that $[\mathbf{x}'_n, \mathbf{y}]$ is the reflected ray of $[\mathbf{x}, \mathbf{x}'_n]$. Note that $\mathbf{x}'_n$ depends on $\mathbf{x}$ and $\mathbf{y}$.

**Proposition 1** Under Hypothesis 1

$$S_{\nu,3}^{E}(\mathbf{x}) := \sum_{n=1}^{M} \frac{1}{4\pi} \int_{\Omega} R_{\nu}(\mathbf{x}'_n, \frac{\mathbf{x}-\mathbf{x}'_n}{|\mathbf{x}-\mathbf{x}'_n|}) \frac{e^{-\int_{[\mathbf{x},\mathbf{x}'_n]\cup[\mathbf{x}_n,\mathbf{y}]} \kappa}}{(|\mathbf{x}-\mathbf{x}'_n| + |\mathbf{x}'_n-\mathbf{y}|)^2} S(\mathbf{y})d\mathbf{y}.$$

*Proof* Let $\mathbf{x}(s) := \mathbf{x}_{\Sigma} - \omega's$. By Lemma 1,

$$\int_{\mathbb{S}^2} \int_0^{\tau_{\mathbf{x}_{\Sigma}, \omega'}} \left[ R_{\nu}(\mathbf{x}_{\Sigma}, \omega) e^{-\int_{[\mathbf{x}, \mathbf{x}'_{\Sigma}] \cup [\mathbf{x}_{\Sigma}, \mathbf{x}(s)]} \kappa} S(\mathbf{x}(s))ds \right] d\omega \\
= \int_{\Omega} R_{\nu}(\mathbf{x}_{\Sigma}, \omega) S(\mathbf{y}) \frac{e^{-\int_{[\mathbf{x}, \mathbf{x}_{\Sigma}] \cup [\mathbf{x}_{\Sigma}, \mathbf{y}]} \kappa}}{(|\mathbf{x}-\mathbf{x}_{\Sigma}| + |\mathbf{x}_{\Sigma}-\mathbf{y}|)^2} d\mathbf{y},$$
provided that $[x_Σ, y]$ is reflected from $[x, x_Σ]$. Now, by hypothesis, if $x$ and $y$ are given in $\Omega$ there are only a finite number of $x_Σ \in \Gamma$ for which $[x_Σ, y]$ is reflected from $[x, x_Σ]$, (see Figure 1). $\square$

**Proposition 2** Let Hypothesis 1 hold. Then the source terms from the boundaries are

$$S_{\nu,1}^{E}(\mathbf{x}) = \frac{1}{4\pi} \int_{\Gamma} Q_{\nu}(\mathbf{y}, \frac{\mathbf{y}-\mathbf{x}}{|\mathbf{y}-\mathbf{x}|}) \frac{[(\mathbf{y}-\mathbf{x}) \cdot \mathbf{n}(\mathbf{y})]_-}{|\mathbf{y}-\mathbf{x}|^3} e^{-\int_{[\mathbf{x},\mathbf{y}]} \kappa} d\Gamma(\mathbf{y}), \quad (10)$$

$$S_{\nu,2}^{E}(\mathbf{x}) = \sum_{n=1}^{M} \frac{1}{4\pi} \int_{\Gamma} R_{\nu}(\mathbf{x}'_{n}, \frac{\mathbf{x}-\mathbf{x}'_{n}}{|\mathbf{x}-\mathbf{x}'_{n}|}) Q_{\nu}(\mathbf{y}, \frac{\mathbf{x}'_{n}-\mathbf{y}}{|\mathbf{x}'_{n}-\mathbf{y}|}) \frac{[(\mathbf{x}'_{n}-\mathbf{y}) \cdot \mathbf{n}(\mathbf{y})]_- e^{-\int_{[\mathbf{x},\mathbf{x}'_{n}] \cup [\mathbf{x}'_{n},\mathbf{y}]} \kappa}}{|\mathbf{x}'_{n}-\mathbf{y}|(|\mathbf{x}-\mathbf{x}'_{n}| + |\mathbf{x}'_{n}-\mathbf{y}|)^2} d\Gamma(\mathbf{y}). \quad (11)$$

Recall that $\mathbf{x}'_n$ depends on $\mathbf{y}$.

*Proof:* Recall that a solid angle integral at $\mathbf{x}$ of a surface $\Sigma$ is

$$\int_{\mathbb{S}2} f(\mathbf{x}, \mathbf{x}') d\omega(\mathbf{x}') = \int_{\Sigma} f(\mathbf{x}, \mathbf{x}') \frac{[(\mathbf{x}-\mathbf{x}') \cdot \mathbf{n}(\mathbf{x}')]_-}{|\mathbf{x}-\mathbf{x}'|} \frac{d\Sigma(\mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|^2}.$$

Therefore, from the definition of $S_{\nu,2}^E$ above we see that (10) holds.

To prove (11) we start from the definition of $S_{\nu,2}^E$ above. For clarity let us assume that $Q_\nu$ and $R_\nu$ do not depend on $\omega$.

Observe that if a ray from $\mathbf{x}$ in the direction $-\omega$ does not hit, after reflection at $\mathbf{x}'$ on some $\Gamma_R$, a boundary $\Gamma_Q$ at $\mathbf{y}$ where $Q_\nu(\mathbf{y})$ is non zero, then $\omega$ does not contribute to $S_{\nu,2}^E$. Thus, we can use the solid angle of $\Gamma_Q$. However the solid angle is not seen from $\mathbf{x}$ but from $\bar{\mathbf{x}}$, the symmetric of $\mathbf{x}$ with respect to the tangent plane to $\Gamma_R$ at $\mathbf{x}'$. As the distance from $\bar{\mathbf{x}}$ to $\mathbf{y}$ is also $|\mathbf{x}-\mathbf{x}'| + |\mathbf{x}'-\mathbf{y}|$, we obtain (11). $\square$

**Corollary 1**

$$J_{\nu}(\mathbf{x}) = \bar{S}_{\nu}^{E}(\mathbf{x}) + \bar{\mathcal{J}}[S_{\nu}](\mathbf{x}), \quad (12)$$

with $\bar{S}_{\nu}^E(\mathbf{x}) := S_{\nu,1}^E(\mathbf{x}) + S_{\nu,2}^E(\mathbf{x})$ given by Proposition 2 and

$$\bar{\mathcal{J}}[S](\mathbf{x}) = \frac{1}{4\pi} \int_{\Omega} \left[ \frac{e^{-\int_{[\mathbf{x},\mathbf{y}]} \kappa}}{|\mathbf{y}-\mathbf{x}|^2} + \sum_{n=1}^{M} \frac{e^{-\int_{[\mathbf{x},\mathbf{x}'_n] \cup [\mathbf{x}'_n,\mathbf{y}]} \kappa}}{(|\mathbf{x}-\mathbf{x}'_n| + |\mathbf{x}'_n-\mathbf{y}|)^2} R_{\nu}(\mathbf{x}'_n, \frac{\mathbf{x}-\mathbf{x}'_n}{|\mathbf{x}-\mathbf{x}'_n|}) \right] S(\mathbf{y}) d\mathbf{y}. \quad (13)$$

## 1.2 Example

Assume that $\Gamma = \Gamma_Q \cup \Gamma_R$ and $Q_\nu(\mathbf{x}, \omega) = [\omega \cdot \mathbf{n}(\mathbf{x})]_- Q^0$ with $Q^0 > 0$ on $\Gamma_Q$ and $0$ on $\Gamma_R$. Assume $R_\nu(\mathbf{x}, \omega) = R^0$ with $R_0 > 0$ on $\Gamma_R$ and $0$ on $\Gamma_Q$. Assume that
there is never more than one reflection point on $\Gamma_R$, i.e. $M = 1$. Then

$$ \bar{S}_{\nu}^{E}(\mathbf{x}) = \frac{Q^0}{4\pi} \int_{\Gamma_Q} \left[ \left( \frac{[(\mathbf{y} - \mathbf{x}) \cdot \mathbf{n}(\mathbf{y})]_-}{|\mathbf{y} - \mathbf{x}|^2} \right)^2 e^{-\int_{[\mathbf{x}, \mathbf{y}]} \kappa} + R^0 \frac{[(\mathbf{x}'_1 - \mathbf{y}) \cdot \mathbf{n}(\mathbf{y})]_-)^2 e^{-\int_{[\mathbf{x}, \mathbf{x}'_1] \cup [\mathbf{x}'_1, \mathbf{y}]} \kappa}}{|\mathbf{x}'_1 - \mathbf{y}|^2 (|\mathbf{x} - \mathbf{x}'_1| + |\mathbf{x}'_1 - \mathbf{y}|)^2} \right] d\Gamma(\mathbf{y}), $$

$$ \bar{\mathcal{J}}[S](\mathbf{x}) = \frac{1}{4\pi} \int_{\Omega} \left[ \frac{e^{-\int_{[\mathbf{x}, \mathbf{y}]} \kappa}}{|\mathbf{y} - \mathbf{x}|^2} + R^0 \frac{e^{-\int_{[\mathbf{x}, \mathbf{x}'_1] \cup [\mathbf{x}'_1, \mathbf{y}]} \kappa}}{(|\mathbf{x} - \mathbf{x}'_1| + |\mathbf{x}'_1 - \mathbf{y}|)^2} \right] S(\mathbf{y})d\mathbf{y}. $$

<figure><img src="image_2.png" /><figcaption>Figure 1: In this configuration the source $\Gamma_Q$ is the upper square. An RC is imposed on the lower plane $\Gamma_R$. $S_{\nu}^E$ has an integral of the solid angle of the upper square seen from $\mathbf{x}$ plus an integral of the solid angle of the upper square seen from $\bar{\mathbf{x}}$, the symmetric of $\mathbf{x}$ with respect to $\Gamma_R$.</figcaption></figure>

## 1.3 Fixed Point Iterations

Consider the fixed point iterations initialized with $T^0$ and $J^0 = 0$.

**Algorithm** For $k = 0, 1, \dots$:

$$ \text{Set } S_{\nu}^{k}(\mathbf{x}) = \kappa(1 - a)B_{\nu}(T^{k}) + \kappa a J_{\nu}^{k}. $$

$$ \text{Set } J_{\nu}^{k+1}(\mathbf{x}) = \bar{S}_{\nu}^{E}(\mathbf{x}) + \bar{\mathcal{J}}[S_{\nu}^{k}](\mathbf{x}). $$

Compute $T^{k+1}$ by solving (using Newton algorithm) for each $\mathbf{x} \in \Omega$ (14)

$$ \int_{0}^{\infty} \kappa_{\nu}(1 - a_{\nu})(J_{\nu}^{k+1} - B_{\nu}(T^{k+1}))d\nu = 0. $$
**Proposition 3** Let {$J_{\nu}^{*}, T^{*}$} be the solution. If $T^{0}(\mathbf{x}) > T^{*}(\mathbf{x})$, $\forall \mathbf{x} \in \Omega$ then the iterations are monotone decreasing: $T^{k}(\mathbf{x}) > T^{k+1} > T^{*}(\mathbf{x})$, $\forall \mathbf{x} \in \Omega$. Conversely if $T^{0}(\mathbf{x}) < T^{*}(\mathbf{x})$, $\forall \mathbf{x} \in \Omega$ then the iterations are monotone increasing: $T^{k}(\mathbf{x}) < T^{k+1} < T^{*}(\mathbf{x})$, $\forall \mathbf{x} \in \Omega$.

*Proof:* Le us prove it for the monotone increasing sequence.

By subtracting the definition $J_{\nu}^{k}$ from that of $J_{\nu}^{k+1}$ and using the linearity of $\bar{\mathcal{J}}$, we obtain

$$J_{\nu}^{k+1}(\mathbf{x}) - J_{\nu}^{k}(\mathbf{x}) = \bar{\mathcal{J}}[S_{\nu}^{k} - S_{\nu}^{k-1}](\mathbf{x}).$$

As $\mathcal{J}$ is a strictly positive operator, if $S_{\nu}^{k} > S_{\nu}^{k-1}$ for all $\mathbf{x}$ then $J_{\nu}^{k+1}(\mathbf{x}) > J_{\nu}^{k}(\mathbf{x})$. The equation for $T^{k+1}$ is also monotone in the sense that

$$J_{\nu}^{k}(\mathbf{x}) > J_{\nu}^{k}(\mathbf{x}) \implies B_{\nu}(T^{k+1}) > B_{\nu}(T^{k}) \implies T^{k+1} > T^{k},$$

because $B_{\nu}$ is increasing in $T$.

Conclusion: if $T^{1} > T^{0}$ and $S^{1} > S^{0}$ then $T^{k+1} > T^{k}$ for all $k$. One sure way to impose it is to choose $T^{0} = 0$ and $J^{0} = 0$.

To prove that $T^{k} < T^{*}$ we observe that

$$J_{\nu}^{k}(\mathbf{x}) - J_{\nu}^{*}(\mathbf{x}) = \bar{\mathcal{J}}[S_{\nu}^{k-1} - S_{\nu}^{*}](\mathbf{x}).$$

$$\text{Hence } S_{\nu}^{k-1} < S_{\nu}^{*} \implies J_{\nu}^{k}(\mathbf{x}) < J_{\nu}^{*}(\mathbf{x}) \implies T^{k} < T^{*}. \quad \square$$

Discretization seems to preserve this property (see Figure 2).

**Remark 3** Henceforth, convergence and uniqueness can probably be proved as in [7], but there are technical difficulties of functional analysis which may not be appropriately discussed here.

## 2 FEM discretization and Compressed H-Matrices

For clarity consider example 1.2. As the values of $Q^{0}$ and $R^{0}$ take different values on $\Gamma_{Q}$ and $\Gamma_{R}$, we write $Q^{0}(\mathbf{x})$ and $R^{0}(\mathbf{x})$.

The domain $\Omega$ is discretized by a tetrahedral mesh; the boundary $\Gamma$ is discretized by a triangular mesh, not necessarily conforming with the volume mesh.

Let {$\mathbf{x}^{j}$}$_{1}^{N}$ be the vertices of the tetraedra of $\Omega$ and {$\tilde{\mathbf{x}}^{l}$}$_{1}^{L}$ the vertices of the triangles of $\Gamma$.

A continuous $P^{1}$ interpolation of $J$ on the tetrahedral mesh is:

$$J(\mathbf{x}) = \sum_{j=1}^{N} J_{j} \hat{w}^{j}(\mathbf{x}) \text{ where } \hat{w}^{j} \text{ is the } P^{1}\text{- Finite Element hat function of vertex } \mathbf{x}^{j}.$$
Then

$$S_{\nu,j} := aJ_{\nu,j} + (1-a)B_{\nu}(T_j), \quad J_{\nu,i} := \bar{S}_{\nu,i}^E + \sum_j G_{\kappa}^{ij} S_{\nu,j} \quad \text{where}$$

$$G_{\kappa}^{ij} = \frac{1}{4\pi} \int_{\Omega} \left[ \kappa \frac{e^{-f_{[\mathbf{x}^i, \mathbf{y}]}^{\kappa}}}{|\mathbf{x}^i - \mathbf{y}|^2} dy + \sum_{n=1}^{M} R^0(\mathbf{x}'_n) \frac{e^{-f_{[\mathbf{x}^i, \mathbf{x}'_n] \cup [\mathbf{x}'_n, \mathbf{y}]}^{\kappa}}{(|\mathbf{x}^i - \mathbf{x}'_n| + |\mathbf{x}'_n - \mathbf{y}|)^2} \right] \hat{w}^j(\mathbf{y}) dy$$

$$\text{and where } \bar{S}_{\nu,i}^E = \frac{1}{4\pi} \int_{\Gamma} Q^0(\mathbf{y}) \left[ \left( \frac{[(\mathbf{x}^i - \mathbf{y}) \cdot \mathbf{n}(\mathbf{y})]_-}{|\mathbf{x}^i - \mathbf{y}|^2} \right)^2 e^{-f_{[\mathbf{x}^i, \mathbf{y}]}^{\kappa}} + \sum_{n=1}^{M} R^0(\mathbf{x}'_n) \frac{[(\mathbf{x}'_n - \mathbf{y}) \cdot \mathbf{n}(\mathbf{y})]_-}{|\mathbf{x}'_n - \mathbf{y}|^2 (|\mathbf{x}^i - \mathbf{x}'_n| + |\mathbf{x}'_n - \mathbf{y}|)^2} \right] d\Gamma(\mathbf{y}).$$

The integrals are approximated with quadrature at points $\{\mathbf{x}_q^j\}_{1}^{M_q}$. The points are inside the elements; consequently $|\mathbf{x}^i - \mathbf{x}_q^j|$ is never zero. A formula of degree 5, with $M_q = 14$, is used when $|\mathbf{x}^i - \mathbf{y}|$ is small and of degree 2, with $M_q = 4$, otherwise; the results do not change when higher degrees are used. Fortunately when $\mathbf{x}^i$ is close to $\Gamma$ an analytical formula can be used [7].

To compute $\mathbf{x}'_n$ such that $[\mathbf{y}, \mathbf{x}'_n]$ is the reflected ray of $[\mathbf{x}'_n, \mathbf{x}^i]$ a loop on all the elements of the reflecting boundaries is necessary. This can be expensive, but in the case of planar reflective boundaries the symmetric point $\bar{\mathbf{x}}^i$ is easy to compute and so is the intersection of $[\bar{\mathbf{x}}^i, \mathbf{y}]$ with the reflective boundary.

Finally, to the vector $\{\bar{S}_{\nu,i}^E\}_{i=1}^{N}$ we associate a matrix $\{\bar{S}_{i,l}^E\}_{i,l=1}^{N,L}$ by replacing $Q^0(\mathbf{y})$ above by $\tilde{w}^l(\mathbf{y})$. Then:

$$Q^0(\mathbf{y}) = \sum_{l=1}^{L} Q_l^0 \tilde{w}^l(\mathbf{y}) \implies \bar{S}_{\nu,i}^E = \sum_{l=1}^{L} \bar{S}_{i,l}^E Q_l^0.$$

## 2.1 Compression

For each $\nu$ we have two large dense matrices, $\{\bar{G}_{i,j}\}_{i,j=1}^{N,N}$ and $\{\bar{S}_{i,l}^E\}_{i,l=1}^{N,L}$.

**Remark 4** Note that for each value of $\nu$ two matrices are needed. However on close inspection it is really two matrices for each value of $\kappa_{\nu}$. Very often, less than ten values are sufficient to represent a general $\kappa_{\nu}$ by a piece-wise constant interpolation on these values.

These matrices can be compressed as $\mathcal{H}$-matrices [2],[15],[16] (and the references therein) so that the matrix-vector product has complexity $O(N \ln N)$.

The method works best when the kernel in the integrals decays with the distance between $\mathbf{x}^i$ and $\mathbf{y}$. In all matrices the kernel decays with the square of the distance. The $\mathcal{H}$-matrix approximation views $\mathbf{G}$ as a hierarchical tree of square blocks. The blocks correspond to interactions between clusters of points near $\mathbf{x}^j$ and near $\mathbf{x}^i$. A far-field interaction block can be approximated by a low-rank matrix because its singular value decomposition (SVD) has fast decaying singular values. We use the Partially Pivoted Adaptive Cross-Approximation
(ACA) [3] to approximate the first terms of the SVD of the blocks, because only $r$ rows and $r$ columns are needed instead of the whole block, where $r$ is the rank of the approximation. The rank is a function of a user defined parameter $\epsilon$ connected to the relative Frobenius norm error. Another criterion must be met: if $R_1$ (resp. $R_2$) is the radius of a cluster of points centered at $\mathbf{x}_1$ (resp. $\mathbf{x}_2$), then one goes down the hierarchical tree until the corresponding block satisfies $\max(R_1, R_2) < \eta|\mathbf{x}_1 - \mathbf{x}_2|$ where $\eta$ is a user defined parameter. If a leaf is reached, the block is not compressed and all the elements are computed.

The precision is not guaranteed if $[(\mathbf{x} - \mathbf{y}) \cdot \mathbf{n}(\mathbf{y})]$- jumps from one triangular face to another is large. A similar singularity caused by normals is analyzed for a double layer potential formulation in [2] (Example 3.38, p.148) and a remedy is proposed. To check whether this remedy is needed here we ran two cases, one without compression and one with 97% compression. No difference was observed.

<figure><img src="image_3.png" /><figcaption>Figure 2: Values of $J$, for the academic test, along the $y$ axis at $x = z = 15$ computed with a RC. Convergence versus iteration number $n$. When the scaled temperature is initialized to $T^0 = 0.001$ at $n = 0$ the convergence is monotonously increasing. When $T^0 = 0.44$ the convergence is monotonously decreasing.</figcaption></figure>

## 3 An Academic Test

In [11] a semi-analytic solution of the RTE is given for a geometry shown on Figure 3. In this test $a = 0$ and $\kappa$ is a function of $\mathbf{x}$ but not of $\nu$. Hence the grey formulation can be used for $\bar{I} = \int_0^\infty I_\nu d\nu$. By averaging (14) in $\nu$ and due
to the Stefan-Boltzmann relation, the following holds:

$$ \int_{0}^{\infty} B_{\nu}(T) d\nu = \sigma T^4 \quad \text{with} \quad \sigma = \frac{\pi^4}{15} \quad \implies \quad \bar{J}^{k+1}(\mathbf{x}) = \bar{S}^E(\mathbf{x}) + \bar{J}[\kappa\sigma(T^k)^4](\mathbf{x}), \quad \kappa\sigma(T^{k+1})^4 = \kappa\bar{J}^{k+1}. \quad (15) $$

## 3.1 The Geometry

The outer container is $D = (0, 60) \times (0, 100) \times (0, 60)$, in cm. A cube $C = [0, 10]^3$, inside $D$, radiates with intensity $Q^0 = 0.1$. A rectangular cylinder prolonging the radiating cube $(0, 10) \times (10, 100) \times (0, 10)$ has a low absorption $\kappa = 10^{-4}$ while the rest has $\kappa = 0.1$. In Kobayashi's test case 1A there is no scattering, and the three planes containing the origin reflects the radiations perfectly ($R_{\nu} = 1$): $(O, x, z)$, $(O, x, y)$, $(O, y, z)$.

Unfortunately the present method cannot handle volumic radiating region. Consequently we have kept the geometry but only the 3 faces of $C$ inside $D$ radiate in all directions $\omega$ with intensity $Q^0[\omega \cdot n]_-$, where $n$ is the normal to the cube's face pointing inside the cube. The domain is $\Omega = D \setminus C$ (see Figure 3). We refer to this case as Test-3.

## 3.2 Results

To assert the precision of the method we consider first only one reflective plane, $\Gamma_R = (0, y, z)$ and a constant $\kappa = 0.1$. We refer to this case as Test-1. Test-2 is Test-1 with $\kappa$ is as in Test-3.

First we verify, on Test-3, that the convergence is monotone increasing if $T^0$ is small and monotone decreasing if $T^0$ is large (Figure 2). Note that the monotone increasing sequence converges faster.

Next, we compare the results, on Test-1, with a computation on a domain $\bar{D} = (-60, 60) \times (0, 100) \times (0, 60)$ which is $D$ plus the symmetric of $D$ with respect to the plane $(0, y, z)$. This is because reflection on a plane is equivalent to extending the domain by symmetry with respect to that plane.

Figure 4 shows, on Test-1, level surfaces of $J$ computed on the symmetrized domain (but restricted to the original domain) and compared with the same level surfaces but computed with the RC. Surfaces with similar colors should be near each other. In fact the difference is not visible except near $z = 0$.

<figure><img src="image_4.png" /><figcaption>Figure 3: A small cube (colored blue and red on the figure) radiates normally to its faces in a medium which has a very small absorption coefficient $\kappa = 10^{-4}$ in the cylinder prolonging the cube and $\kappa = 0.1$ elsewhere.</figcaption></figure>
Figure 5 shows, on Test-2, the level surfaces of $J$ computed with the RC, and the same level surfaces but computed without any RC on the $(O, y, z)$ plane. It is seen that surfaces with similar color are far from each other. By comparing Figure 4 with Figure 5 we see that the RC does almost the same as symmetry and that no RC at all is a non viable approximation for this problem.

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

Figure 4: Level surfaces of $J$ using a log scale computed with $\kappa = 0.1$ and only one reflective plane, $(0, y, z)$ facing us, slightly to the left. Comparison between a computation done with the RC and a computation done on a symmetrized domain, double in size. Surfaces of equal colors are so near each other that it is hard to distinguish them.

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

Figure 5: Same as in Figure 4 but the RC is not used in one computation. Surfaces of equal colors are far from those using the RC, indicating the absolute necessity of a RC. Here $\kappa$ is as in Test-3.

Similarly, Figure 6 shows $x \mapsto J(x, 25, 25)$, computed on the symmetrized domain.
main, or with the radiative condition or without it.

Finally, Figure 7 shows $x \mapsto J(x, 25, 25)$, computed with the RC on 4 meshes. The same first 4 meshes are used in Table 1 where the theoretical complexity $N^{3/\sqrt{N}} \ln N$ is approximately observed. The compressing ratio for the surface and the volume matrices are shown too.

<figure><img src="image_7.png" /><figcaption>Figure 6: Values of $J$ along the $x$ axis at $y = z = 25$ computed by different methods on the coarse mesh for Test-1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Values of $J$ along the $x$ axis at $y = z = 25$ computed on Test-3 with different meshes.</figcaption></figure>

<table><caption>Table 1: CPU time on a MacBookPro M1, and compression level (C.L.) for Test-3.</caption><thead><tr><th>N vertices</th><th>Surface C.L.</th><th>Volume C.L.</th><th>CPU</th><th><sup>10<sup>5</sup></sup>&frasl;<sub>N</sub> CPU<br>N<sup>3/√N</sup> log N</th></tr></thead><tbody><tr><td>2758</td><td>0.43</td><td>0.60</td><td>5.3"</td><td>1.73</td></tr><tr><td>8003</td><td>0.56</td><td>0.77</td><td>16.5"</td><td>1.14</td></tr><tr><td>26189</td><td>0.67</td><td>0.89</td><td>79.4"</td><td>1.00</td></tr><tr><td>84042</td><td>0.75</td><td>0.95</td><td>389"</td><td>0.93</td></tr><tr><td>195974</td><td>0.82</td><td>0.97</td><td>1563"</td><td>1.13</td></tr></tbody></table>

Test 1A of [11], denoted here Test-3, has been computed, i.e. non constant $\kappa$ and 3 reflective planes. The surface levels of $J$ are shown on Figure 10. Convergence versus mesh size on the line $(x, 25, 25)$ is shown on Figure 7.

The comparison with the data in [11] on the line $(5, y, 5)$ is shown on Figure 8. But since the radiative sources are different (volumic in Kobayashi's and surfacic in our case) we have scaled the result with Kobayashi's value at $x = 5, y = 15, z = 5$.

Finally the $L^2$ error is computed by using the finest mesh, $N = 195974$ as a
reference solution. The results are displayed on Figure 9. It shows the $L^2$-error versus $h := \sqrt[3]{N}$, in log-log scales.

<figure><img src="image_9.png" /><figcaption>Figure 8: Values of $J$ versus $y \ge 15$ at $x = z = 5$ and comparison with the values given in [11]. A scaling is applied so that the radiation intensities coincide at $y = 15$ (because [11] is given for volumic source and the present method handles only surface sources).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Log-log plot of $L^2$ error versus the average mesh size $h = \sqrt[3]{N}$ for Test-3. The line $-x + 1.8$ indicates an error $O(h)$. The reference solution is computed on a mesh with $N = 195974$. The plotted points are computed on meshes with $N$ as in Table 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Kobayashi's test: Level surfaces of $J$ using a log scale. The reflective planes are the $(O, x, z)$, $(O, y, z)$, $(O, x, y)$. The origin $O$ is the lower right corner.</figcaption></figure>
# 4 The Chamonix Valley

In [5] the temperature in the Chamonix Valley due to sunlight was studied. With units in 10km, the emitting domain (the ground) is a rectangle [-0.2, 3.32] x [-3.35, 0.163], with the Chamonix city at (1.5, -1.5). The 3D domains is the emitting domain extruded above the ground up to z = 1, i.e. 10km altitude. The Mont Blanc, in the lower left part of the map is 4807m high. The domain is discretized into tetraedra by a 3D automatic mesh generator from a surface mesh.

Naturally the results are affected by the domain truncation because points near the boundaries receive less than half the scattered light. Now RC can be applied to the 4 vertical planes of the truncation. Note that near the corners there is still a light deficiency which could be corrected by allowing multiple reflections (a programming challenge).

## 4.1 Settings

The ground surface radiates, proportionally to the vertical component of the normal $n_z$, the light of a black body at temperature 300°C at all frequencies (but mostly infrared). The intensity was set at $Q^0 = 2.5$ so as to obtain meaningful temperatures, but since the Earth is not in thermal equilibrium with the sunlight it receives, this choice is arbitrary. In any case rescaling is easily done as $J$ is proportional to $Q^0$.

All mountains are covered with snow above 2500m. The snow-covered ground emits only $0.3Q^0$.

The ground surface mesh has 95K vertices and the volume mesh has 786K vertices. In general 12 fixed point iterations are sufficient to find the temperature $T$ from $J_ν$ and these decrease the error by 6 orders of magnitude.

The surface-to-volume matrix compressed to level 0.942. The volume-to-volume matrix compressed to level 0.982.

## 4.2 Test 1: the Grey Case

In this test $κ$ depends on the altitude, but not on $ν$: $κ = \frac{1}{2}(1 - az)$, with $a = \frac{3}{4}$, except in the cloud. The cloud is a layer between altitude $z_m = 0.3$, i.e; 3000m and $z_M = 0.7$, i.e.7000m where $κ$ is multiplied by a Gaussian random number of means 0.2 and variance 0.8. Scattering is only in the cloud with $a = 0.3(z - z_m)_+(z_M - z)_+/(4(z_M - z_m)^2)$.

The program ran on a French national supercomputer in 5'42" using 1920 processors and 12 threads per MPI process. The surface-to-volume matrix was constructed in 38.5" with RC and 13.53 without. The volume-to-volume matrix took 145" with RC and 95.6 without.
<figure><img src="image_12.png" /><figcaption>Figure 11: Ground temperatures (in °C) computed with RC on the 4 vertical boundaries (left) and without them (right). The third figure displays the difference *T* with RC minus *T* without RC.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Ground and vertical temperatures (in °C) computed with RC in the valley of Chamonix. The mesh is shown in blue on the ground and the intersection of the mesh with the vertical plane is shown in white.</figcaption></figure>
Figure 11 shows the computed temperatures (in Celcius) on the ground with and without RC. The difference between the two temperature fields is also displayed: it is noticeably hotter everywhere by a few degrees when computed with RC. Finally, Figure 12 shows the mesh and the temperatures on the ground and on a vertical plane across the Chamonix valley. The parabolic shape of the mountains increases the temperatures above the ground.

## 4.3 Test 2: the General Case

The setting does not change except that now $κ$ varies with altitude as before but with $a = \frac{3}{4}$ and its dependence on $ν$ is read from the Gemini website [12] (see [5] for details). 583 points are needed to discretize $ν$, but only 8 values are retained for a piecewise discontinuous approximation of $κ_ν$ in the exponentials in the matrices. The computing time is roughly 8 times that of the grey case. The temperature versus altitude above the Chamonix city is plotted on Figure 13. The sudden temperature increase just above the ground is persistent with mesh refinement near the ground. The same computation was done in the same domain with the same mesh but on a flat ground $z = 0$ (the domain is parallelipedic). Then there is no sudden increase, implying that the sudden increase is due to the radiation in a U-shaped valley.

<figure><img src="image_14.png" /><figcaption>Figure 13: Values of $T(1.5, -1.5, z)$ versus $z$ for Chamonix and for a flat ground. In both cases, the dotted curves are the results with $κ$ modified by adding $CO_2$.</figcaption></figure>

<caption>Figure 14: $κ_ν$ from the Gemini measurements and $κ_ν$ modified by adding $CO_2$. The solid blue curve is $J(1.5, -1.5, 0.5)$ versus wavelength $c/ν$.</caption>

In a second computation the Gemini values for $κ_ν$ are modified to be 1 in the range $ν = (3/18, 3/14)$ to simulate an increase of $CO_2$ in the atmosphere. On Figure 13, it is seen that, in this configuration, the $CO_2$ increase the temperature near the ground and decrease it at high altitude. The light intensity $J$ at $(1.5, -1.5, 0.5)$ versus wavelength, $c/ν$, ($c \approx 3$ is the scaled speed of light), is plotted on Figure 14. Notice that the computation
captures complex details due to the discontinuities of $\nu \mapsto \kappa_{\nu}$.

## Conclusion

Compressed H-matrices is an ideal tool for RTE in integral form because the complexity of the method is $O(N^{3/\sqrt{N}} \ln N)$ where N is the number of vertices in the 3D mesh and because it can handle frequency dependent absorption and scattering coefficients at the expense of a finite number of compressed matrices and a finite number of matrix-vector products.

In the present study, the integral nonlinear formulation of RTE studied in [5] has been extended to handle reflective boundary conditions. The monotonicity property of the iterative solver is kept. The discretization with the finite element method is the same. However it is much harder to write a general computer code because of the complexity of possible multiple reflections, as in ray tracing. Hence in this article the numerical validation has been done only for a finite number of plane reflective boundaries and with at most one reflection per ray. For the academic test case and for the Chamonix valley it is essential to add reflective conditions for accuracy.

## Acknowledgement

We would like to thank warmly Frédéric Hecht, for his constant good will to adapt FreeFEM to our needs. The FreeFEM++ app can be downloaded from www.freefem.org

Some computations were made on the machine Joliot-Curie of the national computing center TGCC-GENCI under allocation A0120607330.

## References

[1] M. Bebendorf. Approximation of boundary element matrices. 86(4):565–589, 2000.

[2] M. Bebendorf. *Hierarchical Matrices*. Lecture notes in science and engineering. Springer, Heidelberg, 2008.

[3] S. Boerm, L. Grasedyck, and W. Hackbusch. Hybrid cross approximation of integral operators. *Numerische Mathematik*, 10(12):221–249, 2005.

[4] Yuwei Fan, Jing An, and Lexing Ying. Fast algorithms for integral formulations of steady-state radiative transfer equation. *J. Comp. Physics*, 380(1):191–211, 2019.

[5] F. Golse, F. Hecht, O. Pironneau, D. Smetz, and P.-H. Tournier. Radiative transfer for variable 3d atmospheres. *J. Comp. Physics*, 475(111864):1–19, 2023.
[6] F. Golse and O. Pironneau. Radiative transfer in a fluid. *RACSAM*, Springer, Volume dedicated to I. Diaz(doi.org/10.):https://doi.org/10.1007/s13398-022-01362-x, 2022.

[7] F. Golse and O. Pironneau. Stratified radiative transfer in a fluid and numerical applications to earth science. *SIAM Journal on Numerical Analysis*, 60(5):2963-3000, 2022.

[8] W. Hackbusch. A sparse matrix arithmetic based on h-matrices. part i: Introduction to h-matrices. *Computing*, 62(2):89-108, 1999.

[9] F. Hecht. New developments in freefem++. *J. Numer. Math.*, 20:251-265, 2012.

[10] P. Jolivet, M.A. Badri, and Y. Favennec. Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning. *Journal of Computational Physics*, 437:110313, 2021.

[11] K. Kobayashi, N. Sugimura, and Y. Nagaya. 3d radiation transport benchmark problems and results for simple geometries with void region. *Prog. Nucl. Energy*, 39:119-144, 2001.

[12] S.D. Lord. Earth atmosphere transmittance measurements. Technical report, NASA Technical Memorandum 103957, www.gemini.edu/observing/telescopes-and-sites/sites#Transmission, 1992.

[13] Dimitri Mihalas and Barbara Weibel Mihalas. *Foundations of radiation hydrodynamics*. New York, Oxford: Oxford University Press, 1984.

[14] G. Pomraning. *The equations of Radiation Hydrodynamics*. Pergamon Press, NY, 1973.

[15] S. Rjasanow and O. Steinbach. *The Fast Solution of Boundary Integral Equations*. Mathematical and Analytical Techniques with Applications to Engineering. Springer, Heidelberg, 2007.

[16] S. Sauter and C. Schwab. *Boundary Element Methods*, volume 39 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2011.

[17] C. Siewert. On radiative-transfer problems with reflective boundary conditions and internal emission. *Journal of Applied Mathematics and Physics (ZAMP)*, 35:144-155, 1984.

[18] T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. P. Kaipio. Hybrid radiative-transfer-diffusion model for optical tomography. *Applied optics*, 44(6):876-886, 2005.
Declarations

* As to the specific input of each author, P.-H. Tournier provided the *H*-matrices part and the interface between *htool* and *FreeFEM*; the rest of the program and the theory were done together.

* The authors have no relevant financial or non-financial interests to disclose.

* The authors have no conflicts of interest to declare that are relevant to the content of this article.

* All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.

* The authors have no financial or proprietary interests in any material discussed in this article.
Entretien avec Jean-Louis Dufresne le 31/3/2023

* Compute absorption A or transmission T, A, T ∈ (0, 1).

$$1 - A = T = e^{\kappa_l \rho l} = e^{\kappa_m \rho m}$$

where $\kappa_m$ has dimension $m^2/kg$ and $\kappa_l \sim m$ called absorption coefficients.

* With RT only, $T \mapsto Z$ is decreasing. With adiabatic behavior it is a straight line of slope -10K/km (or 6K/km). See Thermal Equilibrium of the Atmosphere with a Given Distribution of Relative Humidity Syukuro Manabe and Richard T. Wetherald Journal of the Atmospheric Sciences, Volume 24: Issue 3, Page(s): 241-259

https://journals.ametsoc.org/view/journals/atsc/24/3/1520-0469_1967_024_0241_teotaw_2_0_co_2.xml?tab_body=pdf

* Vents catabatiques

https://en.wikipedia.org/wiki/Katabatic_wind

* Clouds: Importance of polarization not clear, (see the LOA lab of Lille). Scattering due to microparticles: more important in the direction of the light in front of the particle.

* Integration of function with respect to $\nu$. Use Malkmus Narrow-band model. See

https://doi.org/10.1364/JOSA.57.000323

https://doi.org/10.1016/S0022-4073(97)00214-8