L Penazzi
email: lea.penazzi@mines-albi.fr

S Blanco

C Caliot

C Coustet

M El Hafi

R Fournier

J Gautrais

A Golijanek-Je ¸drzejczyk

M Sans

Path integrals formulations leading to propagator evaluation for coupled linear physics in large geometric models

Keywords: Propagator, Monte Carlo method, Large geometric models, Coupled heat transfer

Reformulating linear physics using second kind Fredholm equations is very standard practice. One of the straightforward consequences is that the resulting integrals can be expanded (when the Neumann expansion converges) and probabilized, leading to path statistics and Monte Carlo estimations. An essential feature of these algorithms is that they also allow to estimate propagators for all types of sources, including initial conditions. The resulting practice is a single Monte Carlo run, for one given set of sources, producing propagators that can later be used with any other set of sources for fast simulations, typically as parts of optimization, inversion, sensitivity analysis and command control algorithms. The present paper illustrates how this practice can be extended to problems involving several interacting physics, provided that their coupling is only at the boundary of the system or at interfaces between sub-parts, and may itself be given the form of a second kind Fredholm equation. A full practical implementation is described as part of the Stardis code, with the example of transfering heat via the coupling of radiation, reaction-diffusion and convection as typically expected in the multidisciplinary context of urban climate modeling. Besides, we show how recent advances in computer graphics indicate that these algorithms can be made numerically extremely efficient when facing large CAD geometries: computing the propagator becomes strictly independent of the geometry refinement, i.e. is identical whatever the number of triangles and tetraedra used to numerize the surface and volume descriptions. To the best of our knowledge this is the first report of propagator computations that remains practical for coupled physics in large CAD geometries.

Introduction

First simple illustration. Starting from G. Green's theory, the propagator concept was introduced by R. Feynman as a way to picture, in integral terms, the solution O(� x, t) of a field physics problem at a location � x and time t, when this physics is linear : O(� x, t) is viewed as an integral over all sources S (� x S , t S ) at all locations � x inside the domain D and all times preceding t (down to initial time t I ), multiplied by a scalar ζ(� x, t, � x S , t S ) 1 :

O(� x, t) = � D d� x S � t t I dt S ζ(� x, t, � x S , t S ) S (� x S , t S ) (1)
The propagator ζ(� x, t, � x S , t S ) indicates how each source impacts the solution, and invites for intuitions of the sources being propagated in space and time throughout the system, toward the considered location � x and time t. Historically, this rewriting of G. Green's formalism is mainly significant with regards to the physical pictures it suggests. Here, we will concentrate on translating these pictures in pure computational terms: since ζ(� x, t, � x S , t S ) is independent of the source values, it can be numerically evaluated on its own. Then, any set of sources values can be plugged into Eq. ( 1) to compute the corresponding O(� x, t).

Let us illustrate this concept with a standard practice in radiative transfer, where the factors associated with the propagative point of view are named "shape factors" or "exchange surfaces" (depending on the context and the chosen formulations). Let us take a simple scene with stationary radiative transfer between a camera and two lamps of respective powers P 1 and P 2 . Let O denote the radiative flux incident on a chosen pixel of the camera. Evaluating the shape factor ξ 1 between the first lamp and the target pixel can be obtained by switching on the first lamp alone and solving the radiative transfer equation to get the corresponding pixel flux O 1 (see Fig. 1). The shape factor associated with this first lamp is then ξ 1 = O 1 P 1 . The same radiative transfer equation can also be solved to evaluate the pixel flux O 2 when only the second lamp is on, and the shape factor of the second lamp is ξ 2 = O 2 P 2 (see Fig. 2). Once these factors are known, since they are independent of the powers of the lamps, they can be used to evaluate the pixel flux Õ for any other set of lamp powers P1 and P2 when the two lamps are on simultaneously:

Õ = ξ 1 P1 + ξ 2 P2
(2)

The required number of shape factors can be huge, for example in the domain of infrared radiation where all surface and volume parts of the system can emit radiation, so the question of the efficient numerical evaluation of large numbers of shape factors has raised numerous technical questions, in particular for multiple scattering and multiple reflection configurations with semi-transparent materials. A classical approach is to use reverse Monte Carlo algorithms in which optical paths are tracked from the receptor (the target pixel in our example) backward to sources (the lamps). The procedure can be duplicated for each target of interest (e.g. the other pixels of the camera). Such algorithms, initially designed to estimate O, can also be used to estimate each shape factor (from each lamp to the target pixel), with no additional computational cost: the paths sampled in one single computation can yield estimates of O as well as all the shape factors. In our example, the reverse Monte Carlo computation can then be made only once for a given set of the two lamps power, and does not need to be repeated for each lamp one by one. Once evaluated, the shape factors can be linearly combined with any sets of lamp powers. In this practice, there is one point of concern: we need to consider statistical correlations when evaluating the estimate uncertainty ("error bars") because the estimates are built using the same set of sampled optical paths over different sets of lamps power. Taking this point into account, the approach is then straightforward: only the shape factors required to estimate O are to be computed using a reverse Monte Carlo.

The question in broader terms. Path-integral statistics allow to address the very same question (a single Monte Carlo run computing entirely the propagator) now for advanced linear physics. Without coupling, i.e. when only one single physical phenomenon is at work, the most common approach to propagator computation consists in translating the corresponding partial differential equation (typically a Boltzmann-type transport equation, a reaction-diffusion equation, Maxwell equation, etc) into a second kind Fredholm equation and developping the solution as a Neumann expansion [START_REF] Heinrich | Monte carlo complexity of global solution of integral equations[END_REF][START_REF] Farnoosh | Monte carlo method for solving fredholm integral equations of the second kind[END_REF]. The expansion is then interpreted in statistical terms to define path-integral statistics, i.e. writting the solution as the expectation of a random variable associated to a stochastic path. The corresponding path-integrals address only the solution for a fixed set of sources, but structurally the functional dependance to each source is linear and expressing the propagator using the very same path-integral statistics is straightforward. This extends the illustration of the above simple example. Even for advanced physics or more complex geometries, Monte Carlo algorithms may be designed so that they sample paths along which propagation information can be stored. Then, the propagation information can be directly used when addressing a new set of sources, instead of re-running the Monte Carlo.

In most applicative contexts, computing the propagator with such a Monte Carlo approach would be of great significance. However, this is not yet feasible because of coupling between different physics. Thus, our question becomes the following: can the path-integral approach be theoretically extended to coupled physics and does it remain computationnaly convenient? The present paper briefly initiates a possible answer to the theoretical part: Thanks to double randomization, the Fredholm approach can be extended by probabilizing the coupling. Then, the question of practicability is restricted to only one particular type of coupling: physics interacting via interfaces only (and not at all locations within the field). This restriction allows to stick to a very active field of research in both the physics and computer graphics communities: the design of recursive grids for acceleration of path tracing in complex geometries, leading essentially to Monte Carlo algorithms that are strictly insensitive to geometry refinement: same computation times for scenes described with hundreds of triangles, hundred of voxels, or billions of them [START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF]. By mentionning voxels, we want here to point out that the restriction to coupling via interfaces does not exclude advanced field descriptions, and even multiple physics at work at the same location in the field: only the coupling is not within the field. When randomizing the coupling at the interfaces, the features of the stochastic paths inside each subpart of the field are strictly preserved and the available schemes for accelerating their construction remain unchanged. This leads to Monte Carlo codes computing propagators inside systems involving coupled physics and displaying the same property of being insensitive to geometrical refinement.

The implementation reported here is made within the Stardis code. Stardis is a code conceived as an evolving ensemble of dedicated thermal simulation routines for complex 2D and 3D environments. It is based on Monte Carlo algorithms constructed from reformulations of the principal heat transfer phenomena, conduction, convection, and radiation, resulting in cross-recursive algorithms with the simulation of "thermal paths" exploring space and time until a boundary or initial condition is met. In this sense, it is natural to consider heat transfer phenomena cou-pled. The Stardis framework consists of two primary components, the stardis-solver2 , which is the core library that simulates the coupled heat transfer using Monte Carlo algorithms by sampling these thermal paths without requiring any volumetric mesh, and the project employing the stardis-solver. Presently, the solver is used along the Stardis command line tools 3 which are intented to be reference implementation of a complete workflow using stardis-solver. In addition, the stardis-solver is implemented in SYRTHES 4 , the general thermal software developed by Electricté de France R&D. Stardis was recently used to explore ways toward the definition of new climate services for analysts and designers anticipating climate change in urban area. In this context, Stardis strength was its ability to deal simultaneously with all the spatial and temporal scales involved in the modeling of energy exchanges, from the milimeter scale of windows and heat seals to the kilometers extensions of cities, and from the minute scale of wind and solar fluctuations to typically fifty-year lifetimes of the ground installations to be planed. Radiative transfer inside and outside buildings could be modeled together with full three-dimension heat diffusion inside the solid structures, as well as convection inside each fluid cell, providing the first reported attempt to model energy transfers in a complete city, without any compromise on its geometrical description [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF]. The corresponding approach was name "teapot in a city" by reference to the "teapot in a stadium" paradigm in computer graphics where it was indeed shown that ray tracing Monte Carlo algorithms could deal with quasi-infinite scale ratios on an industrial basis when producing cinema and video-game images. Retaining Stardis for illustration of the present theoretical discussion of propagator estimation in complex systems, we want to highlight an essential idea to be added to the perspective statements of [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF]: using Feynman-Kac strategy to extend the computer graphics strategies to other physics than light transport allows one to model the heat transfer in geometrically complex systems, but also the systematic computation of all the propagators associated to each energy source. We will concentrate on standard heat transfer physics, but the statement is valid for all the physical contexts where same path-integral Monte Carlo strategies were already reported, either using similitudes with slightly modified versions of Stardis (e.g. linear-Boltzmann equation coupled with reaction-diffusion-advection in porous structures [START_REF] Ibarrart | Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm : getting insensitive to geometrical refinement[END_REF]), or with independent codes using the same path-tracing libraries (e.g. electromagnetism, photosynthesis and molecular spectroscopy, in their linear parts [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF][START_REF] Supplis | Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production[END_REF]).

The article is structured as follows:

• Sec. 2 provides the theoretical background.

• Sec. 3 describes the physics involved in the standard version of Stardis.

• Sec. 4 describes the reverse Monte Carlo path sampling strategy solving this model.

• Sec. 5 describes how Stardis stores the propagation data.

• Sec. 6 describes stardis-solver, an implementation of reference used by Stardis.

• Sec. 7 depicts simulation examples.

• Sec. 8 gives some hints towards a generalisation to non-uniform and time-dependent sources, before concluding remarks in Sec. 9.

Theory

From second-kind Fredholm equations to path-statistics

For the sake of exposition, consider a single model for the field of a quantity θ ≡ θ(� x, t) at location � x and time t within a domain Ω of boundary ∂Ω, where sources S are known, in the volume only (no source at the boundary) and where initial conditions are reported to -∞. Let us further assume that this model can be formulated as a second kind Fredholm equation:

θ(� x, t) = � t -∞ dt b � ∂Ω d� x b G b (� x b , t b |� x, t) θ(� x b , t b ) + � t -∞ dt � � Ω d� x � � S (� x � , t � ) + G(� x � , t � |� x, t) θ(� x � , t � ) � (3)
where the Green functions G and G b are known. To illustrate how such a Fredholm equation leads to a path-integral statistical description, let us further assume that θ(� x, t) is known at the boundary: θ(� x, t) = θ b (� x, t) ∀� x ∈ ∂Ω where θ b is fixed (this condition will translate into a coupling issue in section 2.2). Then, assuming convergence as in a standard Neumann expansion, replacing θ(� x � , t � ) by a recursive call to Eq. 3 yields:

θ(� x, t) = � t -∞ dt b,1 � ∂Ω d� x b,1 G b,1 (� x b,1 , t b,1 |� x, t) θ b (� x b,1 , t b,1 ) + � t -∞ dt 1 � Ω d� x 1 � S (� x 1 , t 1 ) + G(� x 1 , t 1 |� x, t) �� t 1 -∞ dt b,2 � ∂Ω d� x b,2 G b,2 (� x b,2 , t b,2 |� x 1 , t 1 ) θ b (� x b,2 , t b,2 ) + � t 1 -∞ dt 2 � Ω d� x 2 � S (� x 2 , t 2 ) + G(� x 2 , t 2 |� x 1 , t 1 ) ... � �� (4)
This expression can be, in turn, probabilized to define a random path backward in time, starting at (� x, � t) and ending at the boundary at a previous time:

θ(� x, t) = P b (� x, t) � t -∞ dt b,1 � ∂Ω d� x b,1 p b (� x b,1 , t b,1 |� x, t) G b,1 (� x b,1 , t b,1 |� x, t) θ b (� x b,1 , t b,1 ) P b (� x, t) p b (� x b,1 , t b,1 |� x, t) + (1 -P b (� x, t)) � t -∞ dt 1 � Ω d� x 1 p(� x 1 , t 1 |� x, t) � S (� x 1 , t 1 ) (1 -P b (� x, t)) p(� x 1 , t 1 |� x, t) + G(� x 1 , t 1 |� x, t) (1 -P b (� x, t)) p(� x 1 , t 1 |� x, t) � P b (� x 1 , t 1 ) � t 1 -∞ dt b,2 � ∂Ω d� x b,2 p b (� x b,2 , t b,2 |� x 1 , t 1 ) G b,2 (� x b,2 , t b,2 |� x 1 , t 1 ) θ b (� x b,2 , t b,2 ) P b (� x 1 , t 1 ) p b (� x b,2 , t b,2 |� x 1 , t 1 ) + (1 -P b (� x 1 , t 1 )) � t 1 -∞ dt 2 � Ω d� x 2 p(� x 2 , t 2 |� x 1 , t 1 ) � S (� x 2 , t 2 ) (1 -P b (� x 1 , t 1 )) p(� x 2 , t 2 |� x 1 , t 1 ) + G(� x 2 , t 2 |� x, t) (1 -P b (� x 1 , t 1 )) p(� x 2 , t 2 |� x 1 , t 1 )
...

���

(5) In this probabilistic reading of the Neumann expansion, the main sum between the integral over the boundary and the integral over the domain is translated into a Bernoulli test between the two integrals. When the boundary integral is retained, a location is sampled at the boundary and the process is stopped because θ has been set to be known at the boundary. When the domain integral is retained, a location is sampled inside the domain, and since θ is unknown, the process is continued from the last sampled location and time, up to a choice of the boundary integral.

In principle, the probabilities P b (to select the boundary at the i-th step) and probability densities p b and p (of the sampled location and time at the i-th step knowing the location and time of the preceding step) are arbitrary, but we can rely on the Monte Carlo literatures dedicated to each physic to indicate the meaningful choices in terms of variance reduction. Overall, in most linear physics where such converging Neumann expansions are available, the theoretical construction and numerical practice of sampling such paths are already available. In box 2.1, we briefly illustrate this starting point with a famous academic example using Feynman-Kac formula. Let us add some details regarding the two main physics addressed in the following sections: radiation physics and reaction-diffusion in confined domains.

Radiative transfer (or linear transport Bolzmann equation in neutronics, biology, etc) is straightforward. The only specificity is that the domain is in phase space and not only in geometrical space. Apart from this, expressing the linear Bolzmann equation in Fredholm terms, and probabilizing it, is very common and leads to elementary pictures such as multiple scattering or multiple reflection path tracing. In terms of of radiative transfer theory, θ is the specific intensity. Its value at a given location in a given direction can be viewed as an average of radiative energy transported along the line of sight. The probabilities P b to reach the boundary are essentially Beer exponential extinction along the line and the probability density of volume collisions, p, is the spatial derivative of this exponential decrease (see Fig. 3). Translating radiative transfer in terms of path-statistics is straightforward and fully rigorous.

Regarding reaction-diffusion equations, the theoretical background is heavier as it involves Brownian motion in confined spaces for which very little can be done with no numerical approximations. When dealing with thermal diffusion with spatially distributed sources, our approach will therefore be approximate. This implies that Eq. 3 will only hold when a spatial discretization is applied after a Brownian motion. This is typically the case for walk-onsphere algorithms that are among the most efficient available path-sampling approaches to confined diffusion [START_REF] Muller | Some Continuous Monte Carlo Methods for the Dirichlet Problem[END_REF][START_REF] Mascagni | �-shell error analysis for "walk on spheres" algorithms[END_REF][START_REF] Hwang | Off-centered "walk-on-spheres" (wos) algorithm[END_REF]. 3: An heterogeneous, multiple scattering and emitting/absorbing medium D is bounded by a partially reflecting wall ∂D. Physical images associated to the evaluation of a local radiative quantity at the location x 0 (for instance, the monospectral radiative intensity) with a reverse Monte Carlo algorithm are as follows. Photons are followed from the location x 0 in the unit direction -� u 0 until the emission/absorption location, either in the volume D or on the wall ∂D (see x b,13 ). At each step j of the optical path within the medium, a free path l j is sampled according to the exponential Beer law, enabling the computation of the next step location x j = x j-1 -l j � u j-1 . If the location x j is in the medium D, the event may be an absorption, a scattering (see x 1 ,x 3 ,x 4 ,x 6 ,x 7 ,x 9 ,x 11 and x 12 ) or a null-collision (see x 2 and x 10 ). A null-collision event corresponds to a pure-forward scattering event in which � u j = � u j-1 (see box 2.1). If the location x j reaches the wall ∂D, the event may be an absorption (see x b,13 ) or a reflexion (see x b,5 and x b,8 ).

At each step, a location is sampled on the smallest sphere tangenting the boundary and this location is projected on the boundary when its distance to the boundary is lower than a numerical parameter �. The probability P b is therefore the fraction of the sphere satisfying this condition (see Fig. 4) and p is a uniform distribution along the rest of the sphere. A slightly modified version of this algorithm will be used for essentially two reasons: 1) including heterogeneous sources and 2) designing an algorithm using line-boundary intersections (and not line-sphere intersections) so that we benefit from the path-tracing acceleration techniques of computer graphics when dealing with complex geometries. The principle is depicted in Fig. 5 and details will be provided in the following sections.

Figure 4: Illustration of the sampling of a random path using the random walk-on-sphere method to estimate density at location x. In order to end the random walk, the boundary of the domain is thickened of a small value � in which the final position x 3 is projected on the boundary.

Feynman-Kac

The most famous example of expressing the solution of a field physics partial differential equation as a second kind Fredholm equation is Feynman-Kac formula. We use it here to illustrate the class of physical problems addressed in the present article. Let us consider the following reaction-diffusion equation:

∂θ ∂t ���� temporal evolution = -� ∇ • � -D � ∇θ � ������������������ diffusion -� ∇ • � � vθ � �������� advection -ν (θ -θ 0 ) ������������ reaction ( 6
)
where D is the diffusion coefficient, � v an advection speed and ν the reaction frequency. When this frequency is function of time and/or space, an overestimate ν can be introduced [START_REF] Galtier | Integral formulation of null-collision Monte Carlo algorithms[END_REF] to write

∂θ ∂t = -� ∇ • � -D � ∇θ � -� ∇ • � � vθ � - ν � θ - � ν ν θ 0 - � 1 - ν ν � θ �� (7)
and applying Feynman-Kac formula [START_REF] Kac | On distributions of certain wiener functionals[END_REF] gives

θ(x, t) = E                 e -ν(t-T ) θ b (X T , T ) ���������������������������������� boundary + � t T dt � νe -ν(t-t � ) � ν (X t � , t � ) ν θ 0 � X t � , t � � + � 1 - ν (X t � , t � ) ν � θ � X t � , t � � � ���������������������������������������������������������������������������������������������������������������������������������������������������������������� volume                 (8)
where X t � is the associated Weiner process backward in time, starting at location x at time t and first encountering the boundary at T . The solution θ(x, t) is herefore expressed as an expectation of a random expression that includes θ (X t � , t � ), i.e. the solution of the very same problem at another location and another time. This implies that the integral formulation of this expectation is a second kind Fredholm equation of the general type that we used as starting point in Eq. 3:

θ(x, t) = � t -∞ dt b � ∂Ω dx b G b (x b , t b |x, t) θ b (x b , t b ) + � t -∞ p T (t b )dt b � t t b dt � � Ω dx � � S (x � , t � ) + G(x � , t � |x, t) θ(x � , t � ) � (9)
with

G b (x b , t b |x, t) = p X T (x b , t b ) e -ν(t-t b ) G(x � , t � |x, t) = νe -ν(t-t � ) p X t � (x � |T = t b ) � 1 - ν(x � , t � ) ν � S (x � , t � ) = e -ν(t-t � ) p X t � (x � |T = t b ) ν(x � , t � ) θ 0 (x � , t � ) (10)
where p X T (x b , t b ) is the probability density of the location x b and time t b at which the Wiener process encounters the boundary and p X t � (x � |T = t b ) is the probability density that X t � = x � knowing T = t b . The only noticeable difference with Eq. 3 is that the t � integral over ] -∞, t] is expressed as a convolution product, requiring the sampling of the boundary encountering time also for the volume part of the Fredholm equation, which makes no difference as far as the following derivations are concerned.

Coupling via interfaces between sub-parts with distinct physics

The two examples above illustrate the broad variety of available path-statistics dealing with linear field physics (see Fig. 3 and Box 2.1). However, they hold when the addressed quantity θ is unique. Here, our objective is to discuss coupling in the particular case where coupling occurs via internal interfaces. Several physics are involved, defining field quantities θ 1 , θ 2 ... θ N only known at the boundary of the overall system. These quantities are not known at the interfaces between two sub-parts of the system. At these interfaces, the coupling constraint is assumed to express each θ i using linear integral operators over the connected fields, i.e. for each � x b ∈ ∂Ω i and each time t b ,

θ i (� x b , t b ) = N � j=1      � t b -∞ dt � b � ∂Ω j d� x b � F bi j (� x b � , t � b |� x b , t b ) θ j (� x b � , t � b ) + � t b -∞ dt � Ω j d� x F i j (� x, t|� x b , t b ) θ j (� x, t)      (11)
where Ω i is the sub-part of Ω where θ i is defined, ∂Ω i its boundary (a part of which may belong to the boundary of the overall system), j is the index referring to a different sub-part of the system over the total N sub-parts, and F i j and F bi j are known (they are null if θ j is not defined in either of the two sub-parts separated by the considered interface). The approach consists in probabilizing this integral constraint exactly the same way the Fredholm equations are probabilized in their respective fields for each single physics (see Eq. 5):

θ i (� x b , t b ) = N � j=1 P i j      P bi j � t b -∞ dt � b � ∂Ω j d� x b � p bi j (� x b � , t � b |� x b , t b ) F bi j (� x b � , t � b |� x b , t b ) θ j (� x b � , t � b ) P bi j p bi j (� x b � , t � b |� x b , t b ) + (1 -P bi j ) � t b -∞ dt � Ω j d� x p i j (� x, t|� x b , t b ) F i j (� x, t|� x b , t b ) θ j (� x, t) (1 -P bi j )p i j (� x, t|� x b , t b )      ( 12
)
where P i j is the probability of selecting the θ j branch when estimating θ i at the interface (see Fig. 6). P i j , P bi j , p i j and p bi j are then the strict equivalent to P, P b , p and p b in Eq. 5.

Figure 6: Ω domain is compound of N sub-domains such as Ω = {Ω 1 , Ω 2 , ..., Ω N } where θ 1 , θ 2 , ..., θ N are only known at the boundary ∂Ω of the overall system. When estimating θ i at a location � x b and a time t b at the interface ∂Ω j between Ω i and Ω j , selecting the θ j branch is given by probability P i j . If this branch is selected, then θ i (� x b , t b ) either takes the value θ j (� x b , t � b ) with probability P bi j (same location, different time), or the value θ j (� x, t) with probability probability 1 -P bi j (different location within domain Ω j , different time). The path sampling goes on until a known temperature at the boundary ∂Ω is reached.

The parallel with the path-statistics described above for single physics is complete and suffices to recursively define statistical paths for the coupled problem:

• When estimating θ i at a given location within Ω i , the path starts as if the physics of θ i was uncoupled and reaches a location at the boundary of Ω i .

• At this stage, for an uncoupled problem, θ i would be known and the path would end. However, the location reached on ∂Ω i may either be on ∂Ω or located at an interface between Ω i and other sub-parts of the system.

• If this location is on ∂Ω, it means that the boundary of the overall system is reached, thus θ i is known and the path ends.

• If this location is at an interface between Ω i and other sub-parts of the system, then a new physics is sampled among those involved at this interface and the question is transformed into the estimation of θ j either inside Ω j or at the boundary ∂Ω j .

• This question is solved by the same algorithm, now for θ j instead of θ i .

• The path is therefore continued and this alternation between the coupled physics occurs at each encountered interface until ∂Ω is reached.

A plain picture of full paths comprised of successions of coupled sub-paths is therefore established. Each sub-path corresponds to one physics (same statistics as those of the corresponding uncoupled sub-path), the full path ends at the overall boundary ∂Ω with a known value θ bk and the index k corresponds to the last visited physics (that of the last sub-path, reaching the boundary). Fig. 7 illustrates the general pattern for a coupled heat transfer situation in which the objective is to evaluate the probe temperature at a specific time and location 5 .

We will prove practical illustrations in the following sections.

Monte Carlo and the storage of propagators

At a first aim, the equations Eqs. ( 5) and ( 12) have been developed to evaluate the quantity θ with a reverse Monte Carlo algorithm. However, the essential point is this: by fully developing a formulation corresponding to Eq. ( 5) up to the coupling described by Eq. ( 12), the propagator ξ is made explicit and can be built up along the successive path samplings. This is exactly the concept described in introduction by the example with the two lamps. A major consequence is that the procedure to evaluate the propagator applies to a large class of physical problems, in complex geometries, with the same strategy as for evaluating a physics as simple as shape factors.

In other words, whatever the number of physics addressed, whatever the number of sub-zones and thus of sampled sub-paths, whatever the shapes of these sub-zones, the evaluation of the propagator ξ is performed recursively along the sampling of these random paths, successively crossing the various physics as diverse as they might be.

As a matter of fact, in Eq. 5, θ is expressed as an expectation along one path within the volume. In Eq. 12, θ i is also expressed as an expectation in which coupling occurs at boundaries. If θ j (x, t) from Eq. 12 is replaced by θ(x, t) of Eq. 5, recursivity appears in a single expectation. Hence, if both equations are coupled, they define a space of infinite dimension, either recursive in volume or at boundaries. Moreover, in this infinite-dimension expectation expression, sources from Green's theory are made explicit and each source contribution can be formally separated to express the propagator itself as an expectation. Hence, as all the sources are visited explicitly during one simulation (i.e. along the path sampling), all propagators (all source contributions) can be stored all along, without any significant additional computation effort. A simulation is performed "as if" θ was the quantity to be estimated for a given set of sources, and, more importantly, when a source is visited, the pre-factor of the corresponding source contribution is stored for later use with different sets of sources. In short, the propagators expectation needs to be estimated only once.

A main strength of the Stardis code used in the following practical illustration is its ability to deal with huge amounts of physical and geometrical data. In [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF], the city is simulated with interacting buildings, each of them described in full details, room per room, and the main message is that the computation time is fully insensitive to the level of refinement of this description 6 . The same observation can be made concerning the computation of propagators: the computation times required for the estimation of the propagators, and for their use in external codes with new sets of sources, are also both insensitive to the refinement level.

Model

System description

At this stage the focus is set on radiation, diffusion and convection coupled inside systems typical of the heattransfer engineering practice. The system is delimited with a system-boundary surface S that is split into N S subsurfaces S i . The internal volume Ω is split into N Ω sub-volumes Ω i of boundaries ∂Ω i (see Fig. 8).

Each sub-volume is either a uniform opaque solid or a perfectly mixed transparent fluid. The contact between adjacent solid sub-volumes is perfect (although Stardis deals with thermal contact resistances, they will not be described here) and the boundary layers at solid-fluid interfaces are not described explicitly: they are summarized by a convective exchange coefficients. The thermal properties of a solid sub-volume Ω i are the thermal conductivity λ i , the mass density ρ i and the mass thermal capacity c i . For a fluid sub-volume, the fluid volume V i is required with the thermal properties ρ i and c i . A power density ψ i can also be prescribed inside each solid sub-volume. There cannot be two fluid sub-volumes adjacent to each other: a fluid sub-volume is always a fluid cell enclosed by solids.

The ensemble of all solid-fluid interfaces (between adjacent sub-volumes of different types) is noted I. It is split into N I sub-interfaces I i . The surface properties are uniform along each sub-interface: the convective exchange coefficient is noted h i ; the surface of the solid is grey of emissivity � i and reflection is modeled using a fraction α i of specular reflection and a fraction 1α i of diffuse reflection.

On each sub-surface S i , the boundary condition can be of the following types:

• type-1 -S i is along a solid sub-volume and the solid temperature is known at this boundary, noted T B,i .

• type-2 -S i is along a solid sub-volume and the boundary flux density is known, noted ϕ B,i .

• type-3 -S i is along a solid sub-volume, a transparent fluid is facing it, and the fluid temperature is known, noted T BF,i . The boundary flux density is then the sum of the convective flux density and the radiative flux density, with uniform values of the convective exchange coefficient h i , the emissivity � i and the specular/diffuse ratio α i . At such a boundary, for incident directions that come from outside the system, the radiance temperature is known, noted θ BR . • type-4 -S i is at the limit of a fluid sub-volume and the limit temperature is known, noted T B,i . This temperature is to be interpreted as the one of a solid surface enclosing the fluid cell, with uniform values of the convective exchange coefficient h i , the emissivity � i and the specular/diffuse ratio α i .

Radiation

The solids are opaque, the fluids are transparent and photon transport is instantaneous: radiative heat transfer can be summarized to instantaneous exchanges between solid surfaces. At a location � y at the surface of a solid sub-volume D i facing a fluid, the radiative flux density ϕ R (� y, t) is the difference between absorption of radiation in all incident directions � ω and emission by the solid due to its local temperature T i (� y, t) :

ϕ R (� y, t) = -� i � σT i (� y, t) 4 - � H i (� y) |� ω.� n i (� y)| I(� ω, � y, t) d� ω � ( 13
)
where I(� ω, � y, t) is the spectrally integrated intensity at � y in direction � ω, σ is the Stefan-Boltzmann constant, � n i (� y) is the unit normal to the solid at � y and H i (� y) is the hemisphere of all incident directions at � y.

It is assumed that radiative transfer can be linearized with respect to temperature around a given reference temperature T ref , which means that

T 4 i ≈ T 4 ref + 4T 3 ref (T i -T ref ) leading to the expression h R = 4� i σT 3 ref .
We then make the choice of translating the spectrally integrated intensity into a radiance temperature θ R = �

D Γ p γ T (� x γ )dγ, i.
e. a mean radiative temperature seen at the solid/fluid interface due to radiative exchanges through the fluid phase.

Observing that

� H i (� y) |� ω.� n i (� y)| π d� ω = 1, Eq. (13) becomes ϕ R (� y, t) = -h R � T i (� y, t) - � H i (� y) |� ω.� n i (� y)| π θ R (� ω, � y, t) d� ω � (14)
We note � z ≡ � z(� y, -� ω) the location of first intersection with a solid sub-volume Ω j of a straight line starting from � y in direction -� ω. If there is no intersection (� z at infinity), then θ R (� ω, � y, t) equals the incident radiance θ BR (� ω, � y, t) known at the system boundary. Otherwise, θ R (� ω, � y, t) = θ R (� ω,� z, t) (pure transport) and θ R (� ω,� z, t) is modeled as the sum of the emission by the solid at temperature T j (� z, t), the specular reflection of incoming radiation in direction -� ω S where -� ω S is the symmetric of � ω around � n j (� z), and the diffuse reflection of radiation incident in all the directions � ω � of the incident hemisphere H j (� z) at � z. Altogether,

                 If � z at ∞ : θ R (� ω, � y, t) = θ BR (� ω, � y, t) If � z ∈ ∂D j : θ R (� ω, � y, t) = � j T j (� z, t) + (1 -� j )α j θ R (-� ω S ,� z, t) + (1 -� j )(1 -α j ) � H j (� z) |� ω � .� n j (� z)| π θ R (� ω � ,� z, t) d� ω � (15)

Diffusion

At any location � x inside a solid sub-volume D i , at any time t, the solid temperature T i ≡ T i (� x, t) is solution of the following heat equation,

ρ i c i ∂T i ∂t = λ i ΔT i + ψ i ( 16
)
where

ψ i ≡ ψ i (� x, t)
is the local value of the power density. The initial condition at time t I is

T i (� x, t I ) = T I,i (� x) (17)
At any location � y at the boundary of D i (i.e. � y ∈ ∂D i ), at any time t, the modeling of the interface or the boundary condition is one of the following :

• If � y is at an interface with another solid sub-volume Ω j ,

λ i � ∇T i .� n i = λ j � ∇T j .� n i (18)
• If � y is at an interface I k with a fluid sub-volume

D j (with h R = 4� k σT 3 ref ), -λ i � ∇T i .� n i = h k (T j -T i ) -h R � T i - � H i (� y) |� ω.� n i (� y)| π θ R (� ω, � y, t) d� ω � (19)
• If � y is at the boundary of the system, in a sub-surface S j with a type-1 boundary condition,

T i = T B, j (20)
• If � y is at the boundary of the system, in a sub-surface S j with a type-2 boundary condition,

-λ i � ∇T i .� n i = ϕ B, j (21)
• If � y is at the boundary of the system, in a sub-surface S j with a type-3 boundary condition (with h R = 4� j σT 3 ref ),

-λ i � ∇T i .� n i = h j (T BF, j -T i ) -h R � T i - � H i (� y) |� ω.� n i (� y)| π θ R (� ω, � y, t) d� ω � (22)

Convection

Inside a fluid sub-volume D i , at any time t, the fluid temperature T i ≡ T i (t) is uniform and its evolution equation is

ρ i c i V i dT i dt = � ∂D i h(� y)(T S (� y) -T i )d� y (23)
where h(� y) and T S (� y) are respectively the convective exchange coefficient and the surface temperature at � y on one of the solid surfaces delimiting the fluid cell. If � y is at an interface I k with a solid sub-volume D j , then h(� y) = h k and T S (� y) = T j (� y, t). If � y is at the boundary of the system, in a sub-surface S j with a type-4 boundary condition, then h(� y) = h j and T S (� y) = T B, j (� y, t).

Path sampling and propagation

In this section, the reverse Monte Carlo algorithm is exposed, where sampling paths is driven by the model above, to estimate a local temperature at location � x and time t, or a radiance temperature at location � x and time t in direction � ω.

When integrated quantities are required (an average temperature on a volume or a surface, a spatially and angularly integrated radiance for simulation of infrared camera pixels, ...), the only algorithmic change is that, prior to initiating a thermal path, � x and/or � ω are sampled accordingly. Further description of such extensions will not be provided here.

A path sampling Monte Carlo algorithm

Depending on the choice of the quantity of interest, the estimate will yield either θ R (� ω, � x, t) for � x inside a fluid subvolume, either T i (� x, t) for � x inside a solid sub-volume, or T i (., t) for any location inside a given fluid sub-volume. In each case, N thermal paths γ j are sampled and each path is used to produce a Monte Carlo weight w γ j . These weights are then averaged to produce an estimate m of the addressed quantity, together with a standard error s associated to this estimate, that can be interpreted in term of a numerical uncertainty.

T i (� x, t) or T i (., t) or θ R (� ω, � x, t) ≈ m = 1 N N � j=1 w γ j (24) s = 1 √ N         1 N N � j=1 w 2 γ j -m 2         1 2 (25)
The focus is here set on the calculation of w γ for any path γ, highlighting its propagative nature and how the information about the sources is collected along the path.

A thermal path is structured as a succession of diffusive, convective and radiative sub-paths. From this point of view, the only difference between the paths used to evaluate θ R (� ω, � x, t) in a fluid, T i (� x, t) in a solid or T i (., t) in a fluid is that they start with a radiative sub-path, a diffusive sub-path or a convective sub-path respectively. Each sub-path can be therefore considered independently, only keeping in mind that: a) at the beginning of the first sub-path the Monte Carlo weight is initialized to w γ = 0, and b) that the end of each sub-path is either the start of a new sub-path, or the end of the whole path γ. Each path γ ends at a location � x γ,end , either inside the system at the initial time t I or at the boundary at a time t γ,end . When it ends with a known incident radiant temperature at the boundary, the corresponding incident direction is � ω γ,end . As announced in the preceding theoretical part (Sec. 2), we start the description "as if" the only objective was the Monte Carlo estimation of temperatures or radiances, and not their associated propagators. This description is exhaustive as far as the path-sampling algorithms are concerned, but we do not recall the theoretical developments justifying them [START_REF] Villefranque | The "teapot in a city": A paradigm shift in urban climate modeling[END_REF][START_REF] Tregan | Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: a general theoretical framework for linear situations[END_REF]. We will focus successively on the sub-paths associated to each of the three physics considered (Secs. 4.1.1, 4.1.2 and 4.1.3), then on the probabilization of the coupling at the interfaces (Sec. 4.2). Once the whole path-sampling Monte Carlo algorithm is set, i.e. when the weight associated to each path is fully defined (Sec. 4.3), we will explain how propagators are constructed for each source by splitting the weight expression into parts and storing the corresponding data (Sec. 5).

Radiative sub-paths

Radiative sub-paths are constructed using a standard backward tracking multiple-reflection algorithm. Starting from � x with the objective of evaluating θ R (� ω, � x, t), a ray is traced in the scene in direction -� ω, looking for a first intersection � z 1 with a solid surface. If no intersection is found (� z 1 is at infinity), then our radiation model says that θ R (� ω, � x, t) = θ BR (� ω, � x, t) where θ BR is a known incident radiance temperature. In this case, the path γ is ended at location � x γ,end = � x, the time t γ,end = t and the direction � ω γ,end = � ω. The Monte Carlo weight is increased by θ BR :

w γ += θ BR (� ω γ,end , � x γ,end , t γ,end ) (26)
Otherwise � z 1 belongs either to a sub-surface S j or a sub-interface S j where the emissivity � j and the specular/diffuse fraction α j are known. A Bernoulli test of probability � j is made to decide whether absorption occurs. If the test is true, the radiative sub-path is ended at � z 1 . If not, reflection occurs with another Bernoulli test to decide between specular or diffuse reflection, with a Lambertian sampling of the reflection direction in the latter case. The path tracing process is then continued from � z 1 in the direction of reflection -� ω 1 , etc, thus defining a succession of possible reflections at locations � z 1 , � z 2 , � z 3 ... until either no reflection is found or absorption occurs. If no reflection is found, the path is ended with the Monte Carlo weight increment of Eq. ( 26). When absorption occurs as a location � z k , then there are two possible cases :

• If � z k belongs to a sub-surface at the system boundary, the temperature T B is known at this surface (a radiative path travels necessarily in a fluid and only a type-4 boundary condition can be encountered), the path γ is ended and the Monte Carlo weight is increased by T B :

w γ += T B (� x γ,end , t γ,end ) (27)
with � x γ,end = � z k and t γ,end = t.

• If � z k belongs to a sub-interface, then the encountered solid is inside the system and its temperature is unknown.

The path γ must be continued with a new sub-path (with no change of w γ ) and a test is made to decide between the three heat transfer modes occurring at this interface: a radiative sub-path back into the fluid, a convective sub-path also in the fluid, or a diffusive path inside the solid. This test is the object of Sec. 4.2.

Radiative sub-paths : Summary

• A radiative sub-path is an instantaneous backward traced ray in a transparent fluid with multiple reflections at solid surfaces.

• If the sub-path encounters a known incident radiance or a known solid temperature, then γ is ended and the Monte Carlo weight is increased by θ BR or T B .

• Otherwise, the Monte Carlo weight is unchanged and at the absorption location, γ is continued with the start of another sub-path at the corresponding solid-fluid interface.

Diffusive sub-paths

Conductive sub-paths are approximate Brownian motions backward in time and space inside a solid. They are constructed as successions of jumps of arbitrary length δ and in isotropically sampled directions. Convergence towards the exact solution is obtained for δ → 0 (Brownian motion is only exact at the limit δ = 0 [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]). However, as the computational time increases considerably when the value of δ decreases, a compromise is required between computational cost and precision. Hence, δ needs to be set sufficiently low to ensure a satisfactory accuracy on the obtained solution and sufficiently high to provide an appropriate computational time.

Starting from � x with the objective of evaluating T (� x, t), the first algorithmic step is the sampling of a backward time shift δt according to an exponential law of parameter

τ i = δ 2 ρ i c i 6λ i , i.e. δt = -τ i ln(r) ( 28
)
where r is sampled uniformly on [0, 1]. If tδt < t I (the backward shift has crossed the initial time), then γ is ended and the Monte Carlo weight is increased by the initial temperature:

w γ += T I (� x γ,end ) (29)
with � x γ,end = � x. Otherwise a direction � u is sampled isotropically in the unit sphere, a jump is made from � x to � x + δ� u and the Monte Carlo weight is increased to account for the local power density ψ:

w γ += β ψ (� x)ψ(� x, t -δt) (30)
with β ψ = δ 2 6λ i . In the vicinity of a solid surface, δ is adjusted depending on � u so that � x + δ� u may either remain inside the solid or reach the solid surface exactly. While � x + δ� u remains in the solid, say D i , then the diffusive sub-path is repeatedly continued from location to location until reaching the surface ∂D i .

When � x + δ� u reaches ∂D i , say at a location � z and time t z , the diffusive sub-path is stopped. If � z belongs to an interface with another solid sub-volume D j , then the temperature of the interface is unknown and γ must be continued with another diffusive sub-path, initiated either inside D i or inside D j . If � z belongs to an interface with a fluid subvolume, then the temperature of the interface is also unknown and γ must also be continued, and the next sub-path can be either a diffusive one, back into D i , or a convective or a radiative one inside or through the fluid. The corresponding tests are described in Sec. 4.2.

If � z is at the boundary of the system, then the algorithm depends on the boundary condition type:

• For a type-1 boundary condition, the boundary temperature T B is known and γ is stopped and the Monte Carlo weight is increased by T B :

w γ += T B (� x γ,end , t γ,end ) (31)
with � x γ,end = � z and t γ,end = t z .

• For a type-2 boundary condition, the location is shifted back into the solid sub-volume, of a distance δ along the normal, the diffusive sub-path is continued from this new location, and the Monte Carlo weight is increased to account for the value of the local flux density:

w γ += β ϕ (� z)ϕ(� z, t z ) (32)
with β ϕ = δ λ i . • For a type-3 boundary condition, neither the boundary temperature nor the flux density is known and γ is continued exactly the same way as for a solid-fluid interface inside the system (see Sec. 4.2). The only difference is that when the following sub-path is a convective one, then the fluid temperature T BF is known and γ is ended. In this case, the Monte Carlo weight is increased by T BF :

w γ += T BF (� x γ,end , t γ,end ) (33)
with � x γ,end = � z and t γ,end = t z (t z is the time at which the � z location was reached).

Diffusive sub-paths : Summary

• A diffusive sub-path is a Brownian motion backward in time inside a solid sub-volume until it reaches either the initial time or the sub-volume boundary.

• If the initial time is reached, then γ is ended and the Monte Carlo weight is increased by T I .

• If the system boundary is reached at location where the temperature is known, then γ is ended and the Monte Carlo weight is increased by T B .

• If the system boundary is reached at location where the flux density is known, then a new diffusive sub-path is initiated, inside the same sub-volume, and the Monte Carlo weight is increased by β ϕ ϕ.

• In all other cases, the diffusive sub-path has reached a location where neither the temperature nor the flux density is known, and a new sub-path (diffusive, convective or radiative) must be initiated from the corresponding interface.

• Along the path, the Monte Carlo weight is increased to account for the local value of the volume source density ψ.

• As Brownian motion is approximated with discrete jumps of length δ, the continuous effect of the source is replaced by a Monte Carlo weight increment of β ψ ψ at each jump.

Convective sub-paths

Convective sub-paths inside a fluid sub-volume D i are independent of their initial location: the fluid cells are perfectly mixed so T i is only a function of time, and the only required information is the time t at which the convective sub-path was initiated. From t, a backward time shift δt is sampled according to an exponential law of parameter

h i V i ρ i c i S i , i.e. δt = h i V i ρ i c i S i ln(r) (34)
where S i and V i are respectively the surface and volume of the fluid cavity, and r is sampled uniformly on [0, 1]. If tδt < t I (the backward shift has crossed the initial time), then γ is ended and the Monte Carlo weight is increased by the initial temperature:

w γ += T I,i (35)
Otherwise a location � z is sampled on ∂D i according to probability density p � Z proportional to the local value of the convective exchange coefficient:

p � Z (� z) = h(� z) � ∂D i h(� z � )d� z � (36)
and the time is shifted to t z = tδt.

If � z is at the system boundary, then this corresponds necessarily to a type-4 boundary condition and the boundary temperature T B is known, so γ is ended and the Monte Carlo weight is increased by T B :

w γ += T B (� x γ,end , t γ,end ) (37)
with � x γ,end = � z and t γ,end = t z . Otherwise � z is at a solid-fluid interface inside the system and the interface temperature is unknown. γ is then continued with a new diffusive, convective or radiative sub-path as described in Sec. 4.2.

Convective sub-paths : Summary

• A convective sub-path inside a fluid sub-volume is only a backward exponential shift in time.

• If the initial time is reached, then γ is ended and the Monte Carlo weight is increased by T I .

• Otherwise a location is sampled on one of the solid surfaces surrounding the fluid.

• If this location is at the boundary of the system and the surface temperature is known, then γ is ended and the Monte Carlo weight is increased by T B .

• Otherwise the convective path has reached an interface where the temperature is unknown and a new sub-path (diffusive, convective or radiative) must be initiated from this interface.

4.2.

Choosing the next sub-path at an interface 4.2.1. Solid-solid interface When describing diffusive sub-paths, we encountered an algorithmic step where a location � z was reached, at the interface between two solid sub-volumes D i and D j , at time t z . At this interface, the temperature was unknown and γ had to be continued. The same question is raised when γ needs to be started at such a location in order to evaluate T (� x, t) with � x = � z and t = t z . This is achieved by first shifting � z along the normal, either by a distance δ i inside D i or by a distance δ j inside D j , where δ i and δ j are the values of the numerical parameter δ used inside D i and inside D j respectively (depending of their characteristic dimensions). Then a diffusive sub-path is started from this shifted location, still at the same date. Choosing the side is made by retaining D i with probability P cond,i and D j with probability P cond, j :

P cond,i = λ i δ i λ i δ i + λ j δ j P cond, j = 1 -P cond,i (38)

Solid-fluid interface

When describing each of the three sub-path types, we encountered the possibility that a location � z is reached, at an interface I k between a solid sub-volume D i and a fluid sub-volume D j , at time t z . At this interface the temperature is unknown and γ is to be continued. The same question is raised when γ needs to be started at such a location in order to evaluate T (� x, t) with � x = � z and t = t z . The same question is also encountered when a diffusive sub-path reaches a location at the system boundary with type-4 boundary condition. We focus the description on the case of an internal interface between D i and D j .

At such an interface, all three heat transfer modes are present: diffusion inside D i , convection and radiation inside Ω j . A test is therefore made to pick among a diffusive, a convective or a radiative sub-path. Conductive, convective and radiative sub-paths are picked with probabilities P cond , P conv and P ray respectively, with

P cond = λ i δ i λ i δ i + h k + h R P conv = h k λ i δ i + h k + h R P ray = 1 -P cond -P conv (39)
If diffusion is retained, then � z is shifted by a distance δ i along the normal inside D i and the diffusive sub-path is initiated at this shifted location. For convection or radiation, the corresponding sub-path is initiated at � z at time t z .

Choosing the next sub-path at a solid-solid or solid-fluid interface : Summary

• Departing from a solid-solid or a solid-fluid interface is made by initiating a sub-path with a heat transfer mode that is sampled according to probabilities reflecting the flux continuity through the interface.

• When a diffusive sub-path is chosen, the starting location is shifted by a distance δ inside the solid.

• In all cases, there is no increment made to the Monte Carlo weight.

The Monte Carlo weight

As mentioned above, each path γ ends at a location � x γ,end , either inside the system at the initial time t I or at the boundary at a time t γ,end . When it ends with a known incident radiant temperature at the boundary, the corresponding incident direction is � ω γ,end . In each case, the Monte Carlo weight is increased by a temperature value that can be T I (� x γ,end ), T B (� x γ,end , t γ,end ), T BF (� x γ,end , t γ,end ) or θ BR (� ω γ,end , � x γ,end , t γ,end ). Let µ γ,end denote the type of ending, from 0 to 3 in the order of this list.

We can then define T γ,end ≡ T γ,end (� x γ,end , t γ,end , � ω γ,end , µ γ,end ) as :

T γ,end (� x γ,end , t γ,end , � ω γ,end , 0) = T I (� x γ,end ) T γ,end (� x γ,end , t γ,end , � ω γ,end , 1) = T B (� x γ,end , t γ,end ) T γ,end (� x γ,end , t γ,end , � ω γ,end , 2) = T BF (� x γ,end , t γ,end ) T γ,end (� x γ,end , t γ,end , � ω γ,end , 3) = θ BR (� ω γ,end , � x γ,end , t γ,end ) (40)
Along γ, diffusive sub-paths may have crossed solid sub-volumes with a volume power source ψ, and we have seen that the Monte Carlo weight was increased by β ψ ψ at each discrete jump location. Let N ψ denote the number of such locations, and � x γ,ψ (k) and t γ,ψ (k) the location and time of the k-th of these Monte Carlo weight increments. Similarly, diffusive sub-paths may have visited boundary locations where the flux density ϕ is known, and we have seen that the Monte Carlo weight was increased by β ϕ ϕ at each such visit. Let N ϕ denote the number of such visits, and � x γ,ϕ (k) and t γ,ϕ (k) the location and time of the k-th of these Monte Carlo weight increments.

With these notations, the complete expression of the Monte Carlo weight associated to γ is

w γ = N ψ � k=1 β ψ (� x γ,ψ (k)) ψ � � x γ,ψ (k), t γ,ψ (k) � + N ϕ � k=1 β ϕ (� x γ,ϕ (k)) ϕ � � x γ,ϕ (k), t γ,ϕ (k) � + T γ,end (� x γ,end , t γ,end , � ω γ,end , µ γ,end ) (41)
The reading of this weight expression illustrates the relation between such path sampling Monte Carlo algorithms and Green's theory. From Green's point of view, the temperature solution of the problem at (� x * , t * ), or the radiance temperature at (� ω * , � x * , t * ) results from an integral of all the sources, T I , T B , T BF , θ BR , ψ and ϕ, at (� x, t, � ω) multiplied by a propagator density that is independent of the source values. From Monte Carlo point of view, it results from the average of a large number of weights that carry the same sources from (� x γ,end , t γ,end , � ω γ,end ) or from (� x γ,ψ (k), t γ,ψ (k)) and (� x γ,ϕ (k), t γ,ϕ (k)), multiplied by factors that are independent of the source values: the factor is 1 for T I , T B , T BF and θ BR , it is β ψ for ψ and β ϕ for ϕ. The paths sample ζ(� x, t, � x S , t S )S (� x S , t S ) in a backward manner. Let us start by considering a simple configuration akin to the radiative transfer example with two lamps of power P 1 and P 2 viewed from a camera pixel.

Here, the heat transfer mode is diffusion inside a cubic solid with two isothermal faces S 1 and S 2 , facing each other, at temperatures T B,1 and T B,2 , the other four faces being adiabatic (see Fig. 9). The addressed quantity is the stationary temperature T (� x) at a location � x inside the solid. In terms of Green's theory, as the problem is stationary, no propagator is required for the initial condition (reported to -∞). There are no volume sources and the only imposed surface flux is null (at the adiabatic faces). The only sources are therefore T B,1 and T B,2 and we note ζ B,1 (� x) and ζ B,2 (� x) the corresponding propagators:

T (� x) = ζ B,1 (� x)T B,1 + ζ B,2 (� x)T B,2 (42)
Considering the expression of the Monte Carlo weight in Eq. ( 41), how can we provide an estimate for ζ B,1 (� x) and one for ζ B,2 (� x) using the same thermal paths as those used to estimate T (� x)?

In the expression of the Monte Carlo weight of the preceding section, for this simple case, the sums over N ψ and N ϕ vanish (no surface flux, no volume flux), � ω γ,end is not used (we estimate a local temperature and not a radiance temperature), t γ,end is unused (the problem is stationary) and µ γ,end = 1 (the path can only end at S 1 or S 2 , i.e. at a boundary with a known solid temperature). Eq. (41) reduces to

w γ = T γ,end (� x γ,end , 1) (43)
with

T γ,end (� x γ,end , 1) = T B,1 if � x γ,end ∈ S 1 and T γ,end (� x γ,end , 1) = T B,2 if � x γ,end ∈ S 2 .
In the spirit of the example used in introduction, let us address ζ B,1 (� x) by "turning off" the second source, i.e. T B,2 = 0, so that ζ B,1 (� x) = T (� x) T B,1 . This defines the Monte Carlo weight to be used for estimating ζ B,1 (� x) as : w B,1 γ =

w γ T B,1 = T γ,end (� x γ,end ,1) T B,1 with T γ,end (� x γ,end , 1) = T B,1 if � x γ,end ∈ S 1 and T γ,end (� x γ,end , 1) = 0 if � x γ,end ∈ S 2 , namely, w B,1 γ = 1 if � x γ,end ∈ S 1 and w B,1 γ = 0 if � x γ,end ∈ S 2 . This writes w B,1 γ = H(� x γ,end ∈ S 1 ) ( 44
)
where H is a test function, taking the value 1 if the condition is valid and 0 otherwise. Similarly,

w B,2 γ = H(� x γ,end ∈ S 2 ) (45)
In algorithmic terms,

• N paths γ j are sampled;

• w B,1 γ j is computed for each path (1 if S 1 is reached, 0 otherwise); • w B,2 γ j is computed for each path (1 if S 2 is reached, 0 otherwise); • the Monte Carlo estimate of ζ B,1 (� x) is m B,1 = 1 N � N j=1 w B,1 γ j ; • the Monte Carlo estimate of ζ B,2 (� x) is m B,2 = 1 N � N j=1 w B,2
γ j . Then, Eq. (42) can be used to estimate the results for any set of source values.

Implementation in Stardis

In Stardis, we consider scenes where boundary conditions can be split into a set of constant and uniform sources: θ iBR does not depend on location, time and direction, and for each geometrical element i, T B,i , ϕ B,i and ψ i are constant and uniform, T I,i does not depend on location. λ i , ρ i , c i , h i , � i and α i are also considered uniform over element i.

Under these assumptions, ζ(� x, t, � x S , t S ) can be aggregated by geometrical element as done above in Eqs. ( 44) and (45). For this, we use the test function to build one weight expression for each geometrical element i and each type of source, following:

w

I,i γ = H(� x γ,end ∈ Ω i )H(µ γ,end = 0) w B,i γ = H(� x γ,end ∈ S i )H(µ γ,end = 1) w BF,i γ = H(� x γ,end ∈ S i )H(µ γ,end = 2) w BR γ = H(µ γ,end = 3) w ψ,i γ = N ψ � k=1 H(� x γ,ψ (k) ∈ Ω i )β ψ (� x γ,ψ (k)) w ϕ,i γ = N ϕ � k=1 H(� x γ,ϕ (k) ∈ S i )β ϕ (� x γ,ϕ (k)) (46)
These weights are then just a splitting of the weights given by Eq. 41 and only them need to be stored while solving the latter. The Monte Carlo algorithm is therefore modified to estimate:

m I,i = 1 N N � j=1 w I,i γ j m B,i = 1 N N � j=1 w B,i γ j m BF,i = 1 N N � j=1 w BF,i γ j m BR = 1 N N � j=1 w BR γ j m ψ,i = 1 N N � j=1 w ψ,i γ j m ϕ,i = 1 N N � j=1 w ϕ,i γ j (47)
Finally, once these estimates have been constructed, they can be used to estimate the addressed quantity for any set of sources:

T (� x, t) or T (t) or θ R (� ω, � x, t) ≈ m = N ω � i=1 m I,i T I,i + N S � i=1 m B,i T B,i + N S � i=1 m BF,i T BF,i +m BR θ BR + N ω � i=1 m ψ,i ψ i + N S � i=1 m ϕ,i ϕ i (48)

Uncertainty estimation

Since thermal paths are all used to estimate the different propagators, each propagator estimate is correlated to another. Furthermore, surfacic and volumetric sources are encountered along thermal paths and therefore each path does not visit only one single source (the final temperature). Altogether, the question of quantifying the uncertainty associated to T (� x, t) requires some attention as soon as it is computed using the data stored along a Monte Carlo ran with another set of sources.

If, in any applicative context, there was some need for estimating the uncertainty of each propagator (and not only of the resulting uncertainty on T (� x, t)), then, from Eq. ( 47), it appears clearly that a Monte Carlo approach is used for estimating each propagator as an average of dedicated weights. We could therefore compute an uncertainty associated to each of these estimates by computing the standard error of the Monte Carlo weights for each propagator. This would indeed provide a faithful information about the uncertainty with which each propagator is known for a given number of sampled thermal-paths. However, this information would be insufficient to estimate the uncertainty of the finally addressed quantity (T (� x, t) in Eq. 48): since the same set of thermal-paths has been used to estimate all the propagators (in one single Monte Carlo run), the propagator estimates are correlated. Estimating the uncertainty of T (� x, t) then requires to take into account the correlation matrix for all the propagators.

Leaving aside the idea of quantifying the uncertainty of each propagator and computing only the uncertainty on T (� x, t) gives us more freedom on how the stored data can be aggregated. All the required information is stored along each path (i.e. which source and how many times this source is encountered). When evaluating the temperature with the Green function, we could therefore recalculate the Monte Carlo weights and hence, recalculate at the same time the square weight to evaluate the uncertainty as in the case of a regular Monte Carlo computation (see Eq. ( 25)). The quantification would be made exactly as if the Monte Carlo had been re-run, with the same samples but other source values. This would include all the above mentioned correlations. However, operating this way may be cumbersome if a high amount of information has been stored (large number of sources, for instance). But then, we can simply gather the information into squared sums and cross-sums for volumetric/surfacic sources:

ssq I,i = N � j=1 � w I,i γ j � 2 ssq B,i = N � j=1 � w B,i γ j � 2 ssq BF,i = N � j=1 � w BF,i γ j � 2 ssq BR = N � j=1 � w BR γ j � 2 ssq ψ,i = N � j=1 � w ψ,i γ j � 2 ssq ϕ,i = N � j=1 � w ϕ,i γ j � 2 sc I,i,ψ,k = N � j=1 w I,i γ j w ψ,k γ j sc I,i,ϕ,k = N � j=1 w I,i γ j w ϕ,k γ j sc B,i,ψ,k = N � j=1 w B,i γ j w ψ,k γ j sc B,i,ϕ,k = N � j=1 w B,i γ j w ϕ,k γ j sc BF,i,ψ,k = N � j=1 w BF,i γ j w ψ,k γ j sc BF,i,ϕ,k = N � j=1 w BF,i γ j w ϕ,k γ j sc BR,ψ,k = N � j=1 w BR γ j w ψ,k γ j sc BR,ϕ,k = N � j=1 w BR γ j w ϕ,k γ j sc ψ,i,ϕ,k = N � j=1 w ψ,i γ j w ϕ,k γ j (49)
We then compute s as:

s = 1 √ N � 1 N � N ω � i=1 � ssq I,i T 2 I,i + N ω � k=1 sc I,i,ψ,k T I,i ψ k + N S � k=1 sc I,i,ϕ,k T I,i ϕ k � + N S � i=1 � ssq B,i T 2 B,i + N ω � k=1 sc B,i,ψ,k T B,i ψ k + N S � k=1 sc B,i,ϕ,k T B,i ϕ k � + N S � i=1 � ssq BF,i T 2 BF,i + N ω � k=1 sc BF,i,ψ,k T BF,i ψ k + N S � k=1 sc BF,i,ϕ,k T BF,i ϕ k � +ssq BR θ 2 BR + N ω � k=1 sc BR,ψ,k θ BR ψ k + N S � k=1 sc BR,ϕ,k θ BR ϕ k + N ω � i=1 ssq ψ,i ψ 2 i + N S � i=1 ssq ϕ,i ϕ 2 i + N ω � i=1 N S � k=1 sc ψ,i,ϕ,k ψ i ϕ k � -m 2 � 1/2
(50)

Implementation

The implementation of stardis-solver, that is presented here, is a reference implementation suitable for execution with conventional computing resources (low-end personal computer).

The source code of the solver is designed to be easy to understand and suitable for training purposes. Users can then rely on this implementation and make it evolve according to their needs.

The current implementation is a compromise between the different possibilities described in Sec. 5. This compromise consists in:

• Grouping the terms related to volume power densities and heat flux densities, restraining heat flux and power to be uniform over time and space.

• Keeping all positions and times for other sources (initial temperature, ambient radiation temperature, fluid temperature, imposed temperature), allowing these sources to vary, either over time or space or both.

Although most of the propagators that we compute in practice are integrals of the Green function over the system parts where the sources are uniform, we make use of the post-fix green for the corresponding parts of the code. The three main parts of interest and detailed below are:

(1) Data structures used to store the Green function (see Sec. 6.2).

(2) Functions used to fill up these data structures in the construction of the Green function (see Sec. 6.3).

(3) Functions using these data structures to evaluate the temperature for a given set of source values (see Sec. 6.4).

The source code for the data structures and functions described in these three sections are grouped in Appendix A at the end of the document.

Data structures

When building the Green function, the Monte-Carlo weights are not computed, but the data needed to compute them is stored, path by path, for later use. This storage requires two different types: one to store the data collected along individual Green paths, and another one to store the Green function itself, including the data of the sampled Green paths, as well as all the shared data referenced by the paths (materials, interfaces, ...):

<green data structures> = <green path data structure> <green function data structure> Path storage. Green paths are constructed by stardis-solver, following the very same algorithm as when evaluating a temperature. The difference between temperature computation and the construction of the Green function is that when a path is sampled, some of the information is stored in the data structure corresponding to the Green path instead of being used on the fly to compute a temperature (see List. 2). Then, each path sampled by the solver results in a Green path data structure storing the information as follows:

<green path data structure> = struct green_path { <path elapsed time> <flux density terms collection> <power density terms collection> <end of path> <miscellaneous variables> };

Elapsed time is trivially a double: <path elapsed time> = double elapsed_time;

Flux and power terms encountered along the path are partially merged and stored in dynamic arrays. Merging is done by material and interface: all contributions along the path are accumulated and stored as a single term associated with a given material or interface. As flux terms can only appear at interfaces, merged flux terms consist of an interface identifier, the involved side and the corresponding cumulated flux value. On the other hand, power terms appear in media, thus merged power terms consist of a medium identifier and the cumulated power value. The end of the path can be of three types: at a boundary (fragments), in a volume (vertex), or a radiative exchange with the surrounding environment. This end of the path is represented by an union which is interpreted according to the value of the field limit_type, which also allows to interpret limit_id as being an identifier of medium (case in volume) or interface (case at boundary); note that the radiative case requires neither an union member nor a limit_id: <end of path> = union { struct sdis_rwalk_vertex vertex; struct sdis_interface_fragment fragment; } limit; unsigned limit_id; enum sdis_green_path_end_type end_type;

Green function storage. The main structure is used to store everything allowing the later evaluation of a temperature estimator. This includes the description of the sampled paths as well as all the shared data referenced by the sampled paths (see List. 1).

<green function data structure> = struct sdis_green_function { <media collection> <interfaces collection> <paths collection> <miscellaneous variables> };

Collections of media and interfaces accumulate the media and interfaces that have been visited when constructing the Green function. Individual paths can then reference this shared information. These collections are hash tables, i.e. associative containers that favor fast and constant time lookup, to ensure unique storage of only the media and interfaces visited by the paths: <media collection> = struct htable_medium media; <interfaces collection> = struct htable_interf interfaces;

The paths collection is the set of the paths sampled for the construction of the Green function. It is a dynamic array that is well suited for iterative storage and path iteration: <paths collection> = struct darray_green_path paths;

Building propagation information

Various functions are needed to fill in the Green function's data structures. They can be divided in two groups:

<build functions> = <functions building the green function> <functions building green paths>

Functions building the Green function itself do not need further description, as they are limited to collections management. Functions building green paths are divided in two groups: <functions building green paths> = <functions that store path ending> <functions that accumulate data along a path>

Storing path ending

Since there are three different ways to end a path, there are three different functions that can be called to store information about the end of paths (see List. The current path ends at an interface that must be available at the time the Green function is evaluated. We start by making sure that is the case, returning an error if the process fails.

<register interface 'interf' against the green function> = res = ensure_interface_registration(handle->green, interf); if(res != RES_OK) return res;

Then the elapsed time is stored:

<store path duration> = handle->path->elapsed_time = elapsed_time;

Then the information related to the location at interface is stored (including position, normal, parametric coordinates and time):

<store the location at interface> = handle->path->limit.fragment = *frag;

Then the interface identifier is stored:

<store identifier of interface 'interf'> = handle->path->limit_id = interface_get_id(interf);
Finally, the type of the path ending is stored:

<store the path ends up at an interface> = handle->path->end_type = SDIS_GREEN_PATH_END_AT_INTERFACE;

The other two functions are built using the same pattern and are sketched thereafter:

Accumulating data along a path

The data accumulated along a path is volume power density terms and flux density terms.

<functions that accumulate data along a path> = <register a power term> <register a flux term>

It is performed through the following two functions, that share the same pattern: Step number Path position

β ψ (� x γ,ψ (1)) value µ γ i ,end identifier 1 � x γ j ∈ Ω i 0 - 2 � x γ j ∈ Ω 1 δ 2 2 6λ 1 - 3 � x γ j ∈ Ω 1 δ 2 2 6λ 1 + δ 2 3 6λ 1 - 4 � x γ j ,end ∈ Ω i δ 2 2 6λ 1 + δ 2 3 6λ 1 0
Table 1: Successive data stored (framed in green rectangles and referred as green path ) during Monte Carlo calculation for the first Green path.

The first path starts at the probe point � x. Contributions β ψ (� x γ,ψ (1)) related to the power density term Ψ 1 in the medium Ω 1 will accumulate when the path crosses the medium, i.e. during the second and third jump (δ 2 and δ 3 ). Path ends with the initial temperature condition T I associated to the medium Ω i (see Fig. 10).

Monte Carlo weight would be as follows:

w γ j = ψ 1 β ψ (� x γ,ψ (1)) + T γ,end (� x γ,end , t γ,end , � ω γ,end , 0) (51)
Data stored for the Green path (framed in green in Tab. 1) are: successive positions � x γ j , power density term contribution β ψ (� x γ,ψ (1)) along the path and the identifier of the boundary condition encountered µ γ,end = 0 (T i,I : initial condition type-0). Step number Path position

β ϕ (� x γ,ϕ (i)) value µ γ i ,end identifier 1 � x γ j ∈ Ω i 0 - 2 � x γ j ∈ S i 0 - 3 � x γ j ∈ Ω i δ 2 2 λ i - 4 � x γ j ∈ S i δ 2 2 λ i - 5 � x γ j ∈ Ω i δ 2 2 λ i + δ 2 4 λ i - 6 � x γ j ,end ∈ S 2 δ 2 2 λ i + δ 2 4 λ i 1
Table 2: Successive data stored (framed in green rectangles and referred as green path ) during Monte Carlo calculation for the second Green path.

The second path starts at the probe point � x. Contributions β ϕ (� x γ,ϕ (i)) related to flux density ϕ i at S i interface will accumulate when the path hits this interface, i.e. on the second and fourth jump (δ 2 et δ 4 ). The path ends at the interface S 2 with the temperature imposed on this interface T B (see Fig. 11).

Monte Carlo weight would be as follows:

w γ j = ϕ i β ϕ (� x γ,ϕ (i)) + T γ,end (� x γ,end , t γ,end , � ω γ,end , 1) (52)
Data stored for the Green path (framed in green in Tab. 2) are: successive positions � x γ j , density flux contribution β ϕ (� x γ,ϕ (i)) along the path and the identifier of the boundary condition encountered µ γ,end = 1 (T B : type-1 boundary condition).

Using propagation information

Once a set of sampled paths is stored in the Green function, dedicated evaluation functions are needed to apply the Green function to a set of source values (see List. 5). These functions are presented thereafter. Note that when using helper functions which are also used for standard Monte-Carlo computations ( solid_get_volumic_power , interface_side_get_flux ), they have to provide a vertex, even though the computation is time and space independent in the Green function context. The volume power along the path is computed by considering each medium encountered along the path and accumulating the corresponding volume power contribution. Each volume power term power_terms[i] (see Eq. ( 46)),

Simulation examples

This section illustrates a typical use of the storage and use of the propagation information described in the previous sections. The geometrical and physical descriptions of the configurations used for this illustration, as well as the Stardis input files, are available in the enclosed zip file. The objective is essentially to show that a computation performed using the stored propagation information recovers the result that would be obtained with a complete MC run (with a particular attention to the associated statistical errorbars) and to illustrate the benefits in terms of computation times. As far as validation is concerned, we concentrate on the parts of the code constructing and using the propagators, not on the main code itself: Stardis is already validated elsewhere [START_REF] Eymet | Synthèse d'images infrarouges sans calcul préalable du champ de température[END_REF]19]. However, we still reproduce parts of this validation by providing, together with each simulation example, a systematic comparison with the solution computed with a standard deterministic solver [START_REF]Comsol multiphysics ®[END_REF] (referred as COMSOL Multiphysics ® hereafter and "Deterministic" in figures). c)). This benchmark has been already used by Sans et. al [START_REF] Sans | Modeling the Flash Method by using a Conducto-Radiative Monte-Carlo Method: Application to Porous Media[END_REF] for the purpose of validating Monte Carlo simulations of coupled diffusion-radiation heat transfer. The open-porosity configuration corresponds to a heterogeneous 3D honeycomb that can be assimilated to a porous medium with open channels, like a heat exchanger configuration. The closed-porosity configuration has 22 enclosed cavities and may be assimilated to an insulation material.

The physical assumptions are those of Sec. 3.1 with same values for λ, ρ, c and ψ throughout the whole solid phase and same values for h and � along all the solid-fluid interfaces. For open porosity, there is one single imposed fluid temperature T BF (type-3 boundary condition), the same inside the channels and outside the system. For closed porosity, T BF is only imposed outside the system (fluid cells temperatures are free). In both cases, the incoming radiance temperature outside the system θ BR is uniform and isotropic, the solid temperature T B is imposed at the top surface (type-1 boundary condition) and a flux density ϕ B is imposed at the bottom surface (type-2 boundary condition). The initial temperature T I is uniform.

The estimated quantity is the temperature T (� x c , t) at the center � x c of the geometry for a given observation time t as a function depending on the six available sources:

• the initial temperature T I ,

• the top boundary solid temperature T B ,

• the ambient fluid temperature T BF ,

• the ambient radiance temperature θ BR ,

• the flux density at the bottom boundary ϕ,

• the power density throughout the solid phase ψ. 12 gathers all the computation times, illustrating that the benefits of using the stored propagation information, instead of running the Monte Carlo, is a computation time reduction by a factor 10 3 to 10 4 [START_REF] Penazzi | Transfer function estimation with SMC method for combined heat transfer: insensitivity to detail refinement of complex geometries[END_REF].

In closer details, the following comments can be made:

• "Propagator" (using the propagator) and "Monte-Carlo" results (re-running a Monte Carlo for control) are in perfect agreement, as expected, because "Propagator" is statistically rigorously equivalent to re-running the Monte Carlo (see . However, for the validation runs, new random numbers are used to sample the paths, whereas all "Propagator" simulations are based on the same sampled paths. Therefore, although the errorbars associated to "Propagator" can be fully trusted, the simulations made for various values of the sources are all correlated: typically, there are no statistical fluctuations in the errorbars, and when a simulation result for a given source value happens to be below the reference (within the errorbar, otherwise the validation would have failed), it remains below the reference for all other values of the source.

• "Propagator" (or "Monte-Carlo") and "Deterministic" are in perfect agreement as long as the linearization of heat transfer remains relevant. It is well-known that, exchanged radiative heat flux being proportional to T 4 , radiative heat transfer causes non-linear propagation. Moreover, the higher the thermal gradients are, the higher such non-linear effects occur and the larger the bias induced by radiative transfer linearization. Here, without radiation and/or any thermal dependance of the thermal properties, "Propagator" predict the correct temperature for a very wide range of sources values (see Stardis linearizes radiation where COMSOL Multiphysics ® does not. Hence, the gaps observed between the "Propagator" and the "Deterministic" method come from the capacity of a given source to increase thermal gradients, and thus to strengthen non-linearity effects (see working on this issue, starting from the non-linear Monte Carlo approach given by [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF]). The latest available versions of Stardis already include a first representation of the T 4 radiation dependance, by use of branching statitsics, but the complete evaluation of this new potential is still under progress. In any case, when nonlinear effects are significant, computing the propagator is meaningless: the descriptions made in the present article are strictly rooted in linear physics and their application will always remain strictly restricted to linear heat transfer.

• There is a notable distinction between open-porosity and closed-porosity as far as computation time is concerned. In the open-porosity case, a thermal path starting at the center of the system can encounter the flow quite rapidly and then the path is stopped because the fluid temperature is known. In the closed-porosity case, a thermal path starting at the center of the system can also encounter the flow inside a fluid cell quite rapidly, but then the flow temperature is unknown and the path is continued until a source is encountered (either the initial condition at any location, or a known temperature at the boundary). These paths are significantly longer and so is the computation time required for their construction. It can even be extremely long if the number of closed cells is increased: Monte Carlo may encounter difficulties when addressing insulating materials with large numbers of closed cells along all directions. However, this difficulty is not reported in "Propagator": even if the computation cost associated to path-sampling in the reference Monte Carlo run is higher in the closedporosity case, the propagation information stored are similar in the open-porosity and the closed-porosity cases and the computation times associated to the use, by "Propagator", of these propagation information are the same in both cases. The computation benefit of using "Propagator" instead of running a full Monte Carlo is therefore stronger for closed porosities.

Going further: Propagation information storage in case of non-uniform and time-dependent sources

When constructing estimates for the propagators, we have split the initial Monte Carlo weight and gathered all the propagation information corresponding to parts of the system where the sources were uniform and constant (see Sec. 5). As was already mentioned, Stardis makes uniformity and constancy assumptions for each geometrical part but 'stardis-solver' does not. And indeed, in the Monte Carlo weight expression of Eq. ( 41), nothing prevents the initial temperatures, the imposed temperatures at the boundary, the incoming radiance temperature, the surface flux densities imposed at the boundary and the volume power densities imposed inside the solids to be non-uniform and non-constant fields:

ψ � � x γ,ψ (k), t γ,ψ (k) � , ϕ � � x γ,ϕ (k), t γ,ϕ (k)
� and T γ,end (� x γ,end , t γ,end , � ω γ,end , µ γ,end ) can hold this information as it is provided, whatever the geometric structure. This means that even running Stardis with uniform and constant sources in each part, if we store the locations, times and directions used when constructing Monte Carlo weights, then this information can be later used to virtually re-run the Monte Carlo simulation with non-uniform time-dependent sources.

In algorithmic terms:

• N thermal paths are sampled exactly as in the Monte Carlo algorithm of Sec. 4.1;

• Along each path γ j , we store the N ϕ, j locations � x γ j ,ϕ (k) and times t γ j ,ϕ (k) at which surface flux densities were accessed;

• Similarly we store the N ψ, j locations � x γ j ,ψ (k) and times t γ j ,ψ (k) at which volume power densities were accessed;

• The information concerning the end of the path are also stored: � x γ j ,end , t γ j ,end , � ω γ j ,end and µ γ j ,end ;

• When a Monte Carlo estimate is required for any set of non-uniform time-dependent source fields TI , TB , TBF , θBR , ψ and φ, is constructed the following way:

T (� x, t) or T (t) or θR (� ω, � x, t) ≈ m = 1 N N � j=1 w j (53)
with w j = H(µ γ j ,end = 0) TI (� x γ j ,end )

+H(µ γ j ,end = 1) TB (� x γ j ,end , t γ j ,end )

+H(µ γ j ,end = 2) TBF (� x γ j ,end , t γ j ,end )

+H(µ γ j ,end = 3) θBR (� x γ j ,end , t γ j ,end , � ω γ j ,end )

+ N ψ, j � k=1 β ψ (� x γ j ,ψ (k)) ψ � � x γ j ,ψ (k), t γ j ,ψ (k) � + N ϕ, j � k=1 β ϕ (� x γ j ,ϕ (k)) φ � � x γ j ,ϕ (k), t γ j ,ϕ (k) � (54)
For this last strategy, dealing with uncertainties is straightforward. The estimate m in Eq. ( 53) is finally built exactly as if the Monte Carlo was re-run using new sources. The uncertainty s associated to m is therefore the same as for any Monte Carlo simulation (using an estimate of the standard deviation of the Monte Carlo weights) as detailed in Eq. [START_REF] Dunn | Monte Carlo methods for design and analysis of radiation detectors[END_REF].

In terms of code changes, storing unmerged flux and power terms allows code simplification, as a significant part of the merged-storage version of the code is about retrieving the term to which the current contribution needs to be merged. On the other hand, allowing all types of non-uniform and time-dependent sources is still at the cost of storing more propagation since information terms data structures have to store the term's location (time and space). The data structures and functions that implement the unmerged storage of sources' contributions are available in Appendix B.

Conclusion

One of the strengths of Monte Carlo approaches is the ease with which information can be stored, during one run, and then used offline to learn about the physics, preserving all geometric features. The simplest example is the storing of the paths themselves (or of a large enough fraction of the paths). For coupled heat transfer, displaying a selection of the paths and analyzing how they visit the system, both in time and space, switching from one heat transfer mode to the other, is indeed a very practical way to learn how the sources are viewed from one location, how their impact is delayed by the inertial parts of the system, and therefore how a design can be adjusted to achieve a given objective. In the present paper, we have left aside these details about the paths themselves. We focused on the act of quantifying the propagation and not on the analysis of the physical phenomena and the coupling processes responsible for the propagation. But these two practices, computing the propagators and visualizing the processes, are worth being considered sideways in all engineering contexts requiring a close understanding of heat transfer physics at the system scale. Therefore, in addition to the functionalities of Stardis described in Sec. 6, a set of post-treatment tools have also been designed to help visualizing and analyzing thermal paths throughout large scale geometries [19].

For automated engineering practices, e.g. inversion, optimization or command algorithms, analyzing the paths is useless; all is needed is the addressed quantity as a function of the sources: the tools described in the present paper are therefore self-sufficient. However there are numerous questions of direct interest to thermal engineers that cannot be addressed this way. These are all the dependencies on parameters that cannot be considered as sources (in the general sense provided by Green's theory). Typically the dependence on emissivities, convective exchange coefficients, conductivities or capacities rises more complex questions. If only sensitivities were required, i.e. the derivative with respect to each parameter, then the general theory of sensitivity evaluation in Monte Carlo algorithms could be used [START_REF] Lapeyre | Monte-Carlo and domain-deformation sensitivities[END_REF], but we would not build the complete dependence (the function) as we did here with the sources. Addressing the complete non-linear dependence on other parameters than sources is not theoretically unfeasible: it has notably been achieved in the field of radiative transfer under the literature name of "Symbolic Monte Carlo" [START_REF] Dunn | Inverse Monte Carlo analysis[END_REF][START_REF] Dunn | Monte Carlo methods for design and analysis of radiation detectors[END_REF][START_REF] Galtier | A symbolic approach for the identification of radiative properties[END_REF][START_REF] Maanane | Symbolic Monte Carlo method applied to the identification of radiative properties of a heterogeneous material[END_REF] and we have started to work on extending these symbolic techniques to coupled heat transfer, with the objective of implementing them inside stardis-solver [START_REF] Penazzi | Toward the use of Symbolic Monte Carlo for Conduction-Radiation Coupling in Complex Geometries[END_REF][START_REF] Sans | Méthode de Monte-Carlo Symbolique pour la caractérisation des propriétés thermiques : application à la méthode flash[END_REF]. By far more difficult would be the question of addressing the dependence on geometrical parameters. Here also, some attempts have already been made in the field of radiative transfer, but to the best of our knowledge and although large impacts could be expected in terms of applications, there is no report available of any attempt to go beyond the computation of derivatives (geometrical sensitivities). Constructing a thermal heat transfer observable as a symbolic function of a geometric parameter is a question that has not yet been addressed.

Another difficult point associated to strong applicative concerns is the withdraw of the linearization of radiative transfer. This linearization is at the heart of present Monte Carlo approaches to coupled heat transfer. There are convincing perspectives as far as handling non-linearities in the Monte Carlo framework is concerned [START_REF] Dauchet | Addressing nonlinearities in Monte Carlo[END_REF], and some of the corresponding propositions could be used to avoid the linearization of radiation, but then the overall coupled physical problem would be non-linear and the concept of propagation could not be used anymore. All our present proposition would then have to be revisited.

Finally, we have highlighted the fact that storing propagators allows the design of fast external solvers: adressing the same observable for other sources without reruning the Monte Carlo. But another quite significant usage is uncertainty propagation. As the model is linear, any statistical distribution of source-uncertainties is propagated without transformation and the propagator tells us how (exactly the same way it computes the observabale when changing the sources). Typically, if a given source has a normal distribution, the technique described in the present paper will provide the center of the resulting normal distribution of the observable and computing its standard deviation will require nothing more than the square of the propagator. If the Monte Carlo run is very accurate, there will be no need for any consideration of the Monte Carlo uncertainties when studying such source-uncertainty propagations. Of course, if the source uncertainty and the Monte Carlo uncertainty are of the same order of magnitude, then further attention will have to be devoted to the resulting statistical correlations, typically when several sources are uncertain, possibly correlated, their associated propagators being themselves correlated because they were computed using the same Monte Carlo run. But again, all the required information is in the data already stored at the Monte Carlo weight level when running Stardis and only a post treatment is required when such further statistics are required. h a n d l e ->p a t h -> e l a p s e d _ t i m e = e l a p s e d _ t i m e ; h a n d l e ->p a t h -> l i m i t . f r a g m e n t = * f r a g ; h a n d l e ->p a t h -> l i m i t _ i d = i n t e r f a c e _ g e t _ i d ( i n t e r f ) ; h a n d l e ->p a t h -> e n d _ t y p e = SDIS_GREEN_PATH_END_AT_INTERFACE ; r e t u r n RES_OK ; } r e s _ T g r e e n _ p a t h _ s e t _ l i m i t _ v e r t e x ( s t r u c t g r e e n _ p a t h _ h a n d l e * h a n d l e ,

Figure 1 :

 1 Figure 1: P 1 light source on and P 2 off.Figure 2: P 2 light source on and P 1 off.

Figure 2 :

 2 Figure 1: P 1 light source on and P 2 off.Figure 2: P 2 light source on and P 1 off.

Figure

  Figure3: An heterogeneous, multiple scattering and emitting/absorbing medium D is bounded by a partially reflecting wall ∂D. Physical images associated to the evaluation of a local radiative quantity at the location x 0 (for instance, the monospectral radiative intensity) with a reverse Monte Carlo algorithm are as follows. Photons are followed from the location x 0 in the unit direction -� u 0 until the emission/absorption location, either in the volume D or on the wall ∂D (see x b,13 ). At each step j of the optical path within the medium, a free path l j is sampled according to the exponential Beer law, enabling the computation of the next step location x j = x j-1 -l j � u j-1 . If the location x j is in the medium D, the event may be an absorption, a scattering (see x 1 ,x 3 ,x 4 ,x 6 ,x 7 ,x 9 ,x 11 and x 12 ) or a null-collision (see x 2 and x 10 ). A null-collision event corresponds to a pure-forward scattering event in which � u j = � u j-1 (see box 2.1). If the location x j reaches the wall ∂D, the event may be an absorption (see x b,13 ) or a reflexion (see x b,5 and x b,8 ).

Figure 5 :

 5 Figure 5: Illustration of a random path sampling compatible with path-tracing acceleration techniques of computer graphics. Each sphere has the same radius δ s and is adjusted when getting close to the domain boundary.

Figure 7 :

 7 Figure 7: A coupled situation is illustrated in which fluid areas and solid areas, represented respectively in blue and brown surfaces, involve heat transfer by conduction, convection and radiation. The green dot in the central solid area is the probe position from which the paths start. Conductive sub-paths are represented in brown, radiative sub-paths are represented in red with a broken line symbol and convective sub-paths are represented in blue with a symbol indicating a capacity. In this example, the paths can either end in the domain, within the solid or fluid areas, if the initial condition has been reached, or on the boundary of the outer rectangle acting as boundary conditions. Ten paths are being displayed. They have been produced using a simplified version of Stardis that we use in didactic contexts (http://www.edstar.cnrs.fr/prod/fr/training/tool/).

Figure 8 :

 8 Figure 8: Sketch of the general configuration. The system boundary is S. Three internal volumes Ω 1 , Ω 2 and Ω 3 are represented. Inside Ω 1 , there is a fluid inclusion Ω i and a sub-solid Ω i+2 . I i is the solid/fluid boundary between Ω 1 and Ω i . ∂Ω i+2 is the solid boundary between Ω i+2 and Ω 1 . Similarly, Ω 2 contains a fluid inclusion Ω i+1 and I i+1 is the solid/fluid boundary. The last volume Ω 3 contains a solid volume Ω i+3 and the solid boundary is ∂Ω i+3 . Conduction occurs in the solid volumes, radiation and convection occur in the different fluid volumes. h is the convective heat transfer coefficient.

5 .

 5 Stardis : storing the propagation data 5.1. Illustration of the principle in the case of two sources

Figure 9 :

 9 Figure 9: Schemes of the solid cube with two isothermal faces S 1 and S 2 , respectively at temperature T B,1 and T B,2 , with a location � x inside the solid volume.

  <flux density terms collection> = struct darray_flux_term flux_terms; <power density terms collection> = struct darray_power_term power_terms;

  unsigned id; /* Identifier of the medium */ };

6. 3 . 3 .

 33 Examples of Green paths vs. Monte Carlo paths Two examples of paths sampled by the Stardis probe temperature solver are shown in Figs. 10 and 11. For each example, the archived information for the construction of the Green path is shown in Tabs. 1 and 2.

Figure 10 :

 10 Figure 10: First path example.

Figure 11 :

 11 Figure 11: Second path example.

  <for each green path stored into 'green'> { <compute the weight of this green path> <accumulate the resulting weight> } <setup the estimator> <finalize the solve function> }The computation of the weight associated to a path is done with a dedicated function (see List. 4): <compute the weight of this green path> = double w; /* Monte Carlo weight to compute */ res = green_function_solve_path(green, ipath, &w); <handle error code returned in 'res'> The dedicated function simply takes into account the different contributions that have been stored in the current path (power, flux and end of path), and uses them to produce the Monte-Carlo weight of the path: <evaluation functions> += res_T green_function_solve_path (struct sdis_green_function* green, const size_t ipath, double* weight) { <local variables> <check input arguments> <compute the power collected along the path> <compute the flux collected along the path> <fetch the end of path> <compute the overall Monte-Carlo weight> <finalize the solve_path fuction> }

  <compute the overall Monte-Carlo weight> = *weight = power + flux + end_temperature;

  Two academic configurations are considered. They are designed as simplified versions of porous media, one with open porosities (see Fig. C.13 a)), the other with closed porosities (see Fig. C.13

  Tests are conducted first without radiation (� = 0, see Fig. C.14, Fig. C.15 and Fig. C.16 ) and then with radiation (black surfaces, � = 1, see Fig. C.17, Fig. C.18 and Fig. C.19). For each case, the propagation information are stored using a single Monte Carlo run. These propagation information are then used to predict T (� x c , t) (and to estimate its uncertainty) when varying the sources values with factors in the [10 -2 , 10 2 ] range around a fixed reference value for each source: T ref I , T ref B , T ref BF , θ ref BR , ϕ ref and ψ ref (results labeled "Propagator" in the figures). Validation is achieved by comparing with standard Monte Carlo results, labeled "Monte-Carlo" in the figures. The perfect adequacy between the "Propagator" and "Monte-Carlo" results in Figs. C.14, C.15,C.16 C.17, C.18 and C.19 validates the implemented code and the quality of the stored propagation information. Fig.

  Fig. C.14 a) and Fig. C.15 a))

Fig. C. 14 ,

 14 Fig. C.15 and Fig. C.16). Close to the set of reference values for sources, temperature gradients were purposely chosen as low as suitable for the frame of linear heat transfer. Thus, "Propagator" and "Deterministic" are in good agreement (see Fig. C.17 b)). Outside this range of source values, "Propagator" fails to predict the correct temperature field because the linearization of radiation, at the heart of Stardis, becomes meaningless (see Fig. C.17 a)). Here, two different physical models are solved:

  Fig. C.17, Fig. C.18 and Fig. C.19). We are presently Open-porosity geometry without radia�ve transfer Open-porosity geometry with rad a�ve transfer Closed-porosity geometry without radia�ve transfer Closed-porosity geometry with radia�ve transfer

Figure 12 :

 12 Figure 12: Ratio of computation times t Propagator t MC for the corresponding t MC . t Propagator is the computation time for the propagator function and t MC is the computation time for the corresponding Monte Carlo computation.

  s t r u c t s d i s _ m e d i u m * mdm, c o n s t s t r u c t s d i s _ r w a l k _ v e r t e x * v e r t , c o n s t d o u b l e e l a p s e d _ t i m e ) { r e s _ T r e s = RES_OK ; ASSERT ( h a n d l e && mdm && v e r t ) ; ASSERT ( h a n d l e ->p a t h -> e n d _ t y p e == SDIS_GREEN_PATH_END_TYPES_COUNT__ ) ; r e s = e n s u r e _ m e d i u m _ r e g i s t r a t i o n ( h a n d l e ->g r e e n , mdm ) ; i f ( r e s != RES_OK ) r e t u r n r e s ; h a n d l e ->p a t h -> e l a p s e d _ t i m e = e l a p s e d _ t i m e ; h a n d l e ->p a t h -> l i m i t . v e r t e x = * v e r t ; h a n d l e ->p a t h -> l i m i t _ i d = m e d i u m _ g e t _ i d (mdm ) ; h a n d l e ->p a t h -> e n d _ t y p e = SDIS_GREEN_PATH_END_IN_VOLUME ; r e t u r n RES_OK ; } r e s _ T g r e e n _ p a t h _ s e t _ l i m i t _ r a d i a t i v e ( s t r u c t g r e e n _ p a t h _ h a n d l e * h a n d l e , c o n s t d o u b l e e l a p s e d _ t i m e ) { ASSERT ( h a n d l e ) ; ASSERT ( h a n d l e ->p a t h -> e n d _ t y p e == SDIS_GREEN_PATH_END_TYPES_COUNT__ ) ; h a n d l e ->p a t h -> e l a p s e d _ t i m e = e l a p s e d _ t i m e ; h a n d l e ->p a t h -> e n d _ t y p e = SDIS_GREEN_PATH_END_RADIATIVE ; r e t u r n RES_OK ; } r e s _ T g r e e n _ p a t h _ a d d _ p o w e r _ t e r m ( s t r u c t g r e e n _ p a t h _ h a n d l e * h a n d l e , s t r u c t s d i s _ m e d i u m * mdm, c o n s t s t r u c t s d i s _ r w a l k _ v e r t e x * v t x , c o n s t d o u b l e v a l ) { s t r u c t g r e e n _ p a t h * p a t h ; s t r u c t p o w e r _ t e r m * t e r m s ; s i z e _ t n t e r m s ; s i z e _ t i t e r m ; u n s i g n e d i d ; r e s _ T r e s = RES_OK ; ASSERT ( h a n d l e && mdm && v t x ) ; / * Unused p o s i t i o n and t i m e : t h e c u r r e n t i m p l e m e n t a t i o n o f t h e g r e e n f u n c t i o n * a s s u m e s t h a t t h e power i s c o n s t a n t i n s p a c e and t i m e p e r medium . * / ( v o i d ) v t x ;r e s = e n s u r e _ m e d i u m _ r e g i s t r a t i o n ( h a n d l e ->g r e e n , mdm ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; p a t h = h a n d l e -> p a t h ; t e r m s = d a r r a y _ p o w e r _ t e r m _ d a t a _ g e t (& p a t h ->p o w e r _ t e r m s ) ; n t e r m s = d a r r a y _ p o w e r _ t e r m _ s i z e _ g e t (& p a t h ->p o w e r _ t e r m s ) ; i d = m e d i u m _ g e t _ i d (mdm ) ; i t e r m = SIZE_MAX ; / * E a r l y f i n d * / i f ( p a t h -> i l a s t _ m e d i u m < n t e r m s && t e r m s [ p a t h -> i l a s t _ m e d i u m ] . i d == i d ) { i t e r m = p a t h -> i l a s t _ m e d i u m ; } e l s e { / * L i n e a r s e a r c h o f t h e s u b m i t t e d medium * / FOR_EACH( i t e r m , 0 , n t e r m s ) i f ( t e r m s [ i t e r m ] . i d == i d ) break ; } / * Add t h e power t e r m t o t h e p a t h w r t t h e s u b m i t t e d medium * / i f ( i t e r m < n t e r m s ) { t e r m s [ i t e r m ] . t e r m += v a l ; } e l s e { s t r u c t p o w e r _ t e r m t e r m = POWER_TERM_NULL__ ; t e r m . t e r m = v a l ; t e r m . i d = i d ; r e s = d a r r a y _ p o w e r _ t e r m _ p u s h _ b a c k (& h a n d l e ->p a t h ->p o w e r _ t e r m s , &t e r m ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; } / * R e g i s t e r t h e s l o t i n t o w h i c h t h e l a s t a c c e s s e d medium l i e s * / CHK( i t e r m < UINT16_MAX ) ; p a t h -> i l a s t _ m e d i u m = ( u i n t 1 6 _ t ) i t e r m ; e x i t : r e t u r n r e s ; e r r o r : g o t o e x i t ; } r e s _ T g r e e n _ p a t h _ a d d _ f l u x _ t e r m ( s t r u c t g r e e n _ p a t h _ h a n d l e * h a n d l e , s t r u c t s d i s _ i n t e r f a c e * i n t e r f , c o n s t s t r u c t s d i s _ i n t e r f a c e _ f r a g m e n t * f r a g , c o n s t d o u b l e v a l ) { s t r u c t g r e e n _ p a t h * p a t h ; s t r u c t f l u x _ t e r m * t e r m s ; s i z e _ t n t e r m s ; s i z e _ t i t e r m ; u n s i g n e d i d ; r e s _ T r e s = RES_OK ; ASSERT ( h a n d l e && i n t e r f && f r a g && v a l >= 0 ) ; r e s = e n s u r e _ i n t e r f a c e _ r e g i s t r a t i o n ( h a n d l e ->g r e e n , i n t e r f ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; p a t h = h a n d l e -> p a t h ; t e r m s = d a r r a y _ f l u x _ t e r m _ d a t a _ g e t (& p a t h -> f l u x _ t e r m s ) ; n t e r m s = d a r r a y _ f l u x _ t e r m _ s i z e _ g e t (& p a t h -> f l u x _ t e r m s ) ; i d = i n t e r f a c e _ g e t _ i d ( i n t e r f ) ; i t e r m = SIZE_MAX ; / * E a r l y f i n d * / i f ( p a t h -> i l a s t _ i n t e r f < n t e r m s && t e r m s [ p a t h -> i l a s t _ i n t e r f ] . i d == i d && t e r m s [ p a t h -> i l a s t _ i n t e r f ] . s i d e == f r a g -> s i d e ) { i t e r m = p a t h -> i l a s t _ i n t e r f ; } e l s e { / * L i n e a r s e a r c h o f t h e s u b m i t t e d i n t e r f a c e * / FOR_EACH( i t e r m , 0 , n t e r m s ) { i f ( t e r m s [ i t e r m ] . i d == i d && t e r m s [ i t e r m ] . s i d e == f r a g -> s i d e ) { break ; } } } / * Add t h e f l u x t e r m t o t h e p a t h w r t t h e s u b m i t t e d i n t e r f a c e * / i f ( i t e r m < n t e r m s ) { t e r m s [ i t e r m ] . t e r m += v a l ; } e l s e { s t r u c t f l u x _ t e r m t e r m = FLUX_TERM_NULL__ ; t e r m . t e r m = v a l ; t e r m . i d = i d ; t e r m . s i d e = f r a g -> s i d e ; r e s = d a r r a y _ f l u x _ t e r m _ p u s h _ b a c k (& h a n d l e ->p a t h -> f l u x _ t e r m s , &t e r m ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; }/ * R e g i s t e r t h e s l o t i n t o w h i c h t h e l a s t a c c e s s e d i n t e r f a c e l i e s * / CHK( i t e r m < UINT16_MAX ) ; p a t h -> i l a s t _ i n t e r f = ( u i n t 1 6 _ t ) i t e r m ; e x i t : r e t u r n r e s ; e r r o r :g o t o e x i t ; } Listing 4: green_function_solve_path s t a t i c r e s _ T g r e e n _ f u n c t i o n _ s o l v e _ p a t h ( s t r u c t s d i s _ g r e e n _ f u n c t i o n * g r e e n , c o n s t s i z e _ t i p a t h , d o u b l e * w e i g h t ) { c o n s t s t r u c t p o w e r _ t e r m * p o w e r _ t e r m s = NULL ; c o n s t s t r u c t f l u x _ t e r m * f l u x _ t e r m s = NULL ; c o n s t s t r u c t g r e e n _ p a t h * p a t h = NULL ; c o n s t s t r u c t s d i s _ m e d i u m * medium = NULL ; c o n s t s t r u c t s d i s _ i n t e r f a c e * i n t e r f = NULL ; s t r u c t s d i s _ s c e n e * s c n = NULL ; s t r u c t s d i s _ r w a l k _ v e r t e x v t x = SDIS_RWALK_VERTEX_NULL ; s t r u c t s d i s _ i n t e r f a c e _ f r a g m e n t f r a g = SDIS_INTERFACE_FRAGMENT_NULL ; d o u b l e power ; d o u b l e f l u x ; d o u b l e e n d _ t e m p e r a t u r e ; s i z e _ t i , n ; r e s _ T r e s = RES_OK ; ASSERT ( g r e e n && i p a t h < d a r r a y _ g r e e n _ p a t h _ s i z e _ g e t (& g r e e n -> p a t h s ) && w e i g h t ) ; p a t h = d a r r a y _ g r e e n _ p a t h _ c d a t a _ g e t (& g r e e n -> p a t h s ) + i p a t h ; i f ( p a t h -> e n d _ t y p e == SDIS_GREEN_PATH_END_ERROR ) { / * R e j e c t e d p a t h * / r e s = RES_BAD_OP ; g o t o e r r o r ; } / * Compute medium power t e r m s * / power = 0 ; n = d a r r a y _ p o w e r _ t e r m _ s i z e _ g e t (& p a t h ->p o w e r _ t e r m s ) ; p o w e r _ t e r m s = d a r r a y _ p o w e r _ t e r m _ c d a t a _ g e t (& p a t h ->p o w e r _ t e r m s ) ; FOR_EACH( i , 0 , n ) { v t x . t i m e = INF ; medium = g r e e n _ f u n c t i o n _ f e t c h _ m e d i u m ( g r e e n , p o w e r _ t e r m s [ i ] . i d ) ; power += p o w e r _ t e r m s [ i ] . t e r m * s o l i d _ g e t _ v o l u m i c _ p o w e r ( medium , &v t x ) ; } / * Compute i n t e r f a c e f l u x e s * / f l u x = 0 ; n = d a r r a y _ f l u x _ t e r m _ s i z e _ g e t (& p a t h -> f l u x _ t e r m s ) ; f l u x _ t e r m s = d a r r a y _ f l u x _ t e r m _ c d a t a _ g e t (& p a t h -> f l u x _ t e r m s ) ; FOR_EACH( i , 0 , n ) { f r a g . t i m e = INF ; f r a g . s i d e = f l u x _ t e r m s [ i ] . s i d e ; i n t e r f = g r e e n _ f u n c t i o n _ f e t c h _ i n t e r f ( g r e e n , f l u x _ t e r m s [ i ] . i d ) ; f l u x += f l u x _ t e r m s [ i ] . t e r m * i n t e r f a c e _ s i d e _ g e t _ f l u x ( i n t e r f , &f r a g ) ; } / *

  s i z e _ t N = 0 ; / * # r e a l i s a t i o n s * / d o u b l e accum = 0 ; d o u b l e accum2 = 0 ; r e s _ T r e s = RES_OK ; i f ( ! g r e e n | | ! o u t _ e s t i m a t o r ) { r e s = RES_BAD_ARG ; g o t o e r r o r ; } n p a t h s = d a r r a y _ g r e e n _ p a t h _ s i z e _ g e t (& g r e e n -> p a t h s ) ; / * C r e a t e t h e e s t i m a t o r * / r e s = e n s u r e _ m e d i u m _ r e g i s t r a t i o n ( h a n d l e ->g r e e n , mdm ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; p a t h = h a n d l e -> p a t h ; t e r m s = d a r r a y _ p o w e r _ t e r m _ d a t a _ g e t (& p a t h ->p o w e r _ t e r m s ) ; n t e r m s = d a r r a y _ p o w e r _ t e r m _ s i z e _ g e t (& p a t h ->p o w e r _ t e r m s ) ; i d = m e d i u m _ g e t _ i d (mdm ) ; / * s t o r e t e r m * / t e r m . t e r m = v a l ; t e r m . i d = i d ; t e r m . v e r t e x = * v t x ; r e s = d a r r a y _ p o w e r _ t e r m _ p u s h _ b a c k (& h a n d l e ->p a t h ->p o w e r _ t e r m s , &t e r m ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; e x i t : r e t u r n r e s ; e r r o r : g o t o e x i t ; } r e s _ T g r e e n _ p a t h _ a d d _ f l u x _ t e r m ( s t r u c t g r e e n _ p a t h _ h a n d l e * h a n d l e , s t r u c t s d i s _ i n t e r f a c e * i n t e r f , c o n s t s t r u c t s d i s _ i n t e r f a c e _ f r a g m e n t * f r a g , c o n s t d o u b l e v a l ) { s t r u c t g r e e n _ p a t h * p a t h ; s t r u c t f l u x _ t e r m * t e r m s ;s t r u c t u n m e r g e d _ f l u x _ t e r m t e r m = FLUX_TERM_NULL__ ; s i z e _ t n t e r m s ; u n s i g n e d i d ; r e s _ T r e s = RES_OK ; ASSERT ( h a n d l e && i n t e r f && f r a g && v a l >= 0 ) ; r e s = e n s u r e _ i n t e r f a c e _ r e g i s t r a t i o n ( h a n d l e ->g r e e n , i n t e r f ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; p a t h = h a n d l e -> p a t h ; t e r m s = d a r r a y _ f l u x _ t e r m _ d a t a _ g e t (& p a t h -> f l u x _ t e r m s ) ; n t e r m s = d a r r a y _ f l u x _ t e r m _ s i z e _ g e t (& p a t h -> f l u x _ t e r m s ) ; i d = i n t e r f a c e _ g e t _ i d ( i n t e r f ) ; / * s t o r e t e r m * / t e r m . t e r m = v a l ; t e r m . i d = i d ; t e r m . f r a g m e n t = * f r a g ; r e s = d a r r a y _ f l u x _ t e r m _ p u s h _ b a c k (& h a n d l e ->p a t h -> f l u x _ t e r m s , &t e r m ) ; i f ( r e s != RES_OK ) g o t o e r r o r ; e x i t : r e t u r n r e s ; e r r o r : g o t o e x i t ; } Listing 8: green_function_solve_path for unmerged terms s t a t i c r e s _ T g r e e n _ f u n c t i o n _ s o l v e _ p a t h ( s t r u c t s d i s _ g r e e n _ f u n c t i o n * g r e e n , c o n s t s i z e _ t i p a t h , d o u b l e * w e i g h t ) { c o n s t s t r u c t p o w e r _ t e r m * p o w e r _ t e r m s = NULL ; c o n s t s t r u c t f l u x _ t e r m * f l u x _ t e r m s = NULL ; c o n s t s t r u c t g r e e n _ p a t h * p a t h = NULL ; c o n s t s t r u c t s d i s _ m e d i u m * medium = NULL ; c o n s t s t r u c t s d i s _ i n t e r f a c e * i n t e r f = NULL ; s t r u c t s d i s _ s c e n e * s c n = NULL ; d o u b l e power ; d o u b l e f l u x ; d o u b l e e n d _ t e m p e r a t u r e ; s i z e _ t i , n ; r e s _ T r e s = RES_OK ; ASSERT ( g r e e n && i p a t h < d a r r a y _ g r e e n _ p a t h _ s i z e _ g e t (& g r e e n -> p a t h s ) && w e i g h t ) ; p a t h = d a r r a y _ g r e e n _ p a t h _ c d a t a _ g e t (& g r e e n -> p a t h s ) + i p a t h ; i f ( p a t h -> e n d _ t y p e == SDIS_GREEN_PATH_END_ERROR ) { / * R e j e c t e d p a t h * / r e s = RES_BAD_OP ; g o t o e r r o r ; } / * Compute medium power t e r m s * / power = 0 ; n = d a r r a y _ p o w e r _ t e r m _ s i z e _ g e t (& p a t h ->p o w e r _ t e r m s ) ; p o w e r _ t e r m s = d a r r a y _ p o w e r _ t e r m _ c d a t a _ g e t (& p a t h ->p o w e r _ t e r m s ) ; FOR_EACH( i , 0 , n ) { medium = g r e e n _ f u n c t i o n _ f e t c h _ m e d i u m ( g r e e n , p o w e r _ t e r m s [ i ] . i d ) ; power += p o w e r _ t e r m s [ i ] . t e r m * s o l i d _ g e t _ v o l u m i c _ p o w e r ( medium , &p o w e r _ t e r m s [ i ] . v e r t e x ) ; } / * Compute i n t e r f a c e f l u x e s * / f l u x = 0 ; n = d a r r a y _ f l u x _ t e r m _ s i z e _ g e t (& p a t h -> f l u x _ t e r m s ) ; f l u x _ t e r m s = d a r r a y _ f l u x _ t e r m _ c d a t a _ g e t (& p a t h -> f l u x _ t e r m s ) ; FOR_EACH( i , 0 , n ) { i n t e r f = g r e e n _ f u n c t i o n _ f e t c h _ i n t e r f ( g r e e n , f l u x _ t e r m s [ i ] . i d ) ; f l u x += f l u x _ t e r m s [ i ] . t e r m * i n t e r f a c e _ s i d e _ g e t _ f l u x ( i n t e r f , &f l u x _ t e r m s [ i ] . f r a g m e n t ) ; } / * Compute p a t h ' s end t e m p e r a t u r e * / s w i t c h ( p a t h -> e n d _ t y p e ) { c a s e SDIS_GREEN_PATH_END_AT_INTERFACE : i n t e r f = g r e e n _ f u n c t i o n _ f e t c h _ i n t e r f ( g r e e n , p a t h -> l i m i t _ i d ) ; e n d _ t e m p e r a t u r e = i n t e r f a c e _ s i d e _ g e t _ t e m p e r a t u r e ( i n t e r f , &p a t h -> l i m i t . f r a g m e n t ) ; break ; c a s e SDIS_GREEN_PATH_END_IN_VOLUME : medium = g r e e n _ f u n c t i o n _ f e t c h _ m e d i u m ( g r e e n , p a t h -> l i m i t _ i d ) ; e n d _ t e m p e r a t u r e = m e d i u m _ g e t _ t e m p e r a t u r e ( medium , &p a t h -> l i m i t . v e r t e x ) ; break ; c a s e SDIS_GREEN_PATH_END_RADIATIVE : SDIS ( g r e e n _ f u n c t i o n _ g e t _ s c e n e ( g r e e n , &s c n ) ) ; SDIS ( s c e n e _ g e t _ a m b i e n t _ r a d i a t i v e _ t e m p e r a t u r e ( scn , &e n d _ t e m p e r a t u r e ) ) ; i f ( e n d _ t e m p e r a t u r e < 0 ) { / * Cannot h a v e i t n e g a t i v e i f u s e d * / r e s = RES_BAD_ARG ; g o t o e r r o r ; } break ; d e f a u l t : FATAL( " U n r e a c h a b l e � c o d e . \ n " ) ; break ; } / * Compute t h e p a t h w e i g h t * / * w e i g h t = power + f l u x + e n d _ t e m p e r a t u r e ;

Figure C. 14 :(Figure C. 15 :Figure C. 16 :Figure C. 17 :=Figure C. 18 :

 1415161718 Figure C.14: Open-porosity geometry without radiative transfer : a) Ambient fluid temperature b) Power density c) Flux density d) Solid boundary temperature e) Initial temperature. Volume and surface of the geometry are noted V and S and L = 4V/S (L = 1m) is retained as the characteristic size. The probe location � x c = (0.5, 0.5, 0.5) (at the center of the solid). The probe time estimation is t * = λt ρcL 2 =0.89 (t = 1 × 10 6 ). The fluid reference temperature T ref BF = 505K. The reference physical parameters are T ref I -T ref BF T ref BF = -0.01 (reference initial temperature T ref I = 500K), T ref B -T ref BF T ref BF

  3):

	<functions that store path ending> =
	<store path's end at an interface>
	<store path's end in a medium>
	<store radiative path's end>
	Let's start by describing the first function:
	<store path's end at an interface> =
	res_T
	green_path_set_limit_interface_fragment
	(struct green_path_handle* handle,
	struct sdis_interface* interf,
	const struct sdis_interface_fragment* frag,
	const double elapsed_time)
	{
	res_T res = RES_OK;
	<check input arguments>
	<register interface 'interf' against the green function>
	<store path duration>
	<store the location at interface>
	<store identifier of interface 'interf'>
	<store the path ends up at an interface>
	return RES_OK;
	}

  Path weight weight , computed from the source values, is then simply the sum of the different contributions:

	<fetch the end of path> =
	switch(path->end_type) {
	case SDIS_GREEN_PATH_END_AT_INTERFACE:
	interf = green_function_fetch_interf(green, path->limit_id);
	frag = path->limit.fragment;
	end_temperature = interface_side_get_temperature(interf, &frag);
	break;
	case SDIS_GREEN_PATH_END_IN_VOLUME:
	medium = green_function_fetch_medium(green, path->limit_id);
	vtx = path->limit.vertex;
	end_temperature = medium_get_temperature(medium, &vtx);
	break;
	case SDIS_GREEN_PATH_END_RADIATIVE:
	SDIS(green_function_get_scene(green, &scn));
	SDIS(scene_get_ambient_radiative_temperature(scn, &end_temperature));
	if(end_temperature < 0) { /* Cannot be negative if used */
	res = RES_BAD_ARG;
	goto error;
	}
	break;
	default: FATAL("Unreachable code.\n"); break;
	}

power , flux and end_temperature .

In such a formulation all source types are integrated over all times in [t I , t] and all locations in D. Boundary sources are therefore spatial Diracs at ∂D and the initial conditions are translated into temporal Diracs at t I .

https://gitlab.com/meso-star/stardis-solver

https://www.meso-star.com/projects/stardis/stardis.html#cli

https://www.edf.fr/en/the-edf-group/inventing-the-future-of-energy/r-d-global-expertise/our-offers/ simulation-softwares/syrthes

Fig.7was created on an online application[START_REF] Edstar | Thermal simulator training tool[END_REF] solving two-dimensional coupled thermal problems. The implementation relies on the same ray-tracing acceleration libraries used for the complex geometries presented in the present article. The use of this application is described in a Bachelor's level learning scenario[START_REF] Edstar | Introduction scenario for the thermal simulater[END_REF].

A similar illustration is provided at https://www.meso-star.com/projects/stardis/stardis.html using a porous medium of increasing refinement levels.

Acknowledgements

This project has received funding from the "Investissement d'avenir" program of the National Agency for Research of the French state under award number "ANR-10-LBX-22-01-SOLSTICE" and was supported by the ANR HIGH-TUNE, grant ANR-16-CE01-0010 and ANR MCG-RAD, grant ANR-18-CE46-0012 and from the Occitanie Region (Projet CLE-2016 ED-Star).

Stardis 0.CPC Library by Technical Editor) Developer's stardis.html Code Ocean capsule: (to be added by Technical Editor) Licensing provisions: GPLv3 Programming language: ANSI C and Python

link to program files: (to be added repository link: https://www.meso-star.com/projects/stardis/

Code structure

Code structure is briefly presented to help the reader understand the topics related to the Green function. All the data structures and functions described thereafter in a literate programming-inspired way [START_REF] Knuth | Literate Programming[END_REF], are located in the file sdis_green.c . The file is structured as follows:

previously accumulated by the function green_path_add_power_term is multiplied by the given volume power density value solid_get_volumic_power(medium, &vtx) . The total volume power contribution obtained for the path is power : The temperature at the end of the path, depending on the type of end, is the given interface temperature value interface_side_get_temperature(interf, &frag) , the given temperature value of the medium medium_get_temperature or the given ambient radiative temperature value sdis_scene_get_ambient_radiative_temperature(scn, &end_temperature) :