Simulation-Driven Optimization of High-Order Meshes in ALE Hydrodynamics

by   Veselin Dobrev, et al.

In this paper we propose tools for high-order mesh optimization and demonstrate their benefits in the context of multi-material Arbitrary Lagrangian-Eulerian (ALE) compressible shock hydrodynamic applications. The mesh optimization process is driven by information provided by the simulation which uses the optimized mesh, such as shock positions, material regions, known error estimates, etc. These simulation features are usually represented discretely, for instance, as finite element functions on the Lagrangian mesh. The discrete nature of the input is critical for the practical applicability of the algorithms we propose and distinguishes this work from approaches that strictly require analytical information. Our methods are based on node movement through a high-order extension of the Target-Matrix Optimization Paradigm (TMOP) of Knupp. The proposed formulation is fully algebraic and relies only on local Jacobian matrices, so it is applicable to all types of mesh elements, in 2D and 3D, and any order of the mesh. We discuss the notions of constructing adaptive target matrices and obtaining their derivatives, reconstructing discrete data in intermediate meshes, node limiting that enables improvement of global mesh quality while preserving space-dependent local mesh features, and appropriate normalization of the objective function. The adaptivity methods are combined with automatic ALE triggers that can provide robustness of the mesh evolution and avoid excessive remap procedures. The benefits of the new high-order TMOP technology are illustrated on several simulations performed in the high-order ALE application BLAST.



There are no comments yet.


page 19

page 21


hr-adaptivity for nonconforming high-order meshes with the target matrix optimization paradigm

We present an hr-adaptivity framework for optimization of high-order mes...

Non-Conforming Mesh Refinement for High-Order Finite Elements

We propose a general algorithm for non-conforming adaptive mesh refineme...

Adaptive Surface Fitting and Tangential Relaxation for High-Order Mesh Optimization

We propose a new approach for controlling the characteristics of certain...

Matrix optimization on universal unitary photonic devices

Universal unitary photonic devices are capable of applying arbitrary uni...

A simple MATLAB program to compute differentiation matrices for arbitrary meshes via Lagrangian interpolation

A MATLAB program for computing differentiation matrices for arbitrary on...

A variational approach to high-order r-adaptation

A variational framework, initially developed for high-order mesh optimis...

Singularity Structure Simplification of Hexahedral Mesh Weighted Ranking

In this paper, we propose an improved singularity structure simplificati...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Lagrangian methods for compressible multi-material shock hydrodynamics VonNeumann1950 ; Benson1992 ; Loubere2004 ; Scovazzi2007 ; Maire2009 are characterized by a computational mesh that moves with the material velocity. A key advantage of using such methods is that they do not produce any numerical dissipation around contact discontinuities and material interfaces, when these are aligned with the element boundaries of the moving mesh. This is a consequence of the fact that the dynamics equations, when written in a Lagrangian frame, do not include nonlinear convective terms. However, the main disadvantage of Lagrangian methods is that they lead to mesh distortion, which results in small time steps or simulation breakdowns. Many problems of interest, i.e., impact simulations, cannot be run to longer evolution times by a purely Lagrangian method. These disadvantages are addressed by the Arbitrary Lagrangian-Eulerian (ALE) approach Hirt1974 ; Barlow2016 ; Maire2007 ; Galera2010 , which involves the notion of periodic mesh optimization. However, while improving the mesh, ALE methods can produce their own numerical errors around the material interfaces, thus counteracting the main advantage of the Lagrangian approach.

In this paper we propose mesh optimization methods that provide mechanisms to control the amount of numerical error resulting from the ALE approach. Keeping unchanged the method’s underlying solution transfer (remap) and Lagrangian phases, we explore how mesh adaptivity through node movement (a.k.a. -adaptivity), adaptive ALE triggers, and restricted node movement can be used to decrease the errors resulting from the overall ALE procedure. This paper extends our previous work on ALE hydrodynamics in the BLAST code Dobrev2012 ; Dobrev2016 ; Dobrev2018 and mesh optimization by the Target-Matrix Optimization Paradigm Knupp2012 ; ETHOS2018 . We incorporate the TMOP technology in the context of a high-order finite element ALE approach for multimaterial flow and introduce several important notions, including limiting of node movement based on the simulation dynamics, proper normalization of the energy functional, adaptivity to discrete simulation features, and adaptive ALE triggers. We focus on high-order curved meshes which have recently shown several mathematical and computational advantages, including optimal convergence rates on domains with curved boundaries/interfaces, symmetry preservation in radial flow Dobrev2012 ; Dobrev2013

and equivalent simulation quality with a smaller number of degrees of freedom

Dobrev2012 ; Tomov2016 ; Boscheri2016 .

Although in this paper we consider a particular application (BLAST), an important aspect of the TMOP approach is its generality. In particular, our methods are fully algebraic and do not involve any low-level geometric manipulations. They only require information about the local Jacobians, and thus can be applied to meshes of any order (including linear meshes), element type, and dimension. The proposed methods can be applied to both monolithic ALE methods Scovazzi2016 ; Bakosi2017 ; Gaburro2017 to define their mesh velocity, and ALE methods that split the Lagrangian, mesh optimization and remap phases. They can also be applied to diffused interface methods, where ALE errors are decreased by shrinking the transition region of the volume fraction functions, or exact interface representation methods Kucharik2010 , where ALE errors can be decreased by obtaining better mesh resolution in the interface regions.

The algorithms presented in this paper rely strictly on node movement and preserve the original topology of the mesh. Reconnection-based remesh approaches are more powerful in terms of controlling the mesh characteristics, but also require more complicated solution transfer procedures, see Loubere2010 .

Our mesh optimization method is based on the TMOP framework described in ETHOS2018 ; IMR2018 ; knupp2019target . This approach is distinguished from similar methods by its emphasis on target-matrix construction methods that permit a greater degree of control over the optimized mesh. Pointwise mesh quality metrics are defined by utilizing sub-zonal information obtained by sampling element Jacobians at element quadrature points. These metrics are capable of measuring shape, size or alignment of the region around the point of interest. TMOP requires predefined target-matrices as a way for the user to incorporate application-specific physical information into the metric that is being optimized. The combination of targets and quality metrics is used to optimize the node positions, so that they are as close as possible to the shape, size and/or alignment of their targets.

The rest of the paper is organized as follows. In Section 2 we review the multi-material ALE hydrodynamics framework in BLAST and the basic TMOP components. Section 3 describes the TMOP optimization functional in further detail, along with the concepts of limited node motion and proper normalization of the resulting nonlinear functional. Mesh adaptivity to discrete simulation features and its relation to dynamic ALE triggers is discussed in Section 4. Section 5 presents several numerical tests that demonstrate the features of the discussed methods. Conclusions are presented in Section 6.

2 Preliminaries

In this section we provide a brief overview of the multi-material ALE framework (as implemented in the BLAST code) that is used to test the proposed mesh optimization methods. We focus only on aspects that are relevant to our mesh optimization purposes, while a complete description of our ALE hydrodynamics methods can be found in Dobrev2018 . We also provide a summary of the basic TMOP components in the context of BLAST.

2.1 Multi-Material ALE Framework

BLAST solves the multi-material Euler equations by clearly separating the Lagrange, remesh and remap phases. It utilizes a single-fluid multi-material description of the system where all materials share a common velocity, but each has its own density and energy. BLAST supports performing an arbitrary number of Lagrangian steps between two ALE steps (remesh+remap), as opposed to continuous remap algorithms that perform an ALE step after every Lagrangian step.

Space discretization

The space discretization is based on high-order finite elements. Velocity is discretized in the finite element space , with basis , where is the computational mesh corresponding to the domain of interest , and is the space dimension. Specific internal energies are discretized in , with basis , where is the material index. Material densities are evolved at certain quadrature points of interest through the notion of pointwise mass conservation. Throughout this manuscript we refer to the pairs of spaces , by which we denote , i.e., the Cartesian product of the space of continuous finite elements on quadrilateral or hexahedral meshes of degree , and , the companion space of discontinuous finite elements of order one less than the kinematic space. As shown in Dobrev2012 the specific choice of corresponds to traditional staggered grid discretizations under additional simplifying assumptions.

Discrete mesh representation

To obtain a discrete representation of the high-order mesh, we start with the set of kinematic scalar basis functions on a reference element . The shape of any element in the mesh is then fully described by a matrix of size whose columns represent the coordinates of the element control points (a.k.a. nodes or element degrees of freedom). Given , we introduce the map whose image is the high-order element :


where we used to denote the -th column of

. To ensure continuity between mesh elements, we define a global vector

of mesh positions that contains the control points for every element. For any element in the mesh, we can compute the Jacobian of the mapping at any reference point as


We will refer to the initial mesh, which is to be optimized, as the Lagrangian mesh. We denote the Lagrangian mesh by .

Material representation and evolution

Material positions are represented on a discrete level as point-wise volume fractions at certain quadrature points of interest. At any such quadrature point , the volume fraction values satisfy and . During the Lagrange phase, these volume fraction are evolved by a sub-zonal closure model, see Dobrev2016 . Before each remap procedure, these quadrature point values are projected to discontinuous functions, which are remapped, by solving an advection-based PDE, to the new mesh. Once the remap step is complete, the remapped

functions are interpolated back to the quadrature points of interest. The above ALE treatment of the material positions leads to numerical diffusion around the material interfaces, which is often the main source of numerical dissipation error in BLAST.

2.2 Overview of TMOP

Mesh optimization in TMOP is driven by moving the mesh nodes (i.e. changing the values of the vector ) in a way that minimizes a mesh quality metric, , over the computational mesh. In this section we explain the meaning of and . A more detailed discussion about the TMOP optimization functional is presented in Section 3.

The weighted Jacobian matrix represents the transformation from a target configuration to the current physical (or active) mesh positions, so that a perfect mesh would result in . At each quadrature point, the matrix is determined using two Jacobian matrices:

  • The target-matrix , which is the Jacobian of the transformation from the reference to the target coordinates. This matrix is defined according to a user-specified method prior to optimization; it defines the desired local properties of the optimal mesh.

  • The Jacobian matrix of the transformation between reference and physical space coordinates, always computed by (2).

Then we have . Note that both and are space-dependent functions of the nodal coordinates. As any Jacobian matrix, can be written as a decomposition of matrices that represent its main geometric properties:


The concept of target construction allows the user specify some of the above components, and thus control the desired geometric property knupp2019target . We utilize this concept further in Section 4 to achieve simulation-based mesh adaptivity.

Figure 1: Schematic representation of the major TMOP matrices.

Given a matrix , the purpose of the local quality metric is to measure the differences between and only in one or more particular components of (3), while being invariant of the others. The choice of quality metric is based on the geometric aspects that the user wants to control. For example, is a shape metric that controls skew and aspect ratio, but is invariant to orientation and volume. Here, is the Frobenius norm of and . Similarly, and are shapesize metrics that controls volume, skew and aspect ratio, but are invariant to orientation. The reader is referred to ETHOS2018 for a more detailed discussion on quality metrics.

3 Recent Improvements in TMOP

In this section we go deeper into the details of the TMOP’s nonlinear objective function and present several improvements over our previous formulation presented in ETHOS2018 ; IMR2018 . These modifications are motivated by practical concerns related to restricting node movement and proper normalization of the different terms in the objective function.

3.1 Objective Function

TMOP optimizes the mesh quality by minimizing a global objective function that depends on the local quality measure throughout the mesh. This minimization is performed by solving with Newton’s method. The objective function has the following general form:


The right-most term is used to limit the node displacements during optimization. This term is discussed in Section 3.2. The other term represents an explicit combination of mesh quality metrics . Although most problems require only one metric (), in other cases mittal2019mesh it might be necessary to combine metrics with different weights in order to make some of them more dominant in certain regions of the domain. The integrals in (4) are computed as


where is the current mesh, is the target element corresponding to the physical element , are the quadrature weights, and the point is the image of the reference quadrature point in the target element . Note that the right-hand side in (5) depends on the mesh positions through the Jacobian matrices used in the definition of . The integration in (5) is performed over the target elements, enforcing that the integral contribution from a given element , relative to the contributions from other elements, is only based on the difference with its target (which is measured by ), and not on its relative size compared to the other elements. In particular, very small elements will not be neglected by the optimization process due to their size. The existence of a minimum for variational optimization problems like (5) has been explored theoretically in Garanzha2014 ; Garanzha2010 .

In general, the various quality metrics have different magnitudes; we have observed orders of magnitude differences in the resulting integrals. The role of the denominators in the first sum of (4) is to normalize the magnitude of the different metrics. The sum of the metric terms is always equal to at the start of the optimization process and is expected to decrease during the iteration. Furthermore, all metric integrals start at the same value of . This normalization also makes the magnitude of each metric integral invariant of , the level of refinement, the specifics of the problem (e.g. the units system), and the size of the domain.

3.2 Limiting of Node Displacements

One advantage of the Lagrangian approach is that it provides a natural form of mesh adaptivity. For example, the mesh spacing usually gets compressed around shock regions. Thus, it might be beneficial to have the ability to preserve, up to some extent, the positions of the Lagrangian mesh . Another important benefit of staying close to the Lagrangian mesh is minimizing errors due to remap. The purpose of the limiting term (the right-most term in in (4)) is to restrict the movement of nodes during the mesh optimization procedure with respect to their initial positions . Note that the limiting term is space-dependent, allowing the preservation of mesh features (e.g. initial resolution) in specific regions of the domain, which is an important requirement for various applications mittal2019mesh . A simple motivational 2D example is shown in Figure 2, where the goal is to preserve the resolution around the central boundary layer, while achieving uniform mesh resolution in the rest of the domain.

Figure 2: Original mesh (left) and optimized mesh (right) through space-dependent restriction of the node displacements.

In what follows, we rely on the fact that there exists a bijective map between every position , in the current mesh , and its initial position in the Lagrangian mesh , which is the case in methods based on node movement. We define the limiting term from (4) as a general function of distance between the current positions and their corresponding Lagrangian positions as


where is a space-dependent scalar function that specifies the acceptable amount of displacement for each mesh node. Note that moving a given position changes the displacement , but it does not change its image on the initial mesh and the value of . The values of are problem-specific and must be specified by the user. For example, in the remesh phase of ALE simulations, can be some fraction of the displacement that occurs between the last two ALE steps.

To achieve the limiting effect, the function has the following properties:

  1. when , and in particular, the initial value of the limiting integral is negligible.

  2. is an increasing function with respect to , i.e., the limiting term grows for bigger displacements.

  3. . This property is used to give the users a physical intuition of the meaning of , i.e., the limiting term will start to dominate the objective function whenever the mesh nodes move to a distance of , or bigger, from their original positions. Note that, due to their normalization, the metric integrals of (4) are always expected to be in .

  4. has two well-defined derivatives with respect to . This enables us to use a standard Newton-based nonlinear solver.

  5. is unitless, i.e., it is invariant to rescaling of the domain or change of the physics units system.

The rate at which the limiting term starts to influence the objective function depends on the particular form of . Our tests usually utilize one of the following two formulas:


As shown in Figure 3, the exponential option becomes active, and starts to grow much faster, only after mesh displacement approaches the maximum allowed value of .

Figure 3: Graphs of (red) and (black) for .

The constant in equation (6) provides proper normalization. In particular, note that the integrals in (6) are over target elements, which may or may not have the notion of size. For example, it is not necessary to set the local volume of a target element which will be used to perform shape-only optimization. To preserve invariance under mesh refinement, we set to


where is the total number of elements in the mesh, and is the total volume of the domain.

A simple 2D example that demonstrates the desired behavior of the normalization factors under mesh refinement is presented in Figure 4 and Table 1. In this example we optimize a third order mesh towards ideal shape equally-sized targets with the shape+size metric . Node displacements are restricted by using and . We observe similar behavior for the different mesh resolutions.

Figure 4: Perturbed initial mesh (left) and optimized meshes using and with for 0,1,2 refinements of the original mesh
Refinements Final Metric part Limiting part Max displacement
0 0.8423 0.7819 0.0604 0.0083
1 0.8422 0.7819 0.0602 0.0083
2 0.8422 0.7819 0.0602 0.0083
Table 1: Results for final functional value and its components for several refinements of the ICF mesh, using and with .

Local mesh optimization

A standard approach to perform local mesh optimization is to eliminate unknowns on the linear algebra level, i.e., to apply the Newton solver on a submesh. When the region of interest changes dynamically in ALE simulations, however, refactoring the linear algebra problem at every remesh step can be cumbersome. Another application of the space-dependent coefficient is that it can be used to perform local mesh optimization. This is demonstrated by another simple 2D test in Figure 5. Again we use a third-order mesh, , and ideal equally-sized targets, but in this case we vary the function in space through

We observe that the chosen metric is active only in the non-restricted region.

Figure 5: Perturbed initial mesh (left) and locally optimized meshes (0 and 2 refinements of the original mesh) using space-varying coefficient .

Combination with ALE Hydrodynamics

As already mentioned, the above methods to restrict the node displacements, or perform local mesh optimization, are useful when one wants to preserve a certain feature that is available in the Lagrangian mesh. Our experience shows that it is a good practice to set to be some fraction of the displacement that occurs between the last two ALE steps, which we denote by . This ties the amount of mesh relaxation to the amount of Lagrangian motion, and in particular, guarantees that nodes that do not move during the Lagrange phase (), do not move in the remesh phase either. Various modifications can be performed to relax this restriction in certain problems, e.g., diffusing to obtain transition regions.

4 Mesh Adaptivity and ALE Triggers

Mesh adaptivity huang2010adaptive is beneficial when a certain dynamic feature is not resolved by the Lagrangian mesh. In this section we demonstrate how the TMOP-based framework adapts the different geometric properties (3) of the high-order mesh to control the numerical errors in ALE methods. We start by describing the coupling between target construction and discrete simulation data, followed by a discussion of the appropriate triggering mechanisms in the remesh phase of ALE simulations.

4.1 Adaptivity to Discrete Simulation Fields

The objective of the r-adaptivity process is to optimize the mesh using information from a discrete function, e.g., a finite element solution function that is evolved during the Lagrangian phase. In TMOP, r-adaptivity is achieved by incorporating the discrete data into the target Jacobian matrices. This process can be split into five major steps:

  1. Choice of adaptation goal, i.e., which geometric properties of the mesh must be controlled and what is their desired behavior. This can be motivated, for example, by known a-priori error estimates.

  2. Derivation of geometric parameters from relevant simulation data.

  3. Definition of pointwise target Jacobian matrices .

  4. Choice of a quality metric that measures differences in the geometric quantities of interest.

  5. Mesh optimization by solving the final nonlinear optimization problem.

The first two steps are simulation and problem-dependent, but once the geometric parameters are known, one can construct the targets and optimize the mesh automatically. Alternative approaches for linear meshes in the context of ALE hydrodynamics are described in Greene2017 ; Vachal2011 .

Next we provide a simple 2D example to illustrate the above procedure. Suppose we want to adapt the mesh size and aspect ratio to a material interface, and suppose that the material position is prescribed by a volume fraction function , as it is the case in BLAST, see Figure 6. The needed geometric parameters can be computed from the gradient of . The aspect ratio is computed by the ratio of the gradient components (i.e., ) and the local size depends on the magnitude of the gradient (i.e., ). To complete the definition of shape, we also choose the target skewness to be the same as that for an ideal element ( for a quadrilateral element). Using this approach, we construct the target matrix corresponding to (3) as


Next we choose the shape+size metric , which measures differences of size, aspect ratio, and skewness between target and physical configurations. This metric is invariant to orientation (note that sets targets that have no rotation, but this is irrelevant due to the choice of ). The final result for a third order mesh is presented in Figure 6 which shows that the mesh resolves the interface. Further details about target construction algorithms can be found in knupp2019target . Note that the above procedure can be applied to various kinds of simulation data coming from the Lagrangian phase, e.g., shock positions, high temperatures, pressure gradients, etc. Alternative adaptivity approach, for the case when is prescribed analytically, is presented in Roca2018 . We also stress the work presented in marcon2019rp , where the authors utilize an approach similar to TMOP to adapt curved meshes in the context of steady-state Navier Stokes equations.

Figure 6: Material indicator function on the initial uniform mesh (left panel) and adaptively optimized mesh (right panel).

Discrete data on intermediate meshes

Since the discrete solution fields are available only on the Lagrangian mesh , adaptive mesh optimization requires a method to evaluate these fields on the intermediate meshes that are obtained during the Newton iteration. All simulations in this paper use the advection-based transfer as described in Section 4.2 of IMR2018 . Another option is high-order interpolation between meshes, see Section 2.3 of mittal2019nonconforming and Section 4.1 of IMR2018 .

Gradients of adaptive target matrices

As we use gradient-based nonlinear solvers to update mesh positions, we have to compute derivatives of with respect to . Since depends on a discrete function of the mesh positions, the derivatives must be taken into account. The formulas for these derivatives and further discussion are presented in Section 3 of IMR2018 .

Combination of adaptivity and limiting of displacements

Since mesh adaptivity may require significant amounts of mesh displacements to achieve the desired geometric characteristics, including the limiting term (6) in the nonlinear functional can be counterproductive. Nevertheless, the limiting effects (e.g. local mesh optimization) can still be achieved by designing a proper space-dependent limiting function , which must be aware of the adaptivity goals.

4.2 Adaptive ALE Triggers

We have briefly discussed ALE triggers in Section 5 of IMR2018 . In this section we make several important improvements to the work presented in IMR2018 , along with a more detailed discussion.

Since BLAST can perform any number of Lagrangian steps between two ALE steps, an important question is when, or how often, to optimize the mesh. Ideally we want a trigger mechanism that does remesh often enough to avoid small computational time steps due to mesh distortion. At the same time we do not want to remesh too often, as ALE calculations can lead to increased numerical dissipation. In BLAST, this dissipation comes from transitions between primal and conservative variables and monotonicity treatment in the remap phase Dobrev2018 , which is a purely numerical process that is independent of the physics of the Euler equations. Furthermore, in the case of mesh adaptivity, higher number of remesh steps slow down the overall simulations as the mesh adaptivity process involves the extra computations explained in Section 4.1.

The above discussion suggests that the ALE trigger mechanism must be aware of the mesh quality during the Lagrangian phase, and induce an ALE step only when the mesh quality is below a certain threshold. It is important to note that mesh quality must be viewed in terms of the mesh quality metric used in the mesh optimization phase; the mesh optimization phase is expected to produce a mesh that is above the trigger threshold.

The TMOP concepts can be utilized to define such mesh optimization triggers, both in the purely geometric and the simulation-driven case. Similar to constructing target-matrices, the user specifies one or more admissible Jacobian matrix, , which defines the transformation from the reference element to the worst element that can be used during the Lagrangian phase. Note that is initialized in the beginning and stays constant throughout the simulation. Then represents the Jacobian of the transformation between target and admissible coordinates, which is used to calculate , the highest admissible mesh quality metric value for the metric of interest. With that, remesh is triggered whenever at any mesh quadrature point. In the case of mesh adaptivity, the target matrices adapt in time as they depend on some dynamic solution fields. Thus and change too, leading to a trigger that is adapted to the current Lagrangian solution.

For example, suppose we adapt small mesh size to high temperature regions. Then at a fixed point in space , the admissible metric value represents the admissible size deviation, given the current temperature at . The value represents the actual size deviation of the Lagrangian mesh, at the current temperature.

The specific formulation of is, of course, problem dependent, e.g., it can be used to trigger remesh steps whenever particular mesh configurations occur. Just as the target-matrices , the admissible matrices can contain information about the local shape, size or alignment, at any quadrature point. As a simple 2D example, we can define as


When is a shape metric, remesh is triggered based on local aspect ratio. When is a shape+size metric, remesh is triggered based on local size and aspect ratio.

ALE triggers can be combined with the use of the limiting term (6), as long as the limiting distances are large enough to allow the optimization process to achieve at all quadrature points.

5 Numerical tests

In this section we report numerical results from the algorithms in the previous sections as implemented in the BLAST code blast . The baseline TMOP methodology is part of the MFEM finite element library mfem ; MFEMPaper2019 , which is freely available at The application-specific target constructions, however, are only included in the BLAST implementation. All physical units in the following tests are based on the unit system.

In what follows, tests that adapt mesh size always utilize a size indicator function with values in , where values close to represent regions with smaller mesh size. Target matrices are constructed as


where is a user-specified size ratio between big and small local size. Unless otherwise specified, we use for all numerical tests. Such targets specify custom size, ideal shape, and zero rotation; in all tests we use metrics that are invariant of rotation. The small local size is approximated by taking into account the total volume of the domain, the volume occupied by , the number of available elements, and the specified :


The calculation of is performed once on the Lagrangian mesh .

5.1 Triple Point

Our first test is the 3-material Triple Point problem described in Dobrev2012 . This test is used to compare different remesh strategies and demonstrate their advantages and drawbacks in terms of accuracy and computational performance.

We start by performing a purely Lagrangian simulation to time 3.5 and record the final position of one of the materials, see Figure 7. We take this Lagrangian position as the true solution, as there is no diffusion of the material interface. For every other test, we define error by measuring differences with the Lagrangian result, for the same material’s position. These differences are measured by sampling both solutions at equidistant points over the domain. This sampling is a nontrivial procedure by itself, as it involves interpolation in physical space of high-order finite element functions on curved meshes; for this we use the methods from Section 2.3 of mittal2019nonconforming . The reported error is the average difference over all points. All simulations utilize a discretization s.t. the mesh is third order.

Table 2 lists all tested approaches. Most cases use a fixed frequency, so that an ALE step is performed every 50 Lagrangian steps. The limited tests optimize the mesh using the shape metric combined with the quadratic limiting function and , where or , see Section 3.2. The adapted tests optimize the mesh using the shape+size metric , so that material interface regions are assigned smaller mesh sizes. We report the total number of Lagrangian steps, ALE steps, remap advection steps (i.e., the total from within all ALE steps), the relative runtime (all tests use the same computer configuration), and the error in the material position due to numerical diffusion. The final material positions for several of the major options is shown in Figure 8.

Figure 7: Mesh and material position resulting from a purely Lagrangian simulation.

Figure 8: Meshes and material positions resulting from the Eulerian (left panel), limited by (center panel), and adapted with automatic trigger (right panel) simulations.
Approach # Lag # ALE # advect runtime error
Lagrangian 20808 0 0 615 0
Eulerian, period of 50 757 16 90 69 0.033
Limited , period of 50 980 20 59 64 0.022
Limited , period of 50 1505 31 55 86 0.020
Adapted, auto trigger 1666 8 137 174 0.017
Adapted, period of 50 1818 37 154 235 0.022
Table 2: Comparison between different ALE strategies for the Triple Point test.

We observe that the Eulerian approach, which always returns the mesh to its original configuration, is fast as it produces large Lagrangian time steps, but is also the most diffusive. Performing limiting with is a clear improvement over the Eulerian case, both in computational time and accuracy, because less advection remap steps are performed, which compensates the smaller size of the Lagrangian time steps. Decreasing the limiting distances by setting preserves the Lagrangian mesh more closely, resulting in smaller time steps and higher compute time. The error is slightly better, as less advection steps are performed.

Since the error in this test is always in the interface region, adapting the mesh size (with adaptive trigger) to this region is the most accurate approach. The automatic trigger is setup as in (10) with , taking into account the local size and shape of any quadrature point in the mesh. For example, an ALE would be triggered whenever the local size differs by a factor of of the required mesh size, or the aspect ratio differs by a factor of of the ideal aspect ratio. Mesh adaptivity, however, leads to slower computations as it leads to bigger mesh displacements in the remesh phase, which leads to more advection steps in the remap phase.

As discussed in Section 4.2 performing mesh adaptivity with too frequent ALE steps can lead to both slower and more diffusive simulations. This is demonstrated by the last case in Table 2, where an ALE step is performed every 50 steps, independent of the mesh quality during the Lagrangian phase.

5.2 2D ALE Simulation of Gas Impact

This test represents a high velocity impact of gasses that was originally proposed in Barlow14 . It involves three materials that represent an impactor, a wall, and the background. This problem is used to demonstrate the method’s behavior in an impact simulation that cannot be executed in Lagrangian frame to final time as it produces large mesh deformations.

The domain is with boundary conditions. The material regions are and for the impactor, for the wall, and the rest is background. The initial horizontal velocity of the impactor is . The problem is run to a final time of on a second order structured mesh. The complete thermodynamic setup of this problem and additional details about our multi-material finite element discretization and overall ALE method can be found in Dobrev2016 ; Dobrev2018 .

The goal of every ALE remesh step is to adapt the size of the mesh towards the locations of the impactor, the wall, and all material interfaces. The adaptivity field is constructed from the simulation’s discrete volume fractions and , and a reconstructed interface indicator function :

where the range of all functions is . Target Jacobians are constructed by (11), and the shape+size metric is used to optimize the mesh. An automatic ALE trigger is utilized by setting in (10).

The time evolution of the material positions and the corresponding mesh is demonstrated in Figure 9. This calculation required 49400 Lagrange time steps, 38 ALE steps, and 760 advection remap steps in total. We observe that the algorithm adapts well to the moving materials and the interface regions, while preserving good overall shape throughout the domain.

Figure 9: Time evolution of density and mesh position at times 0.5 (top left), 2 (top middle), 4 (top right), 6 (bottom left), 8 (bottom middle) and 10 (bottom right) in the 2D gas impact test case.

5.3 2Drz and 3D Full ALE Simulation of Steel Ball Impact on Unstructured Mesh

As a final example, we consider both 2Drz and full 3D calculations of a high velocity impact of a steel ball against an aluminum plate as described in HowellBall02 ; Dobrev2016 . The problem consists of a spherical steel projectile of radius and initial velocity in the z-direction of impacting a cylindrical plate of aluminum with a radius of and a thickness of . Both the steel and aluminum materials use a Gruneisen equation of state combined with an elastic perfectly plastic “strength” model as described in Dobrev2014 . This test demonstrates the ability of the method to adapt 2D and 3D unstructured meshes in the context of a practical ALE impact simulation.

In Figure 10 we show results of the 2D (axisymmetric) simulation on a sequence of refined 2D, high-order, unstructured meshes. Similarly to the previous example, we adapt to the position of the ball impactor, the aluminum plate, and all material interfaces, using the shape+size metric . The adapted trigger is setup with in (10). For the sequence of mesh refinements, the calculations required a total of , and Lagrange time steps for the coarse, medium and fine meshes, respectively, to reach a final simulation time of . The total number of remesh steps performed for the each 2Drz calculation is , and for the coarse, medium and fine meshes, respectively. The total number of remap advection steps is , and .

Figure 10: Density and mesh positions at times (top) and (bottom) using three different mesh resolutions for the 2Drz ball impact problem.

In Figure 11 we show results of the full 3D simulation on a high-order 3D unstructured mesh with resolution identical to the previous 2D medium case in plane. The adaptivity parameters for this calculation are and for the 3D version of (10). This calculation requires a total of Lagrange time steps, remesh steps, and advection remap steps to reach a final simulation time of .

Figure 11: Density and mesh positions at times (left) and (right) for the 3D ball impact problem.

6 Conclusion

In this paper we have presented algorithms for simulation-driven optimization of high-order curved meshes, with application to high-order arbitrary Lagrangian-Eulerian simulations blast . These allow the user to optimize the shape and size of the mesh, preserve geometric features, adapt the mesh characteristics to discrete fields of interest, and control the ALE frequency by adaptive triggers. These methods provide flexible control over the balance between computational time and amount of numerical dissipation in ALE simulations.

A future area of research will be to combine the existing TMOP technology with adaptive mesh refinement. This would allow to optimize shape, while achieving arbitrary small local size, which is currently constrained by the fact that the mesh topology is fixed. We also plan to add capabilities for tangential relaxation around discrete surfaces, which would be highly beneficial in multi-material ALE simulations.

In memoriam

This paper is dedicated to the memory of Dr. Douglas Nelson Woods (January 11th 1985 - September 11th 2019), promising young scientist and post-doctoral research fellow at Los Alamos National Laboratory. Our thoughts and wishes go to his wife Jessica, to his parents Susan and Tom, to his sister Rebecca and to his brother Chris, whom he left behind.


This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.


  • (1) P. Knupp, Introducing the target-matrix paradigm for mesh optimization by node movement, Engr. with Comptr. 28 (4) (2012) 419–429.
  • (2) BLAST: High-order finite element Lagrangian hydrocode,
  • (3) J. V. Neumann, R. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950).
  • (4) D. J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Comput. Methods Appl. Mech. Eng. 99 (1992) 235–394.
  • (5) R. Loubere, J. Ovadia, R. Abgrall, A Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems, Int J. Num. Meth. Fluids 44 (2004) 645–663.
  • (6) G. Scovazzi, M. Christon, T. Hughes, J. Shadid, Stabilized shock hydrodynamics: I. a Lagrangian method, Comp. Methods in Applied Mech. Eng. 196 (2007) 923–966.
  • (7) P.-H. Maire, A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes, J. Comput. Phys. 228 (7) (2009) 2931–2425.
  • (8) C. W. Hirt, A. A. Amsden, J. L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 277–253.
  • (9) A. J. Barlow, P.-H. Maire, W. J. Rider, R. N. Rieben, M. J. Shashkov, Arbitrary Lagrangian-Eulerian methods for modeling high-speed compressible multimaterial flows, J. Comput. Phys. 322 (2016) 603–665.
  • (10) P.-H. Maire, J. Breil, S. Galera, A cell-centred arbitrary Lagrangian–Eulerian (ALE) method, Int J. Num. Meth. Fluids 56 (2004) 1161–1166.
  • (11) S. Galera, P.-H. Maire, J. Breil, A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction, J. Comput. Phys. 229 (16) (2010) 5755–5787.
  • (12) V. Dobrev, T. Kolev, R. Rieben, High-order curvilinear finite element methods for Lagrangian hydrodynamics, SIAM J. Sci. Comp. 34 (5) (2012) 606–641.
  • (13) V. A. Dobrev, T. V. Kolev, R. N. Rieben, V. Z. Tomov, Multi-material closure model for high-order finite element Lagrangian hydrodynamics, Int. J. Numer. Meth. Fluids 82 (10) (2016) 689–706.
  • (14) R. W. Anderson, V. A. Dobrev, T. V. Kolev, R. N. Rieben, V. Z. Tomov, High-order multi-material ALE hydrodynamics, SIAM J. Sci. Comp. 40 (1) (2018) B32–B58.
  • (15) V. A. Dobrev, P. Knupp, T. V. Kolev, K. Mittal, V. Z. Tomov, The Target-Matrix Optimization Paradigm for high-order meshes, SIAM J. Sci. Comp. 41 (1) (2019) B50–B68.
  • (16) V. Dobrev, T. Ellis, T. Kolev, R. Rieben, High-order curvilinear finite elements for axisymmetric Lagrangian hydrodynamics, Computers & Fluids (2013) 58–69.
  • (17) J.-L. Guermond, B. Popov, V. Tomov, Entropy-viscosity method for the single material Euler equations in Lagrangian frame, Computer Methods in Applied Mechanics and Engineering 300 (2016) 402 – 426.
  • (18) W. Boscheri, M. Dumbser, High order accurate direct arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes on moving curvilinear unstructured meshes, Computers & Fluids 136 (2016) 48 – 66.
  • (19) X. Zeng, G. Scovazzi, A variational multiscale finite element method for monolithic ALE computations of shock hydrodynamics using nodal elements, J. Comput. Phys. 315 (2016) 577–608.
  • (20) J. Bakosi, J. Waltz, N. Morgan, Improved ALE mesh velocities for complex flows, Int J. Num. Meth. Fluids 85 (11) (2017) 662–671.
  • (21) E. Gaburro, M. Dumbser, M. Castro, Direct Arbitrary-Lagrangian-Eulerian finite volume schemes on moving nonconforming unstructured meshes, Computers and Fluids 159 (15) (2017) 254–275.
  • (22) M. Kucharik, R. Garimella, S. Schofield, M. Shashkov, A comparative study of interface reconstruction methods for multi‐material ALE simulations, J. Comput. Phys. 229 (7) (2010) 2432–2452.
  • (23) R. Loubère, P.-H. Maire, M. Shashkov, J. Breil, S. Galera, ReALE: A reconnection-based arbitrary-Lagrangian–Eulerian method, J. Comput. Phys. 229 (12) (2010) 4724–4761.
  • (24) V. A. Dobrev, P. Knupp, T. V. Kolev, V. Z. Tomov, Towards Simulation-Driven Optimization of High-Order Meshes by the Target-Matrix Optimization Paradigm, Springer International Publishing, 2019, pp. 285–302.
  • (25) P. Knupp, Target formulation and construction in mesh quality improvement, Tech. Rep. LLNL-TR-795097, Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States) (2019).
  • (26) K. Mittal, P. Fischer, Mesh smoothing for the spectral element method, Journal of Scientific Computing 78 (2) (2019) 1152–1173.
  • (27) V. Garanzha, L. Kudryavtseva, S. Utyuzhnikov, Variational method for untangling and optimization of spatial meshes, J. Comput. Appl. Math. 269 (2014) 24–41.
  • (28) V. A. Garanzha, Polyconvex potentials, invertible deformations, and thermodynamically consistent formulation of the nonlinear elasticity equations, Comp. Math. Math. Phys. 50 (9) (2010) 1561–1587.
  • (29) W. Huang, R. D. Russell, Adaptive moving mesh methods, Springer Science & Business Media, 2010.
  • (30) P. T. Greene, S. P. Schofield, R. Nourgaliev, Dynamic mesh adaptation for front evolution using discontinuous Galerkin based weighted condition number relaxation, Journal of Computational Physics 335 (2017) 664–687.
  • (31) P. Váchal, P.-H. Maire, Discretizations for weighted condition number smoothing on general unstructured meshes, Computers & Fluids 46 (1) (2011) 479 – 485.
  • (32) G. Aparicio-Estrems, A. Gargallo-Peiró, X. Roca, Defining a stretching and alignment aware quality measure for linear and curved 2D meshes, Springer International Publishing, 2019, pp. 37–55.
  • (33) J. Marcon, G. Castiglioni, D. Moxey, S. J. Sherwin, J. Peiró, -adaptation for compressible flows, arXiv preprint arXiv:1909.10973 (2019).
  • (34) K. Mittal, S. Dutta, P. Fischer, Nonconforming schwarz-spectral element methods for incompressible flow, Computers & Fluids 191 (2019) 104237.
  • (35) MFEM: Modular finite element methods library, doi:10.11578/dc.20171025.1248.
  • (36) R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J. Cerveny, V. Dobrev, Y. Dudouit, A. Fisher, T. Kolev, W. Pazner, M. Stowell, V. Tomov, J. Dahm, D. Medina, S. Zampini, MFEM: A modular finite element methods library, in review (2019). arXiv:1911.09220.
  • (37) A. Barlow, R. Hill, M. J. Shashkov, Constrained optimization framework for interface-aware sub-scale dynamics closure model for multimaterial cells in Lagrangian and arbitrary Lagrangian-Eulerian hydrodynamics, J. Comput. Phys. 276 (0) (2014) 92–135.
  • (38) B. P. Howell, G. J. Ball, A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids, J. Comput. Phys. 175 (2002) 128–167.
  • (39) V. A. Dobrev, T. V. Kolev, R. N. Rieben, High-order curvilinear finite elements for elastic-plastic Lagrangian dynamics, Journal of Computational Physics 257 (B) (2014) 1062–1080.