Franck P Vidal

Pierre-Frédéric Villard

Member, IEEE Évelyne Lutton

Tuning of Patient Specific Deformable Models using an Adaptive Evolutionary Optimization Strategy

Keywords: Evolutionary computation, inverse problems, medical simulation, adaptive algorithm

We present and analyze the behavior of an evolutionary algorithm designed to estimate the parameters of a complex organ behavior model. The model is adaptable to account for patient's specificities. The aim is to finely tune the model to be accurately adapted to various real patient datasets. It can then be embedded, for example, in high fidelity simulations of the human physiology. We present here an application focused on respiration modeling.

The algorithm is automatic and adaptive. A compound fitness function has been designed to take into account for various quantities that have to be minimized. The algorithm efficiency is experimentally analyzed on several real test-cases: i) three patient datasets have been acquired with the "breath hold" protocol, and ii) two datasets corresponds to 4D CT scans. Its performance is compared with two traditional methods (downhill simplex and conjugate gradient descent), a random search and a basic realvalued genetic algorithm. The results show that our evolutionary scheme provides more significantly stable and accurate results.

I. INTRODUCTION

High fidelity simulations of the human physiology require complex bio-mechanical models. Such models need also to be adaptable to account for patient's specificities. They can be used in various medical contexts, e.g. reducing motion artifacts in positron emission tomography (PET) and conebeam computed tomography (CBCT) to improve the image quantification [START_REF] Lamare | List-mode-based reconstruction for respiratory motion correction in PET using non-rigid body transformations[END_REF], [START_REF] Zhang | Correction of motion artifacts in cone-beam CT using a patient-specific respiratory motion model[END_REF], accurate dose calculation in radiotherapy treatment planning [START_REF] Lewis | A theoretical model for respiratory motion artifacts in free-breathing CT scans[END_REF], and high fidelity computer-based training simulators [START_REF] Morgan | Applying theory to practice in undergraduate education using high fidelity simulation[END_REF]. The calibration and parametrization of the models are therefore critical to preserve realism and numerical accuracy. For training simulators, parameters are, however, often manually tuned using trial and error. This task is time consuming and it is not possible to ascertain that the results are optimal. Throughout the paper, we will consider patient specific respiration modeling with the deformation of the diaphragm and liver as the application example [START_REF] Villard | A prototype percutaneous transhepatic cholangiography training simulator with real-time breathing motion[END_REF]. For each patient, 15 parameters need to be finely tuned. The parametrization of the model is considered as an inverse problem. It corresponds to fitting an analytic model with 15 parameters to experimental data. An optimization technique is deployed to solve the problem. Our approach makes use of F. Vidal is with the School of Computer Science, Bangor University, Dean Street, Bangor LL57 1UT, UK (e-mail: f.vidal@bangor.ac.uk).

P.-F. Villard is with LORIA, University of Lorraine, France (e-mail: pierrefrederic.villard@univ-lorraine.fr).

É. Lutton is with AVIZ, INRIA Saclay-Île-de-France, France (e-mail: evelyne.lutton@inria.fr).

Manuscript received March 13, 2012; revised June 12, 2012.

an ad-hoc evolutionary algorithm (EA) that is able to successfully explore a search space with 15 dimensions (15-D).

We choose an evolutionary framework because i) EAs can be used when little is known about the function to optimize, e.g. when no derivative is known, ii) this function does not need to be very smooth, iii) EAs can work with any search space, and iv) EAs are less likely to stop at local optima than classical deterministic optimisation methods. Our objective is to propose a fully automatic and adaptive methodology that leads to significantly better tuning. Our approach is generic and can be generalized to other models when gold truth is available and the discrepancies between the model outputs and gold truth can be numerically assessed.

Section II briefly presents related work. The analytic model of respiration is described in Section III. Details about the evolutionary algorithm are provided in Section IV. In Section V, the performance of the algorithm is analyzed and compared with the performance of other methods. Section VI presents a discussion about the limitations of our approach. Some conclusions are finally presented in Section VII.

II. RELATED WORK

A. Optimization Based on Artificial Evolution

Artificial Evolution is the generic name of a large set of techniques that rely on the computer simulation of natural evolution mechanisms. Since the pioneering works of A. Fraser, H.-J. Bremermann, and after them, J. Holland and I. Rechenberg, Artificial Darwinism techniques have progressively gained a major importance in the domain of stochastic optimization and artificial intelligence [START_REF] Fogel | Evolutionary Computation: Toward a New Philosophy of Machine Intelligence[END_REF].

Their basic idea is to copy, in a very rough manner, the principles of natural evolution, which let a population be adapted to its environment. According to Darwin's theory, adaptation is based on a small set of very simple mechanisms: random variations, and survival/reproduction of the fittest individuals. Computer scientists have transposed this scheme into optimization algorithms. Their major advantage is to make only a few assumptions on the function to be optimized (there is no need to have a continuous or derivable function for instance). In short, Evolutionary Optimization considers a population of potential solutions exactly as a population of real individuals who lives, fights and reproduces in a natural environment. The pressure of this environment is replaced by an "optimization" pressure: the function to be optimized is considered as a measurement of the adaptation of the individual to its environment. In this way, individuals that reproduce are the best ones with respect to the problem to be solved. The reproduction consists in generating new solutions using a variation scheme (the genetic operators) that, by analogy to nature, is called mutation if it involves one individual, or crossover if it involves several individuals.

Evolutionary optimization techniques are particularly well suited to complex problems, where classical methods fail, e.g. due to the irregularity of the function or to the complexity of the search space. The versatility of the evolutionary framework has produced a variety of different optimization techniques suited for various purposes (multi-objective, interactive, cooperative) aimed at exploring different search spaces (discrete, combinatorial, continuous, tree-based, graphbased, grammar-based, constrained, limited or infinite). The major reason of this success is the tunable combination of oriented and random search mechanisms embedded in an evolutionary algorithm: it allows injecting a priori incomplete information in the genetic operators, while letting some other more unpredictable components be randomly searched.

Success stories in various application domains prove the power of these techniques as long as they rely on this "partial" randomness ability: the most successful applications are often based on hybridization of classical and evolutionary techniques. In such cases evolution is used as an "orchestra conductor", that combines a priori information, and classical "local" optimizers in the evolution loop.

Considering evolutionary optimization as a "black box" is not a good strategy in general because one may lose an opportunity to adapt to the problem. Adapting the evolutionary mechanisms to the specificities of the problem usually improves the efficiency of the algorithms and reduces its computation time. In particularly difficult problems it is always helpful to compare evolutionary approaches to a pure random optimization to evaluate the improvement due to the "intelligence" set in the genetic operators. This is what we do in Section V. In our context, the random search algorithm generates a given number of sets of 15 random parameters [START_REF] Villard | A method to compute respiration parameters for patient-based simulators[END_REF]. The set of parameters that provides the lowest fitness is extracted and corresponds to the solution of the optimization problem.

B. Medical Applications of Artificial Evolution

The use of evolutionary techniques in medicine is not a recent trend: J.T. Alander clearly shows the growth of interest for this topic since 1990 in the EA community [START_REF] Alander | Indexed bibliography of genetic algorithms in medicine[END_REF]. In recent reviews, S.L. Smith and S. Cagnoni [START_REF] Smith | Introduction to the special issue on medical applications of genetic and evolutionary computation[END_REF], [START_REF]Genetic and Evolutionary Computation: Medical Applications[END_REF] identified five main trends in this domain:

• medical imaging and signal processing (segmentation, registration, change detection, correction, denoising, 3D reconstruction) • medical data and patient records mining • clinical expert systems and knowledge-based expert systems • therapy (tests, implants tuning) • modeling and simulation It is worth to note that these various applications take advantage of the flexibility of the evolutionary scheme, and rely on many variants of it (cooperation, multi-objective, agentbased approaches). It is also possible in many cases to control the trade-off between computational cost and precision of the results via a simple parameter tuning.

The work presented here is related to the first and last items of the trend list presented above. In the literature, model fitting usually relies on generic models, such as deformable image models in 2D [START_REF] Ballerini | Genetic snakes for color images segmentation[END_REF], [START_REF] Mcintosh | Spinal crawlers: Deformable organisms for spinal cord segmentation and analysis[END_REF], and their extensions in 3D [START_REF]Vessel Crawlers: 3D physically-based deformable organisms for vasculature segmentation and analysis[END_REF], [START_REF] Hill | Model-based interpretation of 3D medical images[END_REF]. Similar work on 3D surface registration [START_REF] Chow | Surface registration using a dynamic genetic algorithm[END_REF] is also dealing with generic models for segmentation purposes.

Once again for a segmentation task, Heimann et al. address the problem of liver surface modeling [START_REF] Heimann | A shape-guided deformable model with evolutionary algorithm initialization for 3D soft tissue segmentation[END_REF]. They use statistical shape models combined with a freely deformable, energybased snake model, and were able to obtain an average surface distance of about 1.6 mm in comparison to manual reference segmentations.

Here, we are relying on an accurate respiration model (see Section III) that necessitates the simultaneous fitting of two surface models, for the diaphragm and the liver. In this paper, the evolutionary engine that has been developed is an attempt to deal with these two objectives by linear combination, with weights that adjust dynamically to balance the relative importance of the two objectives (see Section IV-B). To our best knowledge this problem has never been tackled using evolutionary computing.

C. Respiration Models

Respiratory models can be used in various contexts. The method depends on the application that is targeted. In radiotherapy, deformable models have to be accurate and focused on the targeted anatomical structures (e.g. tumor or organs at risk). Such models are generally based on continuous laws of mechanics, and equations are resolved using the finiteelement technique [START_REF] Al-Mayah | Effect of friction and material compressibility on deformable modeling of human lung[END_REF]. In computer graphics animation, the reality is not the focus. It is rather the impression of realism. Therefore many details and organs should be modeled. In [START_REF] Zordan | Breathe easy: model and control of simulated respiration for animation[END_REF], diaphragm, ribs, intercostal muscles and the external skin are modeled and the deformation is achieved by using the massspring technique [START_REF] Nedel | Real time muscle deformations using mass-spring systems[END_REF]. In medical simulation, the realism as well as the precision are needed to teach the users the real conditions of an operation. However the key challenge is often the computing time, e.g. both the visual and haptic renderings require a high refresh rate. In [START_REF] Hostettler | A real-time predictive simulation of abdominal viscera positions during quiet free breathing[END_REF], the method is entirely based on geometrical constraints and the viscera are only modeled using a single envelope that wraps all the organs. We propose here a method that takes into account the motion of each organ due to the respiration using a geometrical method based on mechanical-based parameters.

III. RESPIRATION MODEL

A. Anatomy and Physiology Analysis

Respiration is a complex process that mainly involves two muscles: the diaphragm and the intercostal muscles [START_REF] Moore | Clinically Oriented Anatomy[END_REF]. Respiration has also an indirect influence on all the organs of the abdominal cavity, in particular the liver. It is attached to the diaphragm by the falciform ligament, and hereby follows its deformation and motion. The diaphragm is composed of two parts: i) the tendon where the lungs and the heart lie, and ii) the muscle part that contracts and relaxes. The diaphragm is attached to the lower ribs and to the spine. In our case, the patient is lying on his/her back. We can therefore assume that the influence of the intercostal muscle is negligible.

B. Soft-Tissue Deformation

We will mainly focus here on i) the liver, which is the treated organ, and ii) the diaphragm, which is the active muscle as previously seen. We choose to compute the deformation using the generalized ChainMail extension [START_REF] Gibson | 3D chainmail: a fast algorithm for deforming volumetric objects[END_REF]. The main advantage of this model is its small computing time. Instead of computing the deformation field based on time integration of the forces as other physically-based methods [START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF], it only uses geometrical equations that are quickly computed. It also preserves the use of parameters that have links with bio-mechanical approaches: the compression (S i ), the stretching (Str i ) and the shearing (Sh i ) where i designates the concerned organ.

To decompose the diaphragm into two parts as pointed out in Section III-A, we choose the Cartesian equation of a plane (a.x + b.y + c.z + d = 0). As boundary conditions, we chose to impose a null displacement when the diaphragm is close to the ribs (defined by a distance D ribs ) and to impose a uniform 3D force to the whole tendon part (defined by (F x , F y , F z )). The muscle part should deform. This deformation is governed here by the muscle intrinsic elasticity. The muscle stretches from the point attached to the ribs to the moving tendon.

The liver attachments to the diaphragm are modeled by a distance (D diaph ) determining which points directly follow the diaphragm. This rigid link simulates the compression of the liver by the diaphragm during the respiration process. The other parts of the liver deforms following the ChainMail algorithm.

C. Parameter Analysis

Various parameters have been extracted from the model previously described (see Fig. 1). These parameters are unique to each patient and need to be individually customized. There are bio-mechanical parameters, anatomy-based geometrical constrains and respiration pattern information. More details on these parameters can be found in [START_REF] Villard | A method to compute respiration parameters for patient-based simulators[END_REF].

D. Model Evaluation

To optimize the model parameters we need datasets of different patients but also metrics to evaluate the accuracy of the respiration model. The simulation always starts at the real inspiration state (Geometry G I R ) extracted from the real data. The aim is to reach the real expiration state (Geometry G E R ) also extracted from the real data. In the 4D CT scan case (Geometry G t R ), going through intermediate states t is also achievable. The metrics assesses the accuracy of a simulated geometry G t S compared to the ground truth geometry extracted at the same time step t.

To compute the difference between two geometries, we choose to analyze for each mesh vertex the point-to-surface distance according to [START_REF] Aspert | MESH: Measuring Errors between Surfaces using the Hausdorff distance[END_REF]. It is based on a distance measure

d (p, G ′ ) = min p ′ ∈G ′ p -p ′ (1)
In this application, the root mean square error

(E RMS (G t S , G t R )
) defined by Eq. ( 2) is used.

E RMS G t S , G t R = 1 G t S p∈G t S d (p, G t R ) 2 dS (2)
IV. OPTIMIZATION ALGORITHM An evolutionary algorithm is used to finely tune the parameters of the deformable model. The core of the genetic engine has been written in C++ using templates. It makes the approach generic, so that it can be used for other optimization problems. However, for each specific optimization problem, it is necessary to implement the corresponding genotype (e.g. the search space), some of the genetic operators (e.g. mutation and crossover), as well as the fitness function. This section describes in detail how the EA has been implemented.

A. Genotype

Fifteen parameters in total have been identified in Sections III-B and III-C and are illustrated in Fig. 1. The evolutionary algorithm explores this 15-D search space. Each parameter corresponds to a real number, whose boundary values are known. The genome of individuals is therefore made of 15 floating point numbers.

B. Adaptive Fitness Function

A metric (E RMS (M 0 , M 1 )) is presented in Section III-D to evaluate the discrepancies between two polygon meshes M 0 and M 1 . It is therefore possible to quantify the difference between the mesh (G t S (i)) simulated using the deformable model with the parameters corresponding to Individual i, and the real mesh (G t R ) extracted from the patient's dataset at state t. For each individual, two metrics are computed (one for the diaphragm, and one for the liver). These values can be used to define the fitness function (f itness) corresponding to i. The optimization consists in minimizing f itness. The simplest function is:

f itness(i) = αE RMS G t S (i, diaph), G t R (diaph) + (1 -α)E RMS G t S (i, liver), G t R (liver) (3
) with 0 ≤ α ≤ 1 to give more or less weight to the diaphragm or the liver. Selecting the value of α is not trivial because the numerical quantity of the error for the diaphragm and the liver can be significantly different. If the same weight is applied to both tissue error measurement (i.e. α is equal to 0.5), the predominant quantity will then have more influence on the optimization process. We would actually expect instead the importance of both tissues to be the same during the minimization. In an application such as the real-time simulation of liver puncture, a higher level of fidelity is required for the liver than the diaphragm. Scaling factors on errors are introduced to give the same relative weight to the diaphragm and liver. Eq. 3 becomes:

f itness(i) = α E diaph RMS E RMS G t S (i, diaph), G t R (diaph) + (1 -α) E liver RMS E RMS G t S (i, liver), G t R (liver) (4)
with E liver RMS and E diaph RMS the error metrics of the best individual provided by the previous generation for the liver and diaphragm respectively. For each iteration of the evolution loop, these metrics are updated. For practical reasons, α is rescaled so that the sum α

E diaph RM S + (1-α) E liver RM S
is equal to one.

C. Genetic Engine

An elitist generational genetic approach is used (see Fig. 2). For each optimization loop performed by the genetic engine, a new population based on the previous generation is created. The population size is made constant. Each gene of an individual is initialized using an uniform distribution to randomly chose a value within its range of validity. During the optimization loop, a genetic operator is randomly chosen to create a new individual. Genetic operators are: i) elitism, ii) mutation, iii) crossover, and iv) new blood. They are described in detail below. Each operator is assigned a probability of occurrence, W , X, Y , and Z respectively. The sum of the operators' probability is 100%. Some operators (e.g. the mutation and the crossover) require one or several individuals of the previous generation to be randomly chosen for the reproduction (i.e. the creation of a new individual based on the chosen individual(s)). This process is performed by the selection operator. The optimization stops when a stopping criteria is reached.

1) Tournament Selection: During the evolution, individuals are selected using tournaments. For each selection, n individuals are randomly chosen, without any bias. The best individual amongst the n chosen ones is the individual who is selected.

2) Elitism: W % best individuals of the previous generation are kept in the new population. In practice, even if W is equal to 0%, the best individual is kept. In this case, if nothing better is found, the best individual is not lost.

3) Mutation: X% of the new individuals are created by mutating individuals of the previous generation. It mimics mutations in a genomic sequence that are observed in molecular biology and genetics. For each mutation, an individual is selected using the selection operator. The genotype of the selected individual is copied into the genotype of the new individual. The numerical value of each gene is then slightly altered to introduce spontaneous and random changes. The range of possible random changes is controlled by the value σ. In classical implementations, this value is fixed. However theoretical studies have shown that a variable, adaptive or selfadaptive mutation is beneficial in optimization [START_REF] Bäck | Optimal mutation rates in genetic search[END_REF]. Adaptive mutation is actually a scheme that has been made experimentally evident in natural population (see for instance [START_REF] Hersh | Adaptive mutation and amplification in Escherichia coli: two pathways of genome adaptation under stress[END_REF] about stress mutation in bacteria populations). Here, we have chosen to use an adaptive strategy for σ, in the style of stress mutation. The adaptation rule produces greater σ for "bad individuals" (i.e. large fitness) and smaller ones for "good individual" (small fitness). The idea is to favor large exploration around bad individuals, whilst performing fine tuning in the vicinity of good individuals. σ is then controlled depending on the fitness value. It makes use of a scaled and stretched cosine function:

σ(f ) =          σ min , f < f it min σ max , f > f it max σ min + (σ max -σ min )× cos π× f -f it min f itmax -f it min +1.0 2.0 , otherwise (5)
f corresponds to the fitness of the individual who will undergo a mutation. σ(f ) smoothly varies between the smaller (f it min ) and the larger (f it max ) fitness thresholds respectively. If f is smaller than f it min , σ is then σ min ; if the individual's fitness is greater than f it max , σ is then σ max (with σ min and σ max two constant values set by the user).

Mutation is performed via the addition of a shift offset to the genome's value, whose value is:

offset = range 2 × σ(F itness(i)) × k (6)
with range the interval size of the valid values of the genomes, and σ(F itness(i)) the amount of mutation corresponding to the fitness of Individual i, and k a random number in the interval [-1, 1].

4) Blend Crossover: Y % of the new individuals are created using BLX (for blend crossover) [START_REF] Herrera | A taxonomy for the crossover operator for real-coded genetic algorithms: An experimental study[END_REF]. For each crossover operation, two individuals (i 1 and i 2 ) are selected using the selection operator. Each component (or chromosome) of the genome of the new individual will be given by a random combination of the corresponding chromosome of i 1 and i 2 as follows:

C = kC 1 + (1 -k)C 2 (7)
with C the numerical value of the chromosome of the new individual, C 1 the numerical value of i 1 's chromosome, C 2 the numerical value of i 2 's chromosome, and k a random number between 0 and 1. For each chromosome, a new k value is chosen.

5) New Blood: Z% of individuals are replaced by new random individuals. Such individuals are created using the same method as during the initialization. This operator allows "new blood" to continually enter the population. It corresponds to a "pure random search" component and maintains a minimal diversity level in the population.

6) Stopping Criteria: Stagnation is the stopping criteria that has been chosen, i.e. the optimization process is stopped when artificial evolution does not improve the result. In practice, when s successive evolution loops provide exactly the same best individual, the optimization process is stopped, and this individual is the solution to the optimization problem.

V. RESULTS AND VALIDATION

This section shows the performance of the evolutionary algorithm in term of accuracy. The results are compared with the outputs of: i) a pure random search to evaluate the improvement provided by the genetic operators, ii) a basic real-valued genetic algorithm (GA) used as a black box evolutionary optimization to assess the efficiency of our new genetic operators, and iii) more traditional methods -Downhill simplex method [START_REF] Nelder | A simplex method for function minimization[END_REF] and Powell's conjugate gradient descent method [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF] -for further comparisons. For each tested case, every stochastic optimization process has been repeated 15 times. For each optimization process, the errors were recorded and taken into account during the statistical tests presented in this section. It allows us to ascertain the effectiveness and usefulness of the evolutionary algorithm, i.e. to demonstrate that it outperforms the brute force algorithm, the black box EA, and the classic methods.

A. Input Data

Abdominal images with respiration information are routinely acquired during radiotherapy treatment planning when knowledge of the tumor displacement is required to perform the dose calculations. Five patient specific datasets have been selected (see Table I). Three datasets (Patients A, B and C) have been acquired with the "breath hold" protocol, i.e. with only two time steps corresponding to the inhale and exhale states. The patients are asked to hold their breath following the "ABC" protocol [START_REF] Wong | The use of active breathing control (ABC) to reduce margin for breathing motion[END_REF]. Two datasets (Patients D and E) correspond to 4D CT scans with ten time steps each. The data was acquired over the respiratory cycle while the patient breathes normally. For Patients A, B and C, a single optimization problem each needs to be solved. For patients with 4D CT scans, one optimization problem per time step needs to be solved. Note that for such patients, the values presented in this section correspond to average values over the 10 time steps.

Every CT scan has been segmented to extract the organs that are required by the simulation model. Medical expert collaborators have manually segmented the data using a graphic tablet. The contouring process has been performed using ITKSnap1 . Polygon meshes were then exported after using the Marching Cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF]. Meshes were decimated and smoothed to have about 2,000 vertices per organ. The manual segmentation process of a CT scan takes on average 2 hours per liver and 3 hours per diaphragm. Fig. 3 shows an input medical image with the manual segmentation of Patient B's diaphragm. In total 26 3D CT scans were segmented.

B. Evolutionary Algorithm Parameters

Table II provides a summary of the algorithm's parameters. Two different sizes of population (50 and 200) are used to assess the stability of the algorithm. α is set to 33% to give more weight to the liver than the diaphragm, without neglecting the diaphragm.

C. Performance Comparison

Box plots have been produced for each patient. Fig. 4 presents the results for Patient B as an illustration. Figures for other patients present the same features. It shows the errors for both the diaphragm and the liver after optimization using the different methods. A population of 200 individuals provides relatively stable results with the Ad hoc EA. This population size has therefore been used for further statistical analysis with other methods. Note that for each patient care is given to use the same amount of computing time for the random search than our ad-hoc evolutionary algorithm (see Table III for the numerical values). In one case (Patient E), the basic realvalued genetic algorithm fails to converge in an acceptable time.

Fig. 5 presents the evolution of the diaphragm weight (α/E diaph RMS ) in the fitness function of Patient A (note that the weight of the liver is complimentary and it corresponds to 1 -α/E diaph RMS . The weights are updated for each iteration of the evolution loop (see Section IV-B). It shows that the final weight are consistent amongst the 15 runs. For each run, the shape of the curve is similar. The same phenomena are observed for every patient. It is therefore acceptable to use the "evolved" fitness function for comparison with other optimization methods when applicable.

Table IV shows the average root mean square error (E RMS ) of the liver and diaphragm for each patient and for each optimization technique. It also presents the result of a nonparametric statistical hypothesis test -the Wilcoxon signedrank test -that is used to compare the performance of our Ad hoc EA with other methods. The size of the samples is 15. If the value of the Wilcoxon signed-rank test (W) is negative, then the Ad hoc EA performs better than the other algorithm. If W is close to 0, then both algorithms are performing similarly. The results show that only our evolutionary algorithm can minimize successfully both the error of the liver and diaphragm. The only positive W value has been obtained for Patient A's diaphragm using the gradient descent method. However the corresponding value for the liver was negative. Classic optimization methods fail to explore the 15-D search space to minimize the two error values. This is particularly true for the downhill simplex method. The basic real-valued genetic algorithm is very slow and fails to converge in some cases (see Patient E in Table III). In our application, artificial evolution can effectively explore a large search space without a priori knowledge.

Table V presents the resulting parameters found by our method for the 15 optimization runs of each patient. The average value and the standard deviation are given for each parameter. Some parameters found over 15 runs are relatively close (e.g. the diaphragm plane (a, b, c and d)) whilst others may significantly differ (e.g. D ribs and Str l ).

VI. DISCUSSION

The comparison with pure random search provides a good assessment of the added value of the genetic engine, i.e. the improvement due to the random search oriented by the genetic operators over a "blind" random search. The pure random search used the "evolved" fitness function for comparison. It represents a disadvantage for the evolutionary approach. It is difficult to evaluate the algorithmic cost of the adaptation of the fitness function in the evolutionary algorithm. It is therefore not considered in the comparisons displayed in Fig. 4 and Tables III andIV

A Patient B Patient C Patient D Patient E a 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 b 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 0.0 ± 0.1 c
-0.9 ± 0.0 -0.9 ± 0.0 -0.9 ± 0.0 -0.9 ± 0.0 0.9 ± 0.0 d 114. a given amount of computation time, the best results are obtained using our evolutionary algorithm.

Our genetic engine is also more efficient than the basic real-valued genetic algorithm used as a black box optimizer. Our method provides more accurate results using less computational power. Our algorithm also outperforms the classic optimization methods that were used for comparison purposes.

The diaphragm errors are larger than liver errors for Patients B, D and E. This is due to the fact that the spatial or time resolution of the initial datasets was lower than in the other test cases. For Patient B, the slice thickness is 5 mm, i.e. much larger than the diaphragm thickness (less than 2 mm). Patients D and E make use of 4D CT data. In such cases, the inaccuracies are due to the motion artifact (known for blurring tissue edges). For Patients B, D and E, the reconstructed surface of the diaphragm is then over-smoothed. The quality of the geometrical models is then less accurate. The diaphragm is a relatively thin tissue. Inaccuracies are attenuated for the liver as it is a much larger organ than the diaphragm.

Let's analyze the value of parameters presented in Table V. The plane coefficients (a and b) are very close to zero. It was to be expected as it roughly corresponds to a transverse plane. The force (F x , F y , F z ) is mainly orthogonal to this plane. The ChainMail parameters are very high for both organs. This is due to the large deformations that occur: the contraction and relaxation of diaphragm muscle as well as the soft-tissue behavior of the liver. Finally, the attachment distance is higher from the diaphragm to the ribs (D ribs ) than from the liver to the diaphragm (D diaph ). This is because D diaph mimics a very thin ligament attaching the liver to the diaphragm.

Finally, Fig. 6 shows 3D plots of surface meshes for all the patients. The printed color depends on a lookup table (LUT) that corresponds to the localized error. Its range varies from blue for no error to red for the maximum error. The motion is fairly well compensated using our genetic algorithm. VII. CONCLUSION AND FUTURE WORK We have presented an artificial evolution strategy to finely tune the parameters of a multidimensional model of respiration with soft tissue deformations. Its efficiency has been validated using five datasets of real patients (that is 23 different optimization problems in total). The advantage of artificial evolution over the downhill simplex, the conjugate gradient descent, the purely random search, and a black box basic realvalued genetic algorithm has also been demonstrated. Results obtained using our artificial evolution framework were both more accurate and more stable.

The proposed evolutionary optimization is adaptive in two ways: i) mutation is adapted based on fitness, and ii) the weight the two-objectives compound fitness is automatically balanced along calculations.

We also demonstrated that this compound fitness function effectively takes into account various properties of the model, e.g. minimizing several error values for the liver and the diaphragm. The approach can be generalized to other models when gold truth is available and the discrepancies between the model outputs and gold truth can be numerically assessed.

The current solution that is to balance the different objectives in a single fitness function can be revisited. We will investigate other composition strategies (e.g. multiplicative instead of additive). Other evolutionary approaches will also be considered. A cooperative-coevolutionary approach [START_REF] Pena-Reyes | Fuzzy CoCo: A cooperativecoevolutionary approach to fuzzy modeling[END_REF] can be used as the problem we presented here includes most of the features that have been identified to be difficult to solve using single-population evolutionary algorithms but it is efficiently addressable using cooperative-coevolution:

• the search space is complex (e.g. 15-D)

• the problem can be split (e.g. minimizing the error for both the liver and the diaphragm) • there is a strong interdependence between the components of the solution (e.g. the diaphragm is influencing the liver) Also a Classical multi-objective evolutionary approach like the famous NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] will be also considered for dealing with multiple objectives. Since these methods are able to perform an optimization without setting any priority between the objectives, they provide a set of potential compromises between the various objectives, the Pareto set. To be applicable in our context it is necessary to select a single acceptable solution from the Pareto set at the end of the process. The superiority of such a strategy in our context is not obvious, as the time spent to obtain a good approximation of the whole Pareto set may be better used to perform calculations focused on a single compound and well chosen objective.

FFig. 1 .

 1 Fig. 1. List of the respiration model parameters.

Fig. 2 .

 2 Fig. 2. Evolutionary loop.

Fig. 3 .

 3 Fig. 3. Diaphragm manual segmentation (here in white).

Fig. 4 .

 4 Fig. 4. Root mean square errors for the Patient B's diaphragm and liver using the Ad hoc EA, pure random search, basic real-valued GA, downhill simplex, and conjugate gradient descent methods.

Fig. 5 .

 5 Fig. 5. Evolution of the diaphragm weight, α/E diaph RM S , for 15 optimization runs for Patient B.

Fig. 6 .

 6 Fig. 6. 3D plots of surface meshes with localized errors. The first row shows the initial difference map between real inhale and real exhale states. The second row shows the difference between real exhale and simulated exhale with our genetic algorithm.

TABLE I PATIENT

 I DATASET PROPERTIES.

TABLE III AVERAGE

 III NUMBER OF FITNESS EVALUATIONS.

  . Despite this structural disadvantage, the evolutionary approach surpasses the random search. For

			Patient A	Patient B	Patient C	Patient D	Patient E
			ERMS W ERMS W ERMS W ERMS W ERMS W
	Ad hoc EA	Diaph Liver	2.05 2.95	N.A. N.A.	4.53 N.A. 2.42 N.A.	2.06 N.A. 2.85 N.A.	3.63 2.16	N.A. N.A.	8.84 4.43	N.A. N.A.
	Pure random search	Diaph Liver	2.27 3.18	-96 -88	4.77 3.02	-84 -118	2.17 2.94	-102 -120	4.06 2.24	-93 -42	9.30 4.87	-45 -34
	Basic real-valued GA	Diaph Liver	2.21 3.13	-102 -44	4.58 2.69	-18 -120	2.10 2.88	-76 -78	3.81 2.45	-12 -35	N.A. N.A. N.A. N.A.
	Downhill simplex	Diaph Liver 21.89 -120 3.67 -120	7.15 6.55	-120 -120 13.04 -120 9.65 -120	5.22 8.39	-119 -120 12.72 -120 9.95 -68
	Diaph Conjugate gradient descent Liver	2.00 5.96	96 -120	5.49 6.55	-120 -120	2.33 5.47	-120 -120	5.33 3.55	-112 -109	9.91 6.41	-33 -50

TABLE IV PERFORMANCE

 IV COMPARISON OF THE Ad hoc EVOLUTIONARY ALGORITHM, PURE RANDOM SEARCH, BASIC REAL-VALUED GENETIC ALGORITHM, DOWNHILL SIMPLEX, AND CONJUGATE GRADIENT DESCENT METHODS. E RM S IS THE AVERAGE ROOT MEAN SQUARE ERROR IN MILLIMETERS. W IS WILCOXON SIGNED-RANK TEST BETWEEN THE Ad hoc EA AND THE OTHER METHODS.

	Patient

  1 ± 7.8 113.7 ± 9.2 110.8 ± 13.5 117.3 ± 10.1 -99.5 ± 14.8

	C d	0.4	± 0.2	0.4	± 0.2	0.6	± 0.2	0.5	± 0.2	0.4	± 0.2
	Str d	5.1	± 1.3	5.2	± 1.8	5.6	± 1.9	5.2	± 2.1	5.5	± 1.7
	Sh d	4.7	± 1.7	3.9	± 2.0	4.9	± 2.0	5.0	± 2.4	4.5	± 1.9
	Fx	0.6	± 0.6	3.9	± 0.4	0.6	± 0.1	2.0	± 2.3	0.3	± 2.7
	Fy	3.5	± 0.2	1.8	± 0.7	6.2	± 0.7	-2.4 ± 2.0	-1.3 ± 2.8
	Fz	-5.8 ± 1.0	-8.1 ± 0.5	-7.9 ± 0.4	1.3	± 1.8	4.2	± 3.9
	D ribs	19.7 ± 14.2 15.9 ± 13.6 19.9 ± 5.1	11.3 ± 7.8	11.7 ± 10.6
	C l	0.5	± 0.1	0.4	± 0.2	0.5	± 0.2	0.5	± 0.2	0.5	± 0.2
	Str l	8.2	± 6.8	5.9	± 3.7	16.7 ± 5.5	5.2	± 2.0	4.9	± 1.9
	Sh l	1.5	± 1.7	4.3	± 2.6	5.9	± 1.9	4.6	± 2.4	2.6	± 1.6
	D diaph	3.3	± 0.7	1.3	± 0.7	0.7	± 0.3	4.1	± 2.2	4.0	± 0.7

TABLE V PARAMETERS

 V AFTER OPTIMIZATION USING THE Ad hoc EA.

http://www.itksnap.org/ (last accessed: 31/05/2012)


Francis Filbet

Chang Yang

INVERSE LAX-WENDROFF METHOD FOR BOUNDARY CONDITIONS OF BOLTZMANN TYPE MODELS

Keywords: Inverse Lax-Wendroff procedure, WENO, Boltzmann type models

, to the kinetic equations. Applications in 1D × 3D and 2D × 3D of this algorithm for Boltzmann type operators (BGK, ES-BGK models) are then presented and numerical results illustrate the accuracy properties of this algorithm.

Introduction

We are interested in the numerical approximation of solutions to kinetic equations set in a complex geometry with different type of boundary conditions. Unfortunately, classical structured or unstructured meshes are not appropriate due to the high dimensional property of kinetic problem. In contrast, the Cartesian mesh makes the numerical method efficient and easy to implement. The difficulty is that obviously grid points are usually not located on the physical boundary when using a Cartesian mesh, thus a suitable numerical method to capture the boundary condition on the complex geometry is required. Several numerical methods based on Cartesian mesh have been developed in computational fluid dynamics in last decade. Among these methods, the immersed boundary method (IBM), first introduced by Peskin [START_REF] Peskin | Flow patterns around heart valves. A numerical method[END_REF] for the study of biological fluid mechanics problems, has attracted considerable attention because of its use of regular Cartesian grid and great simplification of tedious grid generation task. The basic idea of immersed boundary method is that the effect of the immersed boundary on the surrounding fluid is represented through the introduction of forcing terms in the momentum equations. In conservation laws, two major classes immersed boundary like methods can be distinguished on different discretization types. The first class is Cartesian cut-cell method [START_REF] Ingram | Developments in Cartesian cut cell methods[END_REF], which is based on a finite volume method. This conceptually simple approach "cuts" solid bodies out of a background Cartesian mesh. Thus we have several polygons (cut-cells) along the boundary. Then the numerical flux at the boundary of these cut-cells are imposed by using the real boundary conditions. This method satisfies well the conservation laws, however to determine the polygons is still a delicate issue. The second class is based on finite difference method. To achieve a high order interior scheme, several ghost points behind the boundary are added. For instance for solving hyperbolic conservations laws, an inverse Lax-Wendroff type procedure is used to impose some artificial values on the ghost points [START_REF] Tan | Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws[END_REF]. The idea of the inverse Lax-Wendroff procedure (ILW) is to use successively the partial differential equation to write the normal derivatives at the inflow boundary in terms of the tangential and time derivatives of given boundary conditions. From these normal derivatives, we can obtain accurate values of ghost points using a Taylor expansion of the solution from a point located on the boundary.

The goal of this paper is to extend the inverse Lax-Wendroff procedure to kinetic equations together with an efficient time discretization technique [START_REF] Filbet | On deterministic approximation of the Boltzmann equation in a boundary domain[END_REF][START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF] for problems where boundary conditions play a significant role in the long time asymptotic behavior of the solution. In particular, for low speed and low Knudsen flows for which DSMC methods are unsuitable. Therefore, the main issue relies on that the inflow is no longer a given function, while it is determined by the outflow. For this, we proceed in three steps: we first compute the outflow at the ghost points. To maintain high order accuracy and to prevent oscillations caused by shocks, we use a weighted essentially non-oscillatory (WENO) type extrapolation to approximate the ghost points by using the values of interior mesh points. In the same time, we can extrapolate the outflow located at the boundary associated with ghost points. Then, we compute the inflow at the boundary by using the outflow obtained in the first step and Maxwell's boundary conditions. Finally, we perform the inverse Lax-Wendroff procedure to approximate the inflow on the ghost points, where the key point is to replace the normal derivatives by a reformulation of the original kinetic equation.

For simplicity, we only consider simple collision operator as we adapt the ellipsoidal statistics BGK or ES-BGK model introduced by Holway [9]. This model gives the correct transport coefficients for Navier-Stokes approximation, so that Boltzmann or ES-BGK simulations are expected to give the same results for dense gases. Let us emphasize that Direct Simulation Monte-Carlo methods (DSMC) have been performed to the ES-BGK model in complex geometry. However DSMC approach is not computationally efficient for nonstationary or low Mach number flows due to the requirement to perform large amounts of data simpling in order to reduce the statistical noise. In contrast, F. Filbet & S. Jin recently proposed a deterministic asymptotic preserving scheme for the ES-BGK model, where the entire equation can be solved explicitly and it can capture the macroscopic fluid dynamic limit even if the small scale determined by the Knudsen number is not numerically resolved [START_REF] Filbet | An Asymptotic Preserving Scheme for the ES-BGK model of the Boltzmann equation[END_REF]. We will use this scheme to solve ES-BGK model while on the boundary the inverse Lax-Wendroff procedure will be applied.

The outline of the paper is as follows. In Section 2 we describe precisely the inverse Lax-Wendroff procedure to Maxwell's boundary condition in 1D and 2D space dimension. Then in Section 3 we present the ES-BGK model and the application of inverse Lax-Wendroff procedure to this model. In Section 4 a various numerical examples are provided in 1D × 3D and 2D × 3D to demonstrate the interest and the efficiency of our method in term of accuracy and complexity. Finally a conclusion and some perspectives are given in Section 5.

Numerical method to Maxwell's boundary conditions

The fundamental kinetic equation for rarefied gas is the Boltzmann equation

(2.1) ∂f ∂t + v ⋅ ∇ x f = 1 ε Q(f ),
which governs the evolution of the density f (t, x, v) of monoatomic particles in the phase, where

x ∈ Ω x ⊂ R d x , v ∈ R 3 .
The collision operator is either given by the full Boltzmann operator

(2.2) Q(f )(v) = R 3 S 2 B( v -v ⋆ , cos θ) f ′ ⋆ f ′ -f ⋆ f dσ dv ⋆
or by a simplified model as the BGK or ES-BKG operator (see the next section). Boltzmann's type collision operator has the fundamental properties of conserving mass, momentum and energy: at the formal level

R 3 Q(f ) φ(v) dv = 0, φ(v) = 1, v, v 2 .
Moreover, the equilibrium is the local Maxwellian distribution namely:

M[f ](t, x, v) = ρ(t, x) (2π T (t, x)) 3 2 exp - u(t, x) -v 2 2 T (t, x) ,
where ρ, u, T are the density, macroscopic velocity and the temperature of the gas, defined by

(2.3) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ρ(t, x) = R 3 f (t, x, v)dv, u(t, x) = 1 ρ(t, x) R 3 vf (t, x, v)dv, T (t, x) = 1 3 ρ(t, x) R 3 u(t, x) -v 2 f (t, x, v)dv.
In order to define completely the mathematical problem (2.1), suitable boundary conditions on ∂Ω x should be appled. Here we consider wall type boundary conditions introduced by Maxwell [START_REF] Maxwell | [END_REF], which is assumed that the fraction (1α) of the emerging particles has been reflected elastically at the wall, whereas the remaining fraction α is thermalized and leaves the wall in a Maxwellian distribution. The parameter α is called accommodation coefficient [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF].

More precisely, for x ∈ ∂Ω x the smooth boundary ∂Ω x is assumed to have a unit inward normal n(x) and for v ⋅ n(x) ≥ 0, we assume that at the solid boundary a fraction α of particles is absorbed by the wall and then re-emitted with the velocities corresponding to those in a still gas at the temperature of the solid wall, while the remaining portion (1α) is perfectly reflected. This is equivalent to impose for the ingoing velocities

(2.4) f (t, x, v) = (1 -α) R[f (t, x, v)] + α M[f (t, x, v)], x ∈ ∂Ω x , v ⋅ n(x) ≥ 0, with 0 ≤ α ≤ 1 and (2.5) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ R[f (t, x, v)] = f (t, x, v -2(v ⋅ n(x))n(x)), M[f (t, x, v)] = µ(t, x) f w (v).
By denoting T w the temperature of the solid boundary, f w is given by

(2.6) f w (v) ∶= exp - v 2 2 T w ,
and the value of µ(t, x) is determined by mass conservation at the surface of the wall for any

t ∈ R + and x ∈ ∂Ω x (2.7) µ(t, x) v⋅n(x)≥0 f w (v)v ⋅ n(x)dv = - v⋅n(x)<0 f (v)v ⋅ n(x)dv.
This boundary condition (2.4) guarantees the global conservation of mass [START_REF] Filbet | On deterministic approximation of the Boltzmann equation in a boundary domain[END_REF].

In this paper we only apply a second order finite difference method to discretize the transport term of (2.1) but higher order schemes [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] may be applied. Then to keep the order of accuracy of the method, two ghost points should be added in each spatial direction. To impose f at the ghost points, we will apply the inverse Lax-Wendroff procedure proposed in [START_REF] Tan | Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws[END_REF] for conservation laws.

Suppose that the distribution function f at time level t n for all interior points are already known, we now construct f at the ghost points.

2.1. One-dimensional case in space. We start with spatially one-dimensional problem, that is d x = 1. In this case the Boltzmann equation reads:

(2.8) ∂f ∂t + v x ∂f ∂x = 1 ε Q(f ), (x, v) ∈ [x l , x r ] × R 3 ,
where x l and x r are the left and right boundaries respectively, v x is the component of phase space corresponding to x-direction. For the boundary condition in spatially one-dimensional case, the inward normal on the boundary in (2.4) is

n(x l ) = ⎛ ⎜ ⎝ 1 0 0 ⎞ ⎟ ⎠ , n(x r ) = ⎛ ⎜ ⎝ -1 0 0 ⎞ ⎟ ⎠ .
To implement the numerical method, we assume the computational domain is a limited domain

[x min , x max ] × [-V, V ] 3 , where (x l , x r ) ⊂ [x min , x max ]. The computational domain is covered by a uniform Cartesian mesh X h × V h , (2.9) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ X h = {x min = x 0 ≤ ⋯ ≤ x i ≤ ⋯ ≤ x nx = x max } , V h = v j = j ∆v, j = (j 1 , . . . , j 3 ) ∈ Z 3 , j ≤ n v .
with the mesh size ∆x and ∆v for space and velocity respectively. We only consider numerical method of ghost points near the left hand side boundary, since the procedure for right hand side boundary is the same. Figure 1 illustrates a portion of mesh near left boundary x l , which is located between x 0 and x 1 .

x

-1 x 0 x l x 1 x 2 x 3 ◾ ◾ • • • ◾ ◾ • • • x v v j v j ⋆ Figure 1.
A portion of mesh in spatially one dimensional case. • is interior point, ◾ is ghost point, is the left hand side boundary.

We construct f at each ghost point in following three steps: we perform an extrapolation of f to compute a high order approximation of the outflow. Then, we compute an approximation of the distribution function at the boundary using Maxwell's boundary conditions. Finally, we apply the inverse Lax-Wendroff procedure for the inflow.

2.1.1. First step: Extrapolation of f for the outflow. At time t = t n we consider the outflow near the point x l , that is f (t, x l , v j ) where v j 1 < 0. We denote by f i,j an approximation of f at (x i , v j )

A natural idea is to extrapolate f at the left boundary x l or the ghost points x 0 and x -1 using the values of f on interior points. For example from the values f 1,j , f 2,j and f 3,j , we can construct a Lagrange polynomial p 2 (x) ∈ P 2 (R). Then by injecting x l , x 0 or x -1 into p 2 (x), we obtain the approximations of f at the ghost points and left boundary, i.e. f l,j , f 0,j and f -1,j . However, when a shock goes out of the boundary, the high order extrapolation may lead to a severe oscillation near the shock. To prevent this, we would like have a lower order accurate but more robust extrapolation. Therefore, a WENO type extrapolation [START_REF] Tan | Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws[END_REF] will be applied and described below (see subsection 2.3) for this purpose.

Second step: Compute boundary conditions at the boundary.

In the previous step, the outflow at the boundary is obtained by extrapolation. To compute the vlaues of f at the inflow boundary, we apply the Maxwell's boundary condition (2.4), i.e.

(2.10)

f l,j = (1 -α) R[f l,j ] + α M[f l,j ].
On the one hand the specular reflection portion is given straightly by the outflow at the left boundary, which is

R[f l,j ] = f l,j ⋆
, where j ⋆ = (-j 1 , j 2 , . . . , j 3 ).

On the other hand the diffuse one is computed by a half Maxwellian

M[f l,j ] = µ l exp - v j 2 2T l ,
where T l is the given temperature at the left wall and µ l is given by

µ l v j ⋅n(x l )≥0 v j ⋅ n(x l ) exp - v j 2 2 T l = - v j ⋅n(x l )≤0 v j ⋅ n(x l ) f l,j .
2.1.3. Third step: Approximation of f at the inflow boundary. Finally we compute the values of f at the ghost points for the inflow boundary. Here we cannot approximate f by an extrapolation, since the distribution function at interior points cannot predict the inflow.

Then we extend the inverse Lax-Wendroff type procedure recently proposed in [START_REF] Huang | Numerical boundary conditions for the fast sweeping high order WENO methods for solving the Eikonal equation[END_REF][START_REF] Tan | Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws[END_REF][START_REF] Xiong | Fifth order fast sweeping WENO scheme for static Hamilton-Jacobi equations with accurate boundary treatment[END_REF] for solving kinetic equations. At the left boundary x l , a first order Taylor expansion gives

f j (x) = f l,j + (x -x l ) ∂f ∂x x=x l + O(∆x 2 ).
Hence a first order approximation of f at ghost points is

(2.11) f s,j = f l,j + (x s -x l ) ∂f ∂x x=x l , s = -1, 0.
We already have f l,j in the second step, thus it remains to obtain an approximation of the first derivative. By reformulating (2.8), we have

(2.12) ∂f ∂x x=x l = 1 v x - ∂f ∂t + 1 ε Q(f ) x=x l .
Now instead of approximating the first derivative ∂ x f x=x l , we compute the time derivative ∂ t f x=x l and the collision operator Q(f ) x=x l . An approximation of the time derivative can be computed by using several f i,j at previous time levels. Different approximation are obtained either a first order approximation reads

∂f ∂t x=x l = f n l,j -f n-1 l,j ∆t ,
where ∆t is the time step, or one can use a WENO type extrapolation to approximate the time derivative (see subsection 2.3 below).

The last term Q(f ) x=x l can be computed explicitly by using f l,j obtained in previous two steps. Clearly this procedure is independent of the values of f at interior points.

Remark 2.1. Let us observe that when α = 0 we have a pure specular reflection boundary condition. A mirror procedure can be used to approximate f at the ghost points. More precisely, by considering the boundary as a mirror, we approximate the distribution at the ghost points f (x s , v j ) as

f (x s , v j ) = f (2 x l -x s , v j ⋆ ), where j ⋆ = (-j 1 , j 2 , . . . , j 3 )
where 2 x lx s is the mirror image point of x s . Since 2 x lx s is located in interior domain, we can approximate f (2 x lx s , v j ⋆ ) by an interpolation procedure.

2.2. Two-dimensional case in space. The previous approach can be generalized to spatially two-dimensional problem. We assume

d x = 2 in equation (2.1) (2.13) ∂f ∂t + v x ∂f ∂x + v y ∂f ∂y = 1 ε Q(f ),
where the distribution function

f (t, x, v) is defined in (t, x, v) ∈ R + × Ω × R 3 with x = (x, y).
We consider a computational domain

[x min , x max ] × [y min , y max ] × [-V, V ] 3 , such that Ω ⊂ [x min , x max ] × [y min , y max ] and f (t, x, v) ≊ 0, for all v ≥ V .
The computational domain is covered by an uniform Cartesian mesh X h × V h , where X h , V h are defined similarly to (2.9). The mesh steps are respectively ∆x, ∆y and ∆v. In Figure 2, we present a portion of spatial mesh near the boundary. From a ghost point x g , we can find an inward normal n, which crosses the boundary at x p . For the 2D case in space, the numerical approximation of the distribution function f at ghost points is similar to the one dimensional case. However, there are two major differences.

x y n θ ◾ ◾ ◾ ◾ ◾ ◾ ◾ • • x g ◾ • • • • x p • • • • • • • • • • i x -2 i x -1 i x i x + 1 i x + 2 i y -2 i y -1 i y i y + 1 i y + 2
First to compute R[f ] in the second step, the corresponding outflow may not locate on phase space mesh. Secondly to approximate the normal derivative in the third step, besides the time derivative and collision operator we need also the tangential derivative at x p . Once again, we present the method in three steps: 2.2.1. First step: Extrapolation of f for outflow. Let us assume that the values of the distribution function f on the grid points in Ω are given. To approximate f at a ghost point, for instance x s , we first construct a stencil E composed of grid points of Ω for the extrapolation. For instance as it is shown in Figure 2, the inward normal n intersects the grid lines y = y iy , y iy+1 , y iy+2 . Then we choose the three nearest point of the cross point in each line, i.e. marked by a large circle. From these nine points, we can build a Lagrange polynomial q 2 (x) ∈ Q 2 (R 2 ). Therefore we evaluate the polynomial q 2 (x) at x s or x p , and obtain an approximation of f at the boundary and at ghost points. As for the 1D case, a WENO type extrapolation can be used to prevent spurious oscillations, which will be detailed in subsection 2.3.

Second step:

Compute boundary conditions at the boundary. In the previous step, we have obtained the outflow f (x p , v ⋅ n < 0) at the boundary x p . By using (2.4) as we did for the 1D case, we can similarly compute the distribution function f for v ⋅ n ≥ 0. However this time to compute the distribution function for specular reflection

R[f (x p , v)] = f (x p , v -2(v ⋅ n)n), ∀v ∈ V h ,
the vector fields v -2 (v ⋅ n) n may not be located on a grid point. Therefore, we interpolate

f in phase space (x p , v -2(v ⋅ n)n) using the values computed from the outflow f (x p , v) such that v ⋅ n ≥ 0.
2.2.3. Third step: Approximation of f at the inflow boundary. We have obtained the values of f at the boundary points x p for all v ∈ V h in previous two steps. Now we reconstruct the values of f for the velocity grid points such that v ⋅ n ≤ 0 at the ghost point x g by a simple Taylor expansion in the inward normal direction. To this end, we set up a local coordinate system at x p by

x = x ŷ = cos θ sin θ -sin θ cos θ x y ,
where θ is the angle between the inward normal n and the x-axis illustrated in Figure 2. Thus the first order approximation of f (x g , v) reads

f (x g , v) ≊ f (x p , v) + (x g -xp ) ∂ f ∂ x (x p , v), where f (x p , v) = f (x p , v) and ∂ f ∂ x (x p , v)
is the first order normal derivative at the boundary x p . To approximate ∂ f ∂ x (x p , v), we use inverse Lax-Wendroff procedure. Firstly, we rewrite the equation (2.13) in the local coordinate system as

(2.14) ∂ f ∂t + vx ∂ f ∂ x + vy ∂ f ∂ ŷ = 1 ε Q( f ),
where vx = v x cos θ + v y sin θ, vy = -v x sin θ + v y cos θ. Then a reformulation of (2.14) yields

(2.15) ∂ f ∂ x (x p , v) = - 1 vx ∂ f ∂t + vy ∂ f ∂ ŷ - 1 ε Q( f ) x=xp .
Finally instead of approximating 2). We then perform an essentially non-oscillatory (ENO) procedure [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes III[END_REF] for this numerical differentiation to avoid the discontinuity.

WENO type extrapolation.

A WENO type extrapolation [START_REF] Tan | Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws[END_REF] was developed to prevent oscillations and maintain accuracy. The key point of WENO type extrapolation is to define smoothness indicators, which is designed to help us choose automatically between the high order accuracy and the low order but more robust extrapolation. Here we describe this method in spatially 1D and 2D cases. Moreover we will give a slightly modified version of the method such that the smoothness indicators are invariant with respect to the scaling of f . 2.3.1. One-dimensional WENO type extrapolation. Assume that we have a stencil of three points E = {x 1 , x 2 , x 3 } showed in Figure 1 and denote the corresponding distribution function by f 1 , f 2 , f 3 . Instead of extrapolating f at ghost point x g by Lagrange polynomial, we use following Taylor expansion

f g = 2 k=0 (x g -x l ) 2 k! d k f dx k x=x l
We aim to obtain a (3k)-th order approximation of

d k f dx k x=x l denoted by f (k) l , k = 0, 1, 2.
Three candidate substencils are given by S r = {x 1 , . . . , x r+1 }, r = 0, 1, 2.

In each substencil S r , we could construct a Lagrange polynomial p r (x)

∈ P r (R) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ p 0 (x) = f 1 , p 1 (x) = f 1 + f 2 -f 1 ∆x (x -x 1 ), p 2 (x) = f 1 + f 2 -f 1 ∆x (x -x 1 ) + f 3 -2f 2 + f 1 2∆x 2 (x -x 1 )(x -x 2 ).
We now look for the WENO type extrapolation in the form

f (k) l = 2 r=0 w r d k p r (x) dx k (x l ),
where w r are the nonlinear weights depending on f i . We expect that f (k) l has (3k)-order accurate in the case f (x) is smooth in S 2 . The nonlinear weights are given by

w r = α r ∑ 2 s=0 α s , with α r = d r (ε + β s ) 2 ,
where ε = 10 -6 and β r are the new smoothness indicators determined by

β 0 = ∆x 2 , β 1 = 1 ε + f 2 1 + f 2 2 2 l=1 x 1 x 0 ∆x 2l-1 d l dx l p 1 (x) 2 dx = (f 2 -f 1 ) 2 ε + f 2 1 + f 2 2 β 2 = 1 ε + f 2 1 + f 2 2 + f 2 3 2 l=1 x 1 x 0 ∆x 2l-1 d l dx l p 2 (x) 2 dx = 1 12(ε + f 2 1 + f 2 2 + f 2 3 ) (61f 2 1 + 160f 2 2 + 25f 2 3 + 74f 1 f 3 -192f 1 f 2 -124f 2 f 3 ).
We remark that the smoothness indicators β 1 and β 2 have the factors

1 ε+∑ r+1 m=1 f 2 m
, which guarantee that the indicators are invariant of the scaling of f i .

Two-dimensional extrapolation.

The two-dimensional extrapolation is a straightforward expansion of 1D case. The substencils S r , r = 0, 1, 2 for extrapolation are chosen around the inward normal n such that we can construct Lagrange polynomial of degree r. For instance in Figure 2, the three substencils are respectively

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ S 0 = {(x ix , y iy )}, S 1 = {(x ix-1 , y iy ), (x ix , y iy ), (x ix , y iy+1 ), (x ix+1 , y iy+1 )},
S 2 = {(x ix-1 , y iy ), (x ix , y iy ), (x ix+1 , y iy ), (x ix-1 , y iy+1 ), (x ix , y iy+1 ), (x ix+1 , y iy+1 ), (x ix , y iy+2 ), (x ix+1 , y iy+2 ), (x ix+2 , y iy+2 )}.

Once the substencils S r are chosen, we could easily construct the Lagrange polynomials in

Q r (R 2 ) q r (x) = r m=0 r l=0 a l,m x l y m satisfying q r (x) = f (x), x ∈ S r . Then the WENO extrapolation has the form (2.16) f (x) = 2 r=0 w r q r (x), x ∈ S r ,
where w r are the nonlinear weights, which are chosen to be

w r = α r ∑ 2 s=0 α s , with α r = d r (ε + β r ) 2 , where ε = 10 -6 , d 0 = ∆x 2 + ∆y 2 , d 1 = ∆x 2 + ∆y 2 , d 2 = 1 -d 0 -d 1 .
β r are the smoothness indicators determined by

β 0 = ∆x 2 + ∆y 2 , β r = 1 ε + ∑ x∈Sr f (x) 2 1≤ α ≤r K K α -1 (D α q r (x)) 2 dx, r = 1, 2,
where α is a multi-index and K = [x p -∆x 2, x p + ∆x 2]× [y p -∆y 2, y p + ∆y 2], x p = (x p , y p ).

Application to the ES-BGK model

The Boltzmann equation (2.1) governs well the evolution of density f in kinetic regime and also in the continuum regime [START_REF] Filbet | An Asymptotic Preserving Scheme for the ES-BGK model of the Boltzmann equation[END_REF]. However the quadratic collision operator Q(f ) has a rather complex form such that it is very difficult to compute. Hence different simpler models have been introduced. The simplest model is the so-called BGK model [START_REF] Bhatnagar | A model for collision processes in gases. Small amplitude processes in charged and neutral one-component systems[END_REF], which is mainly a relaxation towards a Maxwellian equilibrium state

(3.1) Q(f ) = τ ε (M[f ] -f ),
where τ depends on macroscopic quantities ρ and T .

Although it describes the right hydrodynamical limit, the BGK model does not give the Navier-Stokes equation with correct transport coefficients in the Chapman-Enskog expansion. Holway et al. [START_REF] Holway | Kinetic theory of shock structure using an ellipsoidal distribution function[END_REF] proposed the ES-BGK model, where the Maxwellian M[f ] in the relaxation term of (3.1) is replaced by an anisotropic Gaussian G[f ]. This model has correct conservation laws, yields the Navier-Stokes approximation via the Chapman-Enskog expansion with a Prandtl number less than one, and yet is endowed with the entropy condition [START_REF] Andries | The Gaussian-BGK model of Boltzmann equation with small Prandtl numbers[END_REF]. In order to introduce the Gaussian model, we need further notations. Define the opposite of the stress tensor

(3.2) Θ(t, x) = 1 ρ R 3 (v -u) ⊗ (v -u)f (t, x, v)dv.
Therefore the translational temperature is related to the T = tr(Θ) 3. We finally introduce the corrected tensor

T (t, x) = [(1 -ν)T I + νΘ](t, x),
which can be viewed as a linear combination of the initial stress tensor Θ and of the isotropic stress tensor T I developed by a Maxwellian distribution, where I is the identity matrix.

The ES-BGK model introduces a corrected BGK collision operator by replacing the local equilibrium Maxwellian by the Gaussian G[f ] defined by

G[f ] = ρ det(2πT ) exp - (v -u)T -1 (v -u) 2 .
Thus, the corresponding collision operator is now

(3.3) Q(f ) = τ ε (G[f ] -f ),
where τ depends on ρ and T , the parameter -1 2 ≤ ν < 1 is used to modify the value of the Prandtl number through the formula

2 3 ≤ Pr = 1 1 -ν ≤ +∞.
It follows from the above definitions that

(3.4) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ R 3 f (v) dv = R 3 G[f ](v) dv = ρ, R 3 v f (v) dv = R 3 v G[f ](v) dv = ρ u, R 3 v 2 2 f (v) dv = R 3 v 2 2 G[f ](v) dv = E and ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ R 3 (v -u) ⊗ (v -u) f (v) dv = ρ Θ, R 3 (v -u) ⊗ (v -u) G[f ] dv = ρ T .
This implies that this collision operator does indeed conserve mass, momentum and energy as imposed.

In this section, we will first recall the implicit-explicit (IMEX) scheme to the ES-BGK equation proposed in [START_REF] Bhatnagar | A model for collision processes in gases. Small amplitude processes in charged and neutral one-component systems[END_REF]. Then we apply our ILW procedure to treat the boundary condition for ES-BGK model case.

3.1. An IMEX scheme to the ES-BGK equation. We now introduce the time discretization for the ES-BGK equation (2.1), (3.3)

(3.5) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂f ∂t + v ⋅ ∇ x f = τ ε (G[f ] -f ), x ∈ Ω ⊂ R dx , v ∈ R 3 , f (0, x, v) = f 0 (x, v), x ∈ Ω, v ∈ R 3 ,
where τ depends on ρ, u and T . The time discretization is an IMEX scheme. Since the convection term in (3.5) is not stiff, we will treat it explicitly. The source terms on the right hand side of (3.5) will be handled using an implicit solver. We simply apply a first order IMEX scheme,

(3.6) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ f n+1 -f n ∆t + v ⋅ ∇ x f n = τ n+1 ε (G[f n+1 ] -f n+1 ), f 0 (x, v) = f 0 (x, v) .
This can be written as

f n+1 = ε ε + τ n+1 ∆t [f n -∆t v ⋅ ∇ x f n ] + τ n+1 ∆t ε + τ n+1 ∆t G[f n+1 ], (3.7)
where G(f n+1 ) is the anisotropic Maxwellian distribution computed from f n+1 . Although (3.7) appears nonlinearly implicit, since the computation of f n+1 requires the knowledge of G[f n+1 ], it can be solved explicitly. Specifically, upon multiplying (3.7) by φ(v) defined by

φ(v) ∶= 1, v, v 2 2
and use the conservation properties of Q and the definition of

G[f ] in (2.
3), we define the macroscopic quantity U by U ∶= (ρ, ρ u, E) computed from f and get

U n+1 = ε ε + τ n+1 ∆t R 3 φ(v) (f n -∆t v ⋅ ∇ x f n ) dv + τ n+1 ∆t ε + τ n+1 ∆t R 3 φ(v)G[f n+1 ](v) dv,
or simply

(3.8) U n+1 = R 3 φ(v) (f n -∆tv ⋅ ∇ x f n ) dv .
Thus U n+1 can be obtained explicitly. This gives ρ n+1 , u n+1 and T n+1 . Unfortunately, it is not enough to define G[f n+1 ] for which we need ρ n+1 Θ n+1 . Therefore, we define the tensor Σ by

(3.9) Σ n+1 ∶= R 3 v ⊗ v f n+1 dv = ρ n+1 Θ n+1 + u n+1 ⊗ u n+1
and multiply the scheme (3.7) by v ⊗ v. Using the fact that

R 3 v ⊗ v G[f ](v) dv = ρ (T + u ⊗ u) ,
and (3.9), we get that

Σ n+1 = ε ε + (1 -ν) τ n+1 ∆t Σ n -∆t R 3 v ⊗ v v ⋅ ∇ x f n dv (3.10) + (1 -ν) τ n+1 ∆t ε + (1 -ν) τ n+1 ∆t ρ n+1 T n+1 I + u n+1 ⊗ u n+1
Now G[f n+1 ] can be obtained explicitly from U n+1 and Σ n+1 and then f n+1 from (3.7).

Finally the scheme reads

(3.11) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ U n+1 = R 3 φ(v) (f n -∆tv ⋅ ∇ x f n ) dv, Σ n+1 = ε ε + (1 -ν) τ n+1 ∆t Σ n -∆t R 3 v ⊗ v v ⋅ ∇ x f n dv + (1 -ν) τ n+1 ∆t ε + (1 -ν) τ n+1 ∆t ρ n+1 T n+1 I + u n+1 ⊗ u n+1 . f n+1 = ε ε + τ n+1 ∆t [f n -∆t v ⋅ ∇ x f n ] + τ n+1 ∆t ε + τ n+1 ∆t G[f n+1 ],
The scheme (3.11) is an AP scheme for (3.6). On the one hand, although (3.6) is nonlinearly implicit, is can be solved explicitly. On the other hand, the scheme (3.11) preserves the correct asymptotic [START_REF] Filbet | An Asymptotic Preserving Scheme for the ES-BGK model of the Boltzmann equation[END_REF], which means when holding the mesh size and time step fixed and letting the Knudsen number go to zero, the scheme becomes a suitable scheme for the limiting hydrodynamic models.

3.2.

Inverse Lax-Wendroff procedure for boundary conditions. We have described the numerical method for boundary condition to general kinetic equations in spatially 1D and 2D case. To implement this method, it remains to replace the collision operator Q(f ) in (2.12) or (2.15) by the ES-BGK operator (3.3).

Assume that the approximation to the distribution function at the boundary f (x p , v) is known for all v ∈ V h . Then, the macroscopic quantities ρ, u and T at the boundary point x p can be obtained using (2.3) and (3.4). Therefore, substituting these macroscopic quantities in (3.2), we compute the stress tensor Θ at the boundary point x p , such that the corrected tensor T (x p ). Thus G[f ] is computed for all points (x p , v), where v ∈ V h .

Numerical examples

In this section, we present a large variety of test cases in 1d x and 2d x in space and three dimensional in velocity space showing the effectiveness of our method to get an accurate solution of Boltzmann type equations set in a complex geometry with different boundary conditions. We first give an example on a flow generated by gradients of temperature, which has already been treated by DSMC or other various methods [START_REF] Filbet | On deterministic approximation of the Boltzmann equation in a boundary domain[END_REF].

Finally, we present some numerical results in 2d x .

4.1. Smooth solutions. We consider the ES-BGK equation (2.1)-

(3.3) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂f ∂t + v x ∂f ∂x = Q(f ), , x ∈ (-0.5, 0.5), v ∈ R 3 , f (t = 0) = f 0 ,
with an initial datum f 0 which is a perturbation of the constant state in space and a Maxwellian distribution function in velocity, that is,

f 0 (x, v) = ρ 0 (x) (2π) exp - v 2 2 , x ∈ (-0.5, 0.5), v ∈ R 3
with a density ρ 0 = 1 + 0.1 cos(2π x). We consider purely diffusive boundary conditions with a wall temperature T w = 1. The solution is expected to be smooth for short time and then may develop a discontinuity at the boundary, which may propagate in the physical domain.

We perform several numerical simulations on a time interval [0, t end ] with t end = 1, a computational domain in space I 0 = [-π 6, π 6] such that (-1 2, 1 2) ⊂ I 0 and a domain in velocity V = [-8, 8] 3 . Then, we choose a grid in space for I 0 constituted of n x = n points and a grid V h for the velocity space with n v = n points for each direction with respectively n = 32, 64,..., n = 512. Let us emphasize that the boundary points x = -1 2 and x = 1 2 are not exactly located on a grid point. Since we don't know an exact solution of the problem, we compute relative errors. More precisely, an estimation of the relative error in L 1 norm at time T is given by e 2h = f h (T )f 2h (T ) L 1 , where f h represents the approximation computed from a mesh of size h = (∆x, ∆v). The numerical scheme is said to be k-th order if e 2h ≤ C h k , for all 0 < h ≪ 1.

In Table 1 we compute the order of convergence in L 1 norm of our numerical methods. We can clearly see the expected second order convergence. Moreover, we verify experimentally that our scheme is also second-order accurate at the boundary since the discontinuity occuring at v x = 0 is perfectly located.

n x × n v L 1 error
Order L 1 error at the boundary Order 32 2 8.8833 10 -4 X 3.909 10 -3 X 64 2 2.5221 10 -4 1.94 5.832 10 -3 X 128 2 6.5511 10

(3.3), ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂f ∂t + v x ∂f ∂x = 1 ε Q(f ), x ∈ (-1 2, 1 2), v ∈ R 3 , f (t = 0, x, v) = 1 2π T 0 exp - v 2 2 T 0 ,
with T 0 (x) = 1 and we assume purely diffusive boundary conditions on x = -1 2 and x = 1 2, which can be written as

f (t, x, v) = µ(t, x) f w (v), if (x, v x ) ∈ {-1 2} × R + and (x, v x ) ∈ {1 2} × R -,
where µ is given by (2.7). This problem has already been studied in [START_REF] Rader | Direct simulation Monte Carlo convergence behavior of the hard-sphere-gas thermal conductivity for Fourier heat flow[END_REF] using DSMC for the Boltzmann equation or using deterministic approximation using a BGK model for the Boltzmann equation in [START_REF] Ohwada | Investigation of heat transfer problem of a rarefied gas between parallel plates with different temperatures[END_REF][START_REF] Filbet | On deterministic approximation of the Boltzmann equation in a boundary domain[END_REF]. Here we apply our numerical scheme with the ES-BGK operator (3.3) and choose a computational domain in space I 0 = [-π 6, π 6] such that (-1 2, 1 2) ⊂ I 0 and [-8, 8] 3 for the velocity space with a number grid points n = 32 in each direction and the time step ∆t = 0.001.

The main issue here is to capture the correct steady state for which the pressure is a perturbation of a constant state with a Knudsen layer at the boundary [START_REF] Ohwada | Investigation of heat transfer problem of a rarefied gas between parallel plates with different temperatures[END_REF][START_REF] Rader | Direct simulation Monte Carlo convergence behavior of the hard-sphere-gas thermal conductivity for Fourier heat flow[END_REF].

In Figure 4, we represent the stationary solution (obtained approximately at time t end = 25 for ε = 0.1 up to t end = 75 for ε = 0.025) of the temperature and the pressure profile. The results are in a qualitative good agreement with those already obtained in [START_REF] Rader | Direct simulation Monte Carlo convergence behavior of the hard-sphere-gas thermal conductivity for Fourier heat flow[END_REF] with DSMC. More precisely, the boundary layer (Knudsen layer) appears in the density and temperature as well as the pressure, but it is small for all the quantities. The magnitude in the dimensionless density, temperature, and pressure is of order of ε and the thickness of the layer is, say O(ε). In the density and temperature profiles, we cannot observe it unless we magnify the profile in the vicinity of the boundary (see the zoom in Figure 4). Instead, since the pressure is almost constant in the bulk of the gas, we can observe perfectly the boundary layer by magnifying the entire profile. Let us emphasize that, as it is shown in Figure 4 the Knudsen layer is a kinetic effect, which disappears in the fluid limit (ε → 0).

These results provide strong evidence that the present treatment of boundary conditions using WENO extrapolation and inverse Lax-Wendroff method can be used to determine the state of a gas under highly non-equilibrium conditions. Using deterministic methods, we can investigate the behavior of gases for situations in which molecular diffusion is important e.g., thermal diffusion. Also let us mention that a quantitative comparison between our results (3d v with ES-BGK operator) and [START_REF] Rader | Direct simulation Monte Carlo convergence behavior of the hard-sphere-gas thermal conductivity for Fourier heat flow[END_REF] (3d v Boltzmann with hard sphere potential) or [START_REF] Ohwada | Investigation of heat transfer problem of a rarefied gas between parallel plates with different temperatures[END_REF] (3d v BGK) gives a very good agreement on the values of the Knudsen layer and the values of the pressure inside the domain.

4.3.

High-speed flow through a trapezoidal channel. In this section we deal with spatially two-dimensional ES-BGK model in a trapezoidal domain. We attempt to get some steady state as

v x ∂f ∂x + v y ∂f ∂y = 1 ε Q(f ),
where v ∈ Ω x and v ∈ R 3 . Here we will reproduce a numerical test performed in [START_REF] Rjasanow | Stochastic Numerics for the Boltzmann Equation[END_REF] but with our ILW method. The computational domain is a trapezoid

Ω x = {x = (x, y), 0 < x < a, 0 < y < b + x tan(α)}
as shown in Figure 5 for the parameters a = 2.0, b = 0.4, α = arctan(0.2). Boundary conditions are defined separately for each of the four straight pieces

∂Ω x = Γ l ∪ Γ b ∪ Γ r ∪ Γ t
denoting the left, bottom, right and top parts of the boundary respectively. The bottom part represents the axis of symmetry, so we use specular reflection (2.5) there, i.e.

f (x, v) = (f (x, v -2(v ⋅ n(x))n(x)), x ∈ Γ b , v y > 0.
On the right part we are modeling outflow (particles are permanently absorbed), i.e. On the left part there is an incoming flux of particles, i.e.

(4.2) f (x, v) = f in (x, v) = M in (v), x ∈ Γ l , v x > 0, with an inflow Maxwellian M in (v) = ρ in (2πT in ) 3 2 exp - v -V in 2 2T in .
On the top part of the boundary, we consider a diffuse reflection (2.5) of particles, with a Maxwellian distribution function

M Γt (v) = exp - v 2 2T t .
In the numerical experiments we assume

ρ in = 1, T in = 1, T t = 1.05
and consider the inflow velocity in the form

V in = Mach in γT in ⎛ ⎜ ⎝ 1 0 0 ⎞ ⎟ ⎠ ,
where Mach in = 5 and γ = 1.4.

To start the calculation we take an uniform initial solution equal to the values defined by the left boundary condition:

f 0 (x, v) = ρ in (2πT in ) 3 2 exp - v -V in 2 2T in , x ∈ Ω x , v ∈ R 3 .
We define the Mach number from the macroscopic quantities, computing the moments of the distribution function with respect to v ∈ R 3 , by

Mach = u √ γT ,
where c ∶= √ γT is the sound speed. We apply our inverse Lax-Wendroff method to the boundary conditions at ∂Ω x . More precisely, we extrapolate first the outflow at ghost points corresponding to the four straight pieces. Then we impose directly the inflow at the boundaries Γ l and Γ r by (4.1), (4.2), since they are independent of outflow. While the inflow of Γ b and Γ t is computed by specular and diffuse reflection. Finally we use inverse Lax-Wendroff procedure to compute inflow at ghost points. The significant difference between these two case can be observed in Mach number. In the case Kn = 5, the Mach number reaches its maximum at Γ r while in the case Kn = 0.05 its maximum is at Γ l . We can observe also in the case Kn = 0.05 that there is a clear maximum of the density in the middle of the domain. In the same region the temperature reaches its maximum.

4.4. High-speed flow around an object. In this section, we desire to simulate high-speed airflow around a half airfoil (see Figure 12). The boundary is separated by four parts On the right and left hand sides Γ r , Γ l , we use the same incoming flux (4.1), (4.2). On the top part Γ t , the incoming flow is given by the initial value

∂Ω x = Γ l ∪ Γ b ∪ Γ r ∪ Γ t .
f (x, v) = ρ in (2πT in ) 3 2 exp - v -V in 2 2T in , x ∈ Γ t , v ∈ R 3 .
Finally at the bottom Γ b , we use a purely specular reflection boundary condition. The parameters ρ in , T in , T t , γ, ν have the same values as in the previous test. We use again (4.3) as the initial solution.

We note that on the profile of airfoil we cannot use the neighbor points to approximate the tangential derivative ∂ f ∂ ŷ in (2.15). It is because these neighbor points are not on the same straight. Here we approximate the tangential derivative by using the distribution function of interior domain.

In the following tests, we consider only the situations in hydrodynamic regime, i.e. ε = 0.05, for comparing the ones in literature [START_REF] Graebel | Engineering Fluid Mechanics[END_REF][START_REF] Koren | Multigrid and defect correction for the steady Navier-Stokes equations application to aerodynamics[END_REF]. A 150 × 100 mesh is used in domain Ω x . We use a limit velocity domain [-8, 8] 3 with mesh size 48 × 48 × 12. Two different tests of transonic airflow around this half airfoil are considered: Mach in < 1 and Mach in > 1.

We choose first Mach in = 0.85. So we observe in Figures 13-15 that the flow field around the object includes both sub-(Mach < 1) and supersonic (Mach > 1.2) parts. The transonic 2) period begins when first zones of Mach > 1 flow appear around the object. Supersonic flow can decelerate back to subsonic before the trailing edge. In the case Mach in = 1.2 (see in Figures 161718), the zone of Mach > 1 flow increases towards both leading and trailing edges. There is a normal shock created at trailing edge. The flow decelerates over the shock, but remains supersonic. Moreover a normal shock is created ahead of the object, and the only subsonic zone in the flow field is a small area around the object's leading edge.

Conclusion

In this paper we present an accurate method based on Cartesian mesh to deal with complex geometry boundary for kinetic models set in a complex geometry. We desire to reconstruct the distribution function f on some ghost points for computing transport operator. For this we proceed in three steps: first we extrapolate the distribution function f on ghost points for outflow. Then we use the boundary conditions to compute the inflow at the boundary. Finally we implement an inverse Lax-Wendroff procedure to give an accurate approximation

Figure 2 .

 2 Figure 2. Spatially two-dimensional Cartesian mesh. • is interior point, ◾ is ghost point, is the point at the boundary, is the point for extrapolation, the dashed line is the boundary.

Figure 3 .

 3 Figure 3. Smooth solutions: Experimental order of convergence in L 1 norm (1) in the physical domain (2) at the boundary.

4. 2 .

 2 Flow generated by a gradient of temperature. We consider the ES-BGK equation (2.1)-

Figure 4 .

 4 Figure 4. Flow generated by a gradient of temperature: (1) temperature (2) pressure for various Knudsen numbers ε = 0.025, 0.05 and 0.1.

(4. 1 )Figure 5 .

 15 Figure 5. Trapezoidal domain Ω x .

Figure 6 . 5 0Figure 7 . 5 0Figure 8 .

 65758 Figure 6. Density with ε = 5

  .

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Density with ε = 0.05

Figure 12 .Figure 13 .

 1213 Figure 12. Domain including a half airfoil Ω x .

Figure 14 .

 14 Figure 14. Temperature with Mach in = 0.85

Figure 15 .Figure 16 .Figure 17 .Figure 18 .

 15161718 Figure 15. Mach with Mach in = 0.85

  For this, some neighbor points of x p at the boundary are required (See the empty squares in Figure

∂ f ∂ x (x p , v) directly, we approximate the time derivative ∂ f ∂t , tangential derivative ∂ f ∂ ŷ and collision operator Q( f ). Similarly as in spatially 1D case, we compute ∂ f ∂t and Q( f ). It remains to approximate ∂ f ∂ ŷ .

Table 1 .

 1 Smooth solutions: Experimental order of convergence in L 1 norm.

		-5 1.88	2.341 10 -4	4.1
	256 2	1.7829 10 -5 1.91	5.811 10 -5	2.01
	512 2	4.4571 10 -6	2	1.573 10 -5	1.89

The authors are partially supported by the European Research Council ERC Starting Grant 2009, project 239983-NuSiKiMo.

1