Section: New Results

Numerical developments

Numerical developments

Conformal hexahedral mesh coarsening by agglomeration

Participants : Hervé Guillard, Pierre Cargemel, Youssef Mesri [IFPEN] .

This work has been realized in the framework of a PhD contract with IFPEN that aims to produce a coarsening software for hex-dominant meshes. Reservoir simulation involves to compute dynamic flow of different phases in a porous medium. The initial state of the reservoir is usually pre-computed via geo-statistics methods extrapolating measures of the terrain. Therefore, the input of reservoir simulation is given as a very fine mesh containing heterogeneous data and numerical simulation on this fine mesh is usually non-practical. This work is therefore devoted to the study of an agglomeration strategy, to dynamically coarsen this fine hex-dominant mesh. The adaptivity may be driven by physics and/or geometric estimators. Ideally, the coarsening should be applied locally in low gradient regions, whereas high gradient regions keep the fine mesh. This work has been presented in the 23rd International meshing Roundtable [26] . The planned sequel of this work consists to use the notion of Central Voronoi Tesselation (CVT) to treat the regions where the mesh is not structured and to apply this strategy in different physical contexts from plasma physics to petroleum engineering.

Mapped Fourier Methods for stiff problems in toroidal geometry

Participant : Hervé Guillard.

Due to the particular geometry of tokamaks, a lot of numerical codes developed for their numerical simulations use Fourier methods. Fourier spectral or pseudo-spectral methods are extremely efficient for periodic problems. However this efficiency is lost if the solutions have zones of rapid variations or internal layers. For these cases, a large number of Fourier modes are required and this makes the Fourier method unpractical in many cases. This work investigates the use of mapped Fourier method as a way to circumvent this problem. Mapped Fourier method uses instead of the usual Fourier interpolant the composition of the Fourier interpolant with a mapping in such a way that in the computational space, the functions to represent are not stiff. This work gives some examples of the usefulness of this method and apply it to a simple model of pellet injection in tokamaks as an example of its potential interest for complex multi dimensional problem [34] .

Multislope MUSCL method for general unstructured meshes

Participants : Hervé Guillard, Clément Le touze [ONERA] , Angelo Murrone [ONERA] .

To increase the accuracy in finite volume method, the concept of MUSCL reconstruction has been introduced in the pioneering work of van Leer in the 70'. This technique is still one of the most efficient to deal with the existence of discontinuous solutions in numerical simulations. In the MUSCL technique, a discontinuous linear approximation of the solution is reconstructed on each control volume. The main approximation problem of this method is therefore to reconstruct the slope of the solution.

The multislope concept has been recently introduced in the literature to deal with MUSCL reconstructions on triangular and tetrahedral unstructured meshes in the finite volume cell-centered context. Dedicated scalar slopes are used to compute the interpolations on each face of a given element, in opposition to the monoslope methods in which a unique limited gradient is used. The multislope approach reveals less expensive and potentially more accurate than the classical gradient techniques. Besides, it may also help the robustness when dealing with hyperbolic systems involving complex solutions, with large discontinuities and high density ratios. In this work, we have designed a generalized multislope MUSCL method for cell-centered finite volume discretizations. The method is freed from constraints on the mesh topology, thereby operating on completely general unstructured meshes. Moreover optimal second-order accuracy is reached at the faces centroids and the scheme is $L^{\infty }$ stable. Special attention has also been paid to equip the reconstruction procedure with well-adapted dedicated limiters, potentially CFL-dependent. We have shown in [18] the ability of the method to deal with completely general meshes, while exhibiting second-order accuracy.

Development of a two temperature model

Participants : Hervé Guillard, Afeintou Sangam, Elise Estibals.

A two temperature (ions - electrons) model for non-magnetized plasma has been designed. The numerical scheme is a finite volume method with an approximate Riemann solver using the total energy equation and the electron entropy as main variables. This Riemann solver has been validated against standard shock tube problems and incorporated in the PlaTo platform. The solver has been implementated in toroidal geometry and tested successfully on realistic particular flows encountered in this context. The development of a reduced MHD model based on this two temperature scheme is currently studied.

Entropy Preserving Schemes for Conservation Laws

Participants : Christophe Berthon [University of Nantes] , Bruno Dubroca [CEA/DAM/CESTA and University of Bordeaux 1] , Afeintou Sangam.

A relaxation-type scheme has been proposed to approximate weak solutions of Ten-Moments equations with source terms [2]. These equations model compressible anisotropic flows. Following the technique introduced in [44] , the proposed scheme is proved to be entropy preserving.

Eurofusion WPCD: Free boundary equilibrium code and control

Participants : Cédric Boulbe, Blaise Faugeras, Jean François Artaud [IRFM CEA Cadarache] , Vincent Basiuk [IRFM CEA Cadarache] , Emiliano Fable [Max-Planck-Institut für Plasmaphysik, Garching] , Philippe Huyn [IRFM CEA Cadarache] , Eric Nardon [IRFM CEA Cadarache] , Jakub Urban [IPP, Academy of Sciences of the Czech Republic, Prague] .

Our team is involved in the integrated modelling WPCD (Work Package Code Development) Eurofusion. This project is the continuation of the EFDA-ITM project. The goal of WPCD is to provide a european tool for tokamak simulations. Different physical codes can be coupled using Kepler environment. Machine description and physical data have been described using CPO (Consistent Physical Objet) which are used as standardized inputs and outputs for the codes.

In this project, we participate in the coupling of a free boundary equilibrium solver, the European Transport Solver (ETS) and a plasma shape and position controller. The workflow coupling TCV hybrid Simulink controller and Cedres++ using PF circuit connections has been finalized and tested on the TCV tokamak.

A new workflow coupling Cedres++ with ETS and the TCV controller has been developed and is being tested on a TCV test case. This workflow is an evolution of the coupling CEDRES++ - ETS described in [12] .

A successful benchmark between the three free boundary equilibrium codes CEDRES++([15] [47] ) , FREEBIE ( [43] ), and SPIDER ( [49] ) has been done on static test cases. This activity will be continued to compare the time dependent versions of the three codes.

Optimal control for scenario optimization of discharges in tokamaks

Participants : Jacques Blum, Holger Heumann.

In this project we aim for an automatic determination of optimal voltage evolutions via an optimal control formulation based on a system of partial differential equations that describes the evolution of plasma equilibrium in a tokamak. Optimal voltage evolutions are the one that ensure that the evolution of the plasma runs through predescribed, user-defined states, defined e.g. as desired evolution of shape or position. The system of partial differential equations describing the evolution of the plasma is non-linear and we use a finite element formulation together with implicit time stepping for the discretization [15] . With this approach we end up with a large but finite dimensional optimization problem with non-linear constraints. We are using SQP (sequential quadratic programming), known to be one of the fastest methods of such problems, to solve the finite-dimensional optimization problem. The performance of SQP relies on accurate derivatives of the objective function and the constraints. The derivatives related to the free boundary, derived and implements during H. Heumann’s PostDoc 2011/2012 for a static optimal control problem, appeared here again and are one of the important building blocks for treating the transient case. Both in CEDRES++ and FEEQ.M we have now the capability to solve first test cases to define optimal voltage evolutions. In contrast to the static case, where the linear algebraic systems in the SQP iteration remain reasonable small, the solution of the corresponding linear system in the transient case becomes very time-consuming, which somehow limits the applicability. We are testing variants of SQP, such as BFGS-like updates for the reduced Hessian, to see whether we could speed up and improve robustness of our calculations. Fast iterative solver for large sparse linear systems is another option that we started to investigate. Fast iterative solver for linear system in transient optimal control problems governed by partial differential equations is a very active area of research and we hope to benefit from the latest developments.

Boundary reconstruction for the WEST tokamak with VacTH

Participants : Jacques Blum, Sylvain Bremond [IRFM CEA Cadarache] , Cédric Boulbe, Blaise Faugeras, Holger Heumann, Philippe Moreau [IRFM CEA Cadarache] , Eric Nardon [IRFM CEA Cadarache] , Remy Nouailletas [IRFM CEA Cadarache] , François Saint Laurent [IRFM CEA Cadarache] .

This work is under progress in collaboration with the CEA. The control of the plasma in the future WEST tokamak requires the identification of its boundary in real time during a pulse. The code VacTH under development in the team enables such an identification. Several numerical developments and experiments have been conducted in order to prepare the control of the plasma in the WEST tokamak. The equilibrium code CEDRES++, also developed in the team, is used to simulate a real plasma and to generate synthetic magnetic measurements from which the plasma boundary is reconstructed using the code VacTH. A control algorithm developed by the colleagues from the CEA then uses this knowledge of the plasma shape to adapt the currents flowing in the poloidal field coils in order to achieve a desired evolution of the plasma.

Equilibrium reconstruction for ASDEX UpGrade (AUG) with Vacth-Equinox

Participants : Blaise Faugeras, Rui Cohelo [IPFN, IST, Lisbonne] , Patrick Mccarthy [National University of Ireland University College Cork] .

Within the framework of the WPCD EUROFUSION the code VacTH-Equinox has been adapted to enable equilibrium reconstruction for AUG. The identification of the current density pedestal required the development of a specific regularization scheme allowing weaker regularization close to the plasma boundary and stronger close to the magnetic axis.

Taylor-Galerkin stabilized Finite Element

Participants : José Costa, Boniface Nkonga.

The theoretical part of Taylor-Galerkin/Variational multi-scales (TG/VMS) strategy applied to MHD and reduced MHD modeling has been achieved last year. The final method amounts to add in the finite element formulation, a self-adjoint operator associated to the most critical hyperbolic component of the system to be solved. The design of the critical contours and the identification of associated waves to be stabilized is problem dependent and related to the Jacobian matrix. We have focused this year on the validations of this strategy and the improvement of the linearization used for stabilization. For application to plasma configurations with X-point, we need to reconsider the consistency with equilibrium and the Bohm boundary conditions on open flux walls.

Toward full MHD numerical modeling with ${𝒞}^1$ FE

Participants : José Costa, Giorgio Giorgiani, Hervé Guillard, Boniface Nkonga.

In this context the single fluid full MHD model is considered and the divergence free constraint on the magnetic field is achieved by introduction of a potential vector. The use of the potential vector has the additional advantage that the toroidal component is the magnetic flux of the Grad-Shafranov equilibrium. However, using the potential vector as variable introduces higher order derivatives in the system and classical C0 finite elements cannot be directly applied. This is why our finite element strategies use shape/test functions whose derivatives have global continuity in space (smooth finite elements). The global approach uses cross product shape/test functions between poloidal(2D) and toroidal(1D). In the 2D poloidal plane, discretization uses either quadrangular or triangular elements.

This year we have focused on the numerical analysis associated to the full MHD discretization in configurations with open flux surfaces. In order to derive efficient strategies for the full MHD in the potential vector formulation, the Gauge condition on the potential vector and the boundary conditions have been enforced by penalizations. For the Gauge condition it gives rise to element contributions but also boundary integrals that should be computed on curved surfaces that sometime fitted the magnetic surfaces. Equations are formulated in semi-conservative form such as to apply integration by part. Therefore, boundary conditions can be viewed as evolution of fluxes or variables. Integral formulation on the boundary is very useful for higher order finite elements and also for easier treatment of corners. Indeed in this context the boundary conditions are edges/surfaces oriented and boundary corners are driven by the neighborhood edges penalizations. This strategy is the one that will be used for future developments.

2D Quadrangular Cubic Bezier Finite Elements: This finite element is used for a while for reduced MHD models in the software Jorek. Reduced MHD is used to project the momentum equation in a space orthogonal to the equilibrium. When full MHD models are used, the momentum equation needs to be projected in the equilibrium space and this projection should be consistent with the Grag-Shafranov equilibrium that is used to compute the initial state. This has been achieved by a proper computation of the JxB contribution in the momentum equation, taking into account the poloidal variation of the toroidal component of the magnetic field. Detailed analysis has been performed. The next year will be devoted to implementations and numerical validations.

2D Triangular Powel-Sabin Finite Elements: In order to avoid some mesh singularities when using quadrangular meshes for complex geometries and flux surfaces shapes, triangular elements are a possible option. It is not so easy to derive smooth finite element on triangle with reduced number of degree of freedom. The Bell reduced-quintic finite elements we have considered in the previous years have too much unknowns (6 per vertex). Powell-Sabin splines are piecewise quadratic polynomials with a global ${𝒞}^1$-continuity and 3 unknowns per vertex, they have a local support, they form a convex partition of unity, they are stable, and they have a geometrically intuitive interpretation involving control triangles. Construction of the Powel-Sabin splines needs some geometrical tools that have been developed: Minimum area enclosing triangle of a set of control points (sequential and parallel). This construction is applied to each vertex of the triangular mesh and used to derive the local shape/test functions. These Powel-Sabin splines have been used successfully in the area of computer aided geometric design for the modeling and fitting of surfaces. We have used the Powell-Sabin (PS) splines for the approximation of elliptic partial differential equations (including Grad-Shafranov) in a rectangular domain. In this context have recovered the optimal rate of convergence (order 3). Preliminary result has been obtained for hyperbolic isothermal 2D Euler equations with TG/VMS stabilization. Our aim in the coming years is to apply these PS splines to full MHD in a toroidal geometry.

Genuinely multidimensional Riemann Solver

Participants : Jeaniffer Vides, Boniface Nkonga.

Multidimensional Riemann solvers were pioneered by Abgrall. Abgrall, Maire, Nkonga, Després and Loubere have extensively developed them especially as node-solvers for Lagrangian hydrodynamics. Another strain of work comes from explorations by Wendroff and Balsara who took a space-time approach. In this work, the resolved state is obtained via space-time integration over a wave model, just as was done by Wendroff and Balsara. However, an algebraic approach is used for the development of the fluxes. It is, therefore, shown that the multidimensional fluxes can be obtained by application of jump conditions at the boundaries of the wave model. The problem is of course over determined with the result that the shock jump conditions are only satisfied approximately in a least squares sense. Even so, this work gives us new perspective on multidimensional Riemann solvers. The litteral satisfaction of the shock jump conditions (up to least squares approximation) makes it easier to understand multidimensional Riemann solvers as a natural extension of the one-dimensional Riemann solvers. Contributions have also been made on the development of a minimalist wave model, which might help in reducing dissipation. Further innovations are reported on the assembling of fluxes based on the structure of the wave model, and those innovations are potentially useful. For MHD the CT approach consists of constraining the transport of magnetic field so that the divergence is always kept zero. The method relies on exploiting the dualism between the flux components and the electric field. Since the electric field is needed at the edges of the mesh, the multidimensional Riemann solver can also provide the electric field. By running an extensive set of simulations, it is shown that the multidimensional Riemann solver is robust and can be used to obtain divergence-free formulations for MHD that perform well on several stringent calculations. Future work will improve this strategy by enriching the description of the strongly interaction of waves.

Multi scales approximations of “Shallow water” flows

Participants : Jeaniffer Vides, Boniface Nkonga.

The terminology "Shallow water" is used to characterize thin flows on curved surfaces. It is customary for this type of flows to use the incompressible Navier-Stokes equations to asymptotically derive reduced models for the evolution of the depth integrated speed and the thickness of the flow. Reduced model are mainly hyperbolic and finite volume method are often used for their numerical approximation. Approximations strategies are generally structured as follow:

• Construction of a global coordinate system associated with an assumed analytical surface.

• Reduction of the model relatively to the global coordinate system

• Approximation of the surface by a finite number of elements.

• Approximation of the reduced model using the discrete surface.

In the context of real applications, it is presumptuous to expect an analytical formulation of the surface. From the data provided by observation satellites, we can usually extract a discrete description of the surfaces that drives thin flow. Therefore, it is more practical to use the discrete description as the starting point of the resolution strategy. This is the angle of approach that we have considered. We locally define two mesh scales: the element scale and the cell scale. The discrete mapping and the reduced model are defined at the element scale and the average values that evolve in time are defined at the cell scale. First applications have been successfully performed. We will now continue your investigations and include relevant physics at each scale, including sheared flows. We will also examine the use of multi dimensional Riemann solver in this context.

Computational Magnetohydrodynamics with Discrete Differential Forms

Participants : Holger Heumann, Ralf Hiptmair [SAM, ETH Zürich, Switzerland] , Cecilia Pagliantini [SAM, ETH Zürich, Switzerland] .

Differential forms, or equivalently exterior calculus, are a natural framework for electromagnetics; not only for a better understanding of the theoretical foundation, but also for the development of numerical methods. Keywords are the Hodge decomposition or the de Rham complex that are at the bottom of recent development of efficient multigrid methods or stable mixed finite element methods. Thinking in terms of such co-ordinate free differential forms offers considerable benefits as regards to the construction of structure preserving spatial discretizations.

In the present project, we aim at developing a new approach for the numerical treatment of resistive magnetohydrodynamics where a Galerkin discretization of the electromagnetic part based on finite element exterior calculus (FEEC) will be coupled to advanced finite volume methods for the approximation of the balance laws for the fluid.

The latest developments involved the extension and analysis of the stablized Galerkin schemes for advection of differential forms introduced in [48] to the case of time-dependent and non-regular flow fields.

Hamilton-Jacobi Formulation for Vlasov-Poisson

Participants : Holger Heumann, Eric Sonnendrücker [IPP, Max-Planck-Institute Garching, Germany] , Philip J. Morrison [Institute for Fusion Studies, Austin, USA] .

The phase space mapping induced by the solution of the Vlasov-Poisson problem is a symplectic mapping (or canonical transformation in physics literature) solving Hamilton’s equations. In this project we are developing numerical methods that are based on this formulation. We derived and implemented new finite difference schemes for the corresponding Hamilton-Jacobi equation, that circumvent the projection of the distribution function inherent in Lagrangian methods. First numerical results for standard test problems show the ability of increased resolution of fine-scale effects.

Entropy viscosity technique

Participants : Richard Pasquetti, Jean-Luc Guermond [Texas A & M University] , Boyan Popov [Texas A & M University] .

The entropy viscosity technique allows to address hyperbolic equations by introducing a strongly non linear viscous term where needed, especially at shocks. The basic idea is to set up a viscosity from the residual of the entropy inequality together with a $O\left(h\right)$ upper bound proportional to the local wave speed. In view of addressing situations where vaccum may appear in the tokamak, we have considered the shallow water equations with topography and in situations where dry-wet transitions occur. Using a RK scheme in time and a spectral element method (SEM) in space, we have proposed a variant of the entropy technique, that mainly consists of using the viscosity upper bound in the dry parts, to obtain satisfactory results. This work was presented in [30] , [22] and a publication is under review.

Bohm boundary conditions

Participants : Sébastian Minjeaud, Richard Pasquetti.

In the frame of the ANR project ESPOIR, our partners have proposed a penalty method to enforce the Bohm criterion (Mach number greater than one at the tokamak plates). This approach has been justified by considering a “minimal transport model” that consists of a 1-dimensional non linear hyperbolic system of two equations, that govern the evolutions of the density and velocity. The approach and further developments are described in three recent papers published in the journal of computational physics. Considering the same hyperbolic system, we have proposed a direct way to enforce the Bohm criterion in the frame of an explicit time marching. Using a SVV stabilized SEM it is then possible to resolve the same problem with spectral accuracy. This paper is now in press and will be published as a JCP note.

A numerical scheme for fluid-particules flows

Participants : Florent Berthelin, Thierry Goudon, Sebastian Minjeaud.

We propose a numerical scheme for the simulation of fluid-particles flows with two incompressible phases. The numerical strategy is based on a finite volume discretization on staggered grids, with a flavor of kinetic schemes in the definition of the numerical fluxes. We particularly pay attention to the difficulties related to the volume conservation constraint and to the presence of a close-packing term which imposes a threshold on the volume fraction of the disperse phase. We are able to identify stability conditions on the time step to preserve this threshold and the energy dissipation of the original model. The numerical scheme is validated with the simulation of sedimentation flows.

Identification and forecast of ionospheric disturbances

Participants : Didier Auroux, Sebastian Minjeaud.

In the framework of ANR IODISSEE, in order to identify (and forecast) ionospheric disturbances leading to temporary losses of satellite-to-earth communications (GPS, Galileo), we used Striation software for data assimilation. We obtained the adjoint code thanks to automatic differentiation (Tapenade software from Inria). As the data from Demeter satellite were not available, we extracted synthetic data from a generic model run, and we tried to identified some physical parameters (electronic density, atomic mass, number of particles) of the initial condition from the observations. For a small physical time scale (approximately 1 hour), the identification works very well, and it is possible to retrieve the initial condition from a sparse and noisy observations, allowing us to forecast the evolution of the ionospheric plasma - and then to forecast the disturbances and plasma bubbles that trap GPS and Galileo signals. For longer physical time windows (5 to 10 hours), the identification does not work anymore. We plan to work with real data, if possible, and also with a more complex model (for instance Dynamo software).