1 Introduction
Lagrangian methods for compressible multimaterial 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 LagrangianEulerian (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 TargetMatrix Optimization Paradigm Knupp2012 ; ETHOS2018 . We incorporate the TMOP technology in the context of a highorder 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 highorder 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 lowlevel 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. Reconnectionbased 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 targetmatrix construction methods that permit a greater degree of control over the optimized mesh. Pointwise mesh quality metrics are defined by utilizing subzonal 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 targetmatrices as a way for the user to incorporate applicationspecific 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 multimaterial 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 multimaterial 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 MultiMaterial ALE Framework
BLAST solves the multimaterial Euler equations by clearly separating the Lagrange, remesh and remap phases. It utilizes a singlefluid multimaterial 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 highorder 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 highorder 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 highorder element :
(1) 
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(2) 
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 pointwise 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 subzonal closure model, see Dobrev2016 . Before each remap procedure, these quadrature point values are projected to discontinuous functions, which are remapped, by solving an advectionbased 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 targetmatrix , which is the Jacobian of the transformation from the reference to the target coordinates. This matrix is defined according to a userspecified 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 spacedependent functions of the nodal coordinates. As any Jacobian matrix, can be written as a decomposition of matrices that represent its main geometric properties:
(3) 
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 simulationbased mesh adaptivity.
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:
(4) 
The rightmost 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
(5) 
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 righthand 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 rightmost 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 spacedependent, 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.
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
(6) 
where is a spacedependent 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 problemspecific 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:

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

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

. 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 .

has two welldefined derivatives with respect to . This enables us to use a standard Newtonbased nonlinear solver.

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:
(7) 
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 .
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 shapeonly optimization. To preserve invariance under mesh refinement, we set to
(8) 
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 equallysized targets with the shape+size metric . Node displacements are restricted by using and . We observe similar behavior for the different mesh resolutions.
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 
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 spacedependent 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 thirdorder mesh, , and ideal equallysized targets, but in this case we vary the function in space through
We observe that the chosen metric is active only in the nonrestricted region.
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 TMOPbased framework adapts the different geometric properties (3) of the highorder 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 radaptivity 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, radaptivity is achieved by incorporating the discrete data into the target Jacobian matrices. This process can be split into five major steps:

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 apriori error estimates.

Derivation of geometric parameters from relevant simulation data.

Definition of pointwise target Jacobian matrices .

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

Mesh optimization by solving the final nonlinear optimization problem.
The first two steps are simulation and problemdependent, 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
(9) 
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 steadystate Navier Stokes equations.
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 advectionbased transfer as described in Section 4.2 of IMR2018 . Another option is highorder interpolation between meshes, see Section 2.3 of mittal2019nonconforming and Section 4.1 of IMR2018 .
Gradients of adaptive target matrices
As we use gradientbased 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 spacedependent 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 simulationdriven case. Similar to constructing targetmatrices, 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 targetmatrices , 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
(10) 
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 mfem.org. The applicationspecific 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
(11) 
where is a userspecified 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 :
(12) 
The calculation of is performed once on the Lagrangian mesh .
5.1 Triple Point
Our first test is the 3material 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 highorder 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.
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 
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.
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 multimaterial 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.
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 zdirection 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, highorder, 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 .
In Figure 11 we show results of the full 3D simulation on a highorder 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 .
6 Conclusion
In this paper we have presented algorithms for simulationdriven optimization of highorder curved meshes, with application to highorder arbitrary LagrangianEulerian 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 multimaterial ALE simulations.
In memoriam
This paper is dedicated to the memory of Dr. Douglas Nelson Woods (January 11^{th} 1985  September 11^{th} 2019), promising young scientist and postdoctoral 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.
Disclaimer
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.
References
 (1) P. Knupp, Introducing the targetmatrix paradigm for mesh optimization by node movement, Engr. with Comptr. 28 (4) (2012) 419–429.
 (2) BLAST: Highorder finite element Lagrangian hydrocode, www.llnl.gov/CASC/blast.
 (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 Galerkintype 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 highorder cellcentered Lagrangian scheme for twodimensional 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 LagrangianEulerian 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 LagrangianEulerian methods for modeling highspeed compressible multimaterial flows, J. Comput. Phys. 322 (2016) 603–665.
 (10) P.H. Maire, J. Breil, S. Galera, A cellcentred arbitrary Lagrangian–Eulerian (ALE) method, Int J. Num. Meth. Fluids 56 (2004) 1161–1166.
 (11) S. Galera, P.H. Maire, J. Breil, A twodimensional unstructured cellcentered multimaterial ALE scheme using VOF interface reconstruction, J. Comput. Phys. 229 (16) (2010) 5755–5787.
 (12) V. Dobrev, T. Kolev, R. Rieben, Highorder 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, Multimaterial closure model for highorder 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, Highorder multimaterial 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 TargetMatrix Optimization Paradigm for highorder meshes, SIAM J. Sci. Comp. 41 (1) (2019) B50–B68.
 (16) V. Dobrev, T. Ellis, T. Kolev, R. Rieben, Highorder curvilinear finite elements for axisymmetric Lagrangian hydrodynamics, Computers & Fluids (2013) 58–69.
 (17) J.L. Guermond, B. Popov, V. Tomov, Entropyviscosity 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 arbitraryLagrangianEulerian ADERWENO 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 ArbitraryLagrangianEulerian 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 reconnectionbased arbitraryLagrangian–Eulerian method, J. Comput. Phys. 229 (12) (2010) 4724–4761.
 (24) V. A. Dobrev, P. Knupp, T. V. Kolev, V. Z. Tomov, Towards SimulationDriven Optimization of HighOrder Meshes by the TargetMatrix Optimization Paradigm, Springer International Publishing, 2019, pp. 285–302.
 (25) P. Knupp, Target formulation and construction in mesh quality improvement, Tech. Rep. LLNLTR795097, 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. AparicioEstrems, A. GargalloPeiró, 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 schwarzspectral element methods for incompressible flow, Computers & Fluids 191 (2019) 104237.
 (35) MFEM: Modular finite element methods library, mfem.org. 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 interfaceaware subscale dynamics closure model for multimaterial cells in Lagrangian and arbitrary LagrangianEulerian hydrodynamics, J. Comput. Phys. 276 (0) (2014) 92–135.
 (38) B. P. Howell, G. J. Ball, A freeLagrange augmented Godunov method for the simulation of elasticplastic solids, J. Comput. Phys. 175 (2002) 128–167.
 (39) V. A. Dobrev, T. V. Kolev, R. N. Rieben, Highorder curvilinear finite elements for elasticplastic Lagrangian dynamics, Journal of Computational Physics 257 (B) (2014) 1062–1080.
Comments
There are no comments yet.