Brigitte Bidegaray-Fesquet
email: brigitte.bidegaray@imag.fr

Noise and offset in the IIR filtering of event-based sampled data

We investigate the impact of additive noise and uncertainty on levels (offset) on the filtering of level-crossing sampled data. We therefore suppose that either the signal of the levels are affected by an error following a normal distribution with zero mean. The errors are then analyzed in terms of the standard deviation of the normal distribution. For the illustrations an IIR filter, namely Butterworth, is used, but the results have a wider range of applications.

I. INTRODUCTION

For sporadic signals which are often encountered in mobile applications, a way to reduce the power consumption is to use event-driven sampling. The obtained data are then non-uniform [START_REF] Marvasti | Nonuniform Sampling. Theory and Practice, ser. Information Technology: Transmission, Processing and Storage[END_REF] and signal processing techniques have to be redefined. Analog-to-Digital Converters (ADCs) using a non regular samples can lead to interesting power savings compared to Nyquist ADCs [START_REF] Mark | A nonuniform sampling approach to data compression[END_REF], [START_REF] Sayiner | A level-crossing sampling scheme for A/D conversion[END_REF]. The sample number reduction that can be expected for signals with prescribed Hölder regularity has been investigated in [START_REF] Bidégaray-Fesquet | Data driven sampling of oscillating signals[END_REF] using a wavelet description of data.

Previous works have given equivalents of Finite Impulse Response (FIR) filters [START_REF] Aeschlimann | Asynchronous FIR filters: Towards a new digital processing chain[END_REF], [START_REF] Qaisar | Adaptive rate filtering for a signal driven sampling scheme[END_REF] or Infinite Impulse Response (IIR) filters [START_REF] Fesquet | IIR digital filtering of nonuniformly sampled signals via state representation[END_REF]. The effect of noise and offset has been of course intensively studied for classical regular signals. Here we propose a first analysis of the impact of these effects with filters specifically designed for non-uniform data.

We first describe precisely the type of data produced by the level-crossing sampling procedure and explain why we are specially interested in offset impact. Section II is devoted to the description of the implementation of filters for non-uniform data. Additive noise and offset are then addressed in Sections III and IV.

A. Level-crossing sampling

Our purpose is to write algorithms that can be used in asynchronous systems and for which the samples are obtained via an Asynchronous ADC (A-ADC) based on level-crossing sampling [START_REF] Akopyan | A level-crossing flash asynchronous analog-to-digital converter[END_REF], [START_REF] Allier | A new class of asynchronous A/D converters based on time quantization[END_REF]. This produces non-uniform samples in time with quantized amplitudes. Suppose we have 2 M predefined levels. We use here equi-spaced levels, but the algorithms do not take this feature into account. The design of specific sets of levels for specific applications and type of signals is not addressed here.

Each sample is the couple of an amplitude x j and a time delay δt j , computed with a local clock (see Figure 1). No global time is necessary for the implementation of the algorithms, but, for sake of clarity, we will describe them using sample times t j . We have δt j = t j -t j-1 . time

x j-1 x j • • • • • • • • • • • • • • • • • • • • • • • δt j
Fig. 1. Level-crossing sampling.

For data processing, the input analog signal only known at times t j and other information about the signal is lost. We however need an analog signal to derive the filtering procedures and consider an (at most) linearly interpolated signal from the non-uniform samples:

x(t) ∼ n j=1 (a j + b j t)χ [tj-1,tj ] , (1)
where χ I is the characteristic function of interval I.

B. Jitter vs. offset

There is clearly an amplitude-time duality between classical regular sampling and level-crossing sampling. In classical sampling, samples are supposed to be taken at regular time intervals. Times are supposed to be perfectly known and multiples of some basis sampling time t s . Amplitudes are captured at these times and have the precision of the ADC. On the contrary, in the level-crossing paradigm, levels and hence amplitudes are perfectly known, and times (or more precisely delays) are captured when a level is crossed with the precision of the local clock. Classical jitter consists in having an incertitude on the capture time. The amplitude is not captured at the right time, but data processing uses the theoretical regular structure, and computes as if amplitude x k is captured at time kt s , inducing a bias.

The equivalent phenomenon in level-crossing sampling is therefore an incertitude, an offset, on the amplitudes, i.e. the levels. For some reason the level which is recognized by the hardware can be slightly different from its theoretical design value. The time at which this level is crossed is slightly different from the time at which the theoretical level is crossed but the data processing is done with the theoretical amplitude. This induces a dual bias to that for classical jitter.

Except for may be one sample at the beginning or the end of the data, classical jitter does not induce a change in the number of samples. In level-crossing, an error on the level can induce changes in the number of samples if the value of the level is close to local extrema of the data (see Figure 2). In the situation shown in Figure 2 wrong values of levels can induce in particular in peak clipping.

II. NON-UNIFORM IIR FILTERS

In this paper, we will use Infinite Impulse Response (IIR) filters. The only reason is that we have to able to compare in the offset case two filtered signals at different input times, and that we have already proposed an algorithm that allows to compte the filtering result at any time, independently of the input samples times. In the non-uniform context, the ztransform formulation has no meaning. Two context have been studied and will be used here: the state representation of IIR filters [START_REF] Fesquet | IIR digital filtering of nonuniformly sampled signals via state representation[END_REF] or the formulation in the frequency domain [START_REF] Bidégaray-Fesquet | A fully non-uniform approach to FIR filtering[END_REF], [START_REF]Non-uniform filter design in the log-scale[END_REF].

A. Uniform filters for non-uniform data

The state representation consists in considering the usual filter coefficients α k and β k that occur in the Laplace expression of the filter transfer function

ĥ(p) = N k=0 α k p k / N k=0 β k p k . ( 2
)
Here N is the filter order. Then the filtering process consists in computing the state vector S(t), solution to the N-th order linear system

dS(t) dt = AS(t) + Bx(t), (3)
where

A =       0 1 • • • 0 0 . . . . . . . . . . . . . . . 0 0 • • • 1 0 0 0 • • • 0 1 -β 0 -β 1 • • • -β K-2 -β K-1      
is the N ×N state matrix, and B = (0 • • • 0 1) t is the command vector, and then the filtered signal y(t) is simply computed as

y(t) = CS(t) + Dx(t), (4)
where

C = (α 0 -α N β 0 • • • α N -1 -α N β N -1
) is the observation vector and D = α N the direct link coefficient.

Numerical algorithms follow from then the discretization of the differential equation (3) or its integral form

S(t) = e At S(0) + t t0 e A(t-τ ) Bx(τ ) dτ (5)
(see [START_REF] Fesquet | IIR digital filtering of nonuniformly sampled signals via state representation[END_REF] for more details). The non-uniformity of data is handled by the fact that the discretization uses irregular timesteps stemming from the input data times.

B. Non-uniform filters in the frequency domain

We have also developed filters which are based on the interpolation in the frequency domain of the analog filter transfer function H. The filtered signal is then the convolution of the input signal x(t)

y(t) = +∞ -∞ h(t -τ )x(τ ) dτ (6)
with the inverse Fourier transform of H:

h(t) = 1 2π +∞ -∞ H(ω)e iωt dω. (7)
Then aim of this approach is that it is possible to describe with very few filter samples equivalents of very high order filters, and even the ideal low-pass filter

H(ω) = χ [-ωc,ωc] [10].
Here we sample in the frequency domain existing filter. The sampling can be done non-uniformly, i.e. with non equispaced frequencies. In practice we define the filter for positive values of the frequency and suppose it is symmetric. Levelcrossing sampling is used to define the sampling frequencies and we obtain filter samples (ω k , H k ). We describe here linear interpolation, but log-scale interpolation has also been studied [START_REF]Non-uniform filter design in the log-scale[END_REF]. The filter amplitudes H k are complex and it has been showed in [START_REF]Non-uniform filter interpolation in the frequency domain[END_REF] that the best way to perform the interpolation is to treat the absolute values and the phase of the amplitudes separately. This leads to less attenuation in the desired pass domain. For a low-pass filter linear interpolation yields a continuous approximation of H(ω):

H(ω) ∼ K k=1 (ρ 0 k + ρ 1 k ω)e i(θ 0 k +θ 1 k ω) χ J k + (ρ 0 k -ρ 1 k ω) -e i(θ 0 k -θ 1 k ω) χ J - k . (8)
Here

J k = [ω k-1 , ω k ] (ω 0 = 0), J -k = [-ω k , -ω k-1 ]
, and ρ 0 k , ρ 1 k , θ 0 k , and θ 1 k are computed in a straightforward way from the filter samples (ω k , H k ).

Inserting (1) and ( 7) in ( 6) the filtered signal can be cast as

y(t) = n j=0 x j K k=1 h jk (t), (9)
which more or less looks like a FIR filter formula, although we will interpolate IIR filter transfer functions. Details can be found in [START_REF]Non-uniform filter interpolation in the frequency domain[END_REF] and [START_REF] Bidégaray-Fesquet | A fully non-uniform approach to FIR filtering[END_REF]. One advantage of this approach its that we will be able to evaluated y at any time t, possibly different from times t j .

III. ADDITIVE NOISE

The goal is to compare at each data time (not changed here by the addition of noise) filtered signal y(t) obtained from the discretization of x(t) and the filtered signal ỹ(t) associated to x(t) = x(t) + ε(t) via the computation of the mean square error, weighted bu the length of each sample:

E = 1 T j (y(t j ) -ỹ(t j )) 2 δt j , ( 10
)
where T is the total time duration of the signal.

An additive noise ε(t) is added, with zero mean and a standard deviation σ. The ε(t j ) are independent identically distributed. It is therefore expected that the mean square error will be of the order of σ 2 .

The SPASS MATLAB toolbox is used for numerical tests [START_REF]SPASS 2.0: Signal Processing for ASynchronous Systems[END_REF]. It implements a series of routines specific for the numerical processing of non-uniform data. In particular it implements the A-ADC and filters described above.

For the numerical tests we use

x(t) = .5 + .25 cos(2πt) + .25 cos(4 • 2πt) (11)
as input signal. The regular samples are taken at a 1000 Hz rate. The level-crossing sampling of x(t) is performed defining 2 M (M = 3 or 4) equi-spaced levels ranging from .05 to .95. There are 10.000 regular samples, 320 non-uniform samples for M = 3 and 640 non-uniform samples for M = 4.

Filtering is performed with a 5th order low-pass Butterworth filter, with cut-off frequency 3•2π. The resulting filtered signal y(t) is very similar for the classical Butterworth filter and for both implementations of non-uniform filtering, state representation and frequency-domain interpolation.

The standard deviation σ of the additive noise ε(t) is chosen to be comparable with the level quantum q. More precisely we test σ = q/4, q/2, q, and 2q. Note that q depends on M . The filtered signal ỹ is computed for the regular filter, and via state representation for M = 3 and 4. The mean square error is displayed in Figure 3 in log-log scale. The experiment is performed 10 times and a mean error is computed.

The linear behavior in log-log scale is well demonstrated by the test case. The slope of all the curves is close to 2, confirming that the mean square error is proportional to σ 2 . The proportion constant is dependent on the sampling and is of course much lower for the regular case which is less affected by interpolation errors in the computation of the mean square error.

IV. OFFSET

We suppose that there is an offset that is uniform over time, which means that for each realization of the normal distribution used, the set of levels is replaced by a perturbed one. This behavior can be explained for example by an experimental setting where the comparators of the A-ADC have slightly different properties as prescribed by the design. As for the additive noise case we are interested in the mean square error between the filtered signal with an without offset. More realizations of the normal distribution will be necessary to obtain convergence, since only 2 M random numbers are computed for each realization. Another feature is that the number of non-uniform samples vary from one realization to another. We therefore have to compare the signals at the same times, and use the frequency-domain algorithm for this purpose. Beside we can study the dispersion of the number of samples and compute the variance of this number.

We use the same input signal as in Section III, with M = 3 and σ = q/16, q/8, q/4, q/2, and q. The results are displayed in Figure 4. Here no comparison with the regular case is possible. As expected, if the offset standard deviation is of the order of the quantum, the signal and hence its filtering is too much perturbed. For smaller values of the behavior is almost linear in log-log scale (more realization should be used to have better results) and the slope is 1.8, close to 2. Now the variance of the number of samples is computed and displayed in Figure 5. Here again a higher number of realizations should be done to have better results, but the slope is 1. The variance of the number of samples behaves proportionally to the standard deviation.

V. CONCLUSION

We have investigated the impact of additive noise and uncertainty on levels (offset) on the IRR filtering of levelcrossing sampled data. The number of levels that we use is very low and is quite sufficient to obtain filtered signals which are precise enough for the targeted applications. Besides the test cases used here are practical for filter testing since it is easy to check whether the higher frequency (4 in our case) has been filtered out or not. Of course level-crossing sampling is especially efficient in the case of sporadic signal (typically not sum of sine functions!) because of the very low number of non-linear samples compared to Nyqvist sampling. In both cases additive noise and offset, we have been able to connect the standard deviation of the applied noise (on the signal or the levels) to the mean square error between the perturbed and non perturbed signals.

Further results have to be obtained in the case when both the signal and the levels are affected by noise.

Fig. 2 .

 2 Fig. 2. Variation of the number of samples.

Fig. 3 .

 3 Fig. 3. Error in presence of an additive noise; asterisk: regular case, diamond: M = 3, square M = 4.

Fig. 4 .

 4 Fig. 4. Error in presence of offset.

Fig. 5 .

 5 Fig. 5. Variance of the number of samples in presence of offset.

ACKNOWLEDGMENT

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).


Delphine Ramalingom
email: delphine.ramalingom@univ-reunion.fr

Pierre-Henri Cocquet

Alain Bastide

A new interpolation technique to deal with fluid-porous media interfaces for topology optimization of heat transfer

Keywords: Topology optimization, Heat transfer, Interpolation function, Sigmoid function, Continuous adjoint method, Finite volume method

This paper proposes a new interpolation technique based on density approach to solve topology optimization problems for heat transfer. Problems are modeled under the assumptions of steady-state laminar flow using the incompressible Navier-Stokes equations coupled to the convection-diffusion equation through the Boussinesq approximation. The governing equations are discretized using finite volume elements and topology optimization is performed using adjoint sensitivity analysis. Material distribution and effective conductivity are interpolated by two sigmoid functions respectively h τ (α) and k τ (α) in order to provide a continuous transition between the solid and the fluid domains. Comparison with standard interpolation function of the literature (RAMP function) shows a smaller transition zone between the fluid and the solid thereby, avoiding some regularization techniques. In order to validate the new method, numerical applications are investigated on some cases from the literature, namely the single pipe and the bend pipe. Lastly, as

Introduction

Since its introduction by Bendsoe and Sigmund [START_REF] Bendsφe | Topology optimization: theory, methods, and applications[END_REF] for solid mechanics problems, topology optimization has become a powerful and increasingly popular tool for designers and engineers for design process. Topology optimization is a material distribution method for finding the optimal structure, for a given problem subject to design constraints. Contrarily to shape optimization where the topology (i.e. the number of boundaries and connectivity) is predetermined, topology optimization allows introduction of new boundaries during the design process.

Topology optimization was pioneered for Stokes flow by Borrvall and

Petersson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF]. They introduced a friction term yielding the generalized Stokes equations. Gersborg-Hansen [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF] and Olesen et al. [START_REF] Olesen | A highlevel programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow[END_REF] extended topology optimization for fluid flow problems to the Navier-Stokes equations.

In topology optimization, the material distribution is parametrized by defining a design variable α ∈ {0; 1}. This variable is discrete and should either represent solid material (α = 1) or fluid (α = 0). A common approach to solve the topology optimization problem with this discrete value as optimization parameter, is to change it into a continuous one by introducing a porous media with a continuous permeability variable for each element. This method, known as the Brinkman penalization, leads to a problem where flow and (almost) non-flow regions are developed by allowing interpolation between the lower and upper value of permeability. Generally, authors used the density interpolation function proposed by Borrvall and Petersson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] or a reformulated version of their convex and q-parametrized interpolation function. The parameter q > 0 is a penalty parameter that is used to control the level of 'gray' in the optimal design. However, authors had also experienced problems with locally optimal solutions. Therefore, they considered a two-steps solution procedure where the problem was first solved with a small penalty value of q = 0.01 for example and then the result is used as initial case for the problem with a penalty value of q = 0.1 [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF][START_REF] Olesen | A highlevel programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow[END_REF] or q = 1 [START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF]. The mathematical foundation of the interpolation of α was further investigated by Evgrafov [START_REF] Evgrafov | The limits of porous materials in the topology optimization of Stokes flows[END_REF] where the limiting cases of pure fluid and solid were included. Brinkman approach has since been used for several problems as transport problem [START_REF] Andreasen | Topology optimization of microfluidic mixers[END_REF], reactive [START_REF] Okkels | Scaling behavior of optimally structured catalytic microfluidic reactors[END_REF] and transient flows [START_REF] Deng | Topology optimization of unsteady incompressible NavierStokes flows[END_REF][START_REF] Kreissl | Topology optimization for unsteady flow[END_REF], fluid-structure interaction [START_REF] Yoon | Topology optimization for stationary fluidstructure interaction problems using a new monolithic formulation[END_REF] and also flows driven by body forces [START_REF] Deng | Topology optimization of steady and unsteady incompressible NavierStokes flows driven by body forces[END_REF].

A variation of the approach is presented by Guest and Prevost [START_REF] Guest | Topology optimization of creep-ing fluid flows using a DarcyStokes finite element[END_REF]. They proposed to regularize the solid-fluid structure by treating the material phase as a porous medium where fluid flow is governed by Darcy's law. In their approach, flows through voids are governed by Stokes flow and, when the solid phase is impermeable, discrete no-slip condition is simulated by assigning a low permeability to the solid phase. There exists other alternatives to Brinkman penalization in the literature. The level set approach to topology optimization has been applied to fluid flows problems [START_REF] Challis | Level set topology optimization of fluids in Stokes flow[END_REF][START_REF] Zhou | A variational level set method for the topology optimization of steady-state NavierStokes flow[END_REF][START_REF] Kubo | A level set-based topology optimization method for optimal manifold designs with flow uniformity in plate-type microchannel reactors[END_REF], and recently the level set approach was combined with the extended finite element method (XFEM) by Kreissl and Maute [START_REF] Kreissl | Levelset based fluid topology optimization using the extended finite element method[END_REF] and by Jenkins et al. [START_REF] Jenkins | An immersed boundary approach for shape and topology optimization of stationary fluid-structure interaction problems[END_REF].

The main drawback of the level set appraoch is the constraint of remeshing throughout the optimization process.

The second difficulty in topology optimization is to deal with the difference in thermal conductivity in the solid and fluid domains. Most publications interpolate the conductivity using the SIMP method (Solid Isotropic Material with penalization). This method allows to deal with the discrete nature of conductivity material distribution. So, authors [START_REF] Rozvany | Generalized shape optimization without homogenization[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF][START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Yoon | Topology optimization for stationary fluidstructure interaction problems using a new monolithic formulation[END_REF] considered a continuous local thermal conductivity controlled by the design parameter α ranging from 0 to 1. Thanks to this function, the optimization algorithm is able to reallocate thermal conductivity material, or creating holes in its structure to reach the objective function. Moreover, the convex and p-parametrized function interpolation is similar to the density interpolation function of Borrvall and Petersson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] or the RAMP (Rational Approximation of Material Properties)-style function as introduced by Stolpe [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF][START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Lee | Topology Optimization of Convective Cooling System Designs[END_REF] .

Other methods have also been investigated. Matsumori et al. [START_REF] Matsumori | Topology optimization for fluidthermal interaction problems[END_REF] presented results with a linear-interpolated design-dependant volumetric heat generation. Dede [START_REF] Dede | Multiphysics optimization, synthesis, and application of jet impingement target surfaces[END_REF] used a linear interpolation for thermal conductivity. Thus, the main issue is to deal with intermediate design variables and nonphysical flow solutions.

There are three main categories of algorithms to solve topology optimization problems : gradient-free, gradient-based and hessian-based algorithms.

In topology optimization problems with large number of design variables, gradient-based algorithms are used to find accurate solutions efficiently. One of the advantages of the interpolation functions described above is the possibility of using gradient-based continuous optimization methods. These methods are based on derivatives in order to find extrema, and are the so-called sensitivity analysis. It aims at evaluating the derivative of objective function with respect of α. Gradient-based algorithm is widely used by several authors [START_REF] Othmer | Adjoint methods for car aerodynamics[END_REF][START_REF] Dede | Multiphysics topology optimization of heat transfer and fluid flow systems[END_REF][START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Yoon | Topology optimization for stationary fluidstructure interaction problems using a new monolithic formulation[END_REF][START_REF] Nrgaard | Topology optimization of unsteady flow problems using the lattice Boltzmann method[END_REF]. Moreover, since most of the topology optimization problems involve a huge number of design variable, specific gradient-based optimization algorithms must be chosen to handle this difficulty. A famous algorithm from the literature is the MMA (Method of Moving Asymptotes) developed initially by Svanberg [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF]. In order to reduce the computational costs, adjoint approach consisting to calculate the sensitivities of the objective function by an adjoint state has been adopted. Other methods have been explored to reduce the computation cost: the multigrid preconditioned conjugate gradients (MGCC) by Amir et al. [START_REF] Amir | On reducing computational effort in topology optimization: how far can we go?[END_REF], multi-resolution multi-scale topology optimization technique by Kim et al. [START_REF] Kim | Multi-resolution multi-scale topology optimizationa new paradigm[END_REF], the technique of using adaptive design variable fields by Guest et al. [START_REF] Guest | Reducing dimensionality in topology optimization using adaptive design variable fields[END_REF].

Moreover, various regularization techniques based on filtering of either the design variable α or the sensitivity ∂f ∂α [START_REF] Bendsφe | Topology optimization: theory, methods, and applications[END_REF][START_REF] Gersborg | An explicit parameterization for casting constraints in gradient driven topology optimization[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF][START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Koga | Development of heat sink device by using topology optimization[END_REF] exist to ensure well-posed topology optimization problems. The regularization works by defining a certain length scale r 0 below which any features in α or ∂f ∂α are smeared out by the filter; that results in optimized structures with a minimal feature size r 0 independent of the mesh refinement. As mentioned by some authors [START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF][START_REF] Alexandersen | Topology optimisation for coupled convection problems[END_REF], these regularization techniques allow to avoid checkerboard problems.

This paper proposes a new interpolation technique in order to solve a heat transfer topology optimization problem. Design material and effective conductivity are interpolated respectively by a function h τ (α) and another function k τ (α) in order to provide a continuous transition between the solid and the fluid domains. These interpolation functions avoid the use of some regularization techniques because the problem can be solved in one-step without a new value of the convexity parameter. Moreover, these interpolation functions allow a smaller transition zone between the fluid and solid regions.

To prove this claim and get more qualitative results, the size of these transition zones is explicitly computed and comparison with standard RAMP interpolation is lead. In order to validate the new method, some numerical applications are investigated on the single pipe and the bend pipe cases of the literature. Lastly, as two new parameters are introduced thanks to the interpolation functions, we study their impact on the results of a optimization problem.

Governing equation

The main goal of this paper is to solve topology optimization problems for heat transfer in fluid flow. The latter can be written in the general form below:

Minimize J (u, p, θ) = Ω J Ω (u, p, θ)dΩ + Γ J Γ (u, p, θ)dΓ
Subject to Governing equation for (u, p, θ)

Boundary conditions on Γ , (1)
where u, p and θ are respectively the dimensionless velocity, pressure and temperature, Ω is a bounded open set of R d , d = 2, 3, with boundary Γ = ∂Ω. The function J Ω and J Γ are some cost functional modeling a physical effect one wish to minimize. We assume that Γ can be decomposed as Γ = Γ 1 ∪ Γ 2 ∪ Γ in ∪ Γ out where Γ in is the inlet, Γ out the outlet and Γ 1 , Γ 2 are going to be considered as walls.

For the governing equation, the flows considered in this study are assumed Newtonian and incompressible, steady and laminar. The inverse permeability field is introduced in the steady-state Navier-Stokes equations as a source term h τ (α)u yielding a Brinkman model with a convection term. The set of dimensionless equations governing the conservation of momentum, mass and energy for incompressible steady-state fluid flow are the following:

(u • ∇) u = -∇p + Re -1 ∆u -h τ (α) u + Ri θ - → e y in Ω, ∇ • u = 0 in Ω, ∇ • (u θ) = div (Re -1 P r -1 k τ (α) ∇θ) in Ω. (2)
The reduced dimensionless temperature is θ = (T -T 0 )/∆T and the uniform heat flux Φ = ∂ n θ on some part of Γ 1 is equal to λ f ∆T /l, with λ f the thermal conductivity of the fluid. Parameters governing the flow are the Reynolds number defined as Re = U 0 l/ν, with U the reference velocity (i.e mean velocity here) and l is the hydraulic diameter, the Richardson number Ri = Gr/Re 2 where Gr is the Grashof number defined as Gr = g β ∆T l 3 /ν 2 .

h τ (α) corresponds to the ratio between a kinematic viscosity and a permeability and α is the spatially varying design variable field determined by the optimization algorithm. Regions with low permeability can be considered as solid regions since (at least formally) the velocity of the fluid vanishes in such region, and those with very high permeability regions are interpreted as pure fluid. The interpolation function for the adimensional thermal conductivity is given by k τ (α).

We consider the following set of boundary conditions:

u = 0, ∂ n θ = -1 on Γ 1 , u = 0, ∂ n θ = 0 on Γ 2 , u n = 1, u t = 0 θ = 0 on Γ in , ∂ n u = 0, ∂ n θ = 0, p = 0 on Γ out , (3)
where u n and u t are the normal and tangential components of primal velocities, respectively. To summarize, we require a constant horizontal velocity and a constant temperature θ 0 at the inlet, vanishing gradient for both velocity and temperature of fluid at the outlet. Homogeneous Dirichlet for the velocity and Neumann boundary condition are prescribed for the temperature on the walls. It is worth noting that boundary conditions on the outlet have been considered in [START_REF] Ramalingom | Numerical study of natural convection in asymmetrically heated channel considering thermal stratification and surface radiation[END_REF] and ensure that the fluid does not re-enter in the domain.

Interpolation

The goal of topology optimization is to end up with binary designs, i.e avoid that the design variables take other value than those representing the fluid or the solid. The most important thing is to be sure that the intermediate regions (the transition zones) are unattractive with respect to the optimization problem. This is usually carried out by penalizing the intermediate densities with respect to the material parameters, such as inverse permeability and effective conductivity. A standard approach is to use some convex interpolation (RAMP) function [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF][START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF]. In this section, we propose to use another interpolation function hence interpolate the inverse permeability and the effective conductivity using sigmoid functions. After giving the definition of our interpolation function, we show that the transition zones of the sigmoid are much smaller than those of the standard RAMP function therefore motivating the use of such method in topology optimization.

Interpolation of porosity using a sigmoid function

Inverse permeability is interpolated with a sigmoid function

h τ (α, α 0 ) = α max   1 1 + exp (-τ (α -α 0 )) - 1 1 + exp (τ α 0 )   (4)
where α ∈ [0, α max ] with α max being the maximal value h τ can reach. Direct computations show that h τ (0, α 0 ) = 0 and that the following point-wise convergence holds

lim τ →+∞ h τ (α, α 0 ) =          0 if α < α 0 , α 0 /2 if α = α 0 , α max if α > α 0 . (5)
This shows that h τ is a smooth regularization of a Heaviside step function.

From [START_REF] Kreissl | Levelset based fluid topology optimization using the extended finite element method[END_REF], one can see that α 0 can be tuned to control the size of the fluid part in the computational domain. Finally, note that the definition of our interpolation function has to be changed when α 0 = 0 since we require α ≥ 0 and (4) would lead to a regularization of a step function satisfying h τ (0, 0) = 0 and, for any α > 0, h τ (α, 0) → α max /2 as τ → +∞. In order to cover this case, one can use the following interpolation function

h τ (α, 0) = h τ (α) = 2 α max   1 1 + exp (-τ α) - 1 2   ,
which satisfies h τ (0) = 0 and, if α > 0, h τ (α) → α max as τ → +∞.

Comparison with standard RAMP interpolation function

The RAMP function has been introduced in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] and is defined as follows:

h q,R (s) = α max + (α min -α max )(1 -s) 1 + q 1 -s + q , (6)
where s ∈ [0, 1] is the design parameter and α min , α max are respectively the minimal and maximal value the RAMP function can take. Again, some direct computations show that h q,R is a smooth and convex regularization of a Heaviside step function since h q,R (0) = α min and, for any s > 0, lim q→0 h q,R (s) = α max . In the sequel, we use α min = 0 since this does not add any difficulties nor plus to the discussion. We now compare the performances of our sigmoid (4) with the standard RAMP function ( 6) by focusing on the size of their so-called transition zones. First, note that they are both defined on Ω but the transition zone only depends on the value of the design parameter s(x) for x ∈ Ω. As a result, we are going to consider that both functions have real domain and act on the interval [0, 1]. Also, to compare functions having the same domain and being both smooth regularization of step function as q → 0, we rescale the sigmoids (4) thanks to the formula

h q,S (s, α 0 ) = h 1/q (α max s, α 0 ), h q,S,0 (s) = h 1/q (sα max ), s ∈ [0, 1].
By definition, the transition zones are the undesirable values of the design variable s ∈ [0, 1] for which the velocity of the speed is not enough penalized to vanish or the value of the interpolation function is not small enough for considering these zones as fluid. As one can see from Figure 1, the sigmoid function has much smaller transition zone than the RAMP function.

To prove this claim and get more qualitative results we are going to compute explicitly the size of these transition zones. For a given small enough

α max = 2000.
ε > 0 and a large enough M > 0, the latter are defined by

T (f, ε, M ) = {s ∈ [0, 1] | ε ≤ f (s) ≤ M } , f ∈ {h q,R , h q,S , h q,S,0 } .
Since both interpolation function are increasing, they admit an inverse function f -1 and the transition zone is thus given by the interval

T (f, ε, M ) = [f -1 (ε), f -1 (M )
] from which we infer that

|T (f, ε, M ) | = f -1 (M ) -f -1 (ε). ( 7
)
It then only remains to compute the inverse of the sigmoid and RAMP functions. This is actually achieved by solving f (s) = y for a given y ∈ [0, α max ].

These equations can be solved analytically using only direct computations and give

s = (h q,R ) -1 (y) = 1 -q y -α max α max -y + (1 + q)(α min -α max ) = y 1 + q α max q + y (8) s = (h q,S ) -1 (y) = α 0 α max - q α max        -log α max -y(1 + exp(-α 0 /q)) α max exp(-α 0 /q) + y(1 + exp(-α 0 /q)) if s ≤ α 0 /α max log α max -y(1 + exp(-α 0 /q)) α max exp(-α 0 /q) + y(1 + exp(-α 0 /q)) if s ≥ α 0 /α max s = (h q,S,0 ) -1 (y) = -log α max -y α max + y q α max .
The constraint on s in the definition of the inverse of h q,S comes from the fact that we have to solve an equation of the form exp((α 0 -sα max )/q) = g(y)

for a positive function g. Using ( 7) and ( 8), the size of the transition zones are finally

|T (h q,R , ε, M ) | := T R = qα max (1 + q) (α max q + M )(α max q + ε) (M -ε), |T (h q,S , ε, M ) | := T S = q α max log α max -M (1 + exp(-α 0 /q)) α max exp(-α 0 /q) + M (1 + exp(-α 0 /q)) + q α max log α max -ε(1 + exp(-α 0 /q)) α max exp(-α 0 /q) + ε(1 + exp(-α 0 /q)) , |T (h q,S,0 , ε, M ) | := T S,0 = q α max log (α max -M )(α max + ε) (α max + M )(α max -ε) .
For α 0 = 0, one can clearly see that T S,0 < T R . To deal with the case α 0 > 0, note that the parameter q is small since we wish the interpolation function to be close to an ideal step function. Also, we emphasize that α 0 > 0 is going to depend on α max in order to control the percentage of fluid in Ω. Since α max has to be large enough, one gets that z = exp(-α 0 /q) is a small parameter and we can thus expand T S as z → 0. This gives

T S = q α max log (α max -M )(α max -ε) M ε + O(exp(-α 0 /q)),
from which one can see that T S < T R . Thanks to the factor q/α max , we can see that the intermediate zones of the sigmoid are in fact much smaller than those of the RAMP function. Some numerical calculation of the size of transition T R T S T S,0 q = 1 0.6657 0.0380 0.0055 q = 10 -4 0.1665 3.8002 × 10 -6 5.4881 × 10 -7 q = 10 -5 0.0196 3.8002 × 10 -7 5.4881 × 10 -8 zone confirming this fact can be found in Table 1,2 for α max ∈ {200, 2000} and q ∈ {1, 10 -4 , 10 -5 }.

Interpolation of the thermal conductivity

The effective conductivity is interpolated using a sigmoid function similar to (4):

k τ (α, α 0 ) = (k s -k f )   1 1 + exp (-τ (α -α 0 )) - 1 1 + exp (τ α 0 )   + k f (9)
where α ∈ [0, α max ] is the design parameter, k s and k f are respectively the adimensional thermal diffusivity of the solid and of the fluid. Interpolation function ( 9) is again a smooth regularization of a Heaviside step function that satisfies k τ (0, α 0 ) = k f together with the point-wise convergence

lim τ →+∞ k τ (α, α 0 ) =          k f if α < α 0 , (k f + k s )/2 if α = α 0 , k s if α > α 0 .
Similarly to h τ , the case α 0 = 0 needs the following slight modification to still have, when τ → +∞, a fluid-solid step function

k τ (α, 0) = k τ (α) = 2(k s -k f )   1 1 + exp (-τ α) - 1 2   + k f .
The previous function then satisfies k τ (0) = 0 and, for any α > 0, k τ (α) → k s as τ → +∞. A RAMP-like interpolation function is also usually used for the thermal conductivity [START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF]. Computations similar to those of Section 3.2 can be done to prove that the transition zone of the sigmoid (9) are much smaller than those of the RAMP function hence motivating using such technique in topology optimization.

Gradient computation with the adjoint model

We chose to solve the optimization problem (1) with gradient-based optimization algorithm. The latter requires the computation of the gradient of the cost functional with respect to the design variable which needs the sensitivities of the solution to the state equation with respect to α. To do this, we are going to use the adjoint method (see e.g. [START_REF] Herzog | Algorithms for PDEconstrained optimization[END_REF][START_REF] Petra | Model variational inverse problems governed by partial differential equations[END_REF][START_REF] Hinze | Optimization with PDE constraints[END_REF]) since it allows to compute the gradient of J efficiently, only by solving the so-called direct problem (the governing equation) and the adjoint model.

We adopt here the differentiate then discretize approach [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF] which means that we need to compute the continuous adjoint model associated to (1) as done in [START_REF] Qian | Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint[END_REF][START_REF] Othmer | A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows[END_REF]. The Lagrangian associated to the optimization problem [START_REF] Bendsφe | Topology optimization: theory, methods, and applications[END_REF] can then be defined as follows L((u, θ, p), (u * , θ * , p * ), α) := J (u, θ, p) -

Ω p * div u dΩ - Ω u * • (u • ∇) u + ∇p -Re -1 ∆u + h τ (α, α 0 ) u -Ri θ - → e y dΩ - Ω θ * • ∇ • (u θ) -div(Re -1 P r -1 k τ (α, α 0 )∇θ) dΩ - Γ 1 (u • Φ 1 + (∂ n θ -1)ψ 1 ) dΓ - Γ 2 (u • Φ 2 + ∂ n θψ 2 ) dΓ - Γ in ((u -u in ) • Φ in + θψ in ) dΓ - Γout (∂ n u • Φ out + pq out + ∂ n θψ out ) dΓ (10)
where (u * , θ * , p * , Φ l , ψ l , q out ), for l ∈ {1, 2, in, out}, can be seen as Lagrange multiplier (or adjoint variables). We emphasize that the variables (Φ t , ψ t , q out ) are here to enforce the boundary conditions (3).

The adjoint model is defined thanks to the critical point of L with respect to the state variables (u, θ, p). We use below the following notation for the derivative of an application F : x ∈ E → F(x) ∈ F, where E, F are normed spaces

∂F ∂x [δx] := lim ε→0 F(x + εδx) -F(x) ε . (11)
First, since L is linear with respect to the adjoint variables, it is worth noting that we recover the state equation ( 2)-( 3) if, for l ∈ {1, 2, in, out}, we solve ∂L ∂(u * , θ * , p * , Φ l , ψ l , q out ) [δu * , δθ * , δp * , δΦ l , δψ l , δq out ] = 0.

Using now [START_REF] Kangro | Divergence boundary conditions for vector Helmholtz equations with divergence constraints[END_REF] to differentiate the Lagrangian [START_REF] Gersborg | An explicit parameterization for casting constraints in gradient driven topology optimization[END_REF] with respect to the state variables, and integrating by parts to have no terms involving derivative of (δu, δθ, δp), we end up with

∂L ∂(u, θ, p) [δu, δθ, δp] = Ω ∂J Ω ∂(u, θ, p) [δu, δθ, δp] dΩ + Γ ∂J Γ ∂(u, θ, p) [δu, δθ, δp] dΓ+ Ω δu • ∇p * -h τ (α, α 0 )u * -(∇u) T u * + ∇u * u + θ∇θ * + Re -1 ∆u * dΩ + Ω δθ div(Re -1 P r -1 k τ (α, α 0 )∇θ * ) + u • ∇θ * + Ri u * • e y dΩ + Ω δp div u * dΩ - Γ δp (u * • n) dΓ + Γ δu • -n p * -n θ θ * + Re -1 ∂ n u * -(u • n)u * dΓ - Γ Re -1 ∂ n δu • u * dΓ + Γ δθ (-u • n) θ * -Re -1 P r -1 k τ (α, α 0 ) ∂ n θ * dΓ + Γ θ * Re -1 P r -1 k τ (α, α 0 )∂ n δθ dΓ - Γ1 (δu • Φ 1 + ∂ n δθψ 1 ) dΓ - Γ2 (δu • Φ 2 + ∂ n δθψ 2 ) dΓ - Γ in (δu • Φ in + δθψ in ) dΓ - Γout (∂ n δu • Φ out + δpq out + ∂ n δθψ out ) dΓ (12)
Assuming that (δu, δθ, δp, Φ l , ψ l , q out ) = 0 for l ∈ {1, 2, in, out}, that the first derivatives of δu, δθ vanish on Γ and solving (∂L/∂(u, θ, p)) [δu, δθ, δp] = 0 yield the adjoint problem:

∇p * -h τ (α) u * + θ ∇θ * + Re -1 ∆u * + ∇u * u -(∇u) T u * = -∂J Ω ∂u in Ω div u * = -∂J Ω ∂p in Ω Ri u * • - → e y + u • ∇θ * + div (Re -1 P r -1 k τ (α)∇θ * ) = -∂J Ω ∂θ in Ω. ( 13
)
For the adjoint boundary condition, we are going to show how to obtain them on Γ 1 ∪ Γ 2 since the other part of Γ can be done in the same spirit.

On both walls, the primal velocity is fixed to zero (no-slip condition on the walls) and the heat flux either vanishes or is constant. Since div u = 0 and u is prescribed or null value, there will be δu = 0, div δu = 0 and ∂ n δθ = 0.

We recall a formula from [START_REF] Kangro | Divergence boundary conditions for vector Helmholtz equations with divergence constraints[END_REF]Lemma 7] that holds for any vector fields w and reads

∂ n w • n = div w| Γ -div Γ w t -κw • n, (14)
where κ = div n is the curvature of Γ and div Γ is the surface divergence

operator. Since ∂ n δu • u * = (∂ n δu) n • u * n + (∂ n δu) t • u *
t and δu = 0 and div δu = 0, formula [START_REF] Evgrafov | The limits of porous materials in the topology optimization of Stokes flows[END_REF] gives that

∂ n δu • u * = (∂ n δu) t u * t .
The boundary conditions on Γ 1 ∪ Γ 2 are then obtained from ( 12) by considering, for l = 1, 2, the following vanishing terms

Γ l Re -1 ∂ n δu • u * dΓ = Γ Re -1 (∂ n δu) t • u * t dΓ = 0, Γ l δp u * • n - ∂J Γ ∂p dΓ = 0, ( 15
)
Γ l δθ (-u • n) θ * -Re -1 P r -1 k τ (α, α 0 ) ∂ n θ * + ∂J Γ ∂θ dΓ = 0.
The adjoint boundary conditions on Γ in are obtained as above by taking δθ = 0. For the outlet, we use ∂ n δθ = 0, ∂ n δu = 0 and δp = 0. Finally, we end up with the adjoint boundary condition

On Γ 1 ∪ Γ 2 : u * t = 0, u n θ * + Re -1 P r -1 k τ (α, α 0 ) ∂ n θ * = ∂J Ω ∂θ , u * n = ∂J Ω ∂p , On Γ in : u * n = ∂J Γ ∂p u * t = 0, θ * = 0, On Γ out : u n θ * + Re -1 P r -1 k τ (α, α 0 ) ∂ n θ * = ∂J Γ ∂θ , n p * + n θ θ * + Re -1 ∂ n u * + u n u * = ∂J Γ ∂u . ( 16
)
Now solving (∂L/∂α) [δα] = 0 for all δα yields

∂h τ,α 0 ∂α u • u * + Re -1 P r -1 ∂k τ,α 0 ∂α ∇θ • ∇θ * = 0 in Ω ∂k τ,α 0 ∂α θ * = 0 on Γ 1 , (17)
which is the optimality condition. Finally, according to the adjoint method (see e.g. [START_REF] Herzog | Algorithms for PDEconstrained optimization[END_REF][START_REF] Petra | Model variational inverse problems governed by partial differential equations[END_REF][START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]), the gradient of the cost functional J (u, θ, p) at some α is given by

∂J ∂α (α) = ∂h τ,α 0 ∂α u • u * + Re -1 P r -1 ∂k τ,α 0 ∂α ∇θ • ∇θ * , in Ω, ∂J ∂α (α) = ∂k τ,α 0 ∂α θ * on Γ 1 , (18)
where (u, θ, p) satisfy (2-3) and (u * , θ * , p * ) satisfy the adjoint problem [START_REF] Van Dijk | Level-set methods for structural topology optimization: a review[END_REF][START_REF] Evgrafov | The limits of porous materials in the topology optimization of Stokes flows[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF][START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF].

Remark 1. The adjoint variables (Φ l , ψ l , p out ) used to enforce the boundary conditions (3) of the primal problem are not needed to solve the adjoint problem. Nevertheless, they can be determined using [START_REF] Van Dijk | Level-set methods for structural topology optimization: a review[END_REF][START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF] which cancel many terms in [START_REF] Petra | Model variational inverse problems governed by partial differential equations[END_REF] and gives

Φ l = -n p * -n θ θ * + Re -1 ∂ n u * -(u • n)u * + ∂J Γ ∂u on Γ l ψ l = θ * Re -1 P r -1 k τ (α, α 0 ) + ∂J Γ ∂θ on Γ l , l ∈ {1, 2, in}, Φ out = -Re -1 u * , ψ out = θ * Re -1 P r -1 k τ (α, α 0 ), q out = -u * n - ∂J Γ ∂p .
Note finally that these adjoint variables are also not needed to compute the gradient of the cost functional with respect to the design variable since they does not appear in (18).

Numerical example

This section aims to validate the new interpolation technique. To do so, we investigate two objective functions. First objective function is related to the power dissipated by the fluid through the domain Ω and can be evaluated on the basis of total pressure losses as follows:

J 1 (u, p, θ) = - Γ p t u • n dΓ where p t = p + 1 2 u 2 (19)
The second cost function is related to the maximization of the recoverable thermal power from the domain Ω and is given by

J 2 (u, p, θ) = Γ θ u • n dΓ (20)
We finally consider the following objective functional minimize: We notice that the objective function J is a linear combination of both objective functions described above. It's taking place at the inlet and outlet flow boundary conditions and the combination involves the opposite of J 2 since the optimization problem aims at minimizing the combinatory function J . Using the notation from the general optimization problem (1), we have

J (u, p, θ) = J 1 (u, p, θ) -J 2 (u, p, θ), = - Γ (p t + θ) u • n dΓ
J Ω = 0 and J Γ (u, p, θ) = (p + θ) u • n.
The adjoint model is then given by ( 13)-( 16) where all the partial derivative of J Ω are zero and the following boundary conditions hold

u * = 0, ∂ n θ * = 0 on Γ 1 ∪ Γ 2 , u * t = 0, θ * = 0, u * n = -u n on Γ in , θ * u n + Re -1 P r -1 k τ (α) ∂ n θ * = -u n on Γ out , p * + θ * θ + Re -1 ∇ n u * • n + u * n u n = -(θ + p) on Γ out . (22)
As said in section 4, we propose to solve the topology optimization problem with a steepest descent algorithm where the gradient is computed thanks to the adjoint method. The main flow of the algorithm for the topology optimization is described in Table 3. The forward problem and the optimization

Step 0. Initialization: set all the constants Re, Ri, P r

Step 1. Solve the forward problem (2),(3) problem with the Finite Volume Method

Step 2. Compute objective and constraint values

Step 3. Compute sensitivities by adjoint method

Step 4. Evaluate the optimality condition. If a stopping criterion is met, terminate the calculation

Step 5. Update design variables α with α k+1 = -∇J k+1 + β P R k+1 α k and return to step 1 processes are implemented using OpenFOAM [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. In Step 5, the design variables are evaluated by using the conjugated-gradient descent direction method associated to Polack-Ribiere method

β P R k+1 = ∇J T k+1 (∇J k+1 -∇J k ) ∇J T k ∇J k .
After optimization process, several results are compared as optimized designs, velocity magnitude distribution, temperature distribution, and lastly, we focus on fluid-solid interfaces. We also presents the objective function values for each simulation case. Finally, we investigate to vary the slope abscissa and the curve of sigmoid functions in order to evaluate their impact on the optimization problem. That corresponds to consider different values of α 0 and τ . In our study, α 0 will vary in {25, 50, 100} which correspond respectively to {α max /8, α max /4, α max /2} and τ will vary in {0.3, 0.7, 1.0}. In order to perform topology optimization, the α term may take a finite value between 0 and α max = 200.

Studied Cases

Two representative cases of topology optimization in fluid dynamics are investigated : the single pipe case and the bend pipe case. Figure 2 shows schematic illustrations of the two studied cases.

Design of a single pipe.

We present here a benchmark numerical example used recently by Marck [START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF] to illustrate the viability and efficiency of the methodology presented in this paper. The design domain consists in a cavity of side L. The computational domain is discretized using 100 × 100 elements. The inlet flow is located in the center of the west edge for a length of L/5. The outlet flow is lining up with it on the east edge. The rest of west and east edges are assumed to be adiabatic and will be called Γ 2 , whereas the whole north and 5.1.2. Design of a bend pipe.

The second example has been studied by several authors as [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] and [START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF].

The computational domain is also square-shaped, with a adimensional side L = 1, and the design grid is made of 100 × 100 elements. The inlet flow is located at L/5 of the west edge. The outlet flow is located on the south edge, at L/5 from the east edge.The length of both flow boundary conditions is set to 2L/5. The part of the south edge located on the left of the outlet is subjected to a constant flux of temperature and will be called Γ 1 . The rest of the edges are assumed to be adiabatic and will be designated as Γ 2 .

Physical parameters.

For the two representative cases studied in this section, the inlet flow is prescribed with a constant horizontal velocity equal to U 0 = 0.04 m.s [(Aluminum/Air) Diffusivities], k f = 1.

Results

The aim of the problems defined above is to determine the optimal design that connects the inlet to the outlet of the cavity and that minimizes the objective function J subjected to the constraints (2) and [START_REF] Kreissl | Topology optimization for unsteady flow[END_REF].

First of all, we compute J without optimization for the two studied cases and we obtained respectively for the single pipe J = 9.07 10 -3 and for the bend pipe J = 6. 38 10 -3 . Compared to values issued from optimization process and referenced in Figure 8 and 11, they are reduced between 23% and 42% for the single pipe case, and between 18% and 19% for the bend pipe.

So, the new method investigated here allows to both minimize the power dissipated by the fluid through the domain Ω and maximize the recoverable thermal power from the domain to the fluid.

Secondly, we plot velocity magnitude for various α 0 and τ values. After optimization of J in the single pipe, we obtain an optimized design in which fluid is transported through the single pipe in a straight pipe (see Figure 3). It is quite closed to the optimized design obtained by [START_REF] Marck | Topology optimization of heat and mass transfer problems: laminar flow[END_REF] when the mechanical Similarly, fluid circulation forms a half circle whose center would be the left bottom corner of the square domain in the optimized bend pipe (see Figure 4). Contrary to the literature in fluid mechanics, the optimized path obtained for the flow is not as straight as possible. It can be explained by the fact that we consider here a heat transfer problem. So, the optimized pipe is larger and form a real bend because the fluid flow moves away from hot regions.

Thereafter, we compared distribution of the parameter α in the domain for the two studied cases (see Figures 6 and9). We recall that values less than α 0 are considered as fluid regions. We can notice that this method based on α 0 affected some values of α to the fluid regions, which would not be possible with other methods of the literature. Indeed, generally with the in this porous domain, its existence contributes to increasing thermal power recovered by the fluid. That is illustrated in Figure 5.

At last, it can be seen that thermal conductivity is well distributed, continuously, between fluid and solid regions.

Q t = Ω h τ (α) dΩ α max V tot where V tot is the total volume of Ω (23)
Calculation results are referenced in Figures 7 and10. In the single pipe, when α 0 increases, quantity of material diminishes between 27.5% and 28%.

That means fluid domain is more important and that contributes to the increase of J 2 . So, that will influence the value of J . In the bend pipe, when α 0 increases, quantity of material also diminishes between 3.6% and 16%. So, fluid domain is relatively more important than solid domain but not enough to significantly impact on the value of J , hence, the small effect on values of J (see Figure 11). So, parameter α 0 has a significant incidence on the proportion of fluid domain relatively to solid domain and therefore, can influence the value of the functional objective.

Futhermore, we observe a recirculation zone for α 0 = 100 in the example of the optimized single pipe (see Figure 3). That can also be explained by the value of α 0 . Indeed, larger α 0 means that the proportion of fluid elements compared to α max is bigger and thus fluid regions are more important. When α 0 = 100, that means α 0 = α max /2, the algorithm affects some value of α ∈ [0, 100] corresponding to fluid domain and so impacts on the size of fluid regions relatively to the solid regions in the domain. does not seem to impact significantly the optimization results. However, we can observe in Figure 11 that the frontier between the fluid and the solid is sharper when τ = 2.0. In Figure 8, for τ = 0.3, τ = 0.7, we can notice residual material in fluid domain while they do not exist anymore when τ = 2.0. As the same, the porous region identified in the bend pipe case between the main flow and the hot plate is defined mostly with the fluid diffusivity by the interpolation function. This impact does not influence significantly the result of the optimization computation when we compared their value in Figure 11 but it improves the thermal transition between fluid and solid regions. So, parameter τ seems to improve the thermal conductivity distribution during the process optimization.

Conclusion

This study shows that the new interpolation technique based on sigmoid function is a viable method. Design material is well interpolated since intermediate regions between solid and fluid are suppressed hence giving a frontier between the two regions which is well-defined. Besides, the size of transition zones between fluid and solid regions is explicitly computed. The result is that the new interpolation functions reduce their size comparatively to the RAMP interpolation function. Therefore, the transition between fluid and solid is improving thanks to this sigmoid function for the effective conduc- The parameter τ has a small impact on the results of optimization. For the studied cases, the impact was seen after interpolation of thermal conduc-

Figure 1 :

 1 Figure 1: Interpolation functions with q = 1 and α 0 = α max /4. Left: α max = 200. Right:

Table 1 :× 10 - 4 q 9 Table 2 :

 110492 Size of the transition zones for α max = 200, M = α max /2, ε = 0.1 and α 0 = α max /4. = 10 -4 0.6665 4.9517 × 10 -7 5.4926 × 10 -8 q = 10 -5 0.1666 4.9517 × 10 -8 5.4926 × 10 -Size of the transition zones for α max = 2000, M = α max /2, ε = 0.1 and α 0 = α max /4.

  subject to: Governing equation (2), Boundary conditions (3).

Figure 2 :

 2 Figure 2: Single pipe and Bend pipe

- 1 used

 1 in Reynolds number. The inlet temperature is fixed to θ = 0. The outflow temperature condition is fixed to a zero gradient, as well as the velocity gradient. Hot walls are subject to a constant flux of temperature∂ n θ = -1.The Reynolds number based on the characteristic dimension of the inlet length and the average inlet velocity is Re = 500. The Rayleigh number is fixed to Ra = 5 × 10 5 . Prandtl number is set to 0.71. Grashof number has the following expression Gr = Ra/P r. Adimensional thermal parameters are the following: k s = (9.88×10 -5 m 2 .s -1 )/(2.25×10 -5 m 2 .s -1 )

Figure 3 :

 3 Figure 3: Velocity magnitude distribution in the single pipe

Figure 4 :

 4 Figure 4: Velocity magnitude distribution in the bend pipe

Figure 5 :

 5 Figure 5: Temperature distribution in the bend pipe

Figure 8

 8 represents the distribution of material conductivity after the optimization process in the single pipe case. Compared to Figure7, we observe that frontier between the two regions, fluid and solid, is well-established. For the bend pipe case, although intermediate regions exist, k τ interpolation function allows to affect to this porous domain a conductivity closed to the fluid diffusivity (see Figure11). Therefore, the function k τ (α) provides a continuous transition between the solid and the fluid domains.5.2.1. Influence of α 0 .It can clearly be seen from Figures 7 and 10 that different designs are obtained for the different values of α 0 . As detailed in section 3, α 0 is depending of α max in order to control the percentage of fluid in Ω. This parameter has a significant impact on the quantity of material adding in the domain and that can be illustrated by computing the proportion Q t of material added in the domain Ω as follows :

5. 2 . 2 .

 22 Influence of τ . It can clearly be seen from Figures 7 and 10 that different designs are obtained for the different values of τ . The parameter τ in the interpolation functions h τ (α) and k τ indicates the stiffness of the curve of the functions. It

Figure 6 :

 6 Figure 6: α distribution in the single pipe

Figure 7 :

 7 Figure 7: h τ (α) distribution in the single pipe and proportion of material Q t added in the domain after optimization

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: k τ (α) distribution in the single pipe for different values of τ and α 0 and value of functional objectif J after optimization

Figure 11 :

 11 Figure 11: k τ (α) distribution in the bend pipe and value of functional objectif J after optimization

  tivity. Elements of porous domain have been affected to fluid diffusivity or solid conductivity, and therefore, the thermal transition between material is improved. So, this study has examined a new interpolation technique to deal with fluid-porous media interfaces for topology optimization of heat transfer problems. Further research might explore several type of flow regimes in optimization problems using these interpolation functions. Another possible area of future research would be to investigate new cost functions.

Table 3 :

 3 Algorithm of topology optimization


Andrea Bocco
email: andrea.bocco@cea.fr

Yves Durand
email: yves.durand@cea.fr

Florent De Dinechin
email: florent.de-dinechin@insa-lyon.fr

Hardware support for UNUM floating point arithmetic

The Universal NUMber, or UNUM, is a variable length floating-point format conceived to substitute the current one defined in the IEEE 754 standard. UNUM is able, through an internal algebra based on interval arithmetic, to keep track of the precision during operations, offering better result reliability than IEEE 754.

This work discusses the implementation of UNUM arithmetic and reports hardware implementation results of some of the UNUM operators.

Introduction

The IEEE 754 standard [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF] for floating-point (FP) numbers [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] is widely used by the scientific community. However, FP representation intrinsically introduces several problems, such as: accuracy in iterative algorithms (due to cancellation and accumulation errors), memory footprint and energy consumption (as even low-precision applications require long IEEE 754 FP numbers).

UNUM is an alternative representation format introduced by John L. Gustafson [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. It is a variable-length FP format with an interval-based semantics. It is originally meant for scientific computing, where its improved accuracy and interval semantics can be very precious, or for embedded application which may benefit from reduced memory footprint.

However, as far as we know, the only realizations of UNUM arithmetic have been implemented in software [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. It is not clear if an hardware implementation of UNUM can achieve the claimed potential. The aim of this work is therefore to investigate such an hardware implementation of UNUM, and in particular assess its cost compared to a conventional IEEE 754 FPU.

After a presentation of the UNUM format and algebra in Sections 2, implementation choices left open in [START_REF] Gustafson | The End of Error -Unum Computing[END_REF] are discussed in Section 3. The proposed implementation is described in Section 4 and evaluated and compared to classical FP in Section 5.

The UNUM format and algebra

As Figure 1 shows, a UNUM number is composed of several fields. Compared to standard floating points, it has two distinctive features: First, it is a variable-size, self-descriptive format: the es-1. and fs-1 fields contain the bit lengths of the e (exponent) and f (fraction) fields, minus one. The size of these two fields is itself defined by the UNUM environment [START_REF] Gustafson | The End of Error -Unum Computing[END_REF] (UE). For example, when using a (4,6) UE: 4 bits for es-1 and 6 bits for fs-1, the largest UNUM numbers can count up to 16 bits of exponent and 64 bits of fraction (more than double-precision numbers). The smallest ones have 1 bit of fraction and 1 bit of exponent, but are still encoded on 14 bits because of the other fields.

The second distinctive feature is the u field (or ubit, for "uncertainty" bit). When 0, the UNUM represents a scalar number, whose value encoded into the three fields (s, e, f ). When set, it represents an open interval defined between the scalar value obtained when the ubit is 0, and the one with the fraction field incremented by one Unit in the Last Position (ULP). A u-bit of 1 can be interpreted as "the fraction begins with the provided bits, but there are more and they are unknown".

Finally, a real number can also be represented as an interval defined as a tuple of two UNUMs. Such an interval is called a u-bound [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. It can be open or closed at each endpoint, depending on the corresponding ubit. The set of u-bounds is closed under the basic operations. During such operations, the resulting interval endpoints are rounded toward ±∞, to ensure the property that the real result is included in the resulting interval in spite of rounding error.

The longer-term objective of this work is to propose a parametric hardware+software implementation of a complete UNUM system that addresses the dark spots of [START_REF] Gustafson | The End of Error -Unum Computing[END_REF] e.g. the interval size explosion that plagues interval arithmetic. We are interested in studying how UNUM arithmetic can help deploying known (non-automatic) techniques to contain this interval size explosion.

Closer to this work, the promise of lower consumption has to be challenged against 1/ the increased power budget 978-1-5090-6508-0/17/$31.00 ©2017 IEEE PRIME 2017, Giardini Naxos-Taormina, Italy Digital Circuits and Sub-Systems

UNUM hardware implementation choices

Gustafson organizes the computation of UNUMs into three different layers [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. Well-defined conversion functions are used to move numbers among layers. The h-layer, deals with the data format used for human interaction. We will not address it in this work. The u-layer is where the UNUM formats are defined, corresponds to the FP layer in existing systems. The g-layer is the underground layer where the operations are actually carried out (also called the scratchpad layer). Its format is not necessary identical to the u-layer i.e. numbers may be represented with a fixed size in the g-layer, to match the fact that hardware is fixed. Finally, fused operations can be embedded inside the g-layer to compute complex functions with high accuracy [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. Since the UNUM actually implements an algebra of intervals, the g-layer actually handles g-bounds [START_REF] Gustafson | The End of Error -Unum Computing[END_REF] (the g-layer version of u-bounds). What we will call g-numbers in this work are g-bound endpoints.

As Figures 2 and3 show, the u-layer and g-layer operators are composed of different sub-units, interfaced each other using u-bound and g-bound formats.

A first design choice is to expand variable-size UNUMs into a hardware-friendly fixed-size internal UNUM format which retains the es-1 and fs-1 fields, but pads both exponent and mantissa with zeroes to the maximum size allowed by the UE. The conversions between variable-size UNUMs and this internal format involve shifts (computed out of the es-1 and fs-1 fields), but no rounding or special case management.

A second design choice is the g-number format adopted for the g-bound endpoints. Our solution is depicted on Figure 4. It is a large classical FP format, with two additions: 1/ five summary bits encode exceptional cases, in order to avoid area and latency of decoding them from the exponent field. 2/ an endpoint termination flag, true if the interval endpoint is closed, false if it is open. Like for internal UNUMs, the size of this g-number format is chosen large enough to handle the full set of all possible UNUMs for the chosen UE. The leading bit of the mantissa, which is implicit in IEEE754 formats, is here explicit. A signed exponent (instead of the usual biased one [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF]) simplifies the conversion functions algorithm.

Currently, our internal g-numbers can be subnormals, in which case the leading mantissa bit 0. This conservative choice will be discussed in the conclusion.

UNUM units architecture 4.1. Overview

Four u-bound operators (uOP) have been implemented: an u-bound adder u add, an u-bound comparator u cmp, an u-bound multiplier u mul, and an u-bound divider u div.

They all (except u cmp) share the same macro architecture depicted on Figure 2. The U2G sub-unit converts the input u-bounds into g-bounds. The obtained g-bounds are input to the g-bound operators (gOP) detailed below. The G2U sub-unit converts the obtained g-bound back in ubound (not needed in the u cmp since it directly outputs the comparison result flags computed by the g cmp sub-unit).

The uOP differ only by the gOP sub-unit (g add, g cmp, g mul, g div). They share the same conversion functions.

Each gOP unit (except g cmp) has the same interval processing [START_REF] Schulte | A family of variable-precision interval arithmetic processors[END_REF] macro architecture depicted on Figure 3. The endpoints of the two input g-bounds are used to perform the relevant endpoint operations, rounded to the outside of the interval. The maximum and minimum values are then selected as the endpoints of the output g-bound.

The endpoint computations in gOP units follow a classical floating-point scheme [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]: the input mantissas are aligned to a common format if needed, the operation (+,-,*,/) is computed on the aligned operands, the result is normalized and rounded, and the input exceptions [START_REF] Gustafson | The End of Error -Unum Computing[END_REF] are handled.

The next subsections describe the g-bound operators algorithms. The symbols used are: for a gOP unit, G1 and G2 for the two input g-bounds, and Z for the output gbound. For a g-bound A, we note A and A as its lower and upper endpoints. For a real x we also note (x) (resp.

(x)) the rounding of x to the g-format in the -∞ (resp. +∞) direction.

g-bound adder

The g-bound adder has to compute the sum of two gbound provided in input.

Z = Z = G1 + G2 Z = G1 + G2 (1)
Here neither the input multiplexing unit, nor the minimum and maximum endpoint selection units are needed: the g-bound input endpoints are connected directly to the two endpoint computation units, and their output correspond to the final interval endpoints.

Inside each endpoint computation unit, the final endpoint is computed. The algorithm chosen to sum the two gnumbers is the one shown in [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]: the two input mantissa are aligned to a common exponent value; the aligned mantissa are summed (or subtracted, depending on the input signs) together; The final result is normalized and rounded; The input exception are handled.

g-bound comparator

On intervals, the comparison A < B may be true (if A < B), false (if A ≥ B), or "I can't say" (in the other cases). As a side note, for UNUM to supplant FP, all the programming languages have to be modified to support this third option... But this is out of the scope of this paper.

To allow for that, UNUM comparison unit has to compute the six following comparison flags [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]: Lower Than (LT), Greater Than (GT), EQual (EQ), Not Nowhere EQual (NNEQ), Nowhere EQual (NEQ). In order to do that, it compares each g-bound input endpoint with the two endpoints of the other input. Depending how the two input g-bounds intersect (or not) each others, the four output flags are set or unset, also taking care of the exception values.

The hardware implementation of this function therefore consists of six greater-than comparators and six equal comparators. Their results drive a combinational logic that generates the 6 output flags.

g-bound multiplier

The g-bound multiplier has to compute the multiplication of two g-bound provided in input.

Z = min G1G2 , G1G2 , G1G2 , G1G2 Z = max G1G2 , G1G2 , G1G2 , G1G2 (2)
Fortunately, looking only at the signs allows to eliminate two of the 4 inputs to the max and min operations. As a consequence, with some multiplexing of the inputs endpoints, it is enough to instantiate four endpoint multiplication units, one two-input max, and one two-input min.

In addition, in 15 cases out of 16 input endpoint sign combinations, only one multiplication per endpoint is needed. The proposed implementation exploits this to reduce the average power consumption.

Again, the multiplication algorithm itself is a classical FP one [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. In this case, there is no alignment of the mantissas, the two operand are directly multiplied each other. After that, the intermediate result is normalized and rounded.

g-bound divider

The g-bound divider has to compute the division of two g-bound provided in input.

Z = min G1/G2 , G1/G2 , G1/G2 , G1/G2 Z = max G1/G2 , G1/G2 , G1/G2 , G1/G2
(3) However, it is possible to demonstrate here that, looking at the sign of the input g-bound endpoints, at most one endpoint is a good candidate for each output endpoint. Thanks to that, multiplexing correctly the input g-bound endpoints to the endpoint computation units, it is possible to remove the minimum and maximum sub-units and output the endpoints directly.

Synthesis results and comparison between UNUM and IEEE 754

In this work we have described in hardware and validated the four basic floating point UNUM operators (+,-,*,/). The division was validated with the fixed point mantissa division described as behavioral (/ in VHDL). Unfortunately it is unsuitable for synthesis in this context. Therefore there will be no synthesis result for division in this submission.

The designed units were described in VHDL. They were validated using Mentor Questasim against a reference python library [5] on 2 • 10 6 pseudo-random input vectors. Then they were synthesized using Synopsys Design Compiler with the CORE65LPSVT library.

As reference for the designed UNUM units, a 64bit FPU compliant with the IEEE 754 standard is taken from OpenCores [START_REF]FPU 64 bit IEEE 754 compliant[END_REF]. The UE chosen is the (ess = 4, f ss = 6): it is the smallest UE that includes all double-precision IEEE 754 FP numbers. Its g-layer operators have 64-bit significands, and are therefore slightly more accurate than the double-precision ones (53-bit significands).

Table 1 shows the synthesis results of the proposed UNUM units, decomposed in their sub-units as per Section 4. For each unit is showed the timing latency (ns), the area (number of cell/gate instances chosen by the synthesizer and µm 2 ), and the power consumption (mW ). For each operator replicated inside the g-layer (g add and g mul), PRIME 2017, Giardini Naxos-Taormina, Italy Digital Circuits and Sub-Systems we also factor out the replication. This enables comparison of g add and g mul with the corresponding IEEE-754 operators. Note however that the reported g add and g mul also involve input multiplexers and the minimum and maximum units. In this table we have slightly different results for the U2G and G2U units, although these units are identical in each operator. Obviously the synthesizer optimizes them slightly differently in each context.

Comparing gOP and IEEE units synthesis results, the lessons from this experiment are the following. Our g add and g mul are comparable to the IEEE ones (slightly slower and larger as expected, since their fraction size is 64 instead of 53). The real cost here is the x2 and x4 area overheads due to interval arithmetic. The observed power overhead is limited to a factor 2 in the multiplier, as explained in Section 4.4. The lower power consumption of the g add unit is due to its simpler subnormal handling.

As highlighted in Table 1, gOP units individually are not so different with respect to standard IEEE ones. The main differences (and design complexities) are the conversion functions inside each uOP unit. In fact they almost double its latency. If we are to integrate an UNUM unit in a processor, this is not acceptable. It pleads for exposing the g-layer to the instruction set (with a g-layer register file and instructions that operate in the g-layer), and fusing operations in the g-layer as much as possible (as a compiler optimization). This would dilute the overhead of the conversion functions.

But then, we also need explicit conversions G2U instructions to the UNUM format. Program optimization would attempt to use as few as possible of such instructions, typically only when spilling to external memory. Currently, such conversions compute the smallest matching UNUM format. For iterated computations, the fraction size tends to grow and quickly saturates to its maximum. This will negate the potential power saving due to transmitting UNUM numbers stored on a few bits only. To address this problem, Gustafson suggests starting with a small UE and growing it if needed [START_REF] Gustafson | The End of Error -Unum Computing[END_REF]. A more practical alternative in a hardware context would be that the G2U conversion instructions take as argument the number of significand bits to which to round. The issue, of course, is more burden on the programmer. This idea has to be evaluated on actual applications.

Let us now come back to one design choice: subnormal handling in the g-layer. One valid alternative would be to have only normal numbers in the g-layer. To ensure that all UNUM numbers are representable, we would need one more exponent bit than the maximum allowed by the UE, but the corresponding hardware overhead would be negligible. This would simplify and speed-up the g-layer operations, especially the multiplication. On the other hand it would make the conversions more complex since they would have to handle subnormalisation. However, this is consistent with the idea of a g-layer instruction set with as few as possible G2U instructions.

Conclusion

The new UNUM format is presented with two potential benefits. The first is safer numerical computing thanks to its well defined interval semantics. The second is power reduction in the data transfers thanks to a self-descriptive, variable-size format. However, these two benefits come at the cost of more expensive and more power-hungry computing unit. This work presents a first quantitative assessment of this trade-off. It also discusses some of the design choices that must be made.

The next steps are: to implement and synthesize the g-bound divider unit core, to pipeline the designed units, to complement the proposed UNUM units with a matching register file, to implement a full processor, in order to be able to compare the area and power overhead with the data data transfer with respect to IEEE 754 standard operations.

This requires some deeper rethinking, in particular in the control flow management instructions. For example comparisons require three outputs and not only two as in usual processors. This changes the usual control flow. The programming environment must also expose finer precision control than in classical systems.

Finally, the promises of the UNUM format remain to be demonstrated on full-scale applications. This is the long term objective of this research.

Figure 1 .

 1 Figure 1. The UNUM format

Figure 2 .Figure 3 .

 23 Figure 2. The u-layer operator general architecture

Figure 4 .

 4 Figure 4. The g-number format

TABLE 1 .

 1 SYNTHESIS RESULTS OF UNUM AND IEEE 754 UNITS.

	Unit	Timing (ns)	Area (cells)	(µm 2 )	Power (mW )
	u add U2G in add g add 64-bit FP add G2U in add	16.59 3.54 8.98 8.23 4.07	41852 16573 16113 = 2×8057 7183 9166	132398 52429 50973 = 2×25486 23068 28995	1.89 0.46 0.94 = 2× 0.47 1.75 0.49
	u mul U2G in mul g mul 64-bit FP mul G2U in mul	16.55 3.10 9.20 8.51 4.25	83430 13933 63157 = 4×16790 17635 6341	375472 62704 284232 = 4×71058 77086 28536	8.37 0.64 6.86 = 4×1.72 4.53 0.85
	u cmp U2G in cmp g cmp	4.00 3.38 0.62	20273 18651 1622		66512 61191 5321	2.75 2.37 0.37