Titre : MÉTHODES NUMÉRIQUES D'ORDRE ÉLEVÉ POUR LA MODÉLISATION DE PLASMA LASER

Keywords: laser-plasma, d'ordre élevé, ALE, absorption laser laser plasma, high-order, ALE, laser absorption

Cette thèse présente le développement d'une méthode ALE pour la modélisation de l'interaction laser-plasma. La particularité de cette méthode est l'utilisation d'une étape de projection d'ordre élevé. Cette étape de projection consiste en une interpolation conservative des quantités conservatives du maillage Lagrangien sur un maillage régularisé. Afin d'éviter les oscillations numériques non-physiques, les flux numériques d'ordre élevé sont combinés avec des flux numériques d'ordre moins élevé. Ces flux numériques sont obtenu en considérant les quantités conservatives constantes par morceaux. Cette méthode pour la discrétisation cellule-centrée consiste à préserver les maximums locaux pour la densité, la vitesse et l'énergie interne. Aspects particuliers de la méthode sont appliquées pour la projection la quantité de mouvement pour la discrétisation 'staggered'. Nous l'utilisons ici dans le cadre de la projection sous la forme de la méthode Flux Correction Remapping (FCR). Dans cette thèse le volet applicatif concerne la modélisation de l'interaction d'un laser énergétique avec de plasma et des matériaux microstructures. Un intérêt particulier est porté à la modélisation de l'absorption du laser par une mousse de faible densité. L'absorption se fait à deux échelles spatiales simultanément. Ce modèle d'absorption laser à deux échelles est mis en oeuvre dans le code PALE hydrodynamique. Les simulations numériques de la vitesse de pénétration du laser dans une mousse à faible densité sont en bon accord avec les données expérimentales.

At the first place, I would like to thank my advisers: Richard Liska for patient leading and many relevant suggestions and remarks, Vladimir Tikhonchuk for his comments concerning all viewpoints of the presented topic, Jérôme Breil for careful guiding and advices related to numerical methods and Jiří Limpouch for discussions of physics-related issues and many useful comments.

A substantial part of the presented work was done in collaboration with colleagues from Centre Lasers Intenses et Applications (CELIA) and the Los Alamos National Laboratory (LANL). Philippe Nicolai and Sergey Guskov from CELIA deserve my gratitude for discussions and comments concerning the laser propagation in microstructured targets. Special thanks belong to Mikhail Shashkov for his care, advices and presentation of challenging topics during my internship in LANL. I thank Milan Kuchařík for a particular numerical simulation in his staggered code as well as for his useful remarks and comments. I thank the CHLER's developers team in CELIA for providing me access to their cellcentered code. In addition to joint projects with these people, several scientific results have been achieved together with Pavel Váchal and Oldřich Renner.

I appreciate the goodwill of Czech Technical University in Prague, University of Bordeaux, Los Alamos National Laboratory, and the French Embassy in Czech Republic for providing me the necessary support during my study and research.

Introduction

Numerical methods, i.e. the procedures for solving mathematical problem in general, are successfully used for thousands of years. Probably the oldest known evidence of them is the Babylonian table showing an algorithm to calculate the square root of the number two [START_REF] Štoll | Dějiny fyziky[END_REF]. Common numerical methods, which are taught at technical universities, are dated to the period from 17 th to 19 th century, the period of growth of the classical mechanics.

Their names, e.g. the Newton method, Lagrange polynomial, Gaussian elimination or Euler method, are the clear evidence of that. However, the biggest development of the numerical methods started in 1940s, when the first computers became available. Incredibly fast and free-of-mistakes computers replaced the non-effective human factor in the procedures. The work of people was transferred from the execution to the development of novel methods. Numerical methods had become widely applied in the wide range of science and technology. Applicable range of the numerical methods was further broaden by the invention of the high-order class of the methods. Lower requirements of the computational time and memory of the high-order methods enable their application even for complex systems, such as the plasma produced by high-power laser facilities [START_REF] Moses | Ignition on the National Ignition Facility[END_REF][START_REF] Mima | Inertial fusion experiments and theory[END_REF][START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF]. This work is dedicated to the numerical methods for laser-plasma interaction modeling. More precisely, the thesis deals with the plasma created by the interaction of an intense laser pulse with matter. Specifically, we develop the methods for solving the hydrodynamic equations, considering the plasma as a compressible fluid with additional physical processes like absorption of laser light or heat conduction. The hydrodynamic model compromises between a detailed physical description and a computational efficiency. The efficiency is crucial for the simulations of complex processes like the laser-matter interaction experiments. Applicability and limits of the hydrodynamic model are discussed later.

Laser plasma interaction experiments typically cover a wide range of physical parameters. The characteristic density ranges from a vacuum limit to a few times the solid state density and the temperature scales from thousands to hundred millions degrees Celsius 1. Introduction or Kelvin. Research on such high energy density states of the matter has begun about 15 years before the discovery of a laser itself when a release of a huge amount of energy in the nuclear reactions of fission and later fusion has been achieved in explosions. Further down-scaling of the plasma experiments and a possible control of released energy ware enabled by the discovery of lasers as a high power energy source. The energy densities sufficient for strong material modifications and plasma production were achieved in 1960s.

Since that, many practical applications of laser matter interaction have been proposed.

The most challenging among them would be a release of the thermonuclear energy in a hydrogen target compressed by a powerful lasers. This is the goal for the world's top research laboratories. modeling the laser light absorption in a plasma and in microstructured materials in the hydrodynamic code PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF], whereas I performed mainly the numerical analysis of the remapping methods during the time of my staying at CELIA. That is one reason why this thesis is divided into the two complementary parts. However, thanks to the deep scientific background of my supervisors and collaborators, both topics were discussed at both institutes.

Layout of the thesis

The introduction and overview of the recent progress in the field of the laser-plasma modeling is presented in the second chapter, including the description of a hydrodynamic model, Lagrangian and Arbitrary Lagrangian-Eulerian (ALE) numerical methods, a highorder remap and the laser-matter interaction.

The third chapter is devoted to the description of the high-order remapping methods. The remap represents a conservative interpolation between two computational grids.

After the one-dimensional description of the remapping methods and their high-order extensions, various aspects as the bounds-and symmetry-preservation of the remap are applied in a two-dimensional (2D) geometry both for a scalar and a vector variables. At the end of the chapter, a description of the complete remapping algorithm for the Euler 1.2. Aim of the thesis equations of inviscid compressible flows is presented.

Absorption of the laser in plasma is described in the fourth chapter. Physical description of the electromagnetic radiation propagation in a continuous media is followed by a numerical method which model this process. The ray-tracing method is further generalized for its application in the 2D cylindrically symmetric hydrodynamic code PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF].

Special aspects of the interaction of the laser radiation with microstructures materials such as low-density foams are pointed out, and a new efficient method for modeling of the laser absorption in foams is developed.

The fifth chapter is devoted to the demonstration of the numerical properties of various remapping methods and to the particular simulation of the laser light interaction with the low density foams. These simulation results are compared to the experimental data obtained on the PALS [START_REF] Jungwirth | The Prague Asterix Laser System[END_REF] and GEKKO [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF] laser facilities.

Aim of the thesis

The general aim of the thesis is to contribute to the state-of-art knowledge in the highorder ALE numerical methods, remapping in particular, for the Euler equations and application of the ALE methods for modeling of laser plasma interaction processes.

The numerical diffusion during the necessary remap stage of the ALE method is reduced by developed high-order remapping methods. The symmetry preservation of the methods for vectors limitation is investigated in detail in the case of radial flows. These issues are recently largely discussed by the computational fluid dynamics community.

As a particular application, an interaction of a laser radiation with a low-density foam target is investigated. Experiments with microstructured foams having the mean density of a few mg/cm 3 show a significant improvement of a laser beam quality and its smoothing after propagation through a foam. This smoothing is essential for the Inertial Confinement Fusion (ICF) applications as it permits to achieve a much better implosion quality. The exact physical mechanism of the foam ionization by the laser light is of a microscopic nature and standard hydrodynamic codes, which are considering the foam as a continuous media of an equivalent density, overestimate the speed of the laser penetration in foams.

In this thesis, we propose a new model of foam ionization and a modification of a laser energy absorption method. It accounts the microscopic foam structure within the standard continuous hydrodynamic model. This new model matches the experimental results of the laser-foam interaction. These modifications were achieved by the implementation of a 3D ray-tracing algorithm in the 2D cylindrically symmetric code PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF].

Introduction

State of the art

This chapter summarizes the common background knowledge required for hydrodynamic simulations of a laser-generated plasma. The hydrodynamic description of the plasma is presented and followed by appropriate numerical methods. The recent status of highorder hydrodynamic methods for the particular application is reviewed with a special focus on remapping methods. Finally, we introduce a laser absorption as an additional term to the hydrodynamic model and point out the specificity of the absorption in low-density foams.

Hydrodynamic model of a laser produced plasma

In general, the evolution of charged particles in a plasma is well described by the Vlasov-Fokker-Planck equations. These equations for the distribution functions f e ( r, v, t) and f i ( r, v, t) for electrons and ions in a six-dimensional phase space are completed by the Maxwell equations for the electric and magnetic field evolution. However, the direct solution of this system is computationally expensive and it is possible only for special configurations. Particle in cell (PIC) methods represent an example of numerical methods developed for a detail description of the plasma. PIC methods are widely used for modeling of the high-intensity laser-plasma interaction on short time and spatial scales [START_REF] Klimo | Laser-plasma interaction studies in the context of shock ignition: the regime dominated by parametric instabilities[END_REF], where these methods are more adapted to a collisionless plasma and only low densities can be simulated. On the contrary, due to a high computational cost, the PIC methods are not applicable for direct simulations of the interaction of a long laser pulse with a complex target. Hydrodynamic equations, i.e. conservation laws for the mass, momentum and total energy, represent the moments of the Vlasov-Fokker-Planck equations. These hydrodynamic equations are closed by an expression for the pressure as a function of concerned variables. This expression is called an equation of state. For the plasma considered here, the integration of the Vlasov-Fokker-Planck equations over the velocity space brings hydrodynamic equations for two fluids assuming that the velocity distributions functions of If the plasma magnetization can be neglected and the electron and ion temperatures and densities are close to each other, a one-fluid compressible non-viscous hydrodynamics described by the Euler equations is sufficient to model the important physical phenomena in the system. In difference from the neutral gas dynamics, the plasma hydrodynamics includes additional terms describing a heat conductivity, an absorption model of the laser radiation and appropriate equations of state [START_REF] Kuchařík | Arbitrary Lagrangian-Eulerian (ALE) methods in plasma physics[END_REF]. Further possible extensions are represented by a radiation transport, a transport of super-thermal particles or considering two different temperatures for the electrons and ions.

Euler equations in Lagrangian coordinates

The fundamental set of the hydrodynamic Euler equations for compressible non-viscous fluid in Lagrangian coordinates has the form

dρ dt + ρ div u = 0 (2.1)
ρ d u dt + grad p = 0 (2.2)
ρ d dt + p div u = 0 , (2.3)
where ρ stands for the fluid density, u the velocity, p the pressure, = e/ρ -u 2 /2 the specific internal energy (energy per unit mass) and e the total energy density (energy per unit volume). The material derivative d dt along the infinitesimal fluid element path is defined by

d dt = ∂ ∂t + u • grad ,
any infinitesimal fluid element motion is described by an Ordinary Differential Equation (ODE)

d x dt = u ,
and the system is closed by the equation of state p = p(ρ, ) .

Numerical methods

A particular state of the modeled continuum is given by a set of three variables. For this set, more options are available. In the further text, we denote the mass M , the momentum P = M u and the total energy E as the conservative variables with their volume densities 1 ρ, µ, e. Another set of the fluid density ρ, velocity u and pressure p is called the primitive variables.

Numerical methods

To be able to solve the hydrodynamic model of the plasma numerically, we need to divide the investigated space into smaller parts. This process is called discretization.

Generally, there are two approaches for the model of a continuum fluid. We can divide the investigated area into non-moving (spatially-fixed) parts (cells) and investigate a flow and an evolution of the fluid in these cells. This approach is called Eulerian and is advantageously used in many applications, e.g. for steady state airfoil flows.

In the second approach, the fluid itself is divided into pieces (cells) instead of the space. Now, we are interested in the temporal evolution of these fluid cells along the stream lines, which represent the motion of the fluid pieces. This Lagrangian description is more appropriate for the fluids with a moving boundary (both the boundary of the fluid itself or the boundary between materials). Another advantage of the Lagrangian description is its easy application for a large scale compression or expansion. All these situations are typical for the investigated laser-plasma interaction.

High-order Lagrangian methods

If the computational error in a given norm decreases faster than linearly during the grid refinement 2 , we call the method to be of a high-order. In the Eulerian framework with a static computational grid, a very high-order of accuracy is typically achieved by recent numerical methods [START_REF] Clain | A high-order finite volume method for systems of conservation laws-multi-dim. optimal order detection (MOOD)[END_REF] in the smooth regions, while the accuracy on discontinuities is limited to the first order.

For the hydrodynamic Euler equations in Lagrangian coordinates, two types of the spatial discretization are applied, namely staggered and cell-centered.

In the cell-centered discretization, all variables are interpreted as the mean value over the computational cell. This concept originates from the class of finite volume methods. The hydrodynamic code CHLER [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF] is based on the cell-centered numerical scheme 2. State of the art [START_REF] Maire | A cell-centred Lagrangian scheme for two-dimensional compressible flow problems[END_REF] with the acoustic Godunov solver. Thanks to the piecewise linear reconstruction of the primitive variables with the standard slope limiter [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] using the generalized Riemann problem terminology [START_REF] Falcovitz | Generalized Riemann Problems in Computational Fluid Dynamics[END_REF], a high-order extension of the scheme is proposed in [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF]. The second-order numerical scheme is formulated for an unstructured polygonal two-dimensional computational grid. A cylindrical extension of the numerical scheme is presented in [START_REF] Maire | A high-order cell-centered Lagrangian scheme for compressible fluid flows in two-dimensional cylindrical geometry[END_REF]. For a 3D discretization, examples of the cell-centered Godunov-type methods are presented in [START_REF] Maire | Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics[END_REF][START_REF] Carre | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF].

Another method, based on the discontinuous Galerkin numerical scheme [START_REF] Vilar | Cell-centered discontinous Galerkin discretizations for two-dimensional scalar conservation laws on unstructured grids and for one-dimensional Lagrangian hydrodynamics[END_REF][START_REF] Vilar | A high-order Discontinuous Galerking discretization for solving twodimenzional Lagrangian hydrodynamics[END_REF], is able to achieve the third-order of accuracy. However, the deformation of the computational grid and the evolution of the specific variables are solved separately in this scheme. Due to the separation, the direct application of a computational mesh smoothing, which is necessary for the laser-plasma modeling, does not preserve the third-order of accuracy [START_REF] Velechovský | Numerické metody modelování laserového plazmatu[END_REF].

For the staggered discretization, scalar variables (ρ, p, , e) are represented by their mean values in the cells, whereas vector quantities ( u, x) are located at nodes of a computational grid. PALE (Prague ALE) code [START_REF] Liska | ALE Method for Simulations of Laser-Produced Plasmas[END_REF] using 2D quadrilateral, logically orthogonal computational mesh is based on the computation of nodal forces resulting from pressure gradients. In the staggered codes, these forces, together with the artificial viscosity- [START_REF] Caramana | Formulations of artificial viscosity for muti-dimensional shock wave computations[END_REF] and other forces, are used for the velocity calculation and a subsequent mesh motion. In the PALE code, the compatible total energy-conserving algorithm [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] for Lagrangian hydrodynamics is extended by adding the terms modeling the laser-plasma coupling [START_REF] Kuchařík | Arbitrary Lagrangian-Eulerian (ALE) methods in plasma physics[END_REF][START_REF] Váchal | Rezoning and remapping for ALE simulations in fluid dynamics and plasma physics[END_REF].

With my PhD work I contributed to the development and implementation of laser absorption modules in the PALE code.

Computational grids with curvilinear coordinates are investigated to achieve higher order of accuracy in the two dimensional hydrodynamics [START_REF] Vilar | A high-order Discontinuous Galerking discretization for solving twodimenzional Lagrangian hydrodynamics[END_REF][START_REF] Dobrev | High-order curvilinear finite element methods for lagrangian hydrodynamics[END_REF]. These methods represent possible future of the high-order laser plasma modeling. However, current production codes for the laser-prodiced plasma applications rely mainly on the second-order Lagrangian discretizations.

Arbitrary Lagrangian-Eulerian methods

Arbitrary Lagrangian-Eulerian (ALE) methods combine the Lagrangian approach with a moving computational grid with the Eulerian description. For the purpose of plasmagenerated plasma modeling, these methods attract a widespread interest [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF][START_REF] Margolin | Introduction to "An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF][START_REF] Peery | Multi-material ALE methods in unstructured grids[END_REF][START_REF] Anderson | An arbitrary Lagrangian-Eulerian method with adaptive mesh refinement for the solution of the Euler equations[END_REF][START_REF] Morrell | A cell by cell anisotropic adaptive mesh ALE scheme for the numerical solution of the Euler equations[END_REF][START_REF] Loubere | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF][START_REF] Galera | A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction[END_REF][START_REF] Renner | Environmental conditions in near-wall plasmas generated by impact of energetic particle fluxes[END_REF].

High-order remapping methods

Following the original work [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF], the ALE method consists of the three parts. Several time-steps of the Lagrangian calculation (1) are followed by a smoothing and/or untangling of the computational grid [START_REF] Moses | Ignition on the National Ignition Facility[END_REF]. The second step is called rezoning. Finally, a conservative interpolation of the conservative variables (3) from the old Lagrangian grid to the smoothed one is performed. The last step is called remapping or advection, and it represents the Eulerian part of the method.

Although there is no temporal evolution during the remapping step, an exchange of the conservative quantities can be expressed in terms of fluxes.

The main motivation for the ALE methods is to avoid a degeneration of the computational mesh quality. In the Lagrangian approach, the computational grid is attached to the fluid and has to follow its evolution. It is well adapted for a laminar flow, but a vortex or a shear flow appearing in the fluid are resulting in a severe degradation of the computational mesh. Non-convex, self-intersecting or inverted (negative directional-volume)

cells are examples of the degradation. Smoothing and/or untangling of the computational grid are aiming on avoiding of all such situations. A detailed review of the recent rezoning methods is presented in [START_REF] Váchal | Rezoning and remapping for ALE simulations in fluid dynamics and plasma physics[END_REF].

High-order remapping methods

In the previous section, the remapping procedure was introduced as a last part of the ALE algorithm. For ALE transformation where the connectivity of the Lagrangian and rezoned meshes are the same, the conservation of the total mass c M c , components of momentum c P x c , c P y c and the total energy c E c , is easily enforced by the flux form of the remapping method. For each of the conservative quantities in the Lagrangian cell Q c ∈ {M c , P x c , P y c , E c } and the remapped quantity Qc on the rezoned mesh, the flux form can be written as

Qc = Q c + c ∈C(c) F Q c c , F Q c c = -F Q cc , (2.4)
Here F Q c c stands for the numerical flux through the cell interface (edge) (c , c). The set of all neighboring cells to the cell c is denoted by C(c).

A natural existence of discontinuous solutions for the Euler equations brings the next requirement. A remapping method should avoid the development of new local extrema or even numerical oscillations. This requirement is specific for high-order methods, because standard first order remapping methods produce a sufficient quantity of numerical diffu-2. State of the art sion to avoid these problems. A discrete maximum principle (DMP) is adapted in recent remapping methods. According to this principle, the remapped value is bounded by the extreme values in the surrounding cells of the old mesh. The DMP principle can be applied for the both primitive variables [START_REF] Liska | Synchronized flux corrected remapping for ALE methods[END_REF] and conservative variables [START_REF] Kuchařík | Arbitrary Lagrangian-Eulerian (ALE) methods in plasma physics[END_REF][START_REF] Berndt | Two step hybrid remapping (conservative interpolation) for multimaterial arbitrary Lagrangian-Eulerian methods[END_REF]. It has been proved analytically [START_REF] Clain | A high-order finite volume method for systems of conservation laws-multi-dim. optimal order detection (MOOD)[END_REF] that the strict DMP requirement limits the order of remapping methods to the second order, whereas the higher-order is preserved only in sufficiently smooth areas without local extrema. Modifications of the DMP necessary to avoid the loss of convergence around the extrema are described in [START_REF] Clain | A high-order finite volume method for systems of conservation laws-multi-dim. optimal order detection (MOOD)[END_REF][START_REF] Colella | A limiter for PPM that preserves accuracy at smooth extrema[END_REF].

To meet all these requirements, the following general algorithm for the second-order remapping was described in [START_REF] Margolin | Second-order sign-preserving conservative interpolation (remapping) on general grids[END_REF]. The first step is a calculation of a piecewise linear reconstruction of conservative quantities distributions. Limiters [START_REF] Van Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] are applied to reduce slopes of the reconstruction near discontinuities. The second step is a quadrature of the reconstruction. The quadrature can be performed over the exact intersections of the Lagrangian and rezoned grid or over swept regions [START_REF] Kucharik | An efficient linearity-and-boundpreserving remapping method[END_REF]. The latter option costs much less computational time because the swept regions are defined only by the grid movement, avoiding costly calculation of all cell intersections.

Slope limiters for a piecewise polynomial reconstruction

The aim of the limiters is to reduce a slope [START_REF] Van Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] of a piecewise linear reconstruction and of the higher-order term for a general polynomial reconstruction. A wide range of slope limiters even for a very high-order polynomial reconstruction is available for high-order Euler methods with a static, square computational mesh. The limitation is applied in a hierarchical order starting from the highest order terms of the reconstruction. In the remapping context for the ALE method, one needs to handle non-equidistant and nonorthogonal meshes, whereas the main theoretical background is limited to limiters for a second-order piecewise linear reconstruction.

For the piecewise linear reconstruction of a scalar quantity in the remapping context, the bounds preservation is typically enhanced by the classical Barth-Jespersen slope limiter [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF][START_REF] Barth | Numerical methods for gasdynamic systems on unstructured meshes[END_REF]. A comprehensive review and performance of selected slope limiters for a piecewise parabolic3 reconstruction is provided in the first part of the third chapter.

High-order remapping methods

Slope limiters for vectors

The traditional approach for vectors consists in a definition of bounds for the Cartesian vector components. However, it was recognized that such a definition is not coordinate invariant. In particular, the definition gives non-symmetric bounds (for example, different bounds for the radial velocity vector on an equiangular polar mesh for nodes with the same radius but different angles).

An interesting idea, which was suggested in the context of vector field reconstruction, defines the bounds by using the Vector Image Polygon (VIP) [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] constructed as the convex hull of vectors in the neighboring cells. The approach is further extended to the symmetric remapping context in the second part of the third chapter. The original VIP approach has been already applied to remapping in [START_REF] Luttwak | Vector image polygon (VIP) limiters in ale hydrodynamics[END_REF][START_REF] Luttwak | VIP (vector image polygon) multi-dimensional slope limiters for scalar variables[END_REF].

In the framework of cell-centered Lagrangian discretization and vector field reconstruction inside the cell, it was suggested to use bounds related to the projection of the velocity vector to principal axes of a mesh deformation tensor [START_REF] Maire | Contribution to the numerical modeling of Inertial Confinement Fusion[END_REF][START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | Staggered Lagrangian discretization based on cell-centered Riemann solver and associated hydrodynamics scheme[END_REF], that is, to use directions related to the flow. Inspired by the idea, we applied the projection method for a Flux-Corrected Remapping method. The description of the method for a cell-centered discretization is provided in the last part of the third chapter.

Flux-Corrected Remapping

The original Flux-Corrected Transport (FCT) approach was proposed by Boris and Book [START_REF] Boris | Flux-corrected transport I: SHASTA, a fluid transport algorithm that works[END_REF] and further advanced by Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]. The method is often used to solve advection problems and hyperbolic systems of partial differential equations. However, this method can be applied to more general problems. A comprehensive summary of FCT is presented in the book [START_REF]Flux-Corrected Transport. Principles, Algorithms and Applications[END_REF].

The concept of remapping conservative variables for fluid dynamics was suggested in [START_REF] Liska | Synchronized flux corrected remapping for ALE methods[END_REF][START_REF] Váchal | Sequential flux-corrected remapping for ALE methods[END_REF][START_REF] Liska | Optimization-based synchronized flux-corrected conservative interpolation (remapping) of mass and momentum for arbitrary Lagrangian-Eulerian methods[END_REF], and following these papers we also adopt the term Flux-Corrected Remapping (FCR). A similar method, based on the direct optimization process for a high-order remap of a scalar quantity, in presented in [START_REF] Bochev | Fast optimization-based conservative remap of scalar fields through aggregate mass transfer[END_REF][START_REF] Bochev | Optimization-based remap and transport: A divide and conquer strategy for feature-preserving discretizations[END_REF].

The basic idea of FCT/FCR is to avoid overshoots, undershoots and oscillations, which are typical for the high-order methods. It is achieved by a convex combination of some higher-order fluxes F H with low-order fluxes F L , so that the local bounds are preserved.

State of the art

The FCR flux across the interface of two cells has the form

F FCR = CF H + (1 -C)F L = F L + C F H -F L dF , (2.5)
where dF is referred to as the antidiffusive flux and the correction factor 0 ≤ C ≤ 1 controls the amount of the high-order portion of the flux used. The FCR method finds the highest C, i.e. closest to the high-order fluxes, for which the FCR fluxes preserve local bounds. Using the worst case scenario, the global (mesh-wide) optimization problem is decoupled into a set of local problems, wherein the C is computed.

High-energy laser-matter interaction

Laser intensity of a few times 10 14 W/cm 2 on a target can be maintained for the timescale from a few hundreds picoseconds to a few nanoseconds on the present days high-power laser facilities [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF][START_REF] Jungwirth | The Prague Asterix Laser System[END_REF]. For the investigated targets made of low-atomic number materials like plastic or Aluminum, such a laser intensity is sufficient to rapidly ionize the matter and transform it in a plasma. For simplicity in simulations presented in this thesis, a full ionization of the matter is assumed and an ideal gas equation of state is used. For the considered laser wavelengths of 350 nm or 438 nm, the electron and ion distribution functions are close to the Maxwell distribution, satisfying the requirement for the hydrodynamic model.

A deposition of the laser energy into the target is described through the new term [START_REF] Kuchařík | Arbitrary Lagrangian-Eulerian (ALE) methods in plasma physics[END_REF] in (2.3)

ρ d dt + p div u = -div I . (2.6)
This term, which express the divergence of the laser intensity I, is calculated with a separate routine in a hydrodynamic code, taking into account different mechanisms of the laser absorption.

Propagation and absorption of laser radiation in plasma

According to the fundamental plasma theory, the laser radiation can propagate only in a plasma with the free electron density lower than the critical value. This critical electron plasma density n crit e depends only on the laser wavelength λ and the fundamental physical (2.7)

According to this theory, basic numerical methods are tracking the laser beam propagation in the sub-critical plasma. Once the critical density is reached, some amount of the laser energy is absorbed in the target. In these oversimplified models, the amount of the absorbed energy is an adjustable parameter, typically ranging from 0.3 to 1.0.

More advanced models are based on the geometric optics theory. In this case, the laser beam is split in a set of rays each of them carrying a certain amount of energy, and a ray equation is solved for each ray separately. The total absorption in the ray-tracing algorithm is calculated as a sum of contributions of each ray, which are defined by the local plasma parameters. A general three-dimensional ray-tracing algorithm for a laser plasma is described in [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF]. Taking the diffraction of the laser light into account, the recent article [START_REF] Colaitis | Towards modeling of nonlinear laser-plasma interactions with hydrocodes: The thick-ray approach[END_REF] presents another interesting method. In the thick ray approach, the laser beam is tracked and the evolution of a Gaussian beam wave front curvature along the beam is solved additionally.

Basic mechanism of the laser light energy deposition is the inverse bremsstrahlung absorption. In this case, the oscillation energy of an electron in the harmonic laser electric field is transformed in Coulomb collisions with ions into the chaotic movement thus increasing the plasma temperature. This process in dominant for the considered laser parameters.

Moreover, contrary to the non-collisional mechanisms, the inverse bremsstahlung absorption for the investigated laser intensity depends only on the electron density and temperature, which makes this process suitable for hydrodynamic modeling.

Low-density foam target specifics

Low-density foam targets could improve significantly the Inertial Confinement Fusion (ICF) target design. The original proposal [START_REF] Desselberger | Use of X-ray preheated foam layers to reduce beam structure imprint in laser-driven targets[END_REF] assumes smoothing of radiation inhomogeneities by the heat conductivity in the low-density layer encapsuling the inner target. A symmetric implosion of the inner capsule is achieved thanks to the smoothing the energy flux in the preheated foam. Another approach of laser beam smoothing in a subcritical foam was proposed in [START_REF] Kalal | Thermal smoothing by laser-produced plasma of porous matter[END_REF].

Foams are a complex microstructured material, which cannot be described within a con- Then the laser propagation in the ionized foam is described in the approximation of a homogeneous low-density plasma. However, such a model overestimates the shock wave amplitude and the ionization front propagation speed significantly [START_REF] Kalal | Thermal smoothing by laser-produced plasma of porous matter[END_REF].

The same overestimation of the ionization front speed was observed for hydrodynamic simulations of the foams considered as a homogeneous low-density media. A possibility of modeling the microscopic foam structure is investigated in [START_REF] Kapin | Hydrodynamic simulations of laser interactions with low-density foams[END_REF]. This approach provides an appropriate propagation speed of the ionization front. However, its wider application is limited due to numerical difficulties with the structured mesh evolution and the method can not be generalized to 3D.

Another model, applicable for the 1D hydrodynamic simulations is suggested in [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF]. By considering a microscopic model of the foam homogenization process, a time-dependent local absorption coefficient for the laser light is introduced. It opens a possibility to a quantitative account of the microscopic processes in the macroscopic hydrodynamics. A novel method, generalizing this idea of a time dependent laser absorption coefficient, is presented in the fourth chapter of this thesis.

Remapping methods

As stated in the second chapter, remapping represents a special kind of an interpolation.

To be conservative, this interpolation is written in the flux form (2.4). Moreover, numerical fluxes have the form of a definite integral of a reconstructed function over the overlapping region of the Lagrangian and rezoned meshes.

In one spatial dimension, this integral has a unique form and there are reliable methods to perform a piecewise polynomial reconstruction. Therefore, we start with the description of possible reconstructions for the remapping methods.

1D reconstruction

In the one-dimensional case, we denote the nodal quantity, as the position of the node, by a half-step subscript and the cell-based quantity by an integer in a standard way, so that the cell i is surrounded by the nodes i -1/2, i + 1/2. The positions of the old (Lagrangian) grid nodes are denoted by the set {x i-1/2 }, i ∈ 1 . . . N + 1, where N stands for the total number of cells. The cell center is x i = (x i+1/2 + x i-1/2 )/2. The rezoned grid is indicated by the superscript n, e.g. its nodal position is x n i-1/2 . The mean value u i is the volume density of a conservative quantity over cell. It may represent e.g. the cell density, momentum or total energy u ∈ {ρ, µ, e} for the Euler equation. We define the mean

u i = 1 ∆x i x i+1/2 x i-1/2 u R (x) dx , (3.1)
where ∆x i = x i+1/2 -x i-1/2 is the cell volume (length in 1D) and u R (x) the reconstruction of the unknown function u(x). This function is unknown after each Lagrangian timestep. One of the main remapping tasks is to find the appropriate reconstruction u R (x) of unknown function u(x) from the mean values u i . The task is addressed in detail the 1D analysis presented in this subsection. Knowing the reconstruction u R (x), a remapped (new) cell mean is simply given by

Remapping methods

u n i = 1 ∆x n i x n i+1/2 x n i-1/2 u R (x) dx , (3.2)
where

∆x n i = x n i+1/2 -x n i-1/2
is the volume of the new cell, i.e. the cell after the rezone step. The definite integral on the right hand side can be further split, as indicated in Fig. 3.1. One part corresponds to the mean value of the old cell and the rest is a sum of numerical fluxes (2.4), i.e. quadrature of the appropriate reconstructions over the overlapping regions of the Lagrangian and rezoned mesh.

For a remapping method, a natural requirement is the conservation of the total amount of a conservative quantity in a form

N i=1 u i ∆x i = N i=1 u n i ∆x n i . (3.3)
This requirement is fulfilled using the flux form (2.4) of the remap under the assumption of rigid boundary during the rezone, i.e x 1/2 = x n 1/2 and x N +1/2 = x n N +1/2 . After conservation, another requirement is the preservation of local bounds. We start with a simple definition of the bounds given by the local values on the old mesh

u min i = min {u i-1 , u i , u i+1 } u max i = max {u i-1 , u i , u i+1 } . (3.4)
The remapped mean value (and so the remap) is in bounds if ∀i,

u min i ≤ u n i ≤ u max i .
3.1. 1D reconstruction

Piecewise constant reconstruction

The piecewise constant reconstruction has the form u

R i (x) = u D i for x ∈ [x i-1/2 , x i+1/2 ]
, where u D i is an unknown. According to the definition (3.1), we get

u i = 1 ∆x i x i+1/2 x i-1/2 u D i dx = u D i 1 ∆x i x i+1/2 x i-1/2 dx = u D i so, u D i = u i . (3.5)
This low-order reconstruction preserves bounds for the reconstruction as ∀i,

u R i (x) = u D i = u i for x ∈ [x i-1/2 , x i+1/2 ] u R i (x) ≥ min {u i-1 , u i , u i+1 } , u R i (x) ≤ max {u i-1 , u i , u i+1 } . (3.6)
We denote the associated remapping method using this reconstruction as a donor-cell1

or low-order remapping method. The low-order remapping method preserves the bounds (3.4).

Unlimited piecewise linear reconstruction

For the one spatial dimension, the reconstruction in the cell i, i.e. x ∈ [x i-1/2 , x i+1/2 ] has the general form

u R i (x) = u i + u x i (x -x i ) , (3.7)
originating from the Taylor series around the point x i . The condition (3.2) leads to u i = u i .

The two point centered finite difference can by used to estimate the approximation of the first derivative (slope)

u x i = u i+1 -u i-1 x i+1 -x i-1 . ( 3
mm u x i = minmod u x i , β u i+1 -u i x i+1 -x i , β u i -u i-1 x i -x i-1 , (3.9)
where β ∈ (1, 2) is a parameter (the low values of β lead to stronger limitation, whereas the high values lead to better extrema preservation and convergence on smooth solution).

The function minmod is defined as Barth-Jespersen (BJ) limiter [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] is constructed to preserve the bounds (3.4) by definition. The nodal min and max

u min i-1/2 = min (u i-1 , u i ) u max i-1/2 = max (u i-1 , u i )
are combined with the nodal extrapolated values

u u i,i-1/2 = u u i (x i-1/2 ) = u i -u x i ∆x i /2 (3.11) u u i,i+1/2 = u u i (x i+1/2 ) = u i + u x i ∆x i /2 (3.12)
3.1. 1D reconstruction to get the value of the nodal slope limiter

α i±1/2 i =            min 1, u max i±1/2 -u i u u i,i±1/2 -u i for u u i,i±1/2 -u i > 0 1 for u u i,i±1/2 -u i = 0 min 1, u min i±1/2 -u i u u i,i±1/2 -u i for u u i,i±1/2 -u i < 0 (3.13)
The final cell-based value of α i in the cell i is the minimum of the both nodes values

α i = min(α i-1/2 i , α i+1/2 i
). This widely-used BJ limited reconstruction has the final form

BJ u R i (x) = u i + α i u x i (x -x i ) . ( 3
φ i±1/2 i =            1 ∆u i,i±1/2 ((∆u max i ) 2 + 2 )∆ui,i±1/2+2∆ui,i±1/2 2 ∆u max i (∆u max i ) 2 +2∆u i,i±1/2 2 +∆u max i ∆u i,i±1/2 + 2
for ∆u i,i±1/2 > 0

1 ∆u i,i±1/2 ((∆u min i ) 2 + 2 )∆ui,i±1/2+2∆ui,i±1/2 2 ∆u min i (∆u min i ) 2 +2∆u i,i±1/2 2 +∆u min i ∆u i,i±1/2 + 2 for ∆u i,i±1/2 < 0 , (3.15)
where u min i , u max i are given by (3.4)

∆u max i = u max i -u i ∆u min i = u min i -u i ∆u i,i±1/2 = u u i,i±1/2 -u i .
The small parameter , which is of the order of machine precision, is suggested to avoid the division by zero. Again, the cell-based values φ i are computed as

φ i = min φ i-1/2 i , φ i+1/2 i
giving the final reconstruction

Venk u R i (x) = u i + φ i u x i (x -x i ) . (3.16)
Although there is a plenty of other possibilities for general purpose slope limiters, a particular application of the limiters for the scalar remapping methods is most often limited to the BJ limiter or its smooth extensions. This is mainly due to the limiter simplicity, preservation of the local bounds as well as the second order of accuracy on 3. Remapping methods continuum functions and an easy extension for multidimensional problems.

Unlimited piecewise parabolic reconstruction

The piecewise parabolic reconstruction is a natural extension of the previous piecewise linear reconstruction

u R i (x) = u i + u x i (x -x i ) + 1 2 u xx i (x -x i ) 2 . (3.17)
Computation of the unknown coefficients u i , u x i a u xx i is derived below. A least-square procedure together with the conservativity requirement are applied to obtain analytical formulas.

From the conservation of the mean value in a single cell, the first relation between the unknown coefficients is derived

x i+1/2 x i-1/2 u R i (x) dx = x i+1/2 x i-1/2 u i + u x i (x -x i ) + 1 2 u xx i (x -x i ) 2 dx = u i x x i+1/2 x i-1/2 + u x i 1 2 (x -x i ) 2 x i+1/2 x i-1/2 + 1 2 u xx i 1 3 (x -x i ) 3 x i+1/2 x i-1/2 = u i (x i-1/2 -x i+1/2 ∆x i ) + 1 2 u x i (x i+1/2 -x i ∆x i 2 ) 2 -(x i-1/2 -x i - ∆x i 2 ) 2 + + 1 6 u xx i ∆x i 2 3 + ∆x i 2 3 = u i ∆x i + 1 24 u xx i ∆x 3 i = u i ∆x i u i = u i - 1 24 u xx i ∆x 2 i . (3.18)
The remaining coefficients are computed by the minimization of the least-square deviation functional φ(u i , u x i , u xx i ) in the neighboring cells

φ(u i , u x i , u xx i ) = j∈{i-1,i+1}    u j - 1 ∆x j x j+1/2 x j-1/2 u R i (x) dx    2 .

1D reconstruction

Taking (3.17) into account, the deviation functional φ(u i , u x i , u xx i ) becomes

u i-1 -u i -u x i (∆x i /2) 2 -(∆x i-1 + ∆x i /2) 2 2∆x i-1 -u xx i -(∆x i /2) 3 + (∆x i-1 + ∆x i /2) 3 6∆x i-1 2 + + u i+1 -u i -u x i (∆x i+1 + ∆x i /2) 2 -(∆x i /2) 2 2∆x i+1 -u xx i (∆x i+1 + ∆x i /2) 3 -(∆x i /2) 3 6∆x i+1 2 .
In the next step, we replace u i by (3.18)

u i-1 -u i + 1 2 u x i (∆x i-1 + ∆x i ) - 1 12 u xx i (∆x i-1 + ∆x i )(2∆x i-1 + ∆x i ) 2 + + u i+1 -u i - 1 2 u x i (∆x i+1 + ∆x i ) - 1 12 u xx i (∆x i+1 + ∆x i )(2∆x i+1 + ∆x i ) 2 .
To simplify the notation, we further introduce

∆x i,i±1 = ∆x i±1 + ∆x i ∆ 2 x i,i±1 = 2∆x i±1 + ∆x i (3.19)
to get the final formula for the deviation functional

φ(u x i , u xx i ) = u i-1 -u i + 1 2 u x i ∆x i,i-1 - 1 12 u xx i ∆x i,i-1 ∆ 2 x i,i-1 2 + + u i+1 -u i - 1 2 u x i ∆x i,i+1 - 1 12 u xx i ∆x i,i+1 ∆ 2 x i,i+1 2 . (3.20)
To minimize the functional, we express its first derivatives according to the unknown parameters u x i , u xx i and set them equal to zero

0 = ∂φ ∂u x i = u i-1 -u i + 1 2 u x i ∆x i,i-1 - 1 12 u xx i ∆x i,i-1 ∆ 2 x i,i-1 ∆x i,i-1 - -u i+1 -u i - 1 2 u x i ∆x i,i+1 - 1 12 u xx i ∆x i,i+1 ∆ 2 x i,i+1 ∆x i,i+1 0 = ∂φ ∂u xx i = -u i-1 -u i + 1 2 u x i ∆x i,i-1 - 1 12 u xx i ∆x i,i-1 ∆ 2 x i,i-1 1 6 ∆x i,i-1 ∆ 2 x i,i-1 - -u i+1 -u i - 1 2 u x i ∆x i,i+1 - 1 12 u xx i ∆x i,i+1 ∆ 2 x i,i+1 1 6 ∆x i,i+1 ∆ 2 x i,i+1 .
The final formula for the unlimited coefficients is obtained by the solution of this linear 3. Remapping methods system. Here, we introduce a novel labeling u x i = u x i , u xx i = u xx i to distinguish the coefficients of the linear and parabolic reconstructions2 . The unlimited piecewise parabolic reconstruction has the final form

u x i = 2 (u i+1 -u i )∆x i,i-1 ∆ 2 x i,i-1 + (u i -u i-1 )∆x i,i+1 ∆ 2 x i,i+1 ∆x i,i-1 ∆x i,i+1 (∆ 2 x i,i-1 + ∆ 2 x i,i+1 ) (3.21)
u xx i = 12 (u i+1 -u i )∆x i,i-1 -(u i -u i-1 )∆x i,i+1 ∆x i,i-1 ∆x i,i+1 (∆ 2 x i,i-1 + ∆ 2 x i,i+1 ) (3.22) u R i (x) = u i + u x i (x -x i ) + 1 2 u xx i (x -x i ) 2 - 1 12 ∆x 2 i . (3.23)
The last term in the previous formula results from the conservation requirement (3.18).

In the 1D case of the parabolic reconstruction in a cell, the final derivative approximations are calculated using the two neighboring cells. An example of the reconstruction is plotted in Fig. 3.1. 1D reconstruction u x i in a 1D piecewise linear reconstruction can be found

u x i = 2 (u i+1 -u i )∆x i,i+1 + (u i -u i-1 )∆x i,i-1 ∆x 2 i,i-1 + ∆x 2 i,i+1
.

(3.24)

It ii easy to see that the approximations u x i and u x i of the first derivative for the linear and parabolic reconstruction are not equal. The equality occurs only for the special case of the equidistant mesh, where both formulas give the same expression as the central difference approximation (3.8).

Limiters for a piecewise parabolic reconstruction

Keeping the same logic as in the linear case, we present here limiters for the piecewise parabolic reconstruction. Minmod (MM) limiter for the piecewise parabolic reconstruction [START_REF] Yang | A parameter-free generalized moment limiter for high-order methods on unstructured grids[END_REF] is done as sequential application of (3.9) starting with the second derivative

m u xx i = minmod u xx i , β u x i+1 -u x i ∆x i,i+1 /2 , β u x i -u x i-1 ∆x i,i-1 /2 , with β ∈ (1, 2) .
If the limitation of the second derivative is not necessary (i.e. if m u xx i = u xx i ), then we

Remapping methods

set m u x i = u x i , otherwise m u x i = minmod u x i , β u i+1 -u i ∆x i,i+1 /2 , β u i -u i-1 ∆x i,i-1 /2 .
The final formula for the minmod-limited piecewise parabolic reconstruction reconstruction is

u minmod i (x) = u i + m u x i (x -x i ) + 1 2 m u xx i (x -x i ) 2 - 1 12 ∆x 2 i . (3.25)
An example of the limiter behavior is plotted in Fig. 3.3.

Kuzmin-Barth-Jespersen (KBJ) limiter [START_REF] Kuzmin | A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods[END_REF] for the piecewise parabolic reconstruction (Fig. 3.4) is based on the BJ limiter (3.13). The original limiter for the piecewise linear reconstruction can be expressed as a function BJ (u i-1 , u i , u i+1 , u x i , ∆x i ) with the return value α i . An example of the sequential application of the limiter, as it is described bellow,

1D reconstruction

is plotted in Fig. 3.4

α xx i = BJ u x i-1 , u x i , u x i+1 , u xx i , ∆x i α xo i = BJ (u i-1 , u i , u i+1 , u x i , ∆x i ) α x i = max (α xo i , α xx i ) B u xx i = α xx i u xx i (3.26) B u x i = α x i u x i u KBJ i (x) = u i + B u x i (x -x i ) + 1 2 B u xx i (x -x i ) 2 - 1 12 ∆x 2 i .
Nejat limiter [START_REF] Nejat | A Higher-Order Accurate Unstructured Finite Volume Newton-Krylov Algorithm for Inviscid Compressible Flows[END_REF] for the piecewise parabolic reconstruction is expressed as

σ = 1 -tanh(S(φ 0 -φ i )) 2 ,
with the constants φ 0 = 0.8 and S = 20 as in the original thesis [START_REF] Nejat | A Higher-Order Accurate Unstructured Finite Volume Newton-Krylov Algorithm for Inviscid Compressible Flows[END_REF]. The limiter value φ i is given by (3.15), resulting in the final formula

N u xx i = σu xx i N u x i = [(1 -σ)φ i + σ] u x i u Nejat i (x) = u i + N u x i (x -x i ) + 1 2 N u xx i (x -x i ) 2 - 1 12 ∆x 2 i . (3.27)
When used to solve the advection part of a hydrodynamic scheme, the smooth Nejat and the Venkatakrishnan extensions of the original Barth-Jespersen limiter generally show better results, introducing adjustable portion of numerical diffusion. On the contrary, this could be a disadvantage during the remapping stage, as we have no information on how to adjust their parameters. The remapping step does not need to be necessarily repeated after every time-step for the whole mesh, resulting in the different requirements for the limitation.

All the last three limiters applied to a piecewise parabolic reconstruction are suggested mainly to avoid numerical oscillations of the high-order schemes. In the original papers, the preservation of the local bounds is not addressed. Numerical examples of their performance are typically limited to equidistant meshes. In this thesis, we present a comparison of the limiters for a non-equidistant mesh during the remapping stage in the first section of the fifth chapter.

Remapping methods

Piecewise Parabolic Method (PPM) [START_REF] Colella | The piecewise parabolic method (PPM) for gasdynamical simulations[END_REF] is well adapted in many hydrodynamic codes and is extendable to even higher order of accuracy [START_REF] White | A high-order finite volume remapping scheme for nonuniform grids: The piecewise quartic method (PQM)[END_REF]. This limiter satisfies the reconstruction bounds exactly. Another modification [START_REF] Colella | A limiter for PPM that preserves accuracy at smooth extrema[END_REF] have been developed to preserve smooth extrema, which is necessary to maintain its third order of accuracy. Unfortunately, this modification is applicable only for the equidistant meshes.

This method presents another option how to find appropriate reconstruction. Instead of searching second-and first-derivative approximations u xx and u x in the cell center, the parabolic reconstruction can be described by the value u l at the left node x i-1/2 and u r at the right node x i+1/2 of the cell i. A simple relation between these coefficients can by expressed

u x i = u r i -u l i ∆x i u xx i = 12 (u r i + u l i )/2 -u i ∆x 2 i . (3.28)
The nodal approximations are calculated from a higher-order interpolant around the nodes and further modified to get a value bounded by the mean values in the surrounding cells:

u i+1/2 = u i + ∆x i ∆x i + ∆x i+1 (u i+1 -u i ) + 1 ∆x i-1 + ∆x i + ∆x i+1 + ∆x i+1 × × 2∆x i+1 ∆x i ∆x i + ∆x i+1 ∆x i-1 + ∆x i 2∆x i + ∆x i+1 - ∆x i+2 + ∆x i+1 2∆x i+1 + ∆x i (u i+1 -u i )- (3.29) -∆x i ∆x i-1 + ∆x i 2∆x i + ∆x i+1 δu i+1 + ∆x i+1 ∆x i+1 + ∆x i+2 ∆x i + 2∆x i+1 δu i ,
where the average slope δu i of the parabola in the cell i is approximated as

δu i = 0 if (u i+1 -u i )(u i -u i-1 ) < 0 min(|∆x i u x i |, 2|(u i -u i-1 )|, 2|(u i+1 -u i )|) • sign(∆x i u x i ) otherwise, (3.30)
where u x i in the previous expression comes from (3.21). This condition guarantees that u i+1/2 stays in bounds of u i and u i+1 and leads to a steeper representation of the shocks [START_REF] Colella | The piecewise parabolic method (PPM) for gasdynamical simulations[END_REF].

The whole PPM algorithm can be summarized as follows: Initialization of the one-side values at nodes u r i = u l i+1 = u i+1/2 as described above. To reduce the overshoots (the

(u r i -u i )(u i -u l i ) ≤ 0 ,
then we set u r i = u l i = u i to get a constant reconstruction in the cell. The last conditions guarantee the monotonicity of the reconstruction in the cell. If

(u r i -u l i ) u i - 1 2 (u r i + u l i ) > (u r i -u l i ) 2 6 then we set u l i = 3u i -2u r i . Similarly if - (u r i -u l i ) 2 6 > (u r i -u l i ) u i - 1 2 (u r i + u l i )
then u r i = 3u i -2u l i . This gives us the final values. The final reconstruction formula for the piecewise parabolic reconstruction by the PPM method is obtained from (3.23) using the values provided by (3.28).

2D reconstruction of a scalar quantity

The extension of remapping methods to higher (two and three) spatial dimensions is not straightforward. Both the limited reconstruction and quadrature become more complicated and novel approaches and ideas are needed to maintain required properties and computational efficiency. We describe first the reconstruction of a scalar quantity followed by its limitation.

2D reconstruction of a scalar quantity

On the contrary to the one-dimensional case, we do not aim to derive analytical formulas for the piecewise polynomial reconstruction. Therefore, we start directly with the secondorder piecewise quadratic reconstruction [START_REF] Clain | A high-order finite volume method for systems of conservation laws-multi-dim. optimal order detection (MOOD)[END_REF]. Lower-order methods are accessible through nullifying the high-order parts. The mean value (3.1) is extended to the form

u c = 1 V c Ωc u R (y, x) dx dy , (3.31)
where the computational cell c is defined by the polygon3 Ω c with the volume V c .

Unlimited piecewise quadratic reconstruction

Further, deriving formulas for a quadratic reconstruction in a single cell, we omit the cell index c. According to a Taylor expansion, the quadratic reconstruction formula is

u H (x, y) = u -λ 0 + λ 1 (x -x c ) + λ 2 (y -y c ) + λ 3 (x -x c ) 2 + λ 4 (x -x c )(y -y c ) + λ 5 (y -y c ) 2 , (3.32)
where x c and y c (defined below) represent the coordinates of the cell centroid and {λ 0 , . . . , λ 5 } is a set of unknown coefficients of the reconstruction. Similarly to (3.18), the first coefficient is derived from the requirement to preserve the mean value of the conservative quantity over a single cell

λ 0 = λ 3 (x s -x 2 c ) + λ 4 (z c -x c y c ) + λ 5 (y s -y 2 c ) . ( 3
  u k - 1 V k Ω k u H (x, y) dx dy   2 . (3.34)
Here Ω k stands for the area of the cell k with the volume V k and C represents the set of n neighboring cells to the appointed cell. At least 5 neighboring cells are required to get a unique solution for Λ. Here, we choose 8 corner and edge neighbors for quadrilateral and edge neighbors for other polygonal meshes containing at least 5 edge neighbors (e.g.

Voronoi meshes in [START_REF] Loubere | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF]). The minimization process is equivalent to the solution of the following overdetermined linear system in the least square sense

AΛ = B B = (β 1 , . . . , β n ) T A =     α 1,1 • • • α 1,5 . . . . . . α n,1 α n,5     (3.35)
with

β k = u k -u α k,1 = x c k -x c α k,2 = y c k -y c α k,3 = x s k -x s + 2(x 2 c -x c x c k ) α k,4 = z c k -z c + 2x c y c -x c k y c -y c k x c α k,5 = y s k -y s + 2(y 2 c -y c y c k ) .
The right hand side of the system consists of the differences from the mean values. The These boundary integrals take a simple form for the polynomial functions

x c = 1 V Ω x dV = 1 6V q∈P (x 2 q + x q x q-1 + x 2 q-1 )(y q -y q-1 ) (3.36)
y c = 1 V Ω y dV = - 1 6V q∈P (y 2 q + y q y q-1 + y 2 q-1 )(x q -x q-1 ) . (3.37)
Here P is the set of all nodes of the selected cell and we assume (without the loss of generality) an ordering, which allows us to select the previous (q -1) and the next (q + 1)

point in this set in the counter-clockwise direction.

The elements in the remaining sub-matrix in (3.35), which define the second-order terms, can be expressed in the similar form

x s = 1 V Ω x 2 dV = 1 12V q∈P (x q + x q-1 )(x 2 q + x 2 q-1 )(y q -y q-1 ) (3.38)
y s = 1 V Ω y 2 dV = - 1 12V q∈P (y q + y q-1 )(y 2 q + y 2 q-1 )(x q -x q-1 ) (3.39) z c = 1 V Ω xy dV = 1 24V q∈P x 2
q-1 (3y q-1 + y q ) + 2x q-1 x q (y q-1 + y q-1 )+ +x 2 q (y q-1 + 3y q ) (y q -y q-1 ) (3.40)

V = Ω dV = 1 2 q∈P (x q + x q-1 )(y q -y q-1 ) . (3.41)
In the piecewise quadratic case, we use the Singular Value Decomposition (SVD) iterative

Remapping methods

numerical method to solve the least-square system (3.35). An illustrative comparison of the two dimensional piecewise polynomial reconstruction is plotted in Fig. 3.7.

Barth-Jespersen limiter

The Barth-Jespersen limiter for a piecewise linear reconstruction represents the standard limitation during the remapping stage in recent ALE codes. This limiter is used as a benchmark for our methods.

In the case of a linear reconstruction over a convex cell (we still omit the cell index here for simplicity), the extreme values are always achieved at the cell vertices. The reconstruction

u H (x, y) = u + λ 1 (x -x c ) + λ 2 (y -y c )
is therefore, first extrapolated with the unlimited slopes

λ 1 = δu δx unlim , λ 2 = δu δy
unlim for all vertices x n , n ∈ N of the cell (in the same manner as 3.11)

u u n = u u ( x n ) = u + λ 1 (x n -x c ) + λ 2 (y n -y c ) . (3.42)
The local bounds are computed as the min/max of the mean values over the neighboring cells (3 × 3 patch) around the particular cell. The computation of the unlimited slopes has been presented in the previous subsection. The final slopes, which are used for the numerical integration during remap, are then

λ 1 = α δu δx unlim , λ 2 = α δu δy unlim , (3.43)
where the parameter α is calculated according to (3.13), where u u i,i±1/2 is replaced by u u n at each node. Instead of the minimum of the two nodal values α i±1/2 i in 1D, we set α as the minimum of the all nodal values of the 2D computational cell.

This example of the BJ limiter for a piecewise linear reconstruction reveals the first difficulty of a high-order extension. The quadratic reconstructed function can reach extreme value within the computational cell. In the 1D case, the PPM method detects this situation and fix the slope to get a monotone (constant) reconstruction in the cell. In

2D reconstruction of a vector

2D, this process would require using a constrained optimization problem detection. This can be achieved with the Lagrange multiplier method but the solution is not obvious at all. A further motivation to avoid the reconstruction limiting is presented in the next sections.

2D reconstruction of a vector

A straightforward extension of the reconstruction methods to vectors may consist in an application of the methods for each Cartesian component of the vector separately. However, e.g. the polar symmetry of a radial flow can be violated by such an extension [START_REF] Luttwak | Vector image polygon (VIP) limiters in ale hydrodynamics[END_REF].

The radial flows are of special importance for the ICF target design. Preservation of the spherical symmetry in compression of the ICF capsule is crucial to achieve ignition conditions. Therefore, a great care has to be taken to distinguish the physically relevant hydrodynamic instabilities from the numerical instability caused by inappropriate methods used in numerical simulations of the capsule compression. Our particular goal is to design new remapping methods, which are proven to preserve polar symmetry for radial vector fields.

Unlimited piecewise linear reconstruction

Here, we describe a cell-based unlimited piecewise linear reconstruction. An analogical reconstruction on a dual mesh for a staggered discretization is presented in [a2]. The reconstruction of the vector w in the computational cell c has the form the particular gradients can be derived from the mean value of the gradient over the computational cell. After analytical manipulation, a particular gradient component of

w c ( x) = w c + (∇ w) c ( x -x c ) . ( 3

2D reconstruction of a vector

the reconstructed function f ∈ {w x , w y } can be expressed in the form

∂Vc f dy = Vc ∂f ∂x dV ≈ ∂f ∂x Vc dV = δf δx V c δf δx = 1 V c ∂Vc f dy ≈ 1 V c q∈P(c) f q+1 + f q 2 (y q+1 -y q ) = δf δx , (3.46)
where V c is the cell volume (3.41) and the summation (integration path) is done counterclockwise around the cell4 c, as indicated in Fig. 3.9. The derivative approximation with respect to y has the similar form

δf δy = - 1 2V c q∈P(c) (f q+1 + f q )(x q+1 -x q ) . (3.47)
Although the later formula for the gradient approximation is more suitable to analytically prove the symmetry of the calculated gradient than the solution of the two-column submatrix in (3.35), both methods can be used for remap providing comparable results.

Next, we derive analytical constraints for vector limiters, which are necessary to preserve the symmetry. For the particular polar cell indicated in Fig. 3.9 and for a symmetric radial vector field with the magnitude W = | w|, the unlimited slopes are obtained through (3.46) and (3.47) as

δw x δx = 1 r 2 2 -r 2 1 [W (r 2 ) + W (r 1 )](r 2 -r 1 ) + 2[W (r 2 )r 1 -W (r 1 )r 2 ] cos 2 ϕ δw y δy = 1 r 2 2 -r 2 1 [W (r 2 ) -W (r 1 )](r 2 + r 1 ) + 2[W (r 1 )r 2 -W (r 2 )r 1 ] cos 2 ϕ δw x δy = δw y δx = 1 r 2 2 -r 2 1 W (r 2 )r 1 -W (r 1 )r 2 sin 2ϕ (3.48)
The derivative approximations (3.48) in the linear part of the reconstruction can be expressed in a more illustrative way. Multiplying the first terms in curly brackets by

1 = sin 2 ϕ + cos 2 ϕ and introducing the notation ∆W = W (r 2 ) -W (r 1 ), W = W (r 2 )+W (r 1 ) 2 ,
3. Remapping methods ∆r = r 2 -r 1 and r = r 2 +r 1 2 , we get

δw x δx = ∆W ∆r cos 2 ϕ + W r sin 2 ϕ δw y δy = ∆W ∆r sin 2 ϕ + W r cos 2 ϕ δw x δy = δw y δx = ∆W ∆r - W r sin ϕ cos ϕ .
These formulas correspond well to the analytical one, which can be expressed as ∂ ∂r =

∂x

∂ ∂x = cos ϕ ∂ ∂r - sin ϕ r ∂ ∂ϕ → ∂w x ∂x = ∂W ∂r cos 2 ϕ + W r sin 2 ϕ → ∂w y ∂x = ∂W ∂r - W r sin ϕ cos ϕ ∂ ∂y = sin ϕ ∂ ∂r + cos ϕ r ∂ ∂ϕ → ∂w y ∂y = ∂W ∂r sin 2 ϕ + W r cos 2 ϕ → ∂w x ∂y = ∂W ∂r - W r sin ϕ cos ϕ .
To proof the symmetry of the linear part the high-order term (∇ w) c ( x -x c ) of the reconstruction (3.44), we rotate this term clockwise by an angle ϕ by the rotation matrix

R + = cos(ϕ) sin(ϕ) -sin(ϕ) cos(ϕ) . (3.49)
After some algebra, we get

R + (∇ w) c ( x -x c ) = W (r 2 )-W (r 1 ) r 2 -r 1 r cos δ -2 3 r 2 1 +r 1 r 2 +r 2 2 r 1 +r 1 cos ∆ϕ W (r 2 )+W (r 1 ) r 2 +r 1 r sin δ . (3.50)
To make this expression more clear, we introduce local cell coordinates η = η ξ . The first one, η = r cos δ, represents the radial and ξ = r sin δ the angular direction given by the angle δ. The angles δ and ∆ϕ are explained in Fig. 3.9. Finally

R + (∇ w) c ( x -x c ) = (∇ w) η c ( η -η c ) = W (r 2 )-W (r 1 ) r 2 -r 1 0 0 W (r 2 )+W (r 1 ) r 2 +r 1 c η -η c ξ , (3.51)

2D reconstruction of a vector

where η c is a local coordinate of the cell centroid.

The formula (3.51) proves that the unlimited term is symmetric. We started with reconstruction from a general angle ϕ and the result is angle-independent. Symmetry of the reconstruction with respect to the cell axis is obvious e.g. for ϕ = 0. To conclude, the unlimited linear reconstruction preserves the radial symmetry for a radial flow on a equiangular polar mesh. Therefore, we focus on the limiters, whereas the symmetry violation originates.

Limitation of a piecewise linear reconstruction

The limited piecewise linear reconstruction with the unlimited gradient approximation

(∇ w) c (3.
w c ( x) = w c + L c (∇ w) c ( x -x c ) (3.52)
with the limiter matrix

L c = Φ xx c Φ xy c Φ yx c Φ yy c . (3.53)
Setting Φ xy c = Φ yx c = 0 and calculating remaining coefficients independently according to the scalar case is a standard way to limit vectors in many hydrodynamic codes.

To preserve the radial flow symmetry on a polar mesh, we need an aditional condition for the limiter matrix. For the specific cell c with the cell centroid x c on the axis x, the gradient (3.45) has to have only diagonal components for the symmetric radial vector (velocity) field. To explain this, we express each Cartesian component of the linear part of the reconstruction (the linear part (∇ w) c ( x -x c ) is given by the gradient (∇ w) c ) as the third coordinate above the Cartesian space (x, y). Now, the linear part of the reconstruction of the 2D vector is represented by two planes (one for each Cartesian coordinate) in this space. To keep the symmetry, the first plane associated with the ycomponent of the linear part of the velocity reconstruction has to intersect all points on the axis x, e.g. the cell centroid x c , where

(∇ w) c ( x -x c ) = (∇ w) c ( x c -x c ) = 0.
For any other point on the axis x, the reconstructed y-component is anti-symmetric with respect to the axis for symmetric radial fields. On the contrary, the second reconstruction plane associated with the x-component has to be constant in the direction perpendicular to the axis for symmetric radial fields.

These conditions on the non-diagonal elements being equal to zero provide useful con-

Remapping methods

straints for slope limiters. When we apply the limiter (3.53) to the gradient (3.45) (which was originally diagonal) rotated by the matrix R + (3.49) to a general direction, we expect the gradient to remain diagonal after the rotation backwards by the matrix

R -= (R + ) R -L c R + (∇ w) η c = R -Φ xx c Φ xy c Φ yx c Φ yy c R + δw η δη 0 0 δw ξ δξ c . (3.54)
This condition leads to the set of equations

(Φ xx c -Φ yy c ) sin ϕ cos ϕ + (Φ xy c cos 2 ϕ -Φ yx c sin 2 ϕ) = 0 (3.55) (Φ xx c -Φ yy c ) sin ϕ cos ϕ + (Φ yx c cos 2 ϕ -Φ xy c sin 2 ϕ) = 0 (3.56)
with the solution

Φ xy c = Φ yx c = Φ yy c -Φ xx c 2 tan 2ϕ , (3.57)
where ϕ is an angle of rotation. From the necessary condition (3.57) for the symmetry preservation, we immediately see, that the standard limiter form with

Φ xy c = Φ yx c = 0 is symmetric only for Φ xx c = Φ yy c .
To get a sufficient condition for the symmetry preservation, we need a frame-independent 5

2D reconstruction of a vector

Here, we present a few of the frame-invariant definitions. The first method [START_REF] Maire | Contribution to the numerical modeling of Inertial Confinement Fusion[END_REF] is based on a local frame-invariant transformation of the reconstructed vector prior to the limiter application. For example, a rotation into the direction of the vector in the given cell is the frame invariant transformation. We use this transformation because of the close relation between the vector direction and limitation, which is useful to show properties of the limiter, namely the symmetry-and bound-preservation. In this case, the bounds

w ξ,max c = max k∈{C(c),c} R + c w k ξ w η,max c = max k∈{C(c),c} R + c w k η .
Minimization/maximization goes over the set C(c) of all neighboring cells, including the cell c. The rotation is realized thanks to the local rotation matrix

R + c = 1 | w c | w x c w y c -w y c w x c . (3.59)
The same transformation is used to project reconstructed values with the unlimited slopes at the cell vertices (3.42) and at the cell centers following the BJ limiter logic

w u,ξ,η c ( x) = R + c w u c ( x) w ξ,η c = R + c w c ( x) . (3.60)
For each component of w u,ξ,η c ( x), the BJ limiter (3.43) is applied to obtain the limiting coefficients α ξ c and α η c , each of them as in the scalar 2D case. Then, these coefficients are transformed back to the global Cartesian coordinate system. Following this way we obtain the limiter form suitable for the equation (3.52)

L c = R - c α ξ c 0 0 α η c R + c . (3.61)

Remapping methods

The projection of the BJ limiter is simple and it does not introduce any significant computational cost. We clearly recommend an application of this modification whenever the vector limiting is necessary. Figure 3.10 presents an example of the limited vector reconstruction.

However, for a particular application, the limitation in the form (3.61) is not sufficient to preserve a specific bound. For a polar grid and a radial flow, this limiter does not

u (x) u (y) u (x) u (y) u (x) u (y) u (x) u (y) (a) (b) (c) (d)

Vector Image Polygon (VIP) limiter

The Vector Image Polygon (VIP) method [START_REF] Luttwak | Vector image polygon (VIP) limiters in ale hydrodynamics[END_REF] is an example of a more restrictive limiter.

A graphical representation of the VIP limiter in the space of the vector field components The rectangle is bounded by the vector field components in the 3×3 patch around the a given cell c. Now, if we construct a convex hull 7 of these vector field components, the hull always fits the rectangle. Therefore, the VIP limiter is more restrictive than BJ limiter and the slope given by the VIP limiter is always smaller then the slope given by the BJ limiter.

The definition of the bounds is only the first part of the limiting procedure. The second task is the definition of the limitation matrix coefficients (3.53). As it was proven before, the diagonal form of the matrix satisfies the symmetry requirement only for Φ xx c = Φ yy c . Therefore, we investigate the special form of the limiter matrix

L c = φ c I , (3.62)
where I is the identity matrix and φ c is a common scalar value.

To state the VIP limitation more formally, we repeat the limitation procedure for

u nc = (∇ w) c ( x n -x c ) .
This vector is in the limits if it lies in the convex hull (CH) of the set

CH{ w c c , c ∈ C(c)} = P c w c c = w c -w c .
The polygon P c has one to 9 (for quadrilateral mesh) vertices depending on the vector field and the polygon always contains the origin 0. Examples of these polygons are plotted in Fig. 3.12. The zero vector represents the piecewise constant low-order reconstruction w c ( x) = w c . Because of the linear reconstruction inside the cell and because of the scalar limiter φ nc , the achievable reconstructions are located on the line segment ( 0, u nc ). The set of possible limiter values is represented by an intersection of this line segment with the convex hull P c , i.e. it is again a line segment ( 0, u VIP nc ), where u VIP nc is the intersection of the segment ( 0, u nc ) with the P c boundary. The best (least diffusive) limiter value φ c n , given by the bounds in the cell c, is the closest to the u nc and corresponds to the intersection u VIP nc . The construction of the convex hull P c can be avoided by the computation of all intersections of ( 0, u nc ) with all segments ( w cc 1 , w cc 2 ) for all pairs of cells c 1 , c 2 ∈ C(c). To get the best (highest) limiter φ VIP nc (corresponding to the intersection u VIP nc ), the maximal limiter value from the computation of all these intersections is taken.

The limitation based on the piecewise linear reconstruction with the VIP limiter is perfectly symmetric for radial fields on polar meshes. However, there still might be a 3.3. 2D reconstruction of a vector small violation of the lower bound for the radial component, while the preservation of the maximal bound is guaranteed by the definition of the VIP limiter. This undershoot is the main motivation for the suggested correction of vector bounds called Modified VIP (MVIP).

Modified VIP limiter

The Modified VIP limiter is constructed to add a new constraint (see Fig. 3.11 (c))

for the minimal vector magnitude, however, only in situations, where it is suitable. At first, suppose one-dimensional situation with only two vectors in the opposite directions having the same magnitude. In this case, the VIP set of these vectors is a line segment and after a restriction on the minimal magnitude (absolute value in 1D), it would not allow any contribution of the high-order part of piecewise linear reconstruction. However, standard limiters in 1D, e.g. (3.13) are applied for the value itself instead of its magnitude (absolute value), allowing the absolute value of limited reconstruction to be smaller than the absolute value in the neighboring nodes. The similar situation occurs in 2D when e.g. the origin of the coordinate system falls inside the CH of the vectors in the neighboring cell's. In this case, any restriction in on minimal vector magnitude does not make a good sense. Therefore, we choose the following condition for application of the minimal magnitude correction of VIP. If there exist any rotated coordinate system, in which all limiting vectors are located in one quadrant, then these vectors have a small angular discrepancy and the minimal magnitude condition is applied.

To compute the MVIP slope limiter, we extend the previous VIP algorithm. We evaluate the unlimited reconstruction in the cell nodes n and check the bounds for MVIP (instead of VIP (3.64))

∀n ∈ N (c) w lim nc ∈ MVIP c , (3.65)
where

MVIP c = VIP c ∩ R c (0, w min ) for ∀c 1 , c 2 ∈ C(c), w c 1 • w c 2 > 0 VIP c otherwise ,
Here VIP c represents the Convex Hull of cell velocities w c , c ∈ N (c) and R c (0,

w min ) = w, |w| ≥ min k∈{c,C(c)} (| w k |) .
For a particular node n of the cell c, we find the maximal φ VIP nc as described in the previous section. For a vector field with a large angular discrepancy in the 3×3 patch

φ MVIP nc = φ VIP nc for w VIP nc ∈ R c (0, w min ) φ i nc otherwise , (3.66)
where φ i nc is the limiter value given by the intersection of the boundary of R c (0, w min ) (i.e. a circle) with the line segment ( w c , w VIP nc ). The particular reconstruction w VIP nc is now clearly outside the R c (0, w min ) (i.e. inside the circle) and the mean value w c ∈ R c (0, w min ) (i.e. outside the circle). This implies that the intersection exists and is unique, φ i nc ∈ (0, φ VIP nc ) and finally w lim nc ∈ MVIP c . In [a5], we describe the applications of the VIP limiter for the momentum remap in the staggered discretization, derive the MVIP limiter and demonstrate its performance on a set of cyclic remap [START_REF] Margolin | Second-order sign-preserving conservative interpolation (remapping) on general grids[END_REF] tests.

1D remapping methods

In this section, the complete remapping algorithm including the calculation of numerical fluxes is described. Remap is a procedure how from the old integral averages (3.1), which are known, obtain the new integral averages (3.2) on the new rezoned mesh. Density ρ i is an example of the integral average of a conservative quantity. Remap is easier to state in conservative quantities, which is mass in this case. The mass of the cell i in 1D is M i = ρ i ∆x i . For a general distribution integral average, the old value of conservative quantity in the cell i in 1D is Q i = u i ∆x i and the flux form of the remap (2.4) can be written as (see fig. 3.1)

Q n i = Q i + F i+1/2 -F i-1/2 . (3.67)
In the section 3.1, several methods for the reconstruction of the distribution in 1D are described. The second part of the remapping algorithm is the computation of numerical fluxes through cell faces 8 . These fluxes are obtained by the integration of the reconstruction. Assuming a small displacement of the Lagrangian (old) and rezoned grids, the rezoned nodes of the new grid (marked with superscript n) stay inside two neighboring old cells, i.e. x i-3/2 ≤ x n i-1/2 ≤ x i+1/2 , ∀i. In this case, the numerical fluxes F can be computed exactly for polynomial reconstruction. For a particular node motion to the 3.4. 1D remapping methods right, i.e x n i-1/2 ≥ x i-1/2 , we get

F i-1/2 = x n i-1/2 x i-1/2 u R i (x) dx = x n i-1/2 x i-1/2 a i + b i x + c i x 2 dx = (3.68) = x n i-1/2 -x i-1/2 a i + b i x n i-1/2 + x i-1/2 2 + c i 3 x n i-1/2 2 + x n i-1/2 x i-1/2 + x i-1/2 2 ,
where we express the parabolic interpolant in the form a i + b i x + c i x 2 for simplicity. The parabolic term is related to this form by

c i = 1 2 u xx i b i = u x i -u xx i x i a i = u i -u x i x i + 1 2 u xx i x 2 i - 1 12 ∆x 2 i .
For the nodal movement to the left, i.e.

x n i-1/2 ≤ x i-1/2 , the numerical flux is F i-1/2 = x n i-1/2 x i-1/2 u R i-1 (x) dx = x n i-1/2 x i-1/2 a i-1 + b i-1 x + c i-1 x 2 dx = = x n i-1/2 -x i-1/2 × (3.69) × a i-1 + b i-1 x n i-1/2 + x i-1/2 2 + c i-1 3 x n i-1/2 2 + x n i-1/2 x i-1/2 + x i-1/2 2 .
The conservative quantity in the cell i has the form Q i = u i ∆x i . To stay conservative (3.3), the value u n i in a new cell is equal to the sum of the old u i and the difference9 of the numerical fluxes, as indicated in Fig. 3.1. Supposing rigid boundary nodes, i.e.

x 1/2 = x n 1/2 and x N +1/2 = x n N +1/2 , the corresponding boundary numerical fluxes are zero. By a simple rearrangement of the formula, we get the remapped values u n on the new grid

u n i ∆x n i = u i ∆x i + F i+1/2 -F i-1/2 u n i = u i ∆x i ∆x n i + F i+1/2 -F i-1/2 ∆x n i .
(3.70)

Remapping methods

The formula (3.70) together with the numerical fluxes (3.68) and (3.69) provide the final form of a remapping method using a piecewise constant, linear or parabolic reconstruction.

The remapping process can be further extended. Significant properties of the reconstruction are the conservativity and the bound preservation (3.6). Similar properties can be defined for the whole remapping process. The conservativity is then a direct consequence of the numerical fluxes (2.4) form of the remap. Bounds for the remapped function

u n i are defined by values u i-1 , u i and u i+1 (3.4).
The bounds preservation can be enforced by another methods. For example, the so called repair [START_REF] Loubere | A subcell remapping method on staggered polygonal grids for arbitrary-Lagrangian-Eulerian methods[END_REF][START_REF] Shashkov | The repair paradigm and application to conservation laws[END_REF], where the out of bounds quantity is redistributed to the neighboring cells. Another possibility is the FCR method. The method uses a convex combination of a low-(piecewise constant reconstruction based) and high-order (unlimited piecewise linear or parabolic reconstruction based) numerical fluxes to preserve the high order of convergence on a smooth solution and to keep it in bounds near discontinuities. A multidimensional FCR method is described the next section directly for two spatial dimensions.

2D remapping methods for a scalar

An exact integration of the reconstruction over the rezoned mesh displacement is computationally expensive. This is mainly due to the costly calculation of many possible intersections of the mesh cells. In order to avoid the calculation of the intersections of the old and new cells, the swept-area based method [START_REF] Kucharik | An efficient linearity-and-boundpreserving remapping method[END_REF] uses an integration over the areas given by the mesh displacement, as indicated in Fig. 3.13. Thanks to its simplicity, computational efficiency and accuracy preservation, the sweptarea method is widely applied during the remapping step with the piecewise linear BJ 3.5. 2D remapping methods for a scalar limited two dimensional reconstruction. However, a bound-preserving reconstruction does not imply the bound preservation of the remapped quantities in this case. To avoid a non-physical10 repair technique, we prefer application of the FCR (2.5) idea. Contrary to the repair, the FCR method is local and compromises only between the two numerical fluxes obtained by two different methods.

Flux Corrected Remapping

Regardless to a particular numerical flux calculation method, suppose that we know the numerical fluxes F ik from the cell i to the cell k of a particular conservative quantity Q.

The remapped quantity Q has the form (2.4). Or, we can rewrite the formula for the mean of conservative quantity distribution of the u i = Q i /V i and ũi = Qi / Ṽi in the cell 11i with the volume of the Lagrangian cell V i (3.41) and the rezoned cell

Ṽi ũi = 1 Ṽi u i V i + k F ik . (3.71)
According to the particular numerical flux calculation method, the summation goes over all edge-neighboring cells E(i) (swept-based approach) or the node-neighbors N (i) (exact integration).

The FCR bounds are defined for the distribution of the conservative quantity

u min i ≤ ũi ≤ u max i (3.72) u min i = min k∈{i,C(i)} (u k ) u max i = max k∈{i,C(i)} (u k ) C(i) = E(i) ∪ N (i) . ( 3
u max i ≥ ũi u max i Ṽi ≥ ũi Ṽi = u i V i + k F L ik Ũ L i +C ik dF ik .
The anti-diffusive flux dF ik = F H ik -F L ik is given by the low-order flux F L ik using the piecewise constant reconstruction and by the high-order flux F H ik using the piecewise parabolic (or piecewise linear) reconstruction.

u max i Ṽi -Ũ L i ≥ k C ik dF ik u max i Ṽi -Ũ L i ≥ k,dF ik >0 C ik dF ik + k,dF ik <0 C ik dF ik
The key points of the derivation are: the neglecting of the negative contributions k,dF ik <0

C ik dF ik in the sum on the right hand side

u max i Ṽi -Ũ L i ≥ k,dF ik >0
C ik dF ik and the factoring of the coefficients C ik ∈ (0, 1) in front of the sum (both to handle a worst-case FCR scenario) to obtain

C ik ≤ D max i = u max i Ṽi -Ũ L i k,dF ik >0 dF ik . (3.74)
For every cell i, we have now the sufficient condition for its numerical fluxes

C ik ≤ D max i .
The minimal constraint is derived in the very similar way

u min i ≤ ũi u min i Ṽi ≤ ũi Ṽi = u i V i + k F L ik Ũ L i +C ik dF ik 3.5.
2D remapping methods for a scalar

Ũ L i -u min i Ṽi ≥ - k C ik dF ik Ũ L i -u min i Ṽi ≥ - k,dF ik >0 C ik dF ik + k,dF ik <0 C ik dF ik ⇑ Ũ L i -u min i Ṽi ≥ - k,dF ik <0 C ik dF ik (3.75)
giving us the final condition on C ik

C ik ≤ D min i = u min i Ṽi -Ũ L i k,dF ik <0 dF ik . (3.76)
To satisfy all cell-based constraints (3.74) and (3.76) at the interface between two cells i and k, we set

C ik = min D min i , D max i , D min k , D max k . (3.77)
So far, we did not discuss in detail the choice of the low-and high-order numerical fluxes F L ik and F H ik . To get the low-order fluxes, the piecewise constant reconstruction u i (x, y) = u i in the cell i together with the exact integration would be appropriate. To avoid the calculation of cell intersections, we use a different form for the swept-based method

u * ij (x, y) = u i for V ij > 0 u j for V ij < 0 , (3.78)
where u * ij (x, y) represents the constant reconstruction for the particular swept flux F ij = u * ij V i,j through the edge (i, j) and V ij is the oriented volume of the swept area. For both methods with low-order fluxes and a small mesh displacement during the rezone, i.e. that the sum of absolute values of negative volume fluxes does not exceed the original cell volume [START_REF] Vachal | 2D Flux-Corrected Remapping of Momentum with Velocity Bounds Preserved in General Direction[END_REF], the remapped scalar quantity preserves bounds (3.72).

The same options are valid for the high-order numerical fluxes F H ik depending on the integration method. Up to the second-order piecewise quadratic reconstruction (3.32),

Remapping methods

the particular numerical flux has the form

F H ik = S u H dV = (u -λ 0 ) + λ 1 (x S c -x c ) + λ 2 (y S c -y c ) + λ 3 (x S s -2x S c x c + x 2 c )+ + λ 4 (z S c -y S c x c -x S c y c + x c y c ) + λ 5 (y S s -2y S c y c + y 2 c ) V S ,
where the x c , y c are coordinates of the cell centroid and all other terms are given by the integral formulas (3.36 -3.40) over the swept area S corresponding to the edge (i, k), e.g.

x S s = 1 V S S x 2 dV .

2D remapping methods for a vector

Concerning vector remapping, the vector quantities for Euler equations are the velocity and momentum. Vectors are usually treated component by component. For the reconstruction, we showed that this component based treatment preserves the radial symmetry.

The integration of the vector reconstruction is performed for each Cartesian component, introducing no new phenomena. All the methods described for a scalar in the previous section can be used here.

However, there are some differences. The bounds are typically no longer required for the momentum (volume density), but for the rotated velocity u components in the

computational cell i w ξ,min i ≤ ˜ u ξ i ≤ w ξ,max i , w ξ,min / max i = min/max k∈{i,C i } R i u k ξ (3.79) w η,min i ≤ ˜ u η i ≤ w η,max i , w η,min / max i = min/max k∈{i,C i } (R i u k ) η , (3.80)
where the rotation matrix is given by (3.59) with respect to the velocity vector and the bounds are expressed in the local coordinates to preserve the symmetry. Also the FCR has to be performed in a frame-invariant way to preserve the symmetry, similarly to the methods addressed for the vector reconstruction in section 3.3.

Flux Corrected Remapping

Here, we describe the FCR methods for the momentum remap for a special case of the momentum flux given by the previously known mass flux multiplied by a velocity reconstruction at given point. Although the product form of the momentum flux is not necessary to the symmetry preserving FCR remap derivation, we use that form here with 3.6. 2D remapping methods for a vector respect to possible applications. The form is useful e.g. for a particular implementation in the staggered discretization [a2], whereas here, we describe the method for the cell-centered discretization.

The remapped momentum in Cartesian coordinates is

m c w xy c = m c w xy c + s∈S(c) F µ,FCR,xy c,s
, We reorder them and the define available range for the antidiffusive correction of the momentum Q at the lower and upper bound

where
Q ξ,min c = m c w ξ,min c -m c w ξ,L c = m c w ξ,min c -m c w ξ c - s∈S(c) F m c,s R c w L,xy c,s ξ ≤ 0, Q ξ,max c = m c w ξ,max c -m c w ξ,L c = m c w ξ,max c -m c w ξ c - s∈S(c) F m c,s R c w L,xy c,s ξ ≥ 0,
where the inequalities (sign of Q) follow from the fact that the low-order momentum flux preserves the local bounds by construction. Now the constraint for the ξ-component becomes

Q ξ,min c ≤ s∈S(c) C c,s F m c,s   R c w H,xy c,s -w L,xy c,s   ξ ≤ Q ξ,max c (3.84)
and similarly for the η-component

Q η,min c ≤ s∈S(c) C c,s F m c,s   R c w H,xy c,s -w L,xy c,s   η ≤ Q η,max c . (3.85)
Let us denote

dµ ξ c,s = F m c,s   R c w H,xy c,s -w L,xy c,s   ξ , (3.86
)

dµ η c,s = F m c,s   R c w H,xy c,s -w L,xy c,s   η (3.87)
and split the sums according to the sign of the terms dµ ξ c,s , dµ η c,s

Q ξ,min c ≤ s; dµ ξ c,s >0 C c,s dµ ξ c,s + s; dµ ξ c,s <0 C c,s dµ ξ c,s ≤ Q ξ,max c ,
3.6. 2D remapping methods for a vector

Q η,min c ≤ s; dµ η c,s >0 C c,s dµ η c,s + s; dµ η c,s <0 C c,s dµ η c,s ≤ Q η,max c
The sufficient conditions for the bounds preservation, based on the worst case scenario, are

Q ξ,min c ≤ s; dµ ξ c,s <0 C c,s dµ ξ c,s , (3.88a)
Q ξ,max c ≥ s; dµ ξ c,s >0 C c,s dµ ξ c,s , (3.88b)
Q η,min c ≤ s; dµ η c,s >0 C c,s dµ η c,s , (3.88c)
Q η,max c ≥ s; dµ η c,s <0 C c,s dµ η c,s , (3.88d)
where the correction factor C c,s will be kept between 0 (low-order flux) and 1 (high-order flux), so it will not change the sign of the sums. The described approach is symmetric for the radially symmetric velocity field on polar grids and preserves the local velocity bounds in the direction of the flow. However, the preservation of the local velocity bounds in directions ξ and η does not imply the preservation of the bounds in the radial velocity component. If there is a radial velocity field with a constant radial velocity component W ξ > 0, the local bounds use a lower bound on w ξ,min < W ξ , because the projection of the radial velocity from the neighboring cells (located off the axis ξ) to the direction ξ is strictly smaller than W ξ .

Remap of conservative variables for Euler equations

As stated before, an exact integration of the reconstruction over the rezoned mesh displacement is computationally expensive. The second task is the preservation of bounds for density, velocity and specific internal energy. Concerning the conservation during a remapping step, the reconstruction has to be done for the distribution (volume density) of a conservative variable. For the mass, the distribution is density, which is a primitive 3.7. Remap of conservative variables for Euler equations variable. However, this is no more valid neither for the momentum nor for the total energy.

It is not obvious at all how to achieve the bound preservation during the remap concerning the Euler equations by the reconstruction limiting. That is why we focus here on the application of the FCR method to the whole system.

For the Euler hydrodynamic equations, the conservative quantities are the total mass i M i , the total momentum components i P x i , i P y i and the total energy i E i . The conservation is enforced by the flux form of the remapped values

Qi = Q i + k∈C(i) F Q ik . (3.92)
Here Q i ∈ {M i , P x i , P x i , E i } represents the old mean cell quantity and F Q ik the corresponding numerical flux. Other remapped variables (lower-case) are defined as a combination of the conservative ones and the volume V i of the cell i, that is, the mean cell density ρi , velocity ˜ u i , momentum ˜ µ i , total energy ẽi and specific internal energy

˜ i ρi = Mi / Ṽi (3.93) ˜ u i = ˜ µ i /ρ i = ˜ P i / Mi (3.94) ˜ i = ẽi /ρ i -˜ u 2 i /2 = Ẽi / Mi - ˜ P 2 i /(2 M 2 i ) . (3.95)
To be conservative, the reconstruction inside each cell is performed for the distributions of the conservative quantities u i ∈ {ρ i , µ i , e i }.

Flux Corrected Remapping

The application of the Flux Corrected Transport method during the remapping step of the ALE algorithm is reviewed in [START_REF] Liska | Optimization-based synchronized flux-corrected conservative interpolation (remapping) of mass and momentum for arbitrary Lagrangian-Eulerian methods[END_REF]. The remapped vector of the conservative quantities in the cell i has the form

Qi = V i u i + k∈C(i) F Q ik = Q i + k∈C(i) F Q,L ik + C ik (F Q,H ik -F Q,L ik ) dF Q ik = Q L i + k∈C(i) C ik dF Q ik .
(3.96)

The summation goes over the set of all edge neighbors C(i), V i represents the cell volume, Q L i the quantity remapped by a low-order method and dF Q ik the difference between high-

Remapping methods

and low-order numerical fluxes, i.e. so called anti-diffusive flux.

There are different methods to compute the edge-based correction factors

C ik ∈ [0, 1]
for the mass, momentum and total energy remap. In the 1D sequential FCR idea [START_REF] Váchal | Sequential flux-corrected remapping for ALE methods[END_REF], the correction factor for the momentum remap C µ ik is expressed as

C µ ik = C m ik C ik < C m ik i.e
. the correction factor for the mass remap and similarly C E ik = C µ ik C m ik C ik for the energy remap. A local optimization procedure with a constrained minimization problem for each interface (i, k) is performed in the synchronized FCR [START_REF] Liska | Synchronized flux corrected remapping for ALE methods[END_REF][START_REF] Liska | Optimization-based synchronized flux-corrected conservative interpolation (remapping) of mass and momentum for arbitrary Lagrangian-Eulerian methods[END_REF] method to get the optimal value of the independent coefficients (C m ik , C µx ik , C µy ik , C E ik ) for the mass, the Cartesian momentum components and the energy simultaneously.

On the contrary, here we assume the same correction factors

C m ik = C µx ik = C µy ik = C E ik = C ik ,
because of a simple and efficient resolution of all constraints. We require a preservation of the bounds on the remapped density, components of velocity projected in the local flow direction (ξ, η) and specific internal energy

ρ min i ≤ ρi ≤ ρ max i , ρ min / max i = min/max k∈{i,C(i)} (ρ k ) (3.97) w ξ,min i ≤ ˜ u i ξ ≤ w ξ,max i , w ξ,min / max i = min/max k∈{i,C(i)} R i u k ξ (3.98) w η,min i ≤ ˜ u i η ≤ w η,max i , w η,min / max i = min/max k∈{i,C(i)} (R i u k ) η (3.99) min i ≤ ˜ i ≤ max * i , min i = min k∈{i,C(i)} ( k ) (3.100) R i = 1 | u i | ( u i ) x ( u i ) y -( u i ) y ( u i ) x , max * i = max ˜ L i , max k∈{i,C(i)} ( k ) . (3.101)
Because of the form (3.95), even the internal energy ˜ L i remapped by a low-order method can exceed the maximum bound. Therefore, the greater term from the remapped value ˜ L i and the maximum over surrounding cells is taken in (3.101). Thanks to (3.93), the local limiting coefficients for the density are the same as derived for a scalar case (3.74) and (3.76)

C ik ≤ D 1 i = ρ min i Ṽi -M L i k∈C(i) min (0, dF M ik ) (3.102) C ik ≤ D 2 i = ρ max i Ṽi -M L i k∈C(i) max (0, dF M ik ) , (3.103)

Remap of conservative variables for Euler equations

where M L i and dF M ik stand for the low-order remapped mass and the anti-diffusive mass flux.

Concerning the remap of momentum, we need to combine the mass fluxes F M ik and the corresponding remapped mass Mi with the momentum fluxes F P ik . For the minimum bound of the remapped momentum μξ in the ξ-direction and for the ξ-velocity constraint (3.98), we obtain

w ξ,min i ≤ ˜ u i ξ w ξ,min i Mi ≤ P ξ i w ξ,min i   M L i + k∈C(i) C ik dF M ik   ≤ P ξ,L + k∈C(i) C ik dF P ξ ik w ξ,min i M L i -P ξ,L ≤ k∈C(i) C ik dF P ξ ik -w ξ,min i dF M ik .
The sufficient condition for the last equation is

C ik ≤ D 3 i = w ξ,min i M L i -P ξ,L i k∈C(i) min 0, dF P ξ ik -w ξ,min i dF M ik . (3.104)
In a similar way, the next set of constraints is derived for the second inequality (3.98) and the η component (3.99)

C ik ≤ D 4 i = w ξ,max i M L i -P ξ,L i k∈C(i) max 0, dF P ξ ik -w ξ,max i dF M ik (3.105) C ik ≤ D 5 i = w η,min i M L i -P η,L i k∈C(i) min 0, dF P η ik -w η,min i dF M ik (3.106) C ik ≤ D 6 i = w η,max i M L i -P η,L i k∈C(i) max (0, dF P η ik -w η,max i dF M ik ) . ( 3

.107)

These constraints for the momentum are not the same as (3.89) derived for the symmetric FCR vector remap in section 3.6. In the previous method for a momentum vector, the remapped mass was supposed to be known prior to the remap of the momentum. Further, the remapped flux had a special form of the cell mass times the vector quantity. On the

Remapping methods

contrary, here we suppose the same correction factor C ik for the mass and momentum remap, resulting in a different form of the FCR constraints.

Similarly to the velocity, the specific internal energy is not remapped directly. Starting from the required inequality (3.100), we get

min i ≤ ˜ i min i M 2 i ≤ Ẽi Mi - P ξ i 2 + P η i 2 2 .
Further, we use the formula (3.96) for the remapped mass, momentum and total energy

min i M L i + k C ik dF M ik 2 ≤ ẼL + k C ik dF E ik M L + k C ik dF M ik - - 1 2   P ξ,L + k C ik dF P ξ ik 2 + P η,L + k C ik dF P η ik 2  
and after rearrangement, we get

min i ( M L i ) 2 -ẼL i M L i + 1 2 P ξ,L i 2 + P η,L i 2 ≤ ≤ k C ik dF E ik M L i -dF P ξ ik P ξ,L i -dF P η ik P η,L i + dF M ik ẼL i -2 min i mL i dF de,min ik + + l C il dF E ik dF M il - 1 2 dF P ξ ik dF P ξ il - 1 2 dF P η ik dF P η il -min i dF M ik dF M il dF ds,min ikl .
To get a sufficient condition for this inequality, we decrease the negative part (the positive part is neglected later) of the right-hand side by substituting C il = 1 for negative terms in the inner sum and by neglecting the positive terms (in the inner sum).

min i ( M L i ) 2 -ẼL i M L i + 1
C ik ≤ D 7 i = min i ( M L i ) 2 -ẼL i M L i + 1 2 P ξ,L i 2 + P η,L i 2 k∈C(i)
min 0, dF de,min ik + l∈C(i) min 0, dF ds,min ikl .

(3.108) The described FCT-based remapping method for the system of Euler equations is computationally efficient, only a summation of the earlier-calculated terms is required.

C ik ≤ D 8 i = max * i ( M L i ) 2 -ẼL i M L i + 1 2 P ξ,L i 2 + P η,L i 2 k∈C(i) max 0,
Bounds preservation for the density, velocity and specific internal energy is guaranteed by the construction of the constraints. The ideas for the preservation of the radial symmetry are employed in the presented FCR method. The numerical performance of the method is demonstrated on a set of standard hydrodynamic test cases in subsections 5.2.2 -5.2.5. Details of a particular application of our FCR method using piecewise quadratic reconstruction are presented in [a1].

Laser absorption modeling methods

In this chapter, a general description of the laser plasma interaction is supplemented by a presentation of the ray-tracing method for problems of a cylindrical symmetry and a discussion of the low-density foam modeling. These methods are applied to numerical simulations of experiments presented in the next chapter.

We use the CGS (centimeter, gram, second) system of physical units modified for the plasma, where the temperature is expressed in the electron-Volt equivalent instead of Kelvin. In difference from the standard SI units (kilogram, meter, second), the CGS system is widely used in hydrodynamic codes related to plasma physics. A list of all relevant physical constants in given in Tab. [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF] In the final results, we recalculate the energy units to Joule, where 1 J = 1•10 7 erg and the laser intensity to common Watt per centimeter squared.

Laser plasma interaction

We start the description of the modeling of the laser light interaction with the simplest configuration, which is the ideal collisionless plasma. Then, we add the collision term to introduce the inverse-bremsstrahlung absorption, which is the main mechanism of the laser energy dissipation for given experimental conditions.

Laser absorption modeling methods

Propagation in collisionless plasma

In the simplest model, a propagation of the laser radiation with the angular frequency ω in a collisionless plasma is described by the dispersion relation

ω 2 = ω 2 p + k 2 c 2 , (4.1)
with the speed of light c and the wave vector k. The plasma frequency ω p , i.e. the frequency of the electron oscillations in a cold plasma with the electron density n e , mass m e and charge e, is determined by [START_REF] Velechovský | Numerické metody modelování laserového plazmatu[END_REF]

ω p = 4πe 2 n e m e [s -1 ] . (4.2)
For the ω = ω p in (4.1), the wave vector k is equal to zero, defining the critical frequency where r stands for the ray coordinate and d ds is the derivative along this trajectory. According to (4.6), the trajectory is determined by plasma density variations and the angle of incidence of a laser beam on a plasma with respect to the plasma density gradient.

For each ray, the initial spatial position and the direction of the ray on a plasma boundary (represented by the boundary of the computational domain) are given by spatial and angular profiles of the laser beam.

The ray equation (4.6) can be transformed to a system of ordinary differential equations, whereas standard numerical methods such as an explicit Runge-Kutta method can be applied to solve the problem. Let us define a unit vector in the ray direction in the spherical coordinates d r ds = (sin θ cos φ, sin θ sin φ, cos θ) ,

where θ is the polar angle with respect to the vertical axis z and φ is the azimuthal angle in the xy-plane.

Inverse bremsstrahlung absorption

For the investigated plasma parameters, the Coulomb collisions of plasma electrons and ions are the main absorption mechanism of the laser light. The bremsstrahlung is a radiation process induced by an accelerated charged particle. Similarly to the original X-radiation discovered when a linearly accelerated electron is stopped by heavy ions, the electron acceleration in a plasma is due to its Coulomb interaction with the plasma ions. The reverse process of the radiation energy transformation into the kinetic energy of plasma electrons, is called inverse bremsstrahlung.

Starting from the model of a cold plasma, where an electron is placed in a static ion background, we consider its motion in the oscillating laser electric field. By taking into account the deceleration of the electron given by an electron-ion collision frequency ν ei , we can write the equation of the electron motion in the harmonic electric field of laser

E 0 e -iωt d 2 r dt 2 + ν ei d r dt = e m e E 0 e -iωt .
The coherent electron motion in the laser field produces an electric current which modifies 4.1. Laser plasma interaction the laser propagation. The electron collisions contribute to the complex permittivity of the plasma

= 1 - ω 2 p ω 2 + ν 2 ei + i ν ei ω ω 2 p ω 2 + ν 2 ei . (4.9)
Compared to (4.4) the imaginary part of is approximately proportional to the electronion collision frequency ν ei .

The electron-ion collision frequency in an ideal high temperature plasma is described by the Spitzer formula [START_REF] Lebo | Simulation of high-intensity laser-plasma interactions by use of the 2D Lagrangian code ATLANT-HE[END_REF]. For low electron temperatures T e comparable to the Fermi energy E F = 2 2me (3π 2 n e ) 2/3 the electrons are partially degenerate. An approximate expression for the electron-ion collision frequency, which is valid for a larger range of electron temperatures k B T e ≥ E F , reads:

ν ei = 4 3 √ 2πZ 2 e 4 n i ln Λ √ m e (k B T e + E F ) 3/2 [s -1 ] , (4.10)
with the electron density n e = Zn i , the ionization Z, the ion density n i = ρ Amp , the mass density ρ, and the atomic number A. The Coulomb logarithm ln Λ can be expressed in where the absorption coefficient for the inverse bremsstrahlung κ ib has the form

κ ib = 2ω c ( √ ) . (4.12)
4. Laser absorption modeling methods

3D ray-tracing method for 2D cylindrically symmetric codes

The physical relations describing the laser ray propagation and the energy deposition in plasma are derived in the previous section. Here, the issues specific for the energy deposition calculation in a 2D Lagrangian hydrodynamic code are addressed. An overview of a general ray-tracing algorithm is followed by the description of a particular implementation in the cylindrically symmetric hydrodynamic code PALE [START_REF] Liska | ALE Method for Simulations of Laser-Produced Plasmas[END_REF].

A fully 3D ray-tracing scheme [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF] on unstructured grids composed of hexahedra, prisms, pyramids, and tetrahedra for the laser beam evolution is almost uniquely using in Lagrangian codes applied for modeling of a laser-produced plasma. The laser beam is represented by an ensemble of discrete rays propagating in a plasma according to the law of the geometrical-optics. The equations for the ray trajectory are integrated using a first-or a second-order approximation, which is exact for a constant density gradient.

The ray-tracing method is based on the exact solution of the ray equation (4.6) inside each computational cell and the calculation of the ray refraction on the cells boundaries.

The exact solution of the ray trajectory in a cell is a parabola assuming a constant density gradient. Applying the Snell law of refraction, the parabolas are linked together to form the complete ray trajectory inside a computational mesh. The energy deposition along the ray trajectory is calculated according to the equation (4.11). For the 2D Cartesian calculation in the PALE code, the Kaiser's method is directly applicable [START_REF] Velechovský | Numerické metody modelování laserového plazmatu[END_REF].

To initialize a ray-tracing algorithm, the incoming laser beam has to be split into separate rays, whereas each ray is defined by its spatial position and direction on a computational domain boundary.

For the 2D Cartesian case in the PALE code, spatial positions of the rays are initialized at the laser-irradiated mesh boundary, such as a fixed number of rays is equidistantly distributed along each cell boundary. Assuming the target location in the beam focal plane and the beam Rayleigh length much larger than the target size, the direction of each ray can be initialized parallel to the beam axis. For a general target and a single laser beam in the target symmetry axis (e.g. the double-foil target modeling [a4]), we initialize the ray according to the Gaussian beam model in the normal direction with respect to the local beam wavefront. In the latter case, an artificial diffractive was added to each ray in order to avoid non-physical focusing of the beam to a point or a caustic surface (Fig. 4.1) leading to a laser-intensity singularity. For all particular laser-target configurations, the ray-tracing calculation [START_REF] Velechovský | Numerické metody modelování laserového plazmatu[END_REF] can be performed on two-dimensional computational mesh for both Cartesian and cylindrical geometries in the PALE code providing a laser energy absorption on each time step.

Laser energy estimation along the beam axis

However, a special attention is needed for the calculation of the absorbed laser energy along the symmetry axis by 2D ray-tracing algorithms in 2D cylindrically-symmetric codes. This problem is addressed in [START_REF] Van Der Holst | Simulating radiative shock with the CRASH laser package[END_REF] for a ray-tracing algorithm for the 2D Eulerian hydrodynamic code CRASH based on the direct numerical solution of the ODE's (4.8)

in the whole computational domain instead of using the analytical solution inside each computational cell. The direct integration is more suitable for Cartesian rectangular grids because of a simple mapping of the ray trajectory into the regular domain and vice versa, i.e. simple evaluation of the refractive index everywhere inside the domain. The authors of Ref. [START_REF] Van Der Holst | Simulating radiative shock with the CRASH laser package[END_REF] point out that for the laser ray, which is not pointing to the cylindrical axis of symmetry1 , an apparent reflection of the ray calculated with a 3D ray-tracing algorithm appears in the axi-symmetric target plane off the axis whereas the same ray calculated by a 2D ray-tracing always hit the axis.

Although the requirement for the 3D ray-tracing originates from the cylindrical symmetry violation of the laser-target setup, the 3D ray-tracing is useful even for the symmetric configurations of the PALS [START_REF] Jungwirth | The Prague Asterix Laser System[END_REF] laser, which are typically simulated by the 2D cylindrically-symmetric code PALE. Firstly, the artificial diffractive term can be replaced by a random deviation of the laser rays from the wave front normal in order to obtain The hydrodynamic motion remains axi-symmetric. The absorbed energy is averaged over the azimuthal angle so that the overall plasma dynamics is still described by the cylindrically-symmetric model.

3D ray-tracing algorithm

The goal of this subsection is to introduce a computationally-efficient 3D ray-tracing algorithm for a 2D Lagrangian computational mesh supposing a cylindrical symmetry of the hydrodynamic flow. A direct application of the general 3D method [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF] would require an additional discretization in the angular dimension, dramatically increasing the number of computational cells. Instead of this, we assume to have only one computational cell for all azimuthal angles to produce 3D mesh. The supposed 3D cell has a torus-like shape, whereas the original 2D computational cell defines the inner cross section of the torus.

The number of conical faces of the 3D cell is equal to 2D cell's edge number. Supposing are given by the conditions on the cell boundary. Although the piecewise line-trajectory rays are supposed to be inaccurate in the original paper [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF], we adapt a different method [START_REF] Lebo | Simulation of high-intensity laser-plasma interactions by use of the 2D Lagrangian code ATLANT-HE[END_REF] for the computation of the normal of the refractive surface on the cell faces, providing 4.2. 3D ray-tracing method for 2D cylindrically symmetric codes a sufficient precision in the energy deposition. This is illustrated in Fig. 4.2, where the piecewise line method is compared to the direct numerical solution of (4.8). On the cell's boundaries, the rays with the incoming/outgoing direction d in , d out are refracted according to the general Snell law for normalized (unit) vectors [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF]

d out = n in n out n × -n × d in -n 1 - n in n out 2 ( n × d in ) • ( n × d in ) (4.13)
with the refractive index n in in the cell, where the ray comes from and n out in the next cell. Components of the vector of the refractive plane normal n = (n x , n y , n z ) T before normalization are

N x = ∇ R p x p 2 x + p 2 y N y = ∇ R p y p 2 x + p 2 y (4.14) N z = ∇ Z ,
where ∇ = (∇ R , ∇ Z ) T is the refractive index gradient (expressed in the cylindrical coordinates R, Z of the 2D code PALE) in the ray boundary intersection point P = (P R , P Z ) T = ( p 2 x + p 2 y , p z ) T . The normalization is performed in the way to get a correct orientation of the refractive plane

n =    N | N | if N • d in ≤ 0 -N | N | if N • d in > 0 .

Laser-foam interaction

on the interface between the homogenized plasma and a cold non-ionized foam. For the macroscale model, the cell size corresponds to the estimated size of a single pore in the foam.

Time-dependent absorption coefficient in 1D

Here, we briefly describe the model of the laser-absorption in low-density foams, as it was introduced in [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF]. Assuming that no laser radiation is reflected from the low-density foam target, the laser intensity at the plasma boundary is equal to the laser intensity I 0 = I 0 (t) at the given time. Following the notation of the third chapter, the laser intensity at the nodes of the computational mesh for the staggered discretization is defined as

I i+1/2 = I(x i+1/2 ) = I 0 e -κ ef i (t)∆x i , (4.17)
where ∆x i = x i+1/2 -x i-1/2 and the time dependent absorption coefficient κ ef i (t) equals

κ ef i (t) =    2/δ 0 for s e i (t) ≤ δ 0 ρav ρcr α κ ib i for s e i (t) > δ 0 ρav ρcr α . (4.18)
Here, ρ av stands for the average foam mass density and ρ cr for the mass critical density assuming the full ionization

ρ cr = n crit e Am p Z i (4.19)
with the free-electron critical density n crit e given by (4.3). The inverse bremsstrahlung absorption coefficient κ ib i (4.12) is modified with respect to the expanding layer density

κ ib i = ν ei i c ρ e i ρ cr (4.20)
with the effective density

ρ e i = max Am p Z i n e i , δ 0 s e i (t) 1 α ρ av , (4.21)
the electron-ion collision frequency ν ei i (4.10), the proton mass m p , the average mass number A, the average ionization Z i and the free electron density n e i in the cell i. The empirical parameter α < 1 accounts for the fractal structure of the foam. It depends on

Laser absorption modeling methods

the pore size δ 0 , its shape and on the chemical content of the material.

In the model [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], the thickness of a single expanding layer is approximated as

s e i (t) = t n=1 c n i ∆t n (4.22)
with the hydrodynamic time step ∆t n , the recent time index t and the sound speed c n i depending on the current electron temperature T n e i and the initial temperature T 0

c n i = Z i Am p (T n e i -T 0 ). (4.23)
Model parameters corresponding to the given foam are the pore size δ 0 and a foam fractal parameter α. These coefficients are given by the microscopic foam structure according to the scaling law for the solid density ρ s and the characteristic width δ s of the foam element

δ s ≈ δ 0 ρ av ρ s α . (4.24)
Note that α = 1 corresponds to 1D situation with membrane-like structure and α = 1/2 to 2D wire-like foam. The critical expansion factor s cr of a single foam element in the right hand side of the inequalities (4.18) is therefore

ρ cr s 1 α cr ≈ ρ s δ 1 α s s cr ≈ δ s ρ s ρ cr α ≈ δ 0 ρ av ρ s α ρ s ρ cr α ≈ δ 0 ρ av ρ cr α . (4.25)
We have observed that a 1D hydrodynamic simulation employing this method leads to a reduction of the speed of the ionization front in the foam. However, this reduction is not sufficient for the case of the GEKKO experiment [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF].

Microscale foam model and its discretization

We extended the microscopic model [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF] to a more realistic 3D geometry. Here, each computational cell at the hydrodynamic macroscale is supposed to contain one thin layer of the foam. This layer has the initial density of the solid material and it is approximated by a 1D flat layer. The layer is homogeneously heated either by the laser rays, which are entering in the computational cell, or by an electron heat transfer from the surrounding hot plasma. These source terms, i.e. the laser absorption and the heat flux, needed for

n+1 -n ∆t = -p (v n+1 e + v n e ) 2m n - ∆Q c 2ρ c . (4.28)
The source term ∆Q c presents the laser energy absorbed in the cell c. It is calculated according to the formula (4.16) provided by the ray-tracing algorithm and the density ρ c of the macroscale computational cell.

The Spitzer-Harm collisional electron heat conductivity with a flux limiter is used [START_REF] Kuchařík | Arbitrary Lagrangian-Eulerian (ALE) methods in plasma physics[END_REF].

If the internal energy of the macroscale computational cell in the homogenization layer is increased, the microscale cell internal energy is also increased with the same amount 4. Laser absorption modeling methods (4.29).

˜ n+1 = n+1 + max(0, ˜ c -c ) . (4.29)
Here c and ˜ c stands for the macroscale specific internal energy before and after the heat conductivity step. The heat conductivity is calculated on the macroscale level in the PALE code. The increase in the microscale internal energy causes faster expansion of the slab in the following time step. Note that the cell temperature increases as far as the laser is being absorbed. A foam material ionization induced by heat transport starts to effect the computation if the local laser intensity decrease below the heat flux level.

Modifications on the hydrodynamic scale

The foam is modeled as a homogeneous material with the corresponding average density consisting of two states of computational cells. Except of the first layer of cells on the laser-irradiated target side, all other cells are supposed to be non-homogenized at the time t = 0. For the first layer of cells we do not introduce the absorption index modification according to (4.18), instead, we allow the cell to absorb all energy brought with the rays until the cell is expanded sufficiently according to (4.22). This slight modification allows to avoid too strong initial density gradient at the foam surface. Two modifications with respect to the microscale homogenization state of the cells are applied in the ALE hydrodynamic part of the PALE code with a standard heat transfer.

Firstly, the initially cold foam does not respond to pressure gradients and therefore the material stays rigid and does not move before the laser arrives. More precisely, nodal forces responsible for the Lagrangian mesh movement are set to zero for a node with at least one non-homogenized cell attached to the node. Once all the attached cells becomes homogenized, i.e. solid microscale layers associated with the cells expand below the critical density (4.25), we do not apply any restriction for the nodal forces.

Secondly, the foam solid layers are assumed to be opaque and to have a random orientation in the foam. Therefore, the direction of incoming laser rays is randomly altered at the homogenization (ionization) front, i.e. the surface between the cold foam and an ionized homogeneous plasma. These rays have uniform random distribution in the solid angle.

Some rays penetrate the ionization front to be absorbed here. Remaining rays are returned back to homogenized plasma hitting the surface again or escaping from the target.

These modifications take part only in the non-homogenized area of the computational 4.3. Laser-foam interaction grid.

The laser absorption in a homogenized plasma is calculated by the 3D ray-tracing algorithm described above using the bremsstrahlung absorption coefficient (4.12). In the cells of the homogenization layer, the local absorption coefficient κ ib (4. 16) is set to a value high enough to absorb all the laser rays penetrating these cells. Therefore, due to a random reflection of the incoming rays at the irradiated side of the cells, a half of the laser energy, which reached the homogenization layer is absorbed and a half reflected back to the homogenized plasma with the possibility of the future absorption in a different place or by the inverse bremsstrahlung.

Applications

We start this chapter with the numerical tests of the remapping methods, which are discussed in the third chapter. The cyclic remapping and full-hydrodynamic tests are performed to demonstrate the properties of these methods. We focus on the order-of accuracy, bounds-and symmetry-preservation properties of the methods. In the second part of this chapter, applicability of the different absorption methods is investigated for the hydrodynamic modeling of the laser beam interaction with a low-density foam. Finally, we present hydrodynamic simulations of the experiments performed on the PALS [START_REF] Jungwirth | The Prague Asterix Laser System[END_REF] and GEKKO [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF] facilities by taking into account the particular modifications of the hydrodynamic model of the laser-foams interaction, as described in the fourth chapter.

Cyclic remap numerical tests

We start with a cyclic remapping [START_REF] Margolin | Second-order sign-preserving conservative interpolation (remapping) on general grids[END_REF] for a scalar quantity. The cyclic remapping is a numerical test based on a sequence of computational grids. Some analytical function is discretized on the first grid and then remapped to the second one. From the second grid it is again remapped to the next one until the final grid is reached. A comparison of the remapped values with the discretized known analytical function on the final grid provides the error accumulated during all these remapping steps.

1D remapping methods

Here, we provide the comparison of the selected one-dimensional remapping methods described in the third chapter. The order of convergence and the bound-preservation of these methods are numerically demonstrated on cyclic remapping tests.

In this paragraph, we review the notation being used for the particular remapping methods. The remap with the piecewise constant reconstruction is denoted donor, with the piecewise linear reconstruction lin. -in case of unlimited slope described by (3.8) unl., and with the piecewise parabolic reconstruction pp -in case of unlimited coefficients

Numerical error measurement

We assume than the mean u 0 i in all N cells of the initial computational grid x 0 and the grid motion (described later) are known. In n-th remapping pseudo-time step we get the means u n i on the mesh x n and after n t pseudo-timesteps, we get the means u f i on the final grid x f . If the initial and the final grids are the same x 0 = x f (nodes are marked as x i+1/2 and cells are indexed by i according to the section 3.1), then we define the absolute error of the cyclic remap as

E c = n i=1 |u 0 i -u f i |∆x f i ,
where

∆x f i = x f i+1/2 -x f i-1/2
is the volume of the cell i in the grid x f . Generally, if the function u(x) is known and the initial and the final grids x 0 and x f are not the same, the exact values on the final grid can be evaluated as

u e i = 1 ∆x f i x f i+1/2 x f i-1/2 u(x) dx
and we define the overall error of the cyclic1 remap as

E c = n i=1 |u e i -u f i |∆x f i . (5.1)
The bounds definition (3.4) is connected only with the single remapping step. We define

Cyclic remap numerical tests

the out-of-bounds error in n-th pseudo-time step as

E b n = n i=1 max 0, u n i -max u n-1 i ∆x n i + n i=1 max 0, min u n-1 i -u n i ∆x n i ,
where the n-th pseodo-time step bounds are max u n = max{u n i-1 , u n i , u n i+1 } a min u n = min{u n i-1 , u n i , u n i+1 }. Overall out-of-bounds error for the whole cyclic remapping is then the L 1 error in the pseodo-time

E b = nt n=1 E b n ∆t n . (5.2)
Another interesting property of the remap is the preservation of the global bounds of the initial function

E g = n i=1 max 0, u f i -max u ∆x f i + n i=1 max 0, min u -u f i ∆x f i , (5.3)
where the initial bounds are max u = max

i {u 0 i } a min u = min i {u 0 i }.
All the described errors (together with the ratio of E c at different mesh resolutions) for the selected remapping methods are presented in tables 5.1, 5.3, 5.2 and 5.4 for the different initial functions, different computational grid motion and different mesh resolution (meshes with 64-512 cells). Values in the second, third and fourth columns of each table stand for the ratios of remapping errors E c (5.1) during the grid refinement, while the order of convergence (not shown in the tables) are the base-two logarithm of these values. The rest of the columns show the errors E c n 5.1, E b n 5.2 and E g n 5.3 for a mesh with N cells.

Mesh movement

We perform commonly used [START_REF] Margolin | Second-order sign-preserving conservative interpolation (remapping) on general grids[END_REF] computational grid movement, where the grids x n are given by their nodes {x n i-1/2 , i = 1, . . . , n; n = 0, . . . , n t }. The mesh movement (illustrated in Fig. 5.1) is described by the smooth function

x(ξ, t) = (1 -α(t))ξ + α(t)ξ 3 , α(t) = sin(4πt) 2 , 0 ≤ξ ≤ 1 (5.4) 0 ≤t ≤ 1 .
Using this function, the grid sequence x n is given by the grid nodes positions x(ξ i , t n ) where

x n i-1/2 = 5. Applications (a) (b)
t n = n n t , n = 0, . . . , n t ; ξ i = i -1 N , i = 1, . . . , N + 1 . (5.5)
The number of pseudo-time steps n t is set to n t = 5N , where N is the number of cells in given mesh. This condition is related to a standard CFL condition. If we define

CFL = ∆t max i a i-1/2 ∆x i -1 , (5.6)
where the node velocity is

a i-1/2 = dx i-1/2 dt
, then we obtain the maximal CFL number approximately 0.48 for this mesh motion. This condition can be interpreted as the condition on the nodes of the new grid. Each node of the new grid x n+1 has to stay in limits defined by the neighboring cell-centers from the previous grid x n .

Cyclic remap to initial grid of a sinus function

Numerical performance of our methods for the cyclic remap with a different (random) grid movement is demonstrated in [a3]. Here, we use the smooth mesh movement (5. remapped function inside the computational domain. Note that for piecewise parabolic methods we can either reach third order (unl., KBJ and MM) or preserve bounds (pp FCR and PPM), none of the methods is able to reach third order and preserve bounds for functions with a local extrema inside the computational domain. However, the paper [START_REF] Colella | A limiter for PPM that preserves accuracy at smooth extrema[END_REF] might inspire improvements at local extrema.

Sub-cyclic remap of a monotone function

For this test case, the function exp(x) is monotone and has no local extrema in the computational domain x ∈ [0, 1] and thus all remapping methods preserve bounds.

E c 64 E c 512 E b 64 E b 512 E g 64 E g 512 donor 2
.0 2.0 2.0 7.91 • 10 -3 1.01 • 10 -3 0 0 0 0 lin. unl. 3.9 3.9 4.0 1.16 • 10 -5 1.92 • 10 -7 0 0 0 0 lin. BJ 3.8 3.9 4.0 1.13 • 10 -5 1.91 • 10 -7 0 0 0 0 lin. MM 3.9 3.9 4.0 1.16 • 10 -5 1.92 • 10 -7 0 0 0 0 lin. FCR 3.9 3.9 4.0 1.16 • 10 -5 1.92 • 10 -7 0 0 0 0 pp unl. 7.9 7.9 8.0 4.67 • 10 -7 9.32 • 10 -10 0 0 0 0 pp KBJ 7.9 7.9 8.0 4.67 • 10 -7 9.32 • 10 -10 0 0 0 0 pp MM 7.9 7.9 8.0 4.67 • 10 -7 9.32 • 10 -10 0 0 0 0 pp FCR 7.9 7.9 8.

f (x) = 1.4 for x ∈ [1/2, 1].
The convergence is very similar for all methods as expected, i.e. less than first order for the discontinuous function, but the value of errors are different.

The best method is PPM followed by FCR with the piecewise parabolic, respective linear high-order reconstruction. The donor, piecewise linear limited (BJ,MM,FCR) and piecewise parabolic FCR and MM methods preserve bounds as for the non monotone sin test case. The minmod limiter for the piecewise parabolic reconstruction disturbs the bounds much less than the unlimited reconstruction and the method is third-order accurate on a smooth solution. To conclude, we have shown that from the piecewise parabolic methods, only PPM and FCR preserve bounds. A sub-cyclic or a random mesh movement [a3] are required to obtain relevant orders of convergence, i.e. to avoid a super-convergence for piecewise linear methods.

Applications

2D remapping methods for a scalar

In this section, we focus on the preservation of the order of convergence for the two particular numerical flux computation methods from section 3.5, namely: the sweptbased approximate method and the exact integration method with piecewise linear and piecewise quadratic reconstructions. We demonstrate that both methods preserve the third-order of accuracy for the unlimited piecewise quadratic reconstruction. Presented errors of the sub-cyclic remap E c are defined analogically to (5.1).

Pseudo 1D standard advection test analog

As the very first test of this section, we chose an analogue of the standard 1D advection test [START_REF] Vilar | A high-order Discontinuous Galerking discretization for solving twodimenzional Lagrangian hydrodynamics[END_REF] performed on a 2D computational domain, to check the validity of our methods.

The initial condition consists of a Gaussian, rectangular, triangular and parabolic shaped each step the grid is moved by 1/1000 < ∆x = 1/200. This is an analog of an advection test case with a CFL number equal to 0.2, where the grid is static and the solution is moved with a constant velocity. Results in Fig. 5.3 show that the flux corrected remapping method with the piecewise quadratic reconstruction (blue) is superior to other methods.

The piecewise linear method is more affected by the direction of the mesh movement, see e.g. the difference on the top of the parabolic shape. However, due to our definition of the bounds, the piecewise quadratic method is worse than the other 3rd-order smooth extrema-preserving methods for the Eulerian mesh [START_REF] Vilar | A high-order Discontinuous Galerking discretization for solving twodimenzional Lagrangian hydrodynamics[END_REF].

Pseudo 1D test of convergence

The loss of the order of convergence for our method in a presence of a local extrema inside the computational domain is also indicated by the results in Tab. 5.5. We applied the same grid movement as in the previous test with a smooth initial Gaussian function 1 + exp(-(7x -10.5) 2 ). For the third-order method, the error ratio should be equal to 2 3 , compared to 7.02, which is achieved for pp FCR for the finer grids. Similarly lin. BJ is loosing the convergence with the error ratio 3.72 compared to 2 2 . The loss of convergence is due to the strict bounds-preservation requirement. Table 5.5.: Convergence for a 1D smooth problem with a local extrema.

2D test of numerical flux methods

The next 2D remapping test is chosen to verify the order of the remapping method for both, the swept flux [START_REF] Kucharik | An efficient linearity-and-boundpreserving remapping method[END_REF] and the exact intersection-based computation of numerical fluxes.

Note that both methods are equivalent for previous pseudo 1D configurations. We choose the initially equidistant mesh on a rectangular domain [0, 0.5] × [0, 0.5] and we remap a smooth Gaussian function 1+exp (-10[(x -0.25) 2 + (y -0.25) 2 ]). The positions of mesh nodes are given by

x n i-1/2,j-1/2 = x 0 i-1/2,j-1/2 + t sin(2πx 0 i-1/2,j-1/2 ) sin(2πy 0 i-1/2,j-1/2 ) y n i-1/2,j-1/2 = y 0 i-1/2,j-1/2 + t sin(2πx 0 i-1/2,j-1/2 ) sin(2πy 0 i-1/2,j-1/2 ) t = n 10n t
, n = 1, . . . , n t , i, j = 1, . . . , M. The bound-preservation of the FCR method is a direct consequence of the method construction under conditions described in subsection 3.5.1. We do not demonstrate this property for a 2D cyclic remap. However, we checked the bound-preservation in hydrodynamical calculations and our FCR method preserved the bounds in all presented tests during the remap.

2D remapping methods for a vector

In this subsection, we focus on the preservation of the radial symmetry for radial fields of various momentum remapping methods with piecewise linear velocity reconstructions on a dual mesh for a staggered discretization. We describe these methods for the staggered discretization is [a2, a5]. For the cyclic remap on a polar quadrilateral computational domain (r, ϕ) ∈ [0.01, 0.51] × [0, 2π], we apply the grid movement (5.5) for the radial coordinate r(ξ, t), whereas we keep the grid static in the angular direction ϕ

x n i-1/2,j-1/2 = r(ξ i /2 + 0.01, t n ) cos ϕ j

y n i-1/2,j-1/2 = r(ξ i /2 + 0.01, t n ) sin ϕ j , ϕ j = 2πj n j , j = 1, . . . , n j ,
with the n j cells in the angular direction. This grid movement in illustrated in Fig. 5.5. for the unlimited piecewise linear (blue diamonds) and FCRxy (green lines) method (top), respective for the low-order (magenta), unlimited (blue), and symmetric FCRs (red) method (bottom).

Hydrodynamic test cases

Hydrodynamic test cases

We start the section with the 2D Sedov blast wave test case showing the importance of the symmetry-preserving remap of vectors. In this case, the symmetric FCR method with a piecewise linear reconstruction is applied to the momentum remap in a staggered code. The complete FCR with the piecewise quadratic high-order reconstruction for the all conservative quantities is performed in remaining hydrodynamic test cases of this subsection. Remaining tests demonstrate the properties of the complete remapping method from section 3.6.1. These tests were run in the cell-centered hydrodynamic code CHLER [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF] and some of them were presented in [a1]. Two rezone strategies were applied. We refer a mesh smoothing after each Lagrangian time-step as the ALE mode, whereas the rezone back to the original grid after each time-step as the Euler mode. For the first strategy, the particular rezoning method consists of a local minimization problem for the node-based objective functions that define the criteria for grid quality [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF]. Note that the Lagrangian step remains the same when we switch the tested remapping methods.

2D Sedov blast wave on quadrilateral polar mesh

Here, we demonstrate the properties of the symmetry-preserving FCR method (FCRs) applied for the momentum remap for a full hydrodynamic example -the Sedov point explosion test [START_REF] Sedov | Similarity and Dimensional Methods in Mechanics[END_REF]. These results have been already presented in [a2]. We compare the FCRs method to a standard FCR method applied to both momentum components independently (refereed as FCRxy). This FCRxy method may be obtained by setting the local rotation matrix (3.59) to the identity matrix.

The computational domain (r,ϕ) ∈ (0.01, 1.2) × (0, 2π) is covered by an equidistant 100 × 40 polar mesh. At the initial time t = 0, the fluid (ideal gas with γ = 1.4) is static everywhere w = 0 and has an uniform density ρ = 1. The specific internal energy is ε = 10 -8 everywhere except the innermost layer of cells in the center of the domain, where ε = 821.105, which corresponds to the amount of energy in the definition of the Sedov problem presented in [START_REF] Scovazzi | Multi-scale Lagrangian shock hydrodynamics on Q1/P0 finite elements: Theoretical framework and two-dimensional computations[END_REF], transformed to the given polar mesh. The high-energy cells represent a point initial energy generating a circular shock wave spreading from the center. In the final time t = 1, the shock wave reaches the radius r = 1.

This simulation was run in the framework of the staggered research multimaterial ALE (rmALE) code [START_REF] Kucharik | Conservative multi-material remap for staggered discretization[END_REF], employing the compatible discretization of the Euler equations as presented in [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF], and the edge artificial viscosity model from [START_REF] Caramana | Formulations of artificial viscosity for muti-dimensional shock wave computations[END_REF]. For remapping, the

Applications

Applications

flux-based method described in [START_REF] Kucharik | Flux-based approach for conservative remap of multimaterial quantities in 2D arbitrary Lagrangian-Eulerian simulations[END_REF][START_REF] Kucharik | Conservative multi-material remap for staggered discretization[END_REF] is used -all quantities are remapped in a flux form, the remap of fluid momentum is similar to the approach described in section 3.6.1, where the particular application of the method to staggered discretization is described in [a2].

To enhance the effects of remapping, this simulation was run in the Euler mode, where all fluid quantities are remapped back to the initial computational mesh after each Lagrangian step. We present results of four different runs using four velocity reconstruction methods -the low-order, unlimited, component FCRxy in x, y directions, and the symmetric FCRs approach.

Velocity magnitude 2D distributions at the final time moment are shown in Fig. 5.8

for each reconstruction method. In Fig. 5.9, the final velocity distribution in the entire problem is shown as a function of radius. For this particular problem, we can see only a small difference between the low-and high-order approaches at the maximum velocity.

The FCRxy by components causes a severe symmetry violation at the domain center, making this method inapplicable for realistic simulations. After switching to the new symmetric FCR approach, the solution is perfectly symmetric and the velocity profile is very close to the high-order approach.

In Fig. 5.10, we show the final density distribution in the entire problem for the same runs. The differences between the low-and high-order velocity reconstructions is more apparent here. We can observe blurring of the density profile for the component FCRxy approach due to the loss of symmetry in velocity, which propagates to all fluid quantities.

The new method produces a perfectly symmetric solution.

1D Sod shock tube

To demonstrate the bounds preservation of the complete FCR method as it is described in section 3.7, we start with the standard Sod shock tube [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] test shown in Fig. FCR method is slightly more diffusive on the contact discontinuity than the BJ-limited method. This is probably due to fact that the restrictions to numerical fluxes are applied to all fluid quantities simultaneously. The piecewise quadratic FCR method preserves the

Applications

bounds on density, velocity and specific internal energy during the remapping stage of the calculation, whereas the standard BJ limiter does not kill all numerical (dispersion) oscillations.

2D Sod shock tube

The 2D extension of the Sod test addresses the issue of a poor linear approximation on polar meshes by the piecewise linear reconstruction method, i.e. the method with the slope (3.8, 3.24). The computational domain approximating a full circle with a radius 1 is discretized by a polar quadrilateral mesh with 7 cells in the angular direction and 100 cells in the radial one. The initial values of density, velocity and pressure are {1, 0,1} for the cell center radius r < 0.5 and {0.125, 0,0.1} for r > 0.5 with a gas constant γ = 1.4.

The results of the calculations in the pure Lagrangian and Eulerian modes are presented in Fig. 5.12 and in Fig. 5.13 at the final time 0.2.

The piecewise linear remapping method suffers from a strong diffusion caused by very low angular mesh resolution, whereas the piecewise quadratic FCR method is close to the Lagrangian (reference) solution. Note that for the 2D polar Sod test case, the reference density between the tail of the rarefaction wave and the contact discontinuity as well as between the contact and the shock is no longer constant.

-0. 5. Applications

2D Sedov blast wave on polygonal meshes

To demonstrate the applicability of the remapping method for more complex grids, we run the Sedov blast wave similar to conditions described in section 5. [START_REF] Kamm | On efficient generation of numerically robust sedov solutions[END_REF]) in a L1 norm are given Results of the remapping method with a piecewise constant (i.e. the donor method), the FCR with a piecewise linear and the FCR with a piecewise quadratic reconstructions in the ALE mode with a standard rezoning [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF] for an unstructured grid are presented in Fig. 5.14. The figure shows a better preservation of the polar symmetry and a reduction of the mesh imprint by the higher-order reconstruction during the remapping stage.

2D Noh implosion

The better preservation of the symmetry is further demonstrated by results of the Noh implosion test shown in Fig. 5.15 at the final time 0.6. We run the test in the Euler mode.

The computation is initialized with the constant initial density 1 and the radial velocity . The mesh imprint in Fig. 5.15 is reduced and the shock is better captured by the higher-order methods.

The last two hydrodynamic tests demonstrate an advantage of the high-order piecewise quadratic reconstruction being used in the remapping method. An application of the high-order reconstruction improves both the accuracy and the symmetry of the results.

The BJ-limited remapping method with a piecewise linear reconstruction for the density, components of momentum and total energy (without any repair method) fails to pass both these tests (Sedov and Noh) because the negative internal energy is generated.

Our FCR method using the piecewise quadratic reconstruction for the system of Eulerian hydrodynamic equations is bounds-preserving by construction of the constraints.

No repair techniques are needed to maintain the positivity of the internal energy. The performed numerical tests verify the preservation of the radial symmetry and demonstrate the benefits of the piecewise quadratic reconstruction. The numerical results of our piecewise quadratic FCR method are comparable or superior to existing methods based on the limited piecewise linear reconstruction.

Applications

Laser-foam interaction modeling

In this section, we compare different methods of the laser absorption modeling with the aim to evaluate their applicability for the laser-foam interaction simulation. We start the modeling for the GEKKO [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF] laser system parameters [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF]. The target is a 500 µm The main goal is to explain the measured time of the foam burn-through in experiments.

The measured time is significantly longer than the standard hydrodynamic simulations predict [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF]. This is explained by a microscopic foam structure that delays the ionization process. We will demonstrate that an improved model of laser energy absorption in a 5. Applications foam introduced in the chapter 4 allows to significantly improve the quality of simulations, moving the predicted burn-through time closer to the experiment.

Laser absorption models Laser absorption at critical surface

In this crude model the ray trajectories are straight lines and the absorption takes place in a over-critical plasma. Therefore, the rays penetrate through the sub-critical plasma without any refraction. Once the ray hits a computational cell with the electron density larger that the critical value (4.3), all the energy carried with this ray is absorbed in this super-critical cell. For this modeling of the sub-critical foam, we artificially set the critical electron density to 3•10 21 cm -3 below the average electron density of a fully ionized foam (because for this particular method, the foam has to have at least critical density in order to absorb any laser radiation). Therefore, in this crude model, all the laser energy is absorbed near the ionization front (just in the one layer of computational cells). This results in a formation of a weak shock, as illustrated in the first row of Fig. 5.17 for the time 700 ps. This is not what one would expect in an experiment.

For the sub-critical foam modeling, such a model of laser energy absorption on the critical surface is clearly not applicable. The laser energy absorption is strongly localized at the ionization front resulting in a steep increase of the plasma temperature, see Fig. 5.17. Finally, the calculated propagation time of the ionization front through the foam does not agree with the experiment, as indicated in Tab. 5.8.

Inverse bremsstrahlung with parallel straight rays

The laser energy deposition along each ray is calculated according to the local plasma parameters by using the relation (4.16) for the inverse bremsstrahlung absorption. The refraction is neglected, and the rays propagate through the plasma along the straight lines, until they are fully absorbed. We refer this model as parallel rays. The result of this simulation is plotted in the second row of Fig. 5.17. for the time 700 ps.

This method shows a better quantitative agreement with the more sophisticated numerical simulation in [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF], where the parallel ray model has been used. However, it was concluded that it overestimates the ionization wave speed. In the second row in Tab. 5.8, we show the very similar results to [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF], indicationg that a direct simulation of the foam as a homogeneous gas target is not applicable.

Laser-foam interaction modeling

Inverse bremsstrahlung with a 3D ray-tracing

The inverse bremsstrahlung absorption is complemented by the 3D ray-tracing according to (4.14). We refer this method as (3D RT) gas.

The ray-tracing method introduces an additional laser beam filamentatiton. However, the filamentation predicted by ray-tracing methods is a numerical artifact, which is due to the finite distance between the rays. It does not capture such important physical phenomena as the laser beam diffraction and the ponderomotive force. With respect to the ionization front propagation speed, this method (see Tab. A non-homogenized cold computational cell of the foam is heated up by the electron heat flux. The incoming energy is deposited in the solid layer thus inducing its expansion.

Any absorption of the laser light in the structured cell (as in the original 1D model [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF])

would result in the even faster propagation of the ionization front. According to the section 4.3.1 as in [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], the thin layer in each structured cell is initiated the solid state density of 1 g/cm 3 . Once the layer expands to the sub-critical density level according to (4.18), the cell is switched to the homogenized plasma state. Then the laser rays can penetrate the cell and deposit the corresponding amount of their energy inside according to the inverse bremsstrahlung 3D ray-tracing model.

This method shows the strongest reduction of the ionization front propagation in the foam (see the fourth row in Tab. 5.8). This reduction is more significant than the one obtained with the original numerical 1D model [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], which not provide a sufficient reduction of the propagation speed compared to the models presented above in this subsection.

The latter model is referred as (3D RT) foam. Numerical simulations for a set of four models of laser absorption were performed. The calculated electron density, temperature and absorbed power are shown in Fig. 5.17 at the time of 700 ps. Table 5.8 presents the calculated time needed to burn through the foam according to each model. We associate the ionization front position with the cell having the temperature of 100 eV and we follow these cells up to the moment when the front reaches the rear side of the foam.

Laser-foam interaction modeling

Comparing to the experimental data [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF], the ionization front propagation speed is still too high for all laser absorption models. In the experiment, the propagation time of about 1.7 ns was observed with a X-ray streak camera.

For the model of parallel rays (Fig. 5.17), our density and temperature profiles match quite reasonably the results of the advanced laser-plasma code [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF] at the time 500 ps.

For the more detailed absorption models with the ray-tracing, we observe an artificial filamentation in the simulated plasma, which affects the plasma hydrodynamics. Our calculated propagation time of 550 ps with the parallel rays model differs from 700 ps calculated in the advanced code. This difference may be due to the difference in the definition of the ionization front position. We use the temperature-defined position in contrast to a simulated X-ray emission in the advanced code. Moreover, our target contains only a foam while the target in [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF] contains also a plastic foil behind the foam. In order to explore the cause of the difference, we have performed an additional simulation with the foil attached to the rear side of the target (presented below) and calculated the foam X-ray emission.

laser absorption method ionization front propagation time critical surface method 800 ps Inverse bremsstrahlung with parallel rays 550 ps 3D ray tracing with a gas target 550 ps 3D ray tracing with a foam model 900 ps However, the model fails to reproduce the experimental data for the foam with even lower-density [START_REF] Kapin | Hydrodynamic simulations of laser interactions with low-density foams[END_REF]. We observe a precompression of the target, which causes a significant reduction of the front propagation speed. This effect was not observed experimentally [START_REF] Khalenkov | Experience of micro-heterogeneous target fabrication to study energy transport in plasma near critical density[END_REF]. This controversy indicates the necessity of the further improvement of the foam laser absorption model. The physical parameters could be the speed of the expansion of the foil solid elements and the thickness of the homogenization layer.

Concerning the hydrodynamic simulation with the original time-dependent absorption coefficient (4.18) 1D model [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], we observe that the plasma density and temperature profiles are similar to the simulation with the critical surface absorption model along the

Applications

symmetry axis (Fig. 5.17 top). In both simulations, a shock wave is generated, thus precompressing the target and modifying its properties before the homogenization occurs.

Thus these models are not describing correctly the real homogenization process in the foam. It should rather be driven by the laser itself or by the electron heat conductivity.

The main advantage of our new multiscale model is that is includes a detailed microscopic simulation of each foam layer expansion thus providing a more realistic microscopic model compatible with the overall macroscopic simulation.

With the aim of a better comparison with the simulation results presented in [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF], we conducted an additional simulation of the laser foam interaction with a 18 µm polystyrene foil attached to the rear side of the foam. A simple model of the lateral X-ray emission of the foam was included in the PALE code. We calculate the contributions to the X-ray emission of the all computational cells under the assumption of a optically thin plasma. The simulated X-ray streak camera images are presented in Fig. 5.18. Table 5.9

shows estimated light propagation time in the foam. The calculated electron density,

Applications

A comparison of the ionization front propagation time through the foam with and without the plastic foil attached to its rear side is shown in Tab. 5.8 and Tab. 5.9.

There are marginal improvements in the propagation time estimation. The front position was modeled by the foam X-ray self-emission. The propagation time is about 100 -200 ps longer. The delay is explained by the fact that the self-emission occurs behind the ionization front once the polystyrene foil is sufficiently heated up. The ionization time calculated with the 3D ray tracing foam model is closer to the experimental data.

Numerical simulations of the PALS and GEKKO experiments

In this section, we present the numerical simulations of the real experiments of laser interaction with foam targets to demonstrate the performance of our proposed multiscale method for modeling of laser absorption in microstructured targets. The third setup presented in [START_REF] Sunahara | Experimental evidence of foam homogenization[END_REF] consists of a 500 µm plastic foam with the mean density 10 mg/cm 3 and 18 µm polystyrene foil attached to the rear side. It was already discussed in the previous section when comparing different models of laser absorption.

The foam with the effective charge/mass number Z/A = 3.85/7.22 was heated up by the 350 nm GEKKO laser radiation. The 300 J pulse has a triangular temporal shape with 400 ps rise and 3 ns down time. The spatial profile is described by a 6th order super-Gaussian function with a 100 µm radius.

5.4. Numerical simulations of the PALS and GEKKO experiments

Results

The numerical simulation is performed with the PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF] The suggested multiscale model reproduces the laser-front propagation as well as the detailed hydrodynamic simulation directly resolving the microscopic foam structure [START_REF] Kapin | Hydrodynamic simulations of laser interactions with low-density foams[END_REF].

One dimensional structural representation of a foam with several cells in each microscopic dense layer and a low-density pore provides a detail description of the homogenization process. However, this model is numerically expensive and cannot be extended to a real experimental case. On the contrary, our multiscale model needs much less computational cells to reproduce a correct homogenization time and it is more robust providing possibility to run simulations in a realistic nanosecond regime.

A smoothing of large scale inhomogenities in the laser beam intensity distribution during its propagation through sub-critical foam layers can be explained by both the laser-ray refraction and the heat conductivity effects. These effects need to be studied in future research.

Summary

In the main part of this thesis, we have presented the high-order remapping method for the system of Euler equations describing the behavior of an inviscid compressible fluid, applicable e.g. for the laser-produced plasma. Issues related to the symmetryand bound-preservation of the method as well as the high-order of accuracy are investigated. The method is tested on standard cyclic-remapping problems. Incorporated in the complete Arbitrary Lagrangian-Eulerian (ALE) method, the performance of the remap is further demonstrated in the cell-centered code CHLER [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ALE) method[END_REF] for the standard hydrodynamic test cases. Our remapping method using a piecewise quadratic reconstruction and a Symmetric Flux-Corrected Transport ideas is directly applicable to unstructured two-dimensional computational meshes with the fixed mesh topology. Vector limiters preserving a radial symmetry of the radial fields on polar meshes are described in detail in the context of the reconstruction limiting during the remapping step of the ALE method.

Symmetry-and bound-preserving modifications of the limiters are numerically tested in the staggered hydrodynamic code PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF].

Particular aspects of the laser interaction with low-density foam targets are addressed in the remaining part of the thesis. Hydrodynamic simulations of the microstructured target materials require a specific treatment of the laser absorption and the material response.

The method of a local time-dependent absorption coefficient [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF] is generalized to be applicable in the two-dimensional cylindrically-symmetric hydrodynamic codes, as e.g.

PALE [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF]. The generalization includes implementation of the 3D ray-tracing algorithm and a sub-cycled simulation of the expansion of solid foam elements induced by the laser and electron energy fluxes. This new model is tested in the numerical simulations of the laser-foam interaction experiments performed at the PALS [START_REF] Jungwirth | The Prague Asterix Laser System[END_REF] and GEKKO [START_REF] Miyanaga | The GEKKO XII -HIPER (High Intensity Plasma Experimental Research) System Relevant to Ignition Target[END_REF] The detailed numerical simulations of the laser-generated plasmas are out of scope of this thesis. Instead of this, we address the key elements of the particular numerical methods required to perform more efficient and accurate simulations.

The main contribution and the new ideas of the thesis can be summarized to:

• Comparative study of the 1D slope limiters for a piecewise parabolic reconstruction in the remapping context.

• Investigation and development of the slope limiters for vectors and their application in the remapping step of the ALE method.

• Development of a new high-order symmetry-preserving remapping algorithm for the Euler equations enforcing the bounds preservation without artificial repair techniques.

• Development and implementation of the laser-absorption methods in the hydrodynamic code PALE.

• Development of a new reliable numerical method for modeling of the laser interaction with low-density foams.

These particular results are achieved in the context of the research programs of the Department of Physical Electronics of the Czech Technical University in Prague and the Center for Intense Lasers and Applications of the Bordeaux University. They contribute to the development of reliable numerical methods for laser plasma modeling.

Conclusions for further research

In the future research, the symmerty-and bound-preserving high-order remapping method could be extended in 3D and to the computational meshes with changing topology. Thanks to the construction of the method, these extensions are rather straightforward. In particular, the application of the remapping method in the Reconnection ALE (ReALE) scheme [START_REF] Loubere | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF], has a potential to improve the ReALE performance.

The slope limiters for vectors are widely applied in non-linear high-order hydrodynamic methods. The described symmetry-and bounds-preserving modifications of the limiters may find application whenever possible outside the remapping methods.

The new laser-absorption model for low-density foams can be applied in future studies of the laser-beam smoothing and hydrodynamic modeling of other microstructured materials. This opens a way for new more efficient numerical models for the future ICF research.

B. Résumé substantiel en français

B.4. Résultats obtenus

La contribution principale et les nouvelles idées de la thèse peuvent se résumer à:

• Etude comparative des limiteurs de pente 1D pour une reconstruction par morceaux parabolique dans le contexte de remappage.

• Recherche et développement des limiteurs de pente pour les vecteurs et leur application dans l'étape de remappage de la méthode ALE.

• Développement d'une nouvelle symétrie-préservation remappage algorithme d'ordre eleve pour les équations d'Euler l'application de la préservation des limites sans les techniques de réparation artificiels.

• Développement et mise en oeuvre des laser les méthodes d'absorption dans le code PALE hydrodynamique.

• Développement d'une nouvelle méthode numérique fiable pour la modélisation de l'interaction laser avec de faibles mousses de densité.

Ces résultats sont obtenus notamment dans le cadre des programmes du ministère de la Electronics physiques de l'Université technique tchèque de Prague de recherche et le

Center for Lasers et Applications de l'Université de Bordeaux intenses. Ils contribuent au développement de numérique fiable méthodes de modélisation laser-plasma.

B.5. Des recherches supplémentaires

Dans la recherche future, la symmerty et lié haute préserver afin remappage méthode pourrait être étendue en 3D et aux mailles de calcul avec topologie changeante. Merci à la construction de la méthode, ces extensions sont plutôt simples. En particulier, l'application de la méthode de reconfiguration dans le ALE de reconnexion (ReALE) régime, a le potentiel d'améliorer la performance de ReALE.

Les limiteurs de pente pour les vecteurs sont largement appliquées dans non-linéaire méthodes hydrodynamiques ordre eleve. Les symmetry et limites décrites modifications de conservation des limiteurs peuvent trouver une application chaque fois que possible en dehors des méthodes de reconfiguration.

Le nouveau laser modèle d'absorption pour la basse mousses de densité peut être appliqué dans les études futures du laser lissage du faisceau et la modélisation hydrodynamique d'autres matériaux microstructurés. Cela ouvre une voie à de nouveaux modèles numériques plus efficaces pour l'avenir recherche ICF.

  This thesis presents the work performed during my Ph.D. study under the dual supervision of Professors R. Liska and J. Limpouch at the Department of Physical Electronics of the Faculty of Nuclear Sciences and Physical Engineering of the Czech Technical University (CTU) in Prague jointly with Professors V. T. Tikhonchuk and J. Breil at the Center for Intense Lasers and Applications (CELIA) laboratory of the Bordeaux University. My mission in CTU was the development and implementation of new routines for

2 .

 2 State of the art electrons and ions are close to local Maxwellian thermal distributions.

2. 4 .= m e πc 2 e 2 λ 2 .

 42 High-energy laser-matter interaction constants such as the speed of light c, the electron mass m e and charge e n crit e

2 .

 2 State of the art ventional hydrodynamic model. A detailed information about the laser-foam interaction is needed for a reliable utilization of the low-density foam layer. The laser absorption, energy transfer and shock wave propagation in underdense 4 foams have been recently studied experimentally at nanosecond scale lasers PALS and GEKKO. The qualitative theory of the laser foam interaction introduces a fast homogenization stage, where empty spaces inside the foam pores are rapidly filled by the plasma expanding thin solid layers.

Figure 3 . 1 .

 31 Figure 3.1.: Remapping in a single cell. Black nodes of the Lagrangian computational mesh, red nodes of the rezoned mesh and the blue reconstruction function u R (x) are illustrated. Numerical fluxes are represented by the green areas.

  b, c) if a, b and c are positive 0 if a, b and c do not have equal sign -min(|a|, |b|, |c|) if a, b and c are negative (3.10)

3 . 2 .Figure 3 . 2 .

 3232 Figure 3.2.: Single remap of a step (a) and smooth (b) function with the unlimited piecewise parabolic reconstruction. Old grid with the mean values (circles) is in black, new grid is in red. Reconstruction is plotted in blue, its bounds in gray. A zoom is provided in the smaller bottom figures.

Figure 3 . 3 .

 33 Figure 3.3.: Single remap of a step (a) and smooth (b) function with the MM limited piecewise parabolic reconstruction. Old grid with the mean values (circles) is in black, new grid is in red. Reconstruction is plotted in blue, its bounds in gray. A zoom is provided in the smaller bottom figures.

Figure 3 . 4 .

 34 Figure 3.4.: Single remap of a step (a) and smooth (b) function with the KBJ limited piecewise parabolic reconstruction. Old grid with the mean values (circles) is in black, new grid is in red. Reconstruction is plotted in blue, its bounds in gray. A zoom is provided in the smaller bottom figures.

Figure 3

 3 Figure 3.5.: Single remap of a step (a) and smooth (b) function with the PPM piecewise parabolic reconstruction. Old grid with the mean values (circles) is in black, new grid is in red. Reconstruction is plotted in blue, its bounds in gray. A zoom is provided in the smaller bottom figures.

Figure 3 .

 3 [START_REF] Liska | Selected laser plasma simulations by ALE method[END_REF] presents an example of the reconstruction.Comparison of selected reconstructions for different functions is shown in Fig.3.6. On the contrary to the previous methods (unlimited, MM-and KBJ-limited), the PPM method guaranties the bound preservation by its construction.

Figure 3 . 6 .

 36 Figure 3.6.: 1D piecewise parabolic reconstructions with the different limiters for different mean values on different meshes.

. 33 )

 33 Definitions of the terms x s , y s , z c are provided later in the final formulas for the remaining coefficients Λ = (λ 1 , . . . , λ 5 ) T . The coefficients are calculated by minimization of the 3. Remapping methods reconstruction deviation functional in the surrounding cells k∈C

Figure 3 . 7 .

 37 Figure 3.7.: Piecewise constant (a), linear (b) and quadratic (c) reconstruction over a general quadrilateral mesh (solid black line).

. 44 )

 44 Here x c denotes the cell centroid and (∇ w) c an approximation of the gradient. An example of a constant and a linear vector reconstruction in a single cell is plotted in Fig.3.8. The gradient approximation in the given cell c has the form (∇ w) c = least-square fitting formula (3.35) used in the case of the scalar quantity,

Figure 3 . 8 .

 38 Figure 3.8.: Constant (a) and unlimited (b) linear reconstruction (blue) of a constant radial field (red) in a polar cell (black).

Figure 3 . 9 .

 39 Figure 3.9.: Notation of nodes and angles for a polar cell.

  [START_REF] Maire | Contribution to the numerical modeling of Inertial Confinement Fusion[END_REF]) components given by (3.46) and (3.47) has the general form

Figure 3 .

 3 Figure 3.10.: Limited piecewise linear reconstruction (green) of a constant radial field (red) for a polar cell (black).

  for the vector w c = w x c w y c in the local (with respect to the cell c) coordinates ξ, η given by

Figure 3 .

 3 Figure 3.11.: Bounds (from 3×3 stencil around the central cell c) for a vector reconstruction of: scalar extension of the BJ limiter (a), VIP limiter (b), MVIP limiter and an exact spherical limiter (d). The two identical cells with a different angular position (red and blue). Plot axes are the vector components w = u (x) u (y) .

3. 3 .

 3 2D reconstruction of a vector is provided in Fig.3.11 (b). The min/max (bounds) of the BJ limiter, as well as the limiter's projected form, are represented by a rectangle in the space in Fig.3.11 (a).

  the piecewise linear reconstruction(3.52) with the limiter matrix(3.53). The limited reconstruction in the cell vertices has a formw lim nc = w c + φ c (∇ w) c ( x n -x c ) . (3.63)For every cell center x c , we want to find φ c , so that the limited reconstructions at allnodes n ∈ N (c) (N (c)is the set of 4 nodes of the quadrilateral cell c) are kept inside the VIP bounds VIP c of the cell c. We require ∀n ∈ N (c) w lim nc ∈ VIP c (3.64) where VIP c represents the Convex Hull of the neighboring cell's velocities w k , k ∈ C(c) (C(c) is the set of 9 cells of the 3×3 path around the cell c). At first for all nodes n ∈ N (c), we find the maximal φ c ∈ (0, 1), which keeps w lim nc inside the VIP c . Let us call these values φ VIP nc and w VIP nc . Having computed all φ VIP nc corresponding to w VIP nc , we set the final φ c = min n∈N (c) φ VIP nc 7 vector image polygon 3. Remapping methods minimum of the nodal values φ nc which satisfy the condition (3.64). Computation of a particular φ VIP nc can be done in the following way. We denote the unlimited linear part of the reconstruction at the node n

Figure 3 .

 3 Figure 3.12.: Examples of three polygons P c are plotted in yellow (these polygons are shifted from the polygons in Fig. 3.11). Red arrows represent the linear part of min/max velocities w c c and the segments ( w cc 1 , w cc 2 ) are in blue.

  3. Remapping methods around the cell c, MVIP c = VIP c and therefore φ MVIP nc = φ VIP nc . Otherwise, we check if the limited vector w VIP nc lies in MVIP c and set

Figure 3 .

 3 Figure 3.13.: Exact (a) and swept (b) regions for the numerical fluxes quadrature. The dashed line denotes a Lagrangian-and the full line a rezoned-mesh [71].The red arrows mark the movement of the mesh nodes during the rezone.

. 73 )

 73 Starting with the maximal bound and substituting (3.71) and (2.5) into (3.72), we derive 3. Remapping methods the FCR constraints

F 81 ) 83 )

 8183 the flux is a combination of the low-order and high-order fluxes F µ,FCR,xy c,s = F µ,L,xy c,s + C c,s dF µ,xy c,s , dF µ,xy c,s = F µ,H,xy c,s -F µ,L,xy c,s . We express the momentum flux as a product of (scalar) mass flux F m c,s between cells c and s and the reconstructed velocity w c,s in the particular flux position (the centroid of the swept region) c,s + C c,s w H,xy The remapped momentum is rotated to the (ξ, η) coordinates by R c m c w w L,xy c,s + R c C c,s w H,xy of constraints (3.79) and (3.80) for the cell c multiplied by the remapped mass m c m c w ξ,min c w L,xy c,s + C c,s R c w H,xy w L,xy c,s + C c,s R c w H,xy The following point is to find the scalar flux correction factor C c,s for the each flux around the cell c. That is, to find the maximum possible C c,s that satisfies constraints (3.82) and (3.83) for the cell c and analogous constraints for the cell s ∈ C(c). Let us consider the constraints for the ξ-component of the momentum in the cell c (3.82).

,

  Now, we switch from the edge-related correction factors C c,s to the cell-related bounds D c . Let us define D ξ,the correction factors C c,s are chosen so that 0 ≤ C c,s ≤ D ξ,max c for all s such that dµ ξ c,s > 0, ξ c,s >0 C c,s dµ ξ c,s . Carrying out the same considerations for the other constraints, we have four quantities with the following sufficient condition for the bounds preservation (3.88) C c c,s ≤ min D ξ,min c looking at the interface c and s, while and we know all necessary values of R, dµ, Q and D. To compute the final correction factor C c,s , we use the minimum from the two cells c and s C c,s = min C c c,s , C s c,s . (3.91)

2 P ξ,L i 2 + P η,L i 2 ≤kC 3 . 7 .

 22237 ik dF de,min ik + l min 0, dF ds,min ikl , Remap of conservative variables for Euler equations which is translated via the worst case scenario into

8 i( 3 .

 83 100) with respect to (3.93 -3.95) provides a set of cell-based constraints D 1-102 -3.109).

ωλ 2 = 4πe 2 n crit e m e providing the critical density of free electrons n crit e = m e πc 2 e 2 λ 2 [ 4 . 1 .

 2241 k . Laser radiation cannot propagate in a denser plasma because of the electron screening effect. The plasma frequency ω p is a function of the free-electron density n e . Therefore, for a given laser frequency f , the electromagnetic field can propagate in the plasma only for the density smaller the critical value. Concerning laser angular frequency ω = 2πf = 2πc/λ for the laser vacuum wavelength λ, we can write 2πc cm -3 ]. (4.3) No dissipation processes are taken into account in this derivation, resulting in the reflection of the laser radiation. The plasma is a dielectric medium characterized by a for the light refraction index n and the absorption index χ as the real and the imaginary part of the complex permittivity n = ( √ ) and χ = ( √ ) . (4.5) Laser plasma interaction 4.1.2. Ray equation A ray equation describes a propagation of the laser radiation in a plasma in the geometric optics approximation, neglecting the diffraction phenomena. The approximation is valid in a weakly inhomogeneous plasma with the characteristic density scale length much larger than the laser wavelength. The laser beam is presented as a bunch of thin rays each carrying out a small portion of the whole beam power. The propagation of each other separate ray is calculated independently to each other. The ray trajectory is determined by the refractive index n (4.5) according to the ray equation

  the form ln Λ = max[2, ln b 2 max /b 2 min ] , where b max = (k B T e /m e ) 1/2 / max[ω p , ω] and b min = max[Ze 2 /k B T e , /(k b T e /m e ) 1/2 ] are the maximum and the minimum impact parameters in the electron-ion collision. According to (4.5), the permittivity determines both the laser ray propagation via the refraction index n and the laser absorption. The laser energy absorption along the ray trajectory is calculated from dQ ds = -κ ib Q , (4.11)

4. 2 .Figure 4 . 1 .

 241 Figure 4.1.: Rays focusing to a point (a) and caustic surface (b) for the focused laser beam. Gaussian beam distribution at the computational domain boundary (dashed line).

4 .

 4 Laser absorption modeling methods in the focal plane the beam intensity distribution corresponding to experimental conditions. A random ray deviation is better applicable for general situations and it provides a physically-relevant far-field image of an object in focal plane. This method has been introduced in the 3D ray-tracing algorithm of the code PALE. Low-density foams with a random microscopic internal material structure present another example of the configuration, where the 3D ray-tracing method is needed. To simulate the microscopic foam inhomogeneities on macroscopic level, we introduce a random reflection of laser rays at the homogenization surface of the foam. Non-cylindrically-symmetric rays are introduced in both examples, requiring a 3D ray-tracing treatment of the laser energy deposition.

  constant plasma parameters and thus a constant index of refraction n inside the cell, a straight-line solution of the ray-tracing equation 4.6 is obtained instantly dx + p , where x = s/n in the line parameter and constant vectors d = (d x , d y , d z ) T , p = (p x , p y , p z ) T

1 Figure 4 . 2 .

 142 Figure 4.2.: Ray propagation calculated by the piecewise line ray-tracing (solid line) and by the direct numerical method (circle) for the different ray incidence angles θ 0 = {0.1, 0.25, 0.5} rad on a coarse random grid.

4. 3 . 26 )

 326 Laser-foam interaction the calculation of the solid layer expansion at the microscale level, are estimated from the macroscale hydrodynamic model.The solid layer expansion is calculated by a standard hydrodynamic method with several sub-time-steps within each macroscale time-step. When the density of the modeled 1D layer becomes lower then the plasma critical density (4.25), the computational cell is supposed to be homogenized and transparent for the laser radiation.At the microscale level, a symmetric 1D expansion is supposed, so only a half of the layer is modeled. A staggered discretization with the cell centered scalar variables (density ρ, specific internal energy and pressure p of the foam material) and node centered vectors (position s e and velocity v e ) is used. The mass of the layer is divided between two nodes (the central node c and the edge node e) giving the nodal mass m n . The central node is rigid with the local position s o = 0 and velocity v o = 0. A temporal evolution of the velocity v e of the edge node is given by The microscale time step ∆t is obtained from a local CFL condition and the microscale pressure p = p(ρ, ) from the equation of state of a plastic. Here, p ext stands for the external (macroscale) hydrodynamic pressure. No compression of the layer is allowed, and thus an artificial viscosity term may be omitted. Finally, the nodal position s e and the specific internal energy are updated according to

( 3 .

 3 [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] unl.. Abbreviations for the slope and parabolic limiters correspond to BJ -Barth-5. Applications Jespersen (3.14), MM -minmod (3.9), respective (3.25), KBJ -Kuzmin-Barth-Jespersen (3.26) and Nejat -(3.27). In the case of the MM limiter we use β = 1.5 in (3.9) and (3.25) and we use K = 1 for Nejat limiter. The PPM marks the reconstruction according to (3.28) and lin. FCR, resp. pp FCR the FCR methods (2.5, 3.74) with the piecewise linear unlimited, respective piecewise parabolic unlimited reconstructions used for the high-order fluxes.In the cyclic remapping test cases, the first two steps of the ALE method, namely the Lagrangian step and the rezoning step, are replaced by the prescribed motion of the computational grid. The numerical fluxes through the cell interfaces (nodes in 1D) are present only due to the grid motion. The exact solution of the cyclic remap corresponds to the preservation of the initial value if the initial and the final grids are the same.

Figure 5 . 1 .

 51 Figure 5.1.: Computational grid sequence for 5 cells and 25 pseudo-time steps (a) in the pseudo-time interval t ∈ [0, 1] and for 8 cells in 10 steps with t ∈ [-1/8, 1/8] (b). The horizontal axis represent the spatial coordinate x and the vertical one the pseudo-time t.

  4) in all 1D test cases for two different pseudo-time intervals. As the first test case, we choose a smooth initial function sin(2πx)+1 on the computational domain x ∈ [0, 1]. The results indicating the order of convergence, errors and the bound violations for the first test are presented in Tab. 5.1. The sinus function has a local extrema in this domain, leading to the worse convergence of the strictly bounds-preserving methods, such as PPM and FCR. Further, the smooth grid movement described by (5.4) for the pseudo-time t ∈ [0, 1] . 4.1 4.0 4.0 3.99 • 10 -4 6.11 • 10 -6 7.0 • 10 -6 1.4 • 10 -8 9.9 • 10 -6 5.2 • 10 -8 lin. BJ 4.2 4.4 4.4 1.42 • 10 -3 1.72 • 10 -5 0 0 0 lin. MM 4.2 4.5 4.5 2.15 • 10 -3 2.58 • 10 -5 0 0 0 lin. FCR 4.0 3.9 3.9 7.27 • 10 -4 1.19 • 10 -5 0 0 0 pp unl. 8.0 8.0 8.0 9.66 • 10 -5 1.91 • 10 -7 7.0 • 10 -6 1.4 • 10 -8 1.2 • 10 -5 5.5 • 10 -8 pp KBJ 7.9 8.0 8.0 9.66 • 10 -5 1.92 • 10 -7 7.0 • 10 -6 1.4 • 10 -8 1.2 • 10 -5 5.5 • 10 -8 pp MM 7.9 8.0 8.0 9.75 • 10 -5 1.93 • 10 -7 7.0 • 10 -6 1.4 • 10 -8 1.2 • 10 -5 5.5 • 10 -8 pp FCR 5.6 5.6 5.5 4.39 • 10 -4 2.57 • 10 -47 • 10 -4 7.35 • 10 -6 5.4 • 10 -6 1.7 • 10 -8 1.6 • 10 -5 1.3 • 10 -6

Table 5 . 3 . 5 . 1 .

 5351 0 4.67 • 10 -7 9.32 • 10 -: Errors and their ratios for different sub-cyclic remaps of exp(x) on x ∈ [0, 1] for the pseudo-time t ∈ [-1/8, 1/8] on grids with N = 64, 128, 256 and cells. Cyclic remap numerical tests 5.3 shows the convergence of the selected remapping methods for t ∈ [-1/8, 1/8] for the function exp(x). The first-order is reached by the piecewise constant and the secondorder for the piecewise linear methods. For the most of the piecewise parabolic methods, i.e. pp unl, pp KBJ, pp MM, pp FCR and PPM, the fraction 8 indicating the third-order of convergence is reached. The poor convergence (just second order) of the Nejat limiter can be related to the selected constants in the Nejat limiter. Sub-cyclic remap of a discontinuous function Table 5.4 and Figure 5.2 show results for the step function f (x) = 4.6 for x ∈ [0, 1/2],

  . 1.5 1.6 1.5 5.03 • 10 -2 1.32 • 10 -2 7.7 • 10 -5 3.5 • 10 -6 6.6 • 10 -3 1.3 • 10 -3 lin. BJ 1.5 1.6 1.6 3.85 • 10 -2 9.65 • 10 -unl. 1.6 1.6 1.6 4.97 • 10 -2 1.19 • 10 -2 9.1 • 10 -5 4.0 • 10 -6 7.9 • 10 -3 1.7 • 10 -3 pp KBJ 1.6 1.6 1.7 5.05 • 10 -2 1.20 • 10 -2 6.5 • 10 -5 3.2 • 10 -6 7.6 • 10 -3 1.7 • 10 -3 pp MM 1.6 1.6 1.7 4.22 • 10 -2 9.71 • 10 -3 2.7 • 10 -6 1.7 • 10 -8 2.3 • 10 -5 2.7 • 10 -6 pp FCR 1.6 1.6 1.6 3.60 • 10 -2 8.88 • 10 -14 • 10 -2 1.03 • 10 -2 3.7 • 10 -5 1.5 • 10 -6 3.8 • 10 -3 1.1 • 10 -3Table 5.4.: Errors and their ratios for different sub-cyclic remaps of the step function on x ∈ [0, 1] for the pseudo-time t ∈ [-1/8, 1/8] on grids with N = 64, 128, 256 and 512 cells.

Figure 5 . 2 . 5 . 1 .

 5251 Figure 5.2.: Sub-cyclic remapping for the selected remapping methods. The step function for the pseudo-time t ∈ [-1/8, 1/8] on 16 cells. The initial (exact) function marked by the black line.

Figure 5 . 3 .

 53 Figure 5.3.: Advection-like remap test case for a scalar.

Method

  3.96 3.72 1.08•10 -2 2.73•10 -3 7.33•10 -4 pp FCR 7.84 7.02 6.44•10 -3 8.21•10 -4 1.17•10 -4

  81 8.01•10 -5 2.15•10 -5 5.64•10 -6 Parabolic & exact intersections 7.46 7.67 3.25•10 -5 4.36•10 -6 5.69•10 -7 Linear & swept fluxes 3.73 3.82 7.96•10 -5 2.13•10 -5 5.59•10 -6 Parabolic & swept fluxes 7.46 7.67 3.23•10 -5 4.34•10 -6 5.65•10 -7

Figure 5 .Figure 5 . 4 .

 554 Figure 5.4 shows an example of the initial and final mesh for this mesh movement (on a computational domain [0, 1] × [0, 1]), where the corresponding cells on two meshes are plotted with the same color. The initial grid is indicated by a superscript 0, n t = {100, 200, 400} stands for the total number of remapping steps and M = {21, 41, 81} for the number of grid points in each direction. Results in Tab. 5.6 show the expected

Figure 5 . 5 .

 55 Figure 5.5.: 2D cyclic remap: an example of the initial-(a) and final (b) quadrilateral grid and the particular mesh displacement (c) for a polar grid. Colors in (a) and (b) are used for visualization purposes.

5. 1 . 25 1Figure 5 . 6 .Figure 5 . 7 .

 1255657 Figure 5.6.: Contour plots of the velocity magnitude for the cyclic remap of the step function in the radial velocity field on a polar grid. The low-order (a), unlimited (b), BJ-(c) and MVIP-limited (d), and FCRxy (e) method. Out-of-bound values in red.

Figure 5 . 8 .

 58 Figure 5.8.: Velocity magnitude distribution for the Sedov problem at time t = 1 using the low-order (a), unlimited (b), component FCRxy (c), and symmetric FCRs (d) method for velocity reconstruction.

  5.11. For this Riemann problem, the initial condition values of density, velocity and pressure are {1,0,1} for x < 0.5 and {0.125,0,0.1} for x > 0.5 with a gas constant γ = 1.4. The solution of the Sod test case consists of three type of waves, namely the shock, contact discontinuity and rarefaction wave. We use the Eulerian mode, i.e. the remapping at every time step to the initial grid to emphasize the properties of the remap. The final density, velocity and specific internal energy are shown in Fig.5.11 at the time 0.2. The piecewise quadratic

Figure 5 .

 5 Figure 5.11.: 1D Sod shock tube, Eulerian mode, 3 × 100 cells. Density, velocity and specific internal energy for the whole domain (a) and in details (b) for the piecewise linear BJ limited method (red) and for the piecewise quadratic FCR method (blue).

Figure 5 .

 5 Figure 5.13.: Density scatter plot of the whole computational domain (top) and a zoom (bottom) for the 2D polar Sod shock test case in the Eulerian mode for 7 × 100 cells. Piecewise linear BJ limited method in red, piecewise quadratic FCR in blue and a reference Lagrangian solution in green.

2 . 1 .

 21 Because of the different computational domain, the initial internal energy of a central cell is modified to maintain the same amount of the total deposited energy. The first case in the Euler Method 43 4.94•10 -1 3.68•10 -1 2.58•10 -1 lin FCR 1.35 1.51 4.43•10 -1 3.28•10 -1 2.17•10 -1 pp FCR 1.51 1.56 4.13•10 -1 2.73•10 -1 1.75•10 -1

Figure 5 . 5 . 2 .

 552 Figure 5.14.: Scatter plot of density versus radius for unstructured Sedov test case, remapping method with a piecewise constant (a), FCR with a piecewise linear (b) and a piecewise parabolic (c) reconstructions. ALE mode.

- 1

 1 directed towards the origin on the square domain [-0.3,0.3] × [-0.3,0.3] covered by the equidistant rectangular grid with 100 × 100 cells. The initial and boundary pressure is set to 1 • 10 -6 and we apply the exact density condition on the computational domain boundary ρ(r) = (1 -(0.6/r)) 2

Figure 5 .

 5 Figure 5.15.: Scatter plot of density versus radius for Noh test case, remapping method with a piecewise constant (a), FCR with a piecewise linear (b) and a piecewise parabolic (c) reconstructions. Eulerian mode. Exact solution ρ = 16 for r < 0.2 and ρ = 1 -((0.6/r)) 2 for r > 0.2 plotted by a black line.

  Figure 5.16.: Laser penetration through the foam target. The ionization front surface is shown in blue. Central parts of the three laser beams are emphasized. 3D ray-tracing at 700 ps.

5 . 8 )

 58 provides the same inaccurate result as the previous straight ray model.3D ray-tracing with a foam modelAccording to the laser absorption model[START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], the two states of computational cells are defined. The initialization of this homogenization state of the cells follows the logic introduced in subsection 4.3.3, where except of the first layer of cells on the laser-irradiated target side, all other cells are supposed to be non-homogenized at the time t = 0. Concerning the non-homogenized structured cells, the laser rays are fully reflected in a random direction into the half-space defined by the cell position and the direction of incoming laser ray. The cold cell heating in this model is achieved by the electron heat conduction only.

Figure 5 .

 5 Figure 5.17.: Simulated spatial distributions of the foam plasma at 700 ps according to the different laser absorption models. The laser propagates from the right to the left.

Figure 5 . 3 ]Figure 5 .

 535 Figure 5.18.: Modeled X-ray streak camera image of lateral foam self-emission for different laser absorption models. The laser is incident from the left.

5. 4 . 1 .

 41 Target and laser setupThree experimental setups are modeled. The first two correspond to the PALS experiments at the laser wavelength 438 nm[START_REF] Khalenkov | Experience of micro-heterogeneous target fabrication to study energy transport in plasma near critical density[END_REF][START_REF] Limpouch | Laser absorption and energy transfer in foams of various pore structures and chemical compositions[END_REF][START_REF] Borisenko | Regular 3-D networks with clusters for controlled energy transport studies in laser plasma near critical density[END_REF]. The laser pulse with the Gaussian temporal profile and the FWHM (full width at halt maximum) 320 ps delivering energy of 150 -190 J (for the simulation, we use the average value of 170 J). The focal spot has a Gaussian profile with the FWHM 300 µm. Two targets have been modeled: one is the 380 µm TAC foam with the average density 4.5 mg/cm 3 and the second -400 µm TAC with the density 9 mg/cm 3 . The different thickness of the foams is a consequence of a foam fabrication process. The chemical content of the TAC corresponds to the effective charge and mass numbers Z = 4.54 and A = 8.73. The denser foam was supplemented by the 5 µm thick Aluminum foil at the rear foam side. For the comparison with the simulation, we choose the shots no. 28193 and 28205 presented in[START_REF] Limpouch | Laser absorption and energy transfer in foams of various pore structures and chemical compositions[END_REF].

  code. The electron density, temperature and absorbed-energy distributions are shown in Fig. 5.22 for the time 450 ps. The simulated X-ray streak camera images of a lateral foam self-emission are presented in Fig. 5.20. Experimental images are presented in Fig. 5.21. These numerical simulations for the three different parameters show good scaling and qualitative agreement with the experiments [75, 83]. The lowest-density foam in Fig. 5.20 (left) is ionized during the laser pulse duration. The laser penetration through the higher-density foam (center) is slower and the rear-side foil self-emission occurs at around 1.3 ns, after the end of the laser pulse. 5 mg/cm 3 PALS 9 mg/cm 3 GEKKO 10 mg/cm 3

Figure 5 .

 5 Figure 5.20.: Simulated X-ray streak camera image of the lateral foam self-emission.

Figure 5 .Figure 5 . 5 . 4 .

 5554 Figure 5.21.: Observed X-ray streak camera image of the lateral foam self-emission [75, 83].

  laser facilities. These simulations demonstrate the practical interest of the model and of the application part of this thesis. Contrary to the standard hydrodynamic simulations of a foam as a homogeneous material, presented results open a new approach of efficient modeling of microstructured materials. They show a good agreement of the ionization 6. Summary front speed in the sub-critical foam with the experiments.

Newton,B. 2 .

 2 Lagrange polynomiale, élimination de Gauss ou méthode d'Euler, sont la preuve claire de cela. Cependant, le plus grand développement des méthodes numériques a commencé dans les années 1940, lorsque les premiers ordinateurs sont devenus disponibles.Incroyablement rapide et gratuit-de-erreurs ordinateurs ont remplacé le non-facteur humain efficace dans les procédures. Le travail de personnes a été transféré de l'exécution à l'élaboration de nouvelles méthodes. Les méthodes numériques sont devenus largement appliquée dans le large éventail de la science et de la technologie. Gamme applicable des méthodes numériques a été d'élargir encore par l'invention de la haute-classe pour des méthodes. Exigences inférieures du temps de calcul et la mémoire des haut les modes de commande permettent leur application, même pour les systèmes complexes, tels que le plasma produit par de hautes installations laser de puissance.Ce travail est dédié aux méthodes numériques pour laser-plasma modélisation de l'interaction. Plus précisément, la thèse traite le plasma créé par l'interaction d'une impulsion laser intense avec la matière. Plus précisément, on développe les méthodes pour résoudre les équations hydrodynamiques, compte tenu du plasma comme fluide compressible avec d'autres procédés physiques tels que l'absorption de la lumière laser ou la conduction de chaleur. Le modèle hydrodynamique compromis entre une description physique détaillée et une efficacité de calcul. L'efficacité est cruciale pour les simulations de processus complexes comme le laser-importance des expériences d'interaction.Expériences d'interaction laser-plasma couvrent généralement une large gamme de paramètres physiques. La densité caractéristique varie d'une limite de vide pour quelques fois la densité à l'état solide et les échelles de température de quelques milliers à centaines de millions de degrés Celsius ou Kelvin. La recherche sur de tels états de densité d'énergie élevée de l'affaire a commencé environ 15 ans avant la découverte d'un laser luimême lorsque la libération d'une quantité énorme d'énergie dans les réactions nucléaires de fission et de fusion plus tard, a été réalisé dans des explosions. Une réduction supplémentaire des expériences de plasma et un contrôle possible de l'énergie libérée a été activé par la découverte des lasers en tant que source d'énergie à haute puissance. Les densités d'énergie suffisantes pour les modifications de matériau solide et la production de plasma ont été atteints dans les années 1960. Depuis, de nombreuses applications pratiques de la matière d'interaction laser ont été proposés. Le plus difficile parmi eux serait une libération de l'énergie thermonucléaire dans une cible d'hydrogène comprimé par un puissant lasers. Tel est l'objectif pour les laboratoires de recherche haut du monde. Objectif de la thèse L'objectif général de la thèse est de contribuer à l'état de la connaissance de l'art dans les haut-méthodes numériques afin ALE, remappage en particulier, pour les équations d'Euler et de l'application des méthodes de ALE pour la modélisation des processus d'interaction plasma laser. La diffusion numérique lors de l'étape de reconfiguration nécessaire du procédé ALE est réduite par développés methodes remappage d'ordre eleve. La conservation de la symétrie des méthodes pour les vecteurs limitation est étudiée en détail dans le cas des flux radiaux. Ces questions ont été récemment largement discutés par la communauté de calcul de dynamique des fluides. Comme une application particulière, une interaction d'un rayonnement laser avec une cible de mousse à faible densité est étudiée. Des expériences avec des mousses microstructurées ayant la densité moyenne de quelques mg/cm 3 montrent une amélioration significative de la qualité d'un faisceau laser et son lissage après propagation à travers une mousse. Ce lissage est essentiel pour les applications fusion par confinement inertiel (ICF), car elle permet de réaliser une bien meilleure qualité de l'implosion. Le mécanisme physique exact de l'ionisation de mousse par la lumière laser est de nature microscopique et codes hydrodynamiques standard, qui envisagent la mousse comme un média en continu d'une densité équivalente, surestimer la vitesse de la pénétration du laser dans les mousses. Dans cette thèse, nous proposons un nouveau modèle de mousse ionisation et une modification de la méthode d'absorption de l'énergie laser. Il représente la structure de mousse microscopique dans le modèle hydrodynamique continue standard. Ce nouveau modèle correspond aux résultats expérimentaux de l'interaction laser-mousse. Ces modifications ont été atteints par la mise en oeuvre d'un algorithme de ray-tracing 3D en 2D cylindrique symétrique code PALE. B.3. Démarche adoptée Dans la partie principale de cette thèse, nous avons présenté la méthode de reconfiguration de l'ordre eleve pour le système d'équations d'Euler décrivant le comportement d'un fluide non visqueux compressible, par exemple applicable pour le laser-plasma. Les questions liées à l'symmetry et lié la préservation de la méthode ainsi que la haute ordre de précision sont étudiées. La méthode est testée sur cycliques standards -problèmes remappage. Incorporé dans l'arbitraire complet Lagrange-Euler méthode (ALE), la performance de la reconfiguration est également démontré dans la cellule centrée Code CHLER pour les cas de tests hydrodynamiques standard. Notre méthode de reconfiguration en utilisant une reconstruction quadratique par morceaux et un symétrique Flux de transport corrigés idées est directement applicable aux non structurées deux mailles de calcul dimensions avec la topologie maillée fixe. Limiteurs vectorielles en conservant une symétrie radiale des zones radiales sur mailles polaires sont décrits en détail dans le contexte de la reconstruction de limitation lors de l'étape de remappage de la méthode d'ALE. Symmetry et consolidés modifications préservant des limiteurs sont numériquement testés dans le code hydrodynamique quinconce PALE. Des aspects particuliers de l'interaction laser de faible cibles en mousse de densité sont abordées dans la partie restante de la thèse. Simulations hydrodynamiques des matériaux cibles microstructurées nécessitent un traitement spécifique de l'absorption du laser et la réponse du matériau. La méthode d'une heure locale dépend coefficient d'absorption est généralisé pour être applicable dans le deux cylindrique dimensions codes hydrodynamiques symétriques, comme par exemple PALE. La généralisation comprend la mise en oeuvre du rayon 3D -algorithme de traçage et un sous cyclée simulation de l'expansion des éléments de mousse solides induites par le laser et l'énergie des flux d'électrons. Ce nouveau modèle est testé dans les simulations numériques du laserexpériences d'interaction de mousse effectuées à l'PALS et GEKKO installations laser. Ces simulations montrent l'intérêt pratique du modèle et de la partie de l'application de cette thèse. Contrairement aux simulations hydrodynamiques classiques d'une mousse comme un matériau homogène, les résultats présentés ouvrent une nouvelle approche de la modélisation efficace des matériaux microstructurés. Ils montrent un bon accord de la vitesse du front d'ionisation dans la sous -mousse critique avec les expériences.Les simulations numériques détaillées du laser généré plasmas sont hors de portée de cette thèse. Au lieu de cela, nous abordons les principaux éléments des méthodes numériques particulières requises pour effectuer des simulations plus efficaces et plus précises.

  definition of the limiter matrix(3.53).

	1.4						'cell'
							'xy'
	1.2
	1
	0.8
	0.6
	0.4
	0.8	1	1.2	1.4	1.6	1.8	2	2.2

Table 4 .

 4

		.1
	physical constant	symbol value and units
	Boltzmann constant	k B	1.602 • 10 -12 erg/eV
	reduced Planck const.		1.055 • 10 -27 erg•s
	electron charge	e	4.803 • 10 -10 statC
	electron mass	m e	9.109 • 10 -28 g
	speed of light	c	2.998 • 10 10 cm/s

1.: CGS physical constants.

Table 5 .

 5 2.: Errors and their ratios for different sub-cyclic remaps of sin(2πx) + 1 on x ∈ [0, 1] for the pseudo-time t ∈ [-1/8, 1/8] on grids with N = 64, 128, 256 and 512 cells.

Table

	Remap	64 128	128 256	256 512

Table 5 .

 5 6.: 2D cyclic remap of a smooth function. Unlimited piecewise polynomial reconstruction with different methods of numerical flux computation.

Table 5 .

 5 7.: Convergence in density errors for the Sedov problem on a hexagonal grid with 10/20/40 cells per edge [-1.2,1.2]. Three different remapping methods. Eulerian mode. mode on the equidistant hexagonal grids with 20/40/80 cells on the edge [-1.2,1.2] of the computational domain. The errors and convergence ratios for the density error (computed from a comparison with 1D solution on 1000 cells given by

Table 5 . 8

 58

.: Simulated ionization front propagation time through the 500 µm foam (with respect to the 100 eV temperature). The measured propagation time is 1.7 ns.

The proposed 3D ray-tracing foam absorption model demonstrate a certain reduction of the ionization propagation speed moving it closer to the experimental measurement.

  B.1. Contexte de la rechercheLes méthodes numériques, à savoir les procédures de résolution de problèmes mathématiques en général, sont utilisés avec succès pour des milliers d'années. Probablement la plus ancienne preuve connue d'entre elles est la table de Babylone montrant un algorithme pour calculer la racine carrée du nombre deux. Méthodes numériques ordinaires, qui sont enseignées dans les universités techniques, sont datés à la période de 17. à 19. siècle, la période de croissance de la mécanique classique. Leurs noms, par exemple la méthode de

doubling the number of computational cells in each spatial dimension

we use the term parabolic in a 1D context and the term quadratic for more spatial dimensions

the foam with an average density lower than the plasma critical density

or abbreviated donor

but u i = u i , see(3.18)

In most examples, the cells are quadrilaterals, but the method can be applied to other cell shapes.

The integration over cell vertices is useful for a staggered discretization. For a cell-centered discretization, the approximate integration can be performed over centroids of neighboring cells.

angle-independent for the polar case

note that these coordinates differ from ξ, η in(3.51)

nodes in 1D and edges in 2D

In 1D, we denote the fluxes positive if they are in the direction of the x axis.

We call this additional step of a remapping method non-physical because there is no physical reason to introduce additional redistribution of the interpolated quantities.

In 2D case, we omit writing the line over mean values and replace the superscript n by tilde for remapped quantity to distinguish from 1D case.

Remapping methods

These rays are typical for the beams that are not perfectly aligned with the target axis, e.g. the beams of recent multiple-beam ICF laser facilities.

Laser absorption modeling methods

The grid does not cycle back to the initial position, i.e. x f = x 0 , however we keep the name cyclic, or we use the name sub-cyclic to emphasize this difference.

Acknowledgements

This research has been conducted during the years 2011 -2015 and supported in parts by the Czech Technical University grants SGS10/299/OHK4/3T/14 and SGS13/220/OHK4/3T/14, the Czech Science Foundation projects P205/10/0814, P201/12/P554 and 14-21318S, Czech Ministry of Education grants MSM 684077002, and RVO: 68407700.

Parts of this work was performed under the auspices of the National Nuclear Security Administration of the US Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement number 633053.

Laser absorption modeling methods

The case where the square root in (4.13) in negative corresponds to the internal reflection.

Then the ray direction is inverted

(4.15)

To finalize the ray-tracing algorithm, the intersections of the conical cell faces with the rays have to be calculated. The problem leads to the solution of a single quadratic equation at each interface [START_REF] Kaiser | Laser ray tracing and power deposition on an unstructured threedimensional grid[END_REF]. Knowing the ray initial energy Q in and trajectory, the absorption of the laser energy in each particular cell ∆Q is given by

with the κ ib given by (4.12) and the optical path L = | p in -p out | of the ray in the cell.

The 3D ray-tracing algorithm, as described above, is applied for the simulation of a low-density foam target irradiated by an intense laser beam. The results of the simulation are presented in the fifth chapter of this thesis.

Laser-foam interaction

In this section, modifications of a hydrodynamic code are described, which are introduced in order to provide a correct model of the laser absorption in a low density foam. It enables to describe physical phenomena as the ionization front propagation and plasma heating. These modifications allow to simulate microstructured targets with standard hydrodynamic codes without considering a microscopic foam structure on the hydrodynamic level.

An interaction of a low-density foam target with a high-intensity laser pulse is modeled on two scales. At the macroscale level, the Lagrangian hydrodynamic model considers an artificial transient material state in the homogenization layer and a modified laser-plasma interaction, which simulates complex microscopic physics of the laser energy deposition.

The both modifications are governed by an evolution of the detail foam structure calculated at the microscale level.

The microscale physics describes the local homogenization of the foam by an expansion of a thin plastic layer in 1D geometry. This model was firstly introduced in [START_REF] Guskov | Laser-supported ionization wave in under-dense gases and foams[END_REF], where the expansion speed of each thin solid layer was governed by the local plasma sound speed. The microscale structure is introduced only in one layer of computational cells leads to the super-convergence of the piecewise linear unlimited method (errors in the norm (5.1) show the third-order of accuracy corresponding to the errors ratio 8). It seems that errors, accumulated during the mesh movement to the right, are canceled with the errors introduced by the backward motion. To eliminate this super-convergence, we use the same grid motion (5.4) however with the pseudo-time t ∈ [-1/8, 1/8], which we call sub-cyclic remap as the final mesh is different from the initial one (so it is not cyclic).

Sub-cyclic remap of a sinus function

As stated above, the same initial function sin(2πx) + 1 is discretized on the initial mesh, which is however different from the first test case. The grid motion (5. The super-convergence of the piecewise linear unlimited method (lin. unl.) is canceled by the change of the pseudo-time interval. In Tab. 5.1 we see, that the donor is first order, all methods with piecewise linear reconstructions are second order, piecewise parabolic pp unl, pp KBJ and pp MM ate third order, pp FCR and PPM are better than second order and pp Nejat is second order for this test. Note, that lin. unl. method does not satisfy bounds, while limited piecewise linear methods satisfy bounds. From piecewise parabolic methods, only pp FCR and PPM satisfy bounds, when pp unlimited and pp limited by KBJ or MM do not satisfy bounds. This is related to the existence of local extrema of the A.

A.2. Poster presentations

• MultiMat 2013, September 4, 2013, San Francisco, USA.

• Journée de l'Ecole Doctorale, March 5, 2013, Talence, France.

• CCP2012, Conference on Computational Physics, October 16, 2012, Kobe, Japan.