# A space and time fixed point mesh adaptation method

Bastien Sauvage, Frédéric Alauzet, Alain Dervieux

► To cite this version:

Bastien Sauvage, Frédéric Alauzet, Alain Dervieux. A space and time fixed point mesh adaptation method. Journal of Computational Physics, 2024, 519, pp.113389. 10.1016/j.jcp.2024.113389. hal-04715086

HAL Id: hal-04715086

https://inria.hal.science/hal-04715086v1

Submitted on 30 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A space and time fixed point mesh adaptation method

*Preprint, may 27, 2024*

Bastien Sauvage *

* corresponding author

Université Côte d'Azur, INRIA,

2004 Route des lucioles, F-06902 Sophia-Antipolis, bastien.sauvage@inria.fr

Frédéric Alauzet

Université Paris-Saclay- INRIA, 1, rue Honoré d'Estienne d'Orves F-91126

Palaiseau, frederic.alauzet@inria.fr

Alain Dervieux

Société technologique LEMMA, 2000 route des Lucioles, Sophia-Antipolis, France,
and Université Côte d'Azur, INRIA, 2004 Route des lucioles,
F-06902 Sophia-Antipolis, alain.dervieux@inria.fr
**Abstract:** We introduce a space-time metric-based approach for the best set of spatial meshes $M_{opt}(t)$ combined with the best set time step $	au_{opt}(t)$ of space-time complexity $N_{st}$ for the calculation of transient flows with implicit time advancing. Both types of error estimate, feature-based and goal-oriented, are considered for the compressible RANS equations. The case of a mesh which is adapted at each time step, and the case where the mesh is constant during a time subinterval of the whole simulation are theoretically analyzed. These space-time error estimates are then considered inside a Global Transient Fixed Point mesh adaptation algorithm. Applications to flow with vortex shedding past a cylinder are then described.

**Keywords:** compressible flow, mesh adaptation, anisotropic, space-time, implicit
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Space-time continuous mesh</td><td>5</td></tr><tr><td>3</td><td>Space-time error analysis</td><td>7</td></tr><tr><td>3.1</td><td>Feature-based error model for the scalar advection equation</td><td>7</td></tr><tr><td>3.2</td><td>Feature-based error model for CFD</td><td>8</td></tr><tr><td>3.3</td><td>Goal-oriented error model for CFD</td><td>9</td></tr><tr><td>3.4</td><td>Unified error model for CFD</td><td>10</td></tr><tr><td>4</td><td>Analysis for all-time adaptation</td><td>10</td></tr><tr><td>4.1</td><td>Analysis for all-time mesh adaptation at a given time</td><td>10</td></tr><tr><td>4.2</td><td>Analysis for all-time adaptation over time interval</td><td>11</td></tr><tr><td>5</td><td>Analysis with time subintervals</td><td>14</td></tr><tr><td>5.1</td><td>The Global Transient Fixed Point (GTFP) algorithm</td><td>14</td></tr><tr><td>5.2</td><td>Notations for space-time Global Transient Fixed Point</td><td>15</td></tr><tr><td>5.3</td><td>Spatial mesh optimization on a subinterval</td><td>16</td></tr><tr><td>5.4</td><td>Temporal optimization over the time subintervals</td><td>17</td></tr><tr><td>6</td><td>Global Space-Time Transient Fixed Point algorithm</td><td>20</td></tr><tr><td>7</td><td>Numerical experiments</td><td>20</td></tr><tr><td>7.1</td><td>2D flow past a cylinder at Reynolds 3900</td><td>22</td></tr><tr><td>7.2</td><td>2D flow past a cylinder at Reynolds 3900 with multi-mesh adaptation</td><td>27</td></tr><tr><td>7.3</td><td>2D flow past a cylinder at Reynolds 1M</td><td>30</td></tr><tr><td>7.4</td><td>3D flow past a cylinder at Reynolds 3900</td><td>34</td></tr><tr><td>7.5</td><td>3D flow past a cylinder at Reynolds 1M</td><td>37</td></tr><tr><td>8</td><td>Concluding remarks</td><td>40</td></tr><tr><td>9</td><td>Acknowledgements</td><td>41</td></tr></table>
# 1 Introduction

In the context of CFD calculations, the mesh generation for accurate and robust numerical simulations of Reynolds-Averaged Navier-Stokes (RANS) equations is a time-consuming task. This results from the fact that meshes are traditionally considered as an input to the simulation pipeline. The common practices require a manual refinement of the computational domain and an adjustment of the timestep length both based on the a priori knowledge of the solution. This means that mesh generation and timestep length choice mainly rely on the experience and intuition of a skilled engineer to predict the flow and to adapt the discretization to the flow. It is evident that this operation can be manually accomplished only for simple geometries or academic flows where the solution is known. Whilst, this is not the case for complex geometries as well as for flow conditions exhibiting complicated features (e.g., shocks, supersonic shear layers, separation, etc.). In these cases, traditional "good practice guidelines" leads to burden the simulation pipeline or even turn out to be unfeasible.

Metric-based mesh adaptation is an efficient framework to generate adapted anisotropic meshes under time-delivering constraints, when a study at different physical conditions is demanded. It can be done in a feature-based mode, relying on minimizing the interpolation error of one or several sensors (=features) in $L^p$ norm, or in a goal-oriented mode in which the error committed on a scalar output of a PDE is minimized with the use of an adjoint state. Both can be applied to steady calculations and to unsteady ones. We refer to [14, 15] for theoretical statements, and the monograph [9] for detailed descriptions of these methods. We consider the second case, that is the discretization of a PDE in the space-time domain $Q = \Omega \times [0, T]$. We advance in time and want to adapt the discretization to the solution. Two main options differ according to the spatial mesh adaptation.

(1) An important option is *mesh adaptation at each time step*, which consists in building at each time step a new adapted anisotropic mesh $\mathcal{H}(t)$ by defining a metric field $\mathcal{M}(t)$ on the computational domain $\Omega$, which is optimal for this time level. This option may be expensive in terms of computational cost and may be of low accuracy due to errors committed in transferring solutions between too many successive meshes. Therefore this option is interesting from a theoretical point of view but it is of no practical use.

(2) A second strategy consists in *freezing the adapted mesh during several time steps*. Then, it is mandatory that the mesh anticipates the flow behavior during these time steps, in other words during a given time subinterval. In [2, 4, 12] a transient fixed point (TFP) mesh adaptation algorithm has been proposed in order to master this issue. In the first version, the adaptation loop applies successively to each subinterval where the mesh is frozen. The error criterion is of *feature-based* type, measuring the interpolation error of a sensor chosen by the user. The TFP approach has been extended to a *goal-oriented* adaptation in [7] where several analyses were proposed for evaluating the convergence order of the TFP. In this context, the TFP was extended into the Global Transient Fixed Point (GTFP) in which the different meshes take into account a global space-time complexity evaluation and therefore are generated after the complete time resolution. This was mandatory as we had to solve the backward in time adjoint problem. A more complete accuracy and convergence analysis of GTFP, for a feature-based adaptation is proposed in [5] with many numerical examples.

In these works, either the time step is directly specified by the user, or the time step is assumed to be defined via a CFL stability condition related to an explicit time advancing. This is a useful option when an explicit time scheme is applied. In practice, it adapts quite well the time step to the solution as analyzed in [4]. Thanks to the CFL stability condition, the only unknown is the spatial mesh or the set of spatial meshes to perform the simulation.

The choice of a time step is defined in other terms when an implicit time advancing is used. Indeed, the size of the time step is no longer directly related to a stability condition. The time steps which are used can be notably larger than those permitted with an explicit time advancing. Large time steps induce a higher CPU efficiency, but the time approximation accuracy becomes an issue. Too large time
steps degrade the prediction, too small time steps increase the computational cost. In this context, the time step (or time discretization) becomes also an unknown of the mesh adaptation problem.

Many attempts to control the time step size on an accuracy basis can be found in the literature. In [19], the two components of time error, namely truncation and implicit iterative errors are evaluated and controlled. In [18], the authors use an adjoint based output sensibility and a division of elements and time steps, addressing the largest error with the fewest additional space-time elements. A similar approach is combined with a space-time AMR in [13]. Another proposal for combining time adaptation with AMR is presented in [10]. Papers [20, 21] relies on an *a posteriori* analysis and error equidistribution in space and in time. A context closer to our is addressed in [8] where the authors adapt separately time and space. While giving interesting results, choosing a separate adaptation of the time discretization and the space discretization implies limitations in the global accuracy/efficiency compromise in the calculation. For example, a too fine mesh is useless if the time step is too large, a small time step is useless if local mesh size is too large.

The novelty of this work is to present a fully coupled approach where the temporal and the spatial errors are linked, hence the spatial error is impacting the time error and vice versa. Therefore, the goal is to obtain directly the space-time discretization which minimizes an error model under the constraint of a prescribed space-time discretization complexity. This strategy is developed in the framework of the Global Transient Fixed Point (GTFP) mesh adaptation, and extends GTFP to space-time adaptation in a direct manner. For the feature-based approach, the problem of the optimal simultaneous adaptation of the spatial mesh and the time step can be formulated under the form of a severely nonlinear optimization problem where the spatial and the temporal errors are tightly coupled. To solve it, we propose a slight simplification where both errors are loosely coupled in order to design a tractable accurate and efficient time step and mesh adaptive GTFP algorithm. In the case of a goal-oriented criterion, we demonstrate that the proposed formulation applies in a natural way.

In the present paper, the numerical study is restricted to the feature-based approach for which we demonstrate the validity of the new mesh and time step adaptation algorithms. Section 2 sets the notion of space-time continuous mesh which is mandatory to perform the theoretical analysis. Section 3 gives the considered space-time error models: first in the context of the 1D scalar advection and, second, in the context of the compressible Navier-Stokes equation for the feature-based and goal-oriented methods. The space-time error analysis providing the optimal spatial mesh and the optimal temporal mesh is carried out in the case of a space-time adaptation at each time step in Section 4 and in the case of a space-time adaptation for subintervals with the GTFP in Section 5. The Global Space-Time Transient Fixed Point algorithm is given in Section 6 and is applied to several flows with vortex shedding past a cylinder in Section 7. The paper ends, Section 8, with some concluding remarks and perspectives.

## 2 Space-time continuous mesh

The error analysis is based on the continuous mesh framework [14, 15] where the spatial domain discretization, i.e., the mesh, is represented by continuous functions, i.e., metric fields, defined on the computational domain. In this framework, the discrete error model is recast into a continuous one. Then, minimizing the continuous error model gives the optimal continuous mesh, solution of a continuous optimal system. The optimal adapted mesh is obtained by discretizing the continuous metric field, that is to say by generating a unit mesh with respect to the prescribed metric field [11, 17]. This section introduces the notion of space-time continuous mesh, its complexity, and the mechanics to generate the discrete space-time meshes.

A metric field $(M(\mathbf{x}))_{\mathbf{x}\in\Omega}$ is a 3 × 3 matrix field defined on the 3D spatial computational domain $\Omega$, and such that for any $\mathbf{x}$ in $\Omega$, $M(\mathbf{x})$ is symmetric definite positive. A metric field $(M(\mathbf{x}))_{\mathbf{x}\in\Omega}$ is
called (spatial) continuous mesh [14]. A spatial mesh represented by the metric field $(\mathcal{M}(\mathbf{x}))_{\mathbf{x}\in\Omega}$ is any element of the class of meshes which are unit meshes for the metric field $(\mathcal{M}(\mathbf{x}))_{\mathbf{x}\in\Omega}$.<fnref n="1" /> The spatial complexity of a metric field:

$$C_{\text{space}}(\mathcal{M}) = \int_{\Omega} \sqrt{\det(\mathcal{M}(\mathbf{x}))} \mathrm{d}\mathbf{x}, \qquad (1)$$

represents the number of vertices of the spatial discretization modelled by the metric field.

**Definition 2.1 Space-time continuous mesh.** We call space-time continuous mesh $(\mathcal{M}, \tau)$ the knowledge of the following ingredients:

(i) a time step function $\tau: t \in [0, T] \mapsto \tau(t) \in [0, T[$

(ii) for every $t \in [0, T]$ a spatial metric field $(\mathcal{M}(\mathbf{x}, t))_{\mathbf{x}\in\Omega}$ of spatial complexity $\mathcal{N}(t) = C_{\text{space}}(\mathcal{M}(t))$.

**Remark 2.1** Not all space-time continuous mesh allows to derive a space-time mesh. A space-time continuous mesh $(\mathcal{M}, \tau)$ is a valid parametrization of a space-time discretization if, for any $t \in [0, T]$, a unit mesh can be built from $\mathcal{M}(t)$, and if it exits an integer $n_{step} \ge 1$ such that the time step $\tau$ satisfies

$$\int_{0}^{T} (\tau(t))^{-1} \mathrm{d}t \approx n_{step}.$$

**Definition 2.2 Complexity.** The space-time complexity $C(\mathcal{M}, \tau)$ of a space-time continuous mesh $(\mathcal{M}, \tau)$ is:

$$C_{\text{st}}(\mathcal{M}, \tau) = \int_{0}^{T} C_{\text{space}}(\mathcal{M}(t)) (\tau(t))^{-1} \mathrm{d}t. \qquad (2)$$

From a given space-time continuous mesh, we can recover a fully-discrete time-advancing mesh by choosing a unit space-time mesh of it:

**Definition 2.3 Unit space-time mesh.** Given a (valid and sufficiently smooth) space-time continuous mesh $(\mathcal{M}, \tau)$, a discrete space-time mesh $(\mathcal{H}_k, t_k)_k$ is unit with respect to $(\mathcal{M}, \tau)$ if it verifies:

(i) since the time density $(\tau(t))^{-1}$ is positive and satisfies $\int_{0}^{T} (\tau(t))^{-1} \mathrm{d}t \approx n_{step}$ we can successively build time levels $t_k$ by putting:

$$t_k \text{ such that } \int_{t_{k-1}}^{t_k} (\tau(t))^{-1} \mathrm{d}t = 1,$$

stopping when it does not holds, for $n_{step} = \text{integer} \left( \int_{0}^{T} (\tau(t))^{-1} \mathrm{d}t \right)$.

(ii) At any time level $t_k$, $\mathcal{M}_k = \mathcal{M}(t_k)$ is used for generating a unit spatial mesh $\mathcal{H}_k$.

**Remark 2.2** Since each spatial mesh $\mathcal{H}_k$ generated has about $C_{\text{space}}(\mathcal{M}_k)$ vertices, the total degrees of freedom in the time-advancing mesh $((\mathcal{H}_k)_k, (\tau_k)_k)$ is

$$C_{\text{st}}(\mathcal{M}, \tau) = \sum_{k=1}^{n_{\text{step}}} C_{\text{space}}(\mathcal{M}_k).$$

<footnote n="1">¹See [9]. In short, edges of the unit mesh are of length between $\sqrt{1/2}$ and $\sqrt{2}$ for the length induced by the metric $\mathcal{M}$.</footnote>
# 3 Space-time error analysis

In the following, given a space-time continuous mesh, we provide a model of the approximation error estimating the error committed when using a unit discrete space-time mesh with respect to this space-time continuous mesh. First, such a model is proposed for the 1D scalar advection equation, then the case of the compressible Navier-Stokes equations is addressed.

## 3.1 Feature-based error model for the scalar advection equation

Let us consider the scalar advection model

$$u_t + cu_x = 0, \quad x \in \mathbb{R}.$$

We consider the usual continuous $P_1$ FEM approximation on a splitting of $\mathbb{R}$ in intervals $[x_i, x_{i+1}]$:

$$V_h = \{ \varphi_h \in C^0(\mathbb{R}), \text{supp}(\varphi_h) \text{ is compact, } \varphi_h|_{[x_i, x_{i+1}]} \text{ is affine} \},$$

$$u_h \in V_h, \quad \forall \varphi_h \in V_h, \quad (\varphi_h, u_{h,t} + cu_{h,x}) = 0.$$

We choose the simplifying standpoint of a truncation error analysis of each separate term of the equation. We are first interested by the local spatial error:

$$\varepsilon_{\text{space}} = (\varphi_h, u_{h,t} + cu_{h,x} - (u_t + cu_x)) = 0.$$

A rough truncation error estimate writes this error in terms of mesh size $\Delta x$ ($K_x \in \mathbb{R}$):

$$|(\varphi_h, u_{h,t} + cu_{h,x} - (u_t + cu_x))| \le K_x \Delta x (H_{u_t} + H_u) \Delta x,$$

where $H_v$ holds for the absolute value of the Hessian of $v$. In the continuous mesh framework, if the 1D mesh $([x_i, x_{i+1}])_i$ is unit with respect to $(M(x))_{x \in \Omega}$, then we have $\Delta x = M(x)$. In other words, the continuous local error model writes:

$$|(\varphi_h, u_{h,t} + cu_{h,x} - (u_t + cu_x))| \le K_x \text{trace}(M^{-\frac{1}{2}}(H_{u_t} + H_u)M^{-\frac{1}{2}}).$$

with $K_x \in \mathbb{R}$ a constant depending on the dimension $d$. We also have to define, for all $t \in [0, T]$, a time dependent time step function $\Delta t = \tau(t)$, and a discrete time-derivative:

$$u_{h,\tau,t} \approx u_{h,t},$$

and we estimate the local time error ($K_t \in \mathbb{R}$) by:

$$|u_{h,\tau,t} - u_{h,t}| \le K_t \tau^\alpha \left| \frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}} \right|.$$

If the second-order backward differencing formula is used, we have $\alpha = 2$ and $K_t = \frac{1}{3}$.

We can now define a strongly coupled global space-time error model which is based on the $L^p$-norm of the previous local time and space error models:

$$\bar{\mathcal{E}}_{\text{st}}^p(M, \tau) = \int_0^T \int_{\Omega} \left[ K_t \tau^\alpha \left| \frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}} \right| + K_x \Delta x (H_{u_t} + H_u) \Delta x \right]^p dx dt,$$
where $\Delta x$ is the mesh size prescribed by $\mathcal{M}$. The analysis of $\mathcal{E}_{\text{st}}^p$ is rather complex as the temporal and spatial errors are tightly coupled leading to a strongly non-linear model. We therefore propose a slight simplification where the temporal and spatial errors are loosely coupled in order to design a tractable error modeling. In the sequel, the following *loosely coupled global space-time error estimate* is analyzed:

$$ \mathcal{E}_{\text{st}}^{p}(\mathcal{M}, \tau) = \int_{0}^{T} \int_{\Omega} \left[ K_t \tau^{\alpha} \left| \frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}} \right| \right]^{p} + \left[ K_x \Delta x (H_{u_t} + H_u) \Delta x \right]^{p} dxdt, \quad (3) $$

where we have a sum of a temporal error:

$$ \mathcal{E}_{\text{time}}(\mathcal{M}, \tau) = \int_{0}^{T} K_{t}^{p} \tau^{\alpha p} \int_{\Omega} \left| \frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}} \right|^{p} dxdt, $$

and a spatial error:

$$ \mathcal{E}_{\text{space}}(\mathcal{M}, \tau) = \int_{0}^{T} \int_{\Omega} \left( K_x \Delta x (H_{u_t} + H_u) \Delta x \right)^p dxdt. $$

## 3.2 Feature-based error model for CFD

Let $W = (\rho, \rho\mathbf{u}, \rho E)$ be the conservative variables vector where $\rho$ denotes the density (kg/m$^3$), $\mathbf{u}$ the velocity (m/s), $E$ the total energy per mass unit (m$^2$s$^{-2}$). The compressible Navier-Stokes equations reads:

$$ \Psi(W) = 0 \iff \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \\[1.5ex] \displaystyle \frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p - \nabla \cdot \mathcal{T} = 0, \\[1.5ex] \displaystyle \frac{\partial (\rho E)}{\partial t} + \nabla \cdot ((\rho E + p)\mathbf{u}) - \nabla \cdot (\mathcal{T} \cdot \mathbf{u}) - \nabla \cdot (\lambda \nabla \theta) = 0, \\[1.5ex] \text{+Boundary conditions} \end{array} \right. \quad (4) $$

where $p$ is the pressure (N/m$^2$), given by $p = (\gamma - 1)(\rho E - \frac{1}{2}\rho\|\mathbf{u}\|^2)$, where $\gamma$ is constant ($\gamma = 1.4$ in the sequel), $\theta$ the temperature ($K$) such that $\rho C_v \theta = E - \frac{1}{2}\rho\|\mathbf{u}\|^2$ ($C_v$ being the specific heat at constant volume), $\mu$ the laminar dynamic viscosity (kg/(ms)) and $\lambda$ the laminar conductivity. $\mathcal{T}$ is the laminar stress tensor:

$$ \mathcal{T} = \mu \left[ (\nabla \mathbf{u} + \nabla \mathbf{u}^{\mathrm{T}}) - \frac{2}{3} \nabla \cdot \mathbf{u} \mathbf{I} \right]. $$

The variation of nondimensionalized laminar dynamic viscosity and conductivity coefficients $\mu$ and $\lambda$ as function of the dimensional temperature $T$ are defined by Sutherland's law:

$$ \mu = \mu_{\infty} \left( \frac{\theta}{\theta_{\infty}} \right)^{\frac{3}{2}} \left( \frac{\theta_{\infty} + \text{Su}}{\theta + \text{Su}} \right) \quad \text{and} \quad \lambda = \lambda_{\infty} \left( \frac{\theta}{\theta_{\infty}} \right)^{\frac{3}{2}} \left( \frac{\theta_{\infty} + \text{Su}}{\theta + \text{Su}} \right), $$

where $\text{Su} = 110^{\circ}\text{K}$ is the Sutherland temperature and the index $\infty$ denotes reference quantities. The relation linking $\mu$ and $\lambda$ is expressed from the Prandtl laminar number:

$$ \mathrm{Pr} = \frac{\mu C_p}{\lambda} \quad \text{with} \quad \mathrm{Pr} = 0.72 \quad \text{for (dry) air,} $$

where $C_p$ is the specific heat at constant pressure.
For the compressible Navier-Stokes equation, the feature-based (FB) error analysis leads to the following global space-time error model in $L^p$ norm

$$ (\mathcal{E}_{\text{st}}^{\text{FB}})^p(\mathcal{M}, \tau) = \mathcal{E}_{\text{time}}^{\text{FB}}(\mathcal{M}, \tau) + \mathcal{E}_{\text{space}}^{\text{FB}}(\mathcal{M}, \tau), \qquad (5) $$

with

$$ \mathcal{E}_{\text{time}}^{\text{FB}}(\mathcal{M}, \tau) = \int_{0}^{T} K_{t}^{p} \tau^{\alpha p} \int_{\Omega} \left| \frac{\partial^{\alpha+1} W}{\partial t^{\alpha+1}} \right|^{p} \mathrm{d}\mathbf{x}\mathrm{d}t, \qquad (6) $$

and

$$ \mathcal{E}_{\text{space}}^{\text{FB}}(\mathcal{M}, \tau) = \int_{0}^{T} \int_{\Omega} \left( \text{trace} \left( \mathcal{M}^{-\frac{1}{2}}(\mathbf{x}, t) \mathbf{H}(\mathbf{x}, t) \mathcal{M}^{-\frac{1}{2}}(\mathbf{x}, t) \right) \right)^{p} \mathrm{d}\mathbf{x}\mathrm{d}t, \qquad (7) $$

where $\mathbf{H} = |H_{u_t} + H_u|$ depends on sensor $u$ computed from $W$ the solution of the state Equation (4).

### 3.3 Goal-oriented error model for CFD

We consider the goal-oriented formulation as introduced in [1, 7] and keep the notations of these papers. The goal is to minimize the error $(g, W - W_h)$ committed in the approximation of the functional (or scalar output):

$$ j = (g, W), $$

where $W$ is the exact solution of the state Equation (4) and $W_h$ the approximate solution. The novelty with respect to [1, 7] is that the error on the time discretization is also taken into account. This leads to an extra term $\mathcal{E}_{\text{time}}^{\text{GO}}(\mathcal{M}, \tau)$ in the goal-oriented (GO) global space-time error estimate:

$$ |(g, W_h - W)| \approx \mathcal{E}_{\text{st}}^{\text{GO}}(\mathcal{M}, \tau) = \mathcal{E}_{\text{space}}^{\text{GO}}(\mathcal{M}, \tau) + \mathcal{E}_{\text{time}}^{\text{GO}}(\mathcal{M}, \tau). \qquad (8) $$

We remind that the goal-oriented error estimate is in $L^1$ norm thus, here, we have $p = 1$. The global temporal error model is

$$ \mathcal{E}_{\text{time}}^{\text{GO}}(\mathcal{M}, \tau) = \int_{0}^{T} \int_{\Omega} K_t \tau^{\alpha} \left| W^{*} \frac{\partial^{\alpha+1} W}{\partial t^{\alpha+1}} \right| \mathrm{d}\mathbf{x}\mathrm{d}t, \qquad (9) $$

where $W^*$ is the adjoint state, solution of the adjoint system:

$$ -W_t^* + \left(\frac{\partial \Psi}{\partial W}\right)^* W^* = g. $$

We recall the spatial goal-oriented error estimate developed and progressively enriched in [1, 6, 7, 16]. The goal-oriented error model $\mathcal{E}_{\text{space}}^{\text{GO}}(\mathcal{M}, \tau)$ is expressed in terms of Euler fluxes $\mathcal{F}^E$, viscous fluxes $\mathcal{F}^V$, and boundary Euler fluxes $\bar{\mathcal{F}}$:

$$ \begin{aligned} \mathcal{E}_{\text{space}}^{\text{GO}}(\mathcal{M}, \tau) \approx & \int_{0}^{T} \int_{\Omega} |W_t^*| |W - \pi_{\mathcal{M}} W| \mathrm{d}\mathbf{x}\mathrm{d}t \\ & + \int_{0}^{T} \int_{\Omega} \left| \frac{\partial \mathcal{F}^E}{\partial W} \nabla W^* \right| |W - \pi_{\mathcal{M}} W| \mathrm{d}\mathbf{x}\mathrm{d}t \\ & + \int_{0}^{T} \int_{\Omega} \left| \frac{\partial \mathcal{F}^V}{\partial \nabla W} H(W^*) \right| |W - \pi_{\mathcal{M}} W| \mathrm{d}\mathbf{x}\mathrm{d}t \\ & + \int_{0}^{T} \int_{\Gamma} |W^*| |(\bar{\mathcal{F}}(W) - \pi_{\mathcal{M}} \bar{\mathcal{F}}(W)) \cdot \mathbf{n}| \mathrm{d}\Gamma \mathrm{d}t. \end{aligned} $$
Neglecting the boundary term, we note that all the terms have the form of a weighted interpolation error in $L^1$ norm on the conservative variables. We therefore deduce the following global spatial error model:

$$ \mathcal{E}_{\text{space}}^{\text{GO}}(M, \tau) \approx \int_{0}^{T} \int_{\Omega} \text{trace} \left( M^{-\frac{1}{2}}(\mathbf{x}, t) \mathbf{H}(\mathbf{x}, t) M^{-\frac{1}{2}}(\mathbf{x}, t) \right) \mathrm{d}\mathbf{x} \mathrm{d}t, \quad (10) $$

$$ \text{with } \mathbf{H}(\mathbf{x}, t) = \left| W_t^* + \frac{\partial \mathcal{F}^E}{\partial W} \nabla W^* + \frac{\partial \mathcal{F}^V}{\partial \nabla W} H(W^*) \right| |H(W)|, $$

where $H(W)$ (resp. $H(W^*)$) is the Hessian of $W$ (resp. $W^*$).

## 3.4 Unified error model for CFD

By analyzing the feature-based error model (5,6,7) and the goal-oriented error model (8,9,10), we note that their formulations are similar. Consequently, both models can be unified as follows:

$$ \mathcal{E}_{\text{st}}^{p}(M, \tau) = \mathcal{E}_{\text{time}}(M, \tau) + \mathcal{E}_{\text{space}}(M, \tau), \quad (11) $$

with

$$ \mathcal{E}_{\text{time}}(M, \tau) = \int_{0}^{T} K_{t}^{p} \tau^{\alpha p} \int_{\Omega} \left| W^{*} \frac{\partial^{\alpha+1} W}{\partial t^{\alpha+1}} \right|^{p} \mathrm{d}\mathbf{x} \mathrm{d}t, \quad (12) $$

and

$$ \mathcal{E}_{\text{space}}(M, \tau) = \int_{0}^{T} \int_{\Omega} \left( \text{trace} \left( M^{-\frac{1}{2}}(\mathbf{x}, t) \mathbf{H}(\mathbf{x}, t) M^{-\frac{1}{2}}(\mathbf{x}, t) \right) \right)^{p} \mathrm{d}\mathbf{x} \mathrm{d}t, \quad (13) $$

where:

- in the feature-based case: a sensor $u$ is computed from $W$, allowing to replace $W^* \frac{\partial^{\alpha+1} W}{\partial t^{\alpha+1}}$ by a time derivative $\frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}}$ of the sensor, and setting $\mathbf{H} = |H_{u_t} + H_u|$.

- in the goal-oriented case: $p = 1$, $W^*$ is the adjoint state, and $\mathbf{H}$ is defined as in Equation (10).

# 4 Analysis for all-time adaptation

This section defines an optimal space-time adaptive strategy when the mesh is adapted at each time step. Let $\mathcal{N}_{\text{st}}$ be an integer prescribed by the user which represents the space-time mesh complexity, we call mesh adaptation problem the following problem:

Find $(M, \tau)$ which minimizes $\mathcal{E}_{\text{st}}(M, \tau)$ under the constraint $\mathcal{C}_{\text{st}}(M, \tau) = \mathcal{N}_{\text{st}}$.

## 4.1 Analysis for all-time mesh adaptation at a given time

At a fixed time $t$, given a spatial mesh complexity $\mathcal{N}(t)$, we know the optimal metric $\mathcal{M}_{\text{opt}}(\mathbf{x}, t)$ under the constraint $\mathcal{C}_{\text{space}}(\mathcal{M}(t)) = \mathcal{N}(t)$ by minimizing the spatial part of the error model in $L^p$ norm:

$$ \mathcal{M}_{\text{opt}}(t) = \underset{\mathcal{M}}{\operatorname{Arg min}} \int_{\Omega} \left( \text{trace} \left( (\mathcal{M}(\mathbf{x}, t))^{-\frac{1}{2}} \mathbf{H}(\mathbf{x}, t) (\mathcal{M}(\mathbf{x}, t))^{-\frac{1}{2}} \right) \right)^p \mathrm{d}\mathbf{x}. $$
The pointwise optimal metric is given by [14]:

$$ \mathcal{M}_{\text{opt}}(\mathbf{x}, t) = \mathcal{N}(t)^{\frac{2}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} \mathrm{d}\mathbf{x} \right)^{-\frac{2}{d}} (\det \mathbf{H}(\mathbf{x}, t))^{-\frac{1}{2p+d}} \mathbf{H}(\mathbf{x}, t), \quad (14) $$

and the related optimal spatial error is:

$$ \mathcal{E}_{\text{space}}(\mathcal{M}_{\text{opt}}, t) = d^p K_x^p \mathcal{N}(t)^{-\frac{2p}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} \mathrm{d}\mathbf{x} \right)^{\frac{2p+d}{d}}, \quad (15) $$

where $d$ is the spatial domain dimension, and $K_x = \frac{1}{20}$ if $d=3$ or $K_x = \frac{1}{8}$ if $d=2$ (see [9], Corollary 1 of Theorem 4.2.2.).

## 4.2 Analysis for all-time adaptation over time interval

We are now interested in minimizing the space-time error over the time interval $[0, T]$. The unified space-time error model (11,12,13) becomes by using the optimal spatial error term (14):

$$ \mathcal{E}_{\text{st}}^{p}(\mathcal{M}, \tau) = \int_{0}^{T} \left[ K_{t}^{p} \tau(t)^{\alpha p} \int_{\Omega} |W^{*}(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}}|^{p} \mathrm{d}\mathbf{x} + d^{p} K_{x}^{p} \mathcal{N}(t)^{-\frac{2p}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} \mathrm{d}\mathbf{x} \right)^{\frac{2p+d}{d}} \right] dt, $$

under the space-time complexity constraint

$$ C_{\text{st}}(\mathcal{M}, \tau) = \int_{0}^{T} \mathcal{N}(t) (\tau(t))^{-1} dt = \mathcal{N}_{\text{st}}. $$

Analyzing the formulations of $\mathcal{E}_{\text{st}}^{p}(\mathcal{M}, \tau)$ and $C_{\text{st}}(\mathcal{M}, \tau)$, we note that the error model and the constraint are expressed in terms of $\mathcal{N}(t)$ and $\tau$ but no more in term of $\mathcal{M}$. Consequently, we reformulate the optimization problem in terms of $\mathcal{N}(t)$ and $\tau$ by introducing two functions:

$$ \mathfrak{F} : \mathcal{N} \in \mathcal{C}^0[0, T; \mathbb{R}] \mapsto \mathfrak{F}(\mathcal{N}) \in \mathcal{C}^0[0, T; \mathbb{R}], $$

$$ \bar{\mathcal{G}} : (\mathcal{N}, \tau) \in (\mathcal{C}^0[0, T; \mathbb{R}])^2 \mapsto \bar{\mathcal{G}}(\mathcal{N}, \tau) \in \mathcal{C}^0[0, T; \mathbb{R}], $$

such that

$$ \bar{\mathcal{G}}(\mathcal{N}, \tau) = K_t^p \tau(t)^{\alpha p} \int_{\Omega} |W^*(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}}|^p \mathrm{d}\mathbf{x}, $$

$$ \mathfrak{F}(\mathcal{N}) = d^p K_x^p \mathcal{N}(t)^{-\frac{2p}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} \mathrm{d}\mathbf{x} \right)^{\frac{2p+d}{d}}. $$

Then, the space-time mesh adaptation problem becomes :

$$ \text{Find } (\mathcal{N}_{\text{opt}}, \tau_{\text{opt}}) = \underset{\mathcal{N}, \tau}{\operatorname{Arg min}} \int_{0}^{T} (\bar{\mathcal{G}}(\mathcal{N}, \tau) + \mathfrak{F}(\mathcal{N})) \mathrm{d}t, $$

$$ \text{such that } C_{\text{st}}(\mathcal{N}, \tau) = \int_{0}^{T} \mathcal{N}(t) (\tau(t))^{-1} \mathrm{d}t = \mathcal{N}_{\text{st}}. $$

Let us start with a change of variables to avoid nonlinear constraints, namely $(\mathcal{N}, \zeta) = (\mathcal{N}, \mathcal{N} \tau^{-1})$,

$$ \mathcal{G}(\mathcal{N}, \zeta) = K_t^p \zeta(t)^{-\alpha p} \mathcal{N}(t)^{\alpha p} \int_{\Omega} |W^*(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}}|^p \mathrm{d}\mathbf{x}, $$
so the constraint becomes $\int_{0}^{T} \zeta(t) \, dt = \mathcal{N}_{\text{st}}$. Now, by setting

$$ \mathcal{U}(t) = \alpha p K_t^p \int_{\Omega} |W^*(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}}|^p d\mathbf{x}, \\ \mathcal{K}(t) = 2p \, d^{p-1} K_x^p \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}}, $$

the derivatives of functions $\mathcal{F}$ and $\mathcal{G}$ are given by

$$ \begin{aligned} \frac{\partial \mathcal{F}}{\partial \mathcal{N}} \delta \mathcal{N} &= -\mathcal{N}^{-\frac{2p+d}{d}} \mathcal{K} \delta \mathcal{N}, \\ \frac{\partial \mathcal{G}}{\partial \mathcal{N}} \delta \mathcal{N} &= \zeta^{-\alpha p} \mathcal{N}^{\alpha p-1} \mathcal{U} \delta \mathcal{N}, \\ \frac{\partial \mathcal{G}}{\partial \zeta} \delta \zeta &= -\zeta^{-\alpha p-1} \mathcal{N}^{\alpha p} \mathcal{U} \delta \zeta. \end{aligned} \tag{16} $$

The optimal condition writes

$$ \left\{ \begin{array}{l} \displaystyle \int_0^T \left( \frac{\partial \mathcal{F}}{\partial \mathcal{N}} + \frac{\partial \mathcal{G}}{\partial \mathcal{N}} \right) \delta \mathcal{N} \, dt = 0, \quad \forall \delta \mathcal{N}, \\[1em] \displaystyle \int_0^T \frac{\partial \mathcal{G}}{\partial \zeta} \delta \zeta \, dt = 0, \quad \forall \delta \zeta \text{ such that } \int_0^T \delta \zeta \, dt = 0, \end{array} \right. $$

from which we deduce

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathcal{F}}{\partial \mathcal{N}} + \frac{\partial \mathcal{G}}{\partial \mathcal{N}} = 0, \\[1em] \displaystyle \frac{\partial \mathcal{G}}{\partial \zeta} = -C, \end{array} \right. $$

where, according to Relations (16), $C$ is a positive constant not depending in time. The second equation writes

$$ \zeta(t)^{-\alpha p-1} \mathcal{N}(t)^{\alpha p} \mathcal{U}(t) = C, $$

and gives

$$ \mathcal{N}(t) = \left( \frac{C}{\mathcal{U}(t)} \right)^{\frac{1}{\alpha p}} \zeta(t)^{\frac{\alpha p+1}{\alpha p}} \quad \text{or} \quad \zeta(t) = \left( \frac{\mathcal{U}(t)}{C} \right)^{\frac{1}{\alpha p+1}} \mathcal{N}(t)^{\frac{\alpha p}{\alpha p+1}}. \tag{17} $$

Thanks to Equation (17), the first equation gives

$$ \mathcal{N}(t) = C^{-\frac{\alpha d}{\lambda}} \mathcal{K}(t)^{\frac{d(\alpha p+1)}{p\lambda}} \mathcal{U}(t)^{-\frac{d}{p\lambda}}, \quad \text{with} \quad \lambda = 2(\alpha p + 1) + \alpha d. $$

Recalling that $\zeta = \mathcal{N} \tau^{-1}$ we have

$$ \tau(t) = \left( \frac{C}{\mathcal{U}(t)} \right)^{\frac{1}{\alpha p+1}} \mathcal{N}(t)^{\frac{1}{\alpha p+1}}, $$
and the constraint on the complexity gives

$$ \mathcal{N}_{\text{st}} = \int_{0}^{T} \mathcal{N}(t) \tau(t)^{-1} dt = C^{-\frac{\alpha d+2}{\lambda}} \int_{0}^{T} \mathcal{K}(t)^{\frac{\alpha d}{\lambda}} \mathcal{U}(t)^{\frac{2}{\lambda}} dt, $$

from which we deduce

$$ C = \mathcal{N}_{\text{st}}^{-\frac{\lambda}{\alpha d+2}} \left( \int_{0}^{T} \mathcal{K}(t)^{\frac{\alpha d}{\lambda}} \mathcal{U}(t)^{\frac{2}{\lambda}} dt \right)^{\frac{\lambda}{\alpha d+2}}. $$

For clarity, we introduce the global variable $S = \int_{0}^{T} \mathcal{K}(t)^{\frac{\alpha d}{\lambda}} \mathcal{U}(t)^{\frac{2}{\lambda}} dt$ such that $C = \mathcal{N}_{\text{st}}^{-\frac{\lambda}{\alpha d+2}} S^{\frac{\lambda}{\alpha d+2}}$. Finally, using the value of the constant $C$, the solution of the optimization problem is given by

$$ \mathcal{N}_{\text{opt}}(t) = \mathcal{N}_{\text{st}}^{\frac{\alpha d}{\alpha d+2}} S^{-\frac{\alpha d}{\alpha d+2}} \mathcal{K}(t)^{\frac{d(\alpha p+1)}{p\lambda}} \mathcal{U}(t)^{-\frac{d}{p\lambda}}, $$

$$ \tau_{\text{opt}}(t) = \mathcal{N}_{\text{st}}^{-\frac{2}{\alpha d+2}} S^{\frac{2}{\alpha d+2}} \mathcal{K}(t)^{\frac{d}{p\lambda}} \mathcal{U}(t)^{-\frac{2p+d}{p\lambda}}, $$

The optimal spatial continuous mesh at each time step that minimizes the space-time error is simply obtained by using the above solution in Formula (14):

$$ \mathcal{M}_{\text{opt}}(\mathbf{x}, t) = \mathcal{N}_{\text{opt}}(t)^{\frac{2}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} d\mathbf{x} \right)^{-\frac{2}{d}} (\det \mathbf{H}(\mathbf{x}, t))^{-\frac{1}{2p+d}} \mathbf{H}(\mathbf{x}, t). $$

This results shows that the space-time error model optimally distributes the number of vertices at each time step in order to minimize the space-time error.

**Proposition 4.1 (Optimal space-time continuous mesh at each time step)**

We consider the global space-time error model in $L^p$ norm given by Equation (11). We define the two time dependent functions $\mathcal{U}$ and $\mathcal{K}$:

$$ \mathcal{U}(t) = \alpha p K_t^p \int_{\Omega} \left| W^*(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right|^p d\mathbf{x}, $$

$$ \mathcal{K}(t) = 2p d^{p-1} K_x^p \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}}. $$

$K_t$ and $\alpha$ depends on the chosen implicit time integration scheme, for instance for the second-order backward differencing formula we have $\alpha = 2$ and $K_t = \frac{1}{3}$. $K_x$ depends on the spatial domain dimension $d$ with $K_x = \frac{1}{8}$ in 2D and $K_x = \frac{1}{20}$ in 3D. We also introduce the following global variable:

$$ S = \int_{0}^{T} \mathcal{K}(t)^{\frac{\alpha d}{\lambda}} \mathcal{U}(t)^{\frac{2}{\lambda}} dt. $$

Then, for any time $t$, the optimal spatial mesh complexity $\mathcal{N}_{\text{opt}}(t)$ and the optimal time step $\tau_{\text{opt}}(t)$ are given by:

$$ \mathcal{N}_{\text{opt}}(t) = \mathcal{N}_{\text{st}}^{\frac{\alpha d}{\alpha d+2}} S^{-\frac{\alpha d}{\alpha d+2}} \mathcal{K}(t)^{\frac{d(\alpha p+1)}{p\lambda}} \mathcal{U}(t)^{-\frac{d}{p\lambda}}, $$

$$ \tau_{\text{opt}}(t) = \mathcal{N}_{\text{st}}^{-\frac{2}{\alpha d+2}} S^{\frac{2}{\alpha d+2}} \mathcal{K}(t)^{\frac{d}{p\lambda}} \mathcal{U}(t)^{-\frac{2p+d}{p\lambda}}, $$

with $\lambda = 2(\alpha p + 1) + \alpha d$, and the optimal spatial continuous mesh is:

$$ \mathcal{M}_{\text{opt}}(\mathbf{x}, t) = \mathcal{N}_{\text{opt}}(t)^{\frac{2}{d}} \left( \int_{\Omega} (\det \mathbf{H}(\mathbf{x}, t))^{\frac{p}{2p+d}} d\mathbf{x} \right)^{-\frac{2}{d}} (\det \mathbf{H}(\mathbf{x}, t))^{-\frac{1}{2p+d}} \mathbf{H}(\mathbf{x}, t), $$

where for the feature-based case $\mathbf{H} = |H_{u_t} + H_u|$, and for the goal-oriented case $\mathbf{H}$ is given by Relation (10) and $p = 1$.
# 5 Analysis with time subintervals

The analysis of last section considers that the spatial mesh is updated at each time step. As remarked in previous works, this option has two important disadvantages, namely (i) the computational effort consumed in regenerating the mesh at each time step is generally prohibitive for practical applications, and (ii) the loss of accuracy in transferring the solution from one spatial mesh to the next one is also generally too large.

## 5.1 The Global Transient Fixed Point (GTFP) algorithm

To avoid these disadvantages, the Transient Fixed Point (TFP) was introduced in [2, 3]. The idea was to keep the same adapted mesh for a given time subinterval where the flow solver performs several time steps. In this first version, the mesh adaptation loop applies successively to each subinterval. As a result, this approach does not allow an analysis of the space-time error because the error estimation is done subinterval by subinterval. Therefore, we cannot optimize the complete space-time mesh. To solve this issue, the Global Transient Fixed Point (GTFP) was proposed in [5, 7] in which a global space-time error evaluation is done after the complete time resolution of the simulation. This section recalls its main features.

The simulation time frame $[0, T]$ is split into $n_{adap}$ subintervals of same length (Figure 1):

$$[0, T] = [0 = t_0, t_1[ \cup \dots \cup [t_{i-1}, t_i[ \cup \dots \cup [t_{n_{adap}-1}, t_{n_{adap}} = T[ = \bigcup_{i=1}^{n_{adap}} [t_{i-1}, t_i[$

The number $n_{adap}$ of adaptation time subintervals $[t_{i-1}, t_i[$ is a discretization parameter to be specified by the user. Each subinterval $[t_{i-1}, t_i[$ contains a large number of time steps of the flow solver. The GTFP algorithm is schematized in Algorithm 1 where $\mathcal{H}$, $\mathcal{S}$ and $\mathcal{M}$ denote respectively meshes, solutions and metrics. The external loop applies a fixed point on the mesh adaptation process to converge the non-linear space-time mesh adaptation problem. In the internal loop, knowing the spatial mesh, the time subinterval $[t_{i-1}, t_i[$ is divided into $m$ time-integration intervals $[t_i^k, t_i^{k+1}]$, with $k = 0, \dots, m$ and $t_i^0 = t_{i-1}, t_i^m = t_i$. For $t_i^0$ the solution is given by interpolation from the previous mesh, and for any $k = 0, \dots, m-1$, the flow variables are advanced from time level $t_i^k$ to time level $t_i^{k+1}$ by means of the numerical scheme. Before the end of the external loop, a global error analysis produces the complete series of metrics and meshes for the next computation on the $n_{adap}$ subintervals.

Several analyses of the GTFP are proposed in [7] and [9] (Lemma 2.9.3). In particular, for the 3D anisotropic case, a necessary condition for second order spatial convergence is that subintervals $[t_{i-1}, t_i[$

<figure><img src="image_1.png" /><figcaption>Figure 1: Time splitting of the GTFP mesh adaptation algorithm. Subintervals (in green) used for the transient process and time steps (in red).</figcaption></figure>
**Algorithm 1** GTFP: Global Transient Fixed Point for Unsteady Flows [5]

Initial mesh and solution ($\mathcal{H}_0, \mathcal{S}_0^0$) and set targeted space-time complexity $\mathcal{N}_{st}$

# Fixed-point loop to converge the global space-time mesh adaptation problem
For $j = 1, n_{ptfx}$

# Adaptive loop to advance the solution in time on time frame $[0, T]$

1. For $i = 1, n_{adap}$ \# Advance the solution in time in subinterval $[t_{i-1}, t_i[$

(a) $\mathcal{S}_{0,i}^j$ = Interpolate conservatively next subinterval initial solution from $(\mathcal{H}_{i-1}^j, \mathcal{S}_{i-1}^j, \mathcal{H}_i^j)$;

(b) $\mathcal{S}_i^j$ = Compute solution on subinterval from pair $(\mathcal{H}_i^j, \mathcal{S}_{0,i}^j)$;

(c) $|\mathbf{H}_i^j$ = Compute subinterval Hessian-metric from solution sample $(\mathcal{H}_i^j, \{\mathcal{S}_i^j(k)\}_{k=1,nk})$;

EndFor

2. $\mathcal{C}^j$ = Compute space-time complexity from all Hessian-metrics $(\{|\mathbf{H}_i^j\}_{i=1,n_{adap}})$;

3. $\{\mathbf{M}_i^j\}_{i=1,n_{adap}}$ = Compute all subinterval unsteady metrics $(\mathcal{C}^j, \{|\mathbf{H}_i^j\}_{i=1,n_{adap}})$;

4. $\{\mathbf{M}_i^j\}_{i=1,n_{adap}}$ = Metric gradation on all subinterval unsteady metrics $\{\mathbf{M}_i^j\}_{i=1,n_{adap}}$;

5. $\{\mathcal{H}_i^{j+1}\}_{i=1,n_{adap}}$ = Generate all subinterval adapted meshes $(\{\mathcal{H}_i^j, \mathbf{M}_i^j\}_{i=1,n_{adap}})$;

EndFor

are **two times** smaller for a four times spatial smaller error [7]. As concerns space-time convergence in $\mathbf{L}^1([0, T]; \mathbf{L}^1(\Omega))$, which is the convergence in terms of the space-time complexity $\mathcal{N}_{st}$, it can be obtained in many cases at order 8/5 by dividing the time step by a factor 4: $\tau \to \tau/4$ and passing from $(\mathcal{N}_{st}, n_{adap})$ to $(32\mathcal{N}_{st}, 2n_{adap})$ [7]. Further estimates are proposed in [5].

The present work proposes the extension of the GTFP to time-accurate implicit schemes. The central novelty is the simultaneous adaptation of the mesh and the time steps. It relies on the definition of a time-advancing mesh.

## 5.2 Notations for space-time Global Transient Fixed Point

It is useful to specify how the main continuous mesh variables are depending on the time variable $t$. In the context of the GTFP, we have:

1. the mesh complexity $\mathcal{N}(t)$ is a constant over each subinterval:

$$ \forall i = 1, \dots, n_{adap}, \quad \forall t \in [t_{i-1}, t_i[, \quad \mathcal{N}(t) = \mathcal{N}^i \in \mathbb{R}, \qquad (18) $$

2. the continuous mesh (metric field) $(\mathcal{M}(\mathbf{x}, t))_{\mathbf{x} \in \Omega}$ is constant over each subinterval:

$$ \forall i = 1, \dots, n_{adap}, \quad \forall t \in [t_{i-1}, t_i[, \quad (\mathcal{M}(\mathbf{x}, t))_{\mathbf{x} \in \Omega} = (\mathcal{M}^i(\mathbf{x}))_{\mathbf{x} \in \Omega}, \qquad (19) $$

3. the associated adapted mesh is also fixed over each subinterval:

$$ \forall i = 1, \dots, n_{adap}, \quad \forall t \in [t_{i-1}, t_i[, \quad \mathcal{H}(t) = \mathcal{H}^i. \qquad (20) $$
The complexities, continuous meshes and adapted meshes are solely changing when passing from subinterval $[t_{i-1}, t_i]$ to the next one $[t_i, t_{i+1}]$.

This section presents the optimal choice satisfying the above particular properties by solving the related optimization problem. We are again interested in minimizing the *loosely coupled space-time error model*:

$$ \mathcal{E}_{\text{st}}^{p}(\mathcal{M}, \tau) = \int_{0}^{T} \left[ \int_{\Omega} \left( K_t \tau(t)^{\alpha} \left| W^{*}(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right| \right)^{p} d\mathbf{x} + \int_{\Omega} \left( \text{trace} \left( \mathcal{M}^{-\frac{1}{2}}(\mathbf{x}, t) \mathbf{H}(\mathbf{x}, t) \mathcal{M}^{-\frac{1}{2}}(\mathbf{x}, t) \right) \right)^{p} d\mathbf{x} \right] dt, $$

under the space-time complexity constraint

$$ \mathcal{C}_{\text{st}}(\mathcal{M}, \tau) = \int_{0}^{T} \mathcal{N}(t)(\tau(t))^{-1} dt = \mathcal{N}_{\text{st}}. $$

As previously, we denote the first term of our error model $\mathcal{E}_{\text{time}}(\mathcal{M}, \tau)$ and the second term $\mathcal{E}_{\text{space}}(\mathcal{M}, \tau)$. We then rewrite each component of the error model in terms of subintervals using the above notations:

$$ \mathcal{E}_{\text{time}}(\mathcal{M}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \mathcal{E}_{\text{time}}^{i}(\mathcal{M}^{i}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \int_{t_{i-1}}^{t_{i}} K_{t}^{p} \tau(t)^{\alpha p} \left( \int_{\Omega} \left| W^{*}(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right|^{p} d\mathbf{x} \right) dt, $$

$$ \mathcal{E}_{\text{space}}(\mathcal{M}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \mathcal{E}_{\text{space}}^{i}(\mathcal{M}^{i}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \int_{t_{i-1}}^{t_{i}} \int_{\Omega} \left( \text{trace} \left( (\mathcal{M}^{i}(\mathbf{x}))^{-\frac{1}{2}} \mathbf{H}(\mathbf{x}, t) (\mathcal{M}^{i}(\mathbf{x}))^{-\frac{1}{2}} \right) \right)^{p} d\mathbf{x}dt. $$

This defines our space-time error model for a GTFP mesh adaptation with $n_{\text{adap}}$ subintervals.

## 5.3 Spatial mesh optimization on a subinterval

Let $[t_i, t_{i+1}]$ be the considered subinterval. In the spatial error term, only the term $\mathbf{H}(\mathbf{x}, t)$ has a time dependency as the metric field is fixed on the considered subinterval. Its integral over time $\int_{t_{i-1}}^{t_i} \mathbf{H}(\mathbf{x}, t)dt$ has thus to be estimated. Here, we have two choices. Either, we can overestimate this integral using the $L^\infty$ norm in time:

$$ \int_{t_{i-1}}^{t_i} \mathbf{H}(\mathbf{x}, t) dt \le (t_i - t_{i-1}) \max_{t \in [t_i, t_{i+1}]} \mathbf{H}(\mathbf{x}, t) = \mathbf{H}_{L^\infty}^{i}(\mathbf{x}, t). $$

This choice has a better detection of the high temporal variations of the unknowns. Or, the integral can be estimated using the Trapezoidal rule:

$$ \int_{t_{i-1}}^{t_i} \mathbf{H}(\mathbf{x}, t) dt \approx \sum_{k=1}^{n_{\text{step}}} (t_{k+1} - t_k) \frac{\mathbf{H}(\mathbf{x}, t_{k+1}) + \mathbf{H}(\mathbf{x}, t_k)}{2} = \mathbf{H}_{L^1}^{i}(\mathbf{x}, t). $$

This choice is mandatory when considering the goal-oriented approach. In the following, we denote by $\mathbf{H}^i(\mathbf{x}, t)$ one of these two formulations.

The error in space for the considered subinterval becomes

$$ \mathcal{E}_{\text{space}}^{i}(\mathcal{M}^{i}) = \int_{\Omega} \left( \text{trace} \left( (\mathcal{M}^{i}(\mathbf{x}))^{-\frac{1}{2}} \mathbf{H}^{i}(\mathbf{x}) (\mathcal{M}^{i}(\mathbf{x}))^{-\frac{1}{2}} \right) \right)^{p} d\mathbf{x}, $$
and the spatial optimization problem on the subinterval reads

$$ \mathcal{M}_{opt}^{i} = \underset{\mathcal{M}^{i}}{\operatorname{Arg min}} \mathcal{E}_{\text{space}}^{i}(\mathcal{M}^{i}) \quad \text{such that} \quad \mathcal{C}(\mathcal{M}^{i}) = \int_{\Omega} \sqrt{\det \mathcal{M}^{i}(\mathbf{x})} \, d\mathbf{x} = \mathcal{N}^{i}. $$

As previously, the optimal continuous mesh of this problem is

$$ \mathcal{M}_{\text{opt}}^{i}(\mathbf{x}) = (\mathcal{N}^{i})^{\frac{2}{d}} \left( \int_{\Omega} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{-\frac{2}{d}} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{-\frac{1}{2p+d}} \mathbf{H}^{i}(\mathbf{x}), \quad (21) $$

and the corresponding optimal error writes

$$ \mathcal{E}_{\text{space}}^{i}(\mathcal{M}_{\text{opt}}^{i}) = d^{p} K_{x}^{p} (\mathcal{N}^{i})^{-\frac{2p}{d}} \left( \int_{\Omega} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}}. $$

## 5.4 Temporal optimization over the time subintervals

Now, the goal is to performed a global temporal minimization which will provide the optimal continuous mesh for each subinterval (in fact, the optimal spatial complexity for each subinterval as it is the only term which differs) and the optimal temporal mesh, i.e., the optimal time steps. We focus on minimizing the space-time error over the simulation time interval $[0, T]$:

$$ \mathcal{E}_{\text{st}}^{p}((\mathcal{M}^{i})_{i}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \left[ \int_{t_{i-1}}^{t_i} K_{t}^{p} \tau(t)^{\alpha p} \left( \int_{\Omega} \left| W^{*}(\mathbf{x}, t) \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right|^{p} d\mathbf{x} \right) + d^{p} K_{x}^{p} (\mathcal{N}^{i})^{-\frac{2p}{d}} \left( \int_{\Omega} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}} \right] dt, $$

under the space-time complexity constraint

$$ C_{\text{st}}((\mathcal{M}^{i})_{i}, \tau) = \sum_{i=1}^{n_{\text{adap}}} \mathcal{N}^{i} \int_{t_{i-1}}^{t_i} (\tau(t))^{-1} dt = \mathcal{N}_{\text{st}}. $$

As previously, the optimization problem is reformulated in terms of $(\mathcal{N}^i)_i$ and $\tau$ by introducing two functions:

$$ \bar{\mathcal{G}}(\mathcal{N}^i, \tau) = K_t^p \tau(t)^{\alpha p} \int_{\Omega} \left| W^*(\mathbf{x}, t) \cdot \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right|^p d\mathbf{x}, $$

$$ \mathfrak{F}(\mathcal{N}^i) = d^p K_x^p (\mathcal{N}^i)^{-\frac{2p}{d}} \left( \int_{\Omega} \left( \det \mathbf{H}^i(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}}, $$

leading to

$$ \left\{ \begin{array}{l} \text{Find } ((\mathcal{N}_{\text{opt}}^i)_i, \tau_{\text{opt}}) = \underset{(\mathcal{N})_{i, \tau}}{\operatorname{Arg min}} \sum_{i=1}^{n_{\text{adap}}} \left( \int_{t_{i-1}}^{t_i} \bar{\mathcal{G}}(\mathcal{N}^i, \tau) dt + \mathfrak{F}(\mathcal{N}^i) \right), \\ \text{such that } C_{\text{st}}((\mathcal{N}_{\text{opt}}^i)_i, \tau) = \sum_{i=1}^{n_{\text{adap}}} \mathcal{N}^i \int_{t_{i-1}}^{t_i} \tau(t)^{-1} dt = \mathcal{N}_{\text{st}}. \end{array} \right. $$

We consider the change of variables to avoid nonlinear constraints, namely $(\mathcal{N}^i, \zeta^i(t)) = (\mathcal{N}^i, \mathcal{N}^i \tau(t)^{-1})$, and we get

$$ \mathcal{G}(\mathcal{N}^i, \zeta^i) = K_t^p \zeta^i(t)^{-\alpha p} (\mathcal{N}^i)^{\alpha p} \int_{\Omega} \left| W^*(\mathbf{x}, t) \cdot \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}} \right|^p d\mathbf{x}, $$
and the constraint becomes $\sum_{i=0}^{n_{adap}} \int_{t_{i-1}}^{t_i} \zeta^i(t) dt = \mathcal{N}_{st}$. Now, by setting

$$ \mathcal{U}(t) = \alpha p K_t^p \int_{\Omega} |W^*(\mathbf{x}, t) \cdot \frac{\partial^{\alpha+1} W(\mathbf{x}, t)}{\partial t^{\alpha+1}}|^p d\mathbf{x}, \\ \mathcal{K}^i = 2p \, d^{p-1} K_x^p \left( \int_{\Omega} (\det \mathbf{H}^i(\mathbf{x}))^{\frac{p}{2p+d}} \, d\mathbf{x} \right)^{\frac{2p+d}{d}}, $$

the derivatives of functions $\mathfrak{F}$ and $\bar{\mathcal{G}}$ are given by

$$ \frac{\partial \mathfrak{F}}{\partial \mathcal{N}^i} \delta \mathcal{N}^i = -(\mathcal{N}^i)^{-\frac{2p+d}{d}} \mathcal{K}^i \delta \mathcal{N}^i, $$

$$ \frac{\partial \mathcal{G}}{\partial \mathcal{N}^i} \delta \mathcal{N}^i = \zeta^i(t)^{-\alpha p} (\mathcal{N}^i)^{\alpha p-1} \mathcal{U}(t) \delta \mathcal{N}^i, $$

$$ \frac{\partial \mathcal{G}}{\partial \zeta^i} \delta \zeta^i = -\zeta^i(t)^{-\alpha p-1} (\mathcal{N}^i)^{\alpha p} \mathcal{U}(t) \delta \zeta^i. $$

The optimality condition writes

$$ \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^{n_{adap}} \left( \int_{t_{i-1}}^{t_i} \frac{\partial \mathcal{G}}{\partial \mathcal{N}^i} dt + \frac{\partial \mathfrak{F}}{\partial \mathcal{N}^i} \right) \delta \mathcal{N}^i = 0, \quad \forall \delta \mathcal{N}^i, \\[1.5em] \displaystyle \sum_{i=1}^{n_{adap}} \int_{t_{i-1}}^{t_i} \frac{\partial \mathcal{G}}{\partial \zeta^i} \delta \zeta^i dt = 0, \quad \forall \delta \zeta^i \text{ such that } \sum_{i=1}^{n_{adap}} \int_{t_{i-1}}^{t_i} \delta \zeta^i dt = 0. \end{array} \right. $$

From the above relations we deduce

$$ \left\{ \begin{array}{l} \displaystyle \int_{t_{i-1}}^{t_i} \left( \frac{\partial \mathcal{G}}{\partial \mathcal{N}^i} dt + \frac{\partial \mathfrak{F}}{\partial \mathcal{N}^i} \right) = 0, \\[1.5em] \displaystyle \frac{\partial \mathcal{G}}{\partial \zeta^i} = -C. \end{array} \right. $$

where $C > 0$ does not depend on time. The second equation gives

$$ \mathcal{N}^i = \left( \frac{C}{\mathcal{U}(t)} \right)^{\frac{1}{\alpha p}} \zeta^i(t)^{\frac{\alpha p+1}{\alpha p}} \quad \text{or} \quad \zeta^i(t) = \left( \frac{\mathcal{U}(t)}{C} \right)^{\frac{1}{\alpha p+1}} (\mathcal{N}^i)^{\frac{\alpha p}{\alpha p+1}}. \qquad (22) $$

Thanks to Equation (22), the first equation gives

$$ \mathcal{N}^i = C^{-\frac{\alpha d}{\lambda}} (\mathcal{K}^i)^{\frac{d(\alpha p+1)}{\lambda p}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{-\frac{d(\alpha p+1)}{\lambda p}}, \quad \text{with } \lambda = 2(\alpha p + 1) + \alpha d. \qquad (23) $$

Recalling that $\zeta^i(t) = \mathcal{N}^i \tau^{-1}(t)$ we have immediately

$$ \tau^i(t) = \left( \frac{C}{\mathcal{U}(t)} \right)^{\frac{1}{\alpha p+1}} (\mathcal{N}^i)^{\frac{1}{\alpha p+1}}, \qquad (24) $$
then the constraint on the space-time complexity, using Equations (24) then (23) to replace $\tau(t)$ then $\mathcal{N}^i$, gives

$$ \mathcal{N}_{\text{st}} = \sum_{i=1}^{n_{\text{adap}}} \int_{t_{i-1}}^{t_i} \mathcal{N}^i(\tau(t))^{-1} dt = C^{-\frac{1}{\alpha p+1}} \sum_{i=1}^{n_{\text{adap}}} (\mathcal{N}^i)^{\frac{\alpha p}{\alpha p+1}} \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \\ C^{-\frac{\alpha d+2}{\lambda}} \sum_{i=1}^{n_{\text{adap}}} (\mathcal{K}^i)^{\frac{\alpha d}{\lambda}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{\frac{2(\alpha p+1)}{\lambda}}, $$

and we deduce

$$ C = \mathcal{N}_{\text{st}}^{-\frac{\lambda}{\alpha d+2}} \left( \sum_{i=1}^{n_{\text{adap}}} (\mathcal{K}^i)^{\frac{\alpha d}{\lambda}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{\frac{2(\alpha p+1)}{\lambda}} \right)^{\frac{\lambda}{\alpha d+2}}. $$

For readability, we introduce the global variable $S = \sum_{i=1}^{n_{\text{adap}}} (\mathcal{K}^i)^{\frac{\alpha d}{\lambda}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{\frac{2(\alpha p+1)}{\lambda}}$ such that $C = \mathcal{N}_{\text{st}}^{-\frac{\lambda}{\alpha d+2}} S^{\frac{\lambda}{\alpha d+2}}$. Knowing the constant $C$, we can now express the solution of the optimization problem

$$ \mathcal{N}_{\text{opt}}^i = \mathcal{N}_{\text{st}}^{\frac{\alpha d}{\alpha d+2}} S^{-\frac{\alpha d}{\alpha d+2}} (\mathcal{K}^i)^{\frac{d(\alpha p+1)}{\lambda p}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{-\frac{d(\alpha p+1)}{\lambda p}}, $$

$$ \tau_{\text{opt}}^{i}(t) = \mathcal{N}_{\text{st}}^{-\frac{2}{\alpha d+2}} S^{\frac{2}{\alpha d+2}} (\mathcal{K}^{i})^{\frac{d}{\lambda p}} (\mathcal{U}(t))^{-\frac{1}{\alpha p+1}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{-\frac{d}{\lambda p}}. $$

The optimal spatial continuous mesh for a given subinterval $[t_i, t_{i+1}]$ that minimizes the space-time error is simply obtained by using the above solution in Formula (21):

$$ \mathcal{M}_{\text{opt}}^{i}(\mathbf{x}) = (\mathcal{N}_{\text{opt}}^{i})^{\frac{2}{d}} \left( \int_{\Omega} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{-\frac{2}{d}} \left( \det \mathbf{H}^{i}(\mathbf{x}) \right)^{-\frac{1}{2p+d}} \mathbf{H}^{i}(\mathbf{x}). $$

This results shows that the space-time error model optimally distributes the number of vertices for each subinterval in order to minimize the space-time error.

### Proposition 5.1 (Optimal space-time continuous mesh for the GTFP)

We consider the global space-time error model in $L^p$ norm given by Equation (11). We define the time dependent function $\mathcal{U}$ and the variable $\mathcal{K}^i$ defined on each subinterval:

$$ \mathcal{U}(t) = \alpha p K_t^p \int_{\Omega} |W^*(\mathbf{x}, t)|^{\frac{\alpha p+1}{\alpha p+1}} W(\mathbf{x}, t)^{\frac{p}{2p+d}} d\mathbf{x}, \\ \mathcal{K}^i = 2p d^{p-1} K_x^p \left( \int_{\Omega} \left( \det \mathbf{H}^i(\mathbf{x}) \right)^{\frac{p}{2p+d}} d\mathbf{x} \right)^{\frac{2p+d}{d}}, \quad (25) $$

$K_t$ and $\alpha$ depends on the chosen implicit time integration scheme, for instance for the second-order backward differencing formula we have $\alpha = 2$ and $K_t = \frac{1}{3}$. $K_x$ depends on the spatial domain dimension $d$ with $K_x = \frac{1}{8}$ in 2D and $K_x = \frac{1}{20}$ in 3D. We also introduce the following global variable:

$$ S = \sum_{i=1}^{n_{\text{adap}}} (\mathcal{K}^i)^{\frac{\alpha d}{\lambda}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{\frac{2(\alpha p+1)}{\lambda}}, $$
Then, for any subinterval $[t_i, t_{i+1}]$, the optimal spatial mesh complexity $\mathcal{N}_{opt}^i$ and the optimal time step $\tau_{opt}^i(t)$ are given by:

$$ \mathcal{N}_{opt}^i = \mathcal{N}_{st}^{\frac{\alpha d}{\alpha d+2}} S^{-\frac{\alpha d}{\alpha d+2}} (\mathcal{K}^i)^{\frac{d(\alpha p+1)}{\lambda p}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{-\frac{d(\alpha p+1)}{\lambda p}}, \quad (26) $$

$$ \tau_{opt}^{i}(t) = \mathcal{N}_{st}^{-\frac{2}{\alpha d+2}} S^{\frac{2}{\alpha d+2}} (\mathcal{K}^{i})^{\frac{d}{\lambda p}} (\mathcal{U}(t))^{-\frac{1}{\alpha p+1}} \left( \int_{t_{i-1}}^{t_i} \mathcal{U}(t)^{\frac{1}{\alpha p+1}} dt \right)^{-\frac{d}{\lambda p}}. \quad (27) $$

with $\lambda = 2(\alpha p + 1) + \alpha d$, and the optimal spatial continuous mesh is:

$$ \mathcal{M}_{opt}^{i}(\mathbf{x}) = (\mathcal{N}_{opt}^{i})^{\frac{2}{d}} \left( \int_{\Omega} (\det \mathbf{H}^{i}(\mathbf{x}))^{\frac{p}{2p+d}} d\mathbf{x} \right)^{-\frac{2}{d}} (\det \mathbf{H}^{i}(\mathbf{x}))^{-\frac{1}{2p+d}} \mathbf{H}^{i}(\mathbf{x}), \quad (28) $$

where $\mathbf{H}^i$ is defined according to Section 5.3 with for the feature-based case $\mathbf{H} = |H_{u_t} + H_u|$, and for the goal-oriented case $\mathbf{H}$ is given by Relation (10) and $p = 1$.

# 6 Global Space-Time Transient Fixed Point algorithm

The Global Transient Fixed Point (GTFP) algorithm was proposed for specifying automatically a succession of $n_{adap}$ meshes over a decomposition in subintervals used for the transient process, see Figure 1. This algorithm needs to be extended to the space-time error analysis of the previous section where the temporal error and the adapted temporal mesh have to be also managed. This new algorithm, called Global Space-Time Transient Fixed Point (GSTTFP), is presented in Algorithm 2.

From a practical point of view, the flowchart presented in Figure 2 shows when the terms involved in Proposition 5.1 are computed. In a first phase (in purple), the successive time subintervals $[t_{i-1}, t_i[$ are visited for:

- (1a.) interpolating the final solution of the previous subinterval on the new subinterval mesh (if necessary),

- (1b.) computing the flow solution,

- (1c.) according to Equation (25), computing the spatial error $\mathcal{K}$

- (1d.) according to Equation (25), computing the temporal error $\mathcal{U}$.

In a second phase (in green), when all the time subintervals have been computed, the global variables are evaluated: (2.) the optimal spatial complexities $\mathcal{N}_{opt}$ (Eq. (26)), (3.) the optimal time steps $\tau_{opt}$ (Eq. (27)), and (4.) the optimal continuous meshes $\mathcal{M}_{opt}$ (Eq. (28)). Finally, (5.) a new discrete space-time mesh is generated.

# 7 Numerical experiments

This paper focuses on unsteady flows that can be considered as quasi-steady where the unsteady turbulence evolves in a fixed local region. For instance, this is the case of the turbulence behind a cylinder. For this type of flow, it is interesting to consider the specific case of a single spatial adapted mesh for
**Algorithm 2** GST: Global Space-Time Transient Fixed Point for Unsteady Flows

Initial mesh, time step and solution ($\mathcal{H}^0, \tau^0, S_0^0$) and set targeted space-time complexity $\mathcal{N}_{\text{st}}$

\# Fixed-point loop to converge the global space-time mesh adaptation problem
For $j = 1, n_{\text{ptfx}}$

\# Adaptive loop to advance the solution in time on time frame $[0, T]$

1. For $i = 1, n_{\text{adap}}$ \# Advance the solution in time in subinterval $[t_{i-1}, t_i[$

(a) $S_{0,i}^j$ = Interpolate conservatively next subinterval initial solution from ($\mathcal{H}_{i-1}^j, S_{i-1}^j, \mathcal{H}_i^j$);

(b) $S_i^j$ = Compute solution on subinterval from triple ($\mathcal{H}_i^j, \tau_i^j, S_{0,i}^j$);

(c) $|\mathbf{H}_i^j$ = Compute subinterval Hessian-metric from solution sample ($\mathcal{H}_i^j, \{S_i^j(k)\}_{k=1,nk}$);

(d) $(\mathcal{E}_{\text{time}})_i^j$ = Compute time error ($\mathcal{H}_i^j, \tau_i^j, \{S_i^j(k)\}_{k=1,nk}$);

EndFor

2. $\{\mathcal{N}_{\text{opt}}^{i,j+1}\}_{i=1,n_{\text{adap}}} = \text{Compute space complexity from } (S_i^j, \{|\mathbf{H}_i^j\}_{i=1,n_{\text{adap}}})$;

3. $\{\tau_{\text{opt}}^{i,j+1}\}_{i=1,n_{\text{adap}}} = \text{Compute time step from } (S_i^j, \{|\mathbf{H}_i^j\}_{i=1,n_{\text{adap}}})$;

4. $\{\mathbf{M}_i^{j+1}\}_{i=1,n_{\text{adap}}} = \text{Compute all metrics + gradation } (\{\mathcal{N}_{\text{opt}}^{i,j+1}\}_{i=1,n_{\text{adap}}}, \{|\mathbf{H}_i^j\}_{i=1,n_{\text{adap}}})$;

5. $\mathcal{H}_{\text{st}}^{j+1} = (\{\mathcal{H}_i^{j+1}, (t_k)_i^{j+1}\}_{i=1,n_{\text{adap}}}) = \text{Generate all subinterval adapted meshes and time steps } (\{\mathcal{H}_i^j, \mathbf{M}_i^j, \tau_{\text{opt}}^{i,j+1}\}_{i=1,n_{\text{adap}}})$;

EndFor

the entire simulation time frame. This type of approach is also very suitable for Large Eddy Simulation (LES) applications. This must be contrasted with true unsteady flows involving fast dynamics, such as blast wave applications, where considering many sub-intervals is mandatory to optimize the space-time mesh [5].

As regards the error estimate, we only consider the feature-based method with the local Mach number as sensor.

The Navier-Stokes system is discretized in space using a vertex-centered mixed finite volume/finite element numerical scheme applied to unstructured meshes composed of triangles/tetrahedra. The diffusive terms are discretized using P1 Galerkin finite-elements on the triangle/tetrahedra, whereas finite-volumes are used for the convective terms. The numerical approximation of the convective fluxes at the interface of neighboring finite volume cells is based on the HLLC approximate Riemann solver. To obtain second-order accuracy in space, the Monotone Upwind Scheme for Conservation Laws reconstruction method (MUSCL) is used. Time advancing is carried out through an implicit linearized method, based on a second-order accurate backward difference scheme, which means that in our theory $\alpha = 2$ and $K_t = \frac{1}{3}$ in Proposition 5.1.

As concerns mesh adaptation, we keep the ingredients of [2], except that the temporal error needs to be evaluated. It is a third derivative and, as any truncation error, most approximations are highly oscillating. We apply a strong moving average filtering. For all the simulations, the Spalart-Allmaras turbulence model is used [22]. More details can be found in [1, 9].

For all the presented simulations, the physical conditions are: $|\mathbf{u}| = 1 \text{ m/s}$, $\rho = 1 \text{ kg/m}^3$, $\mu = \frac{1}{Re}$
<figure><img src="image_2.png" /><figcaption>Figure 2: Global Space-Time Transient Fixed Point (GSTTFP) flowchart. The “?” holds for testing whether the GSTTFP fixed point is converged or not.</figcaption></figure>

where $Re$ is the Reynolds number. The Mach number is set to 0.3, so the pressure is $p \approx 7.9365$ Pa. Our experience with vortex shedding flows is that using a unique adapted mesh for the several vortex shedding cycles lands to easier flow and error statistics.

## 7.1 2D flow past a cylinder at Reynolds 3900

The first test case is the 2D computation of a flow around a cylinder at Reynolds number 3900. A circular computational domain of radius 20 diameter of the cylinder is considered. The solution is initialized by running the solution for 200 seconds in physical time on an initial radial mesh composed of 12K vertices. Then, the simulation consists in running a period of 10 seconds which corresponds to two vortex shedding cycles. The space-time mesh adaptation is carried out for this 10 seconds time frame by using a single adapted spatial mesh, that is $n_{adap} = 1$. Four space-time complexity values were considered, namely $N_{st}$ equal to 5M, 10M, 20M and 40M.

Figure 3 shows the obtained adapted spatial mesh and the associated final solution for a space-time complexity of 5M and 40M. We note that the mesh is highly refined in the boundary layer region and in the turbulent wake of the cylinder. The adapted mesh is anisotropic in the boundary layer while it is almost isotropic in the wake.
Figure 4 shows the computed time steps as functions of the physical time, i.e., the computed adapted temporal mesh, for the successive iterations of the fixed point, and this, for the four space-time complexities. We note that, after the second fixed point iteration, the curves are more or less identical. This points out the convergence of the GSTTFP algorithm. When we compare the computed time steps for the four complexities, we notice that the curves are very similar with a smaller amplitude for the largest space-time complexities corresponding to a higher time accuracy. It is interesting to note that we observe the quasi-periodicity of the flow in the temporal mesh - four periods are clearly visible which correspond to two vortex shedding cycles - despite the fact that no specific physical criteria have been introduced into the error estimation. This demonstrates that the error estimate is capture the physics of the flow.

By comparison with the previous version of the GTFP algorithm which does not adapt the time step, the extra computations are negligible. Conversely, the optimal choice of time step may induce an important saving in CPU time. As concerns the CPU time, each of the 15 iterations of the fixed point adaptive algorithm for the 40M space-time complexity takes six minutes in serial on a laptop Dell Precision of 2.3 GHz.

Figure 5 presents the evolution of the spatial $\mathcal{E}_{\text{space}}$ and the time $\mathcal{E}_{\text{time}}$ error functionals for the successive adaptation fixed point iterations. When the algorithm starts, we note a quick convergence of the process as only three to four fixed point iterations are sufficient to reach the optimal errors values. This result leads to two remarks.

First, whatever the space-time complexity, both errors converge toward similar values which is in agreement with our theory where the aim is to balance both errors. Indeed, if it is not the case, for instance if the temporal error is larger than the spatial one, this means that we could reduce slightly the number of vertices, increasing slightly the spatial error, and increase the number of time steps (decreasing the temporal error) with the result of making these two errors closer. However, in the proposed method, the optimization is performed in a continuous context, producing an optimality condition,

<figure><img src="image_3.png" /><figcaption>Figure 3: 2D cylinder flow at Reynolds 3900. Adapted mesh (left) and the associated Mach number (right) for a space-time complexity of 5M (top) and 40M (bottom).</figcaption></figure>
which we in turn discretize. The final difference between spatial and temporal errors is therefore due to
the discretization (in space and time). The result obtained confirms this argument because we observe
that when the space-time complexity increases from 5M to 40M, the final gap between the two errors
decreases quickly.

Second, we can analyze the convergence order of the spatial and temporal errors between the space-
time complexity of 5M and 40M. Between both complexities, the size of the space-time mesh, which is of
dimension 3 (2D+t) is increased by a factor 8. For the errors, the reduction factor is 0.017/0.005 = 3.4
from which we deduce a convergence order<fnref n="2" /> of 1.76, close to second order.

<footnote n="2">²We have $\mathcal{E} \sim C N^{-\frac{\alpha}{d}}$ where $\alpha$ is the order of convergence, thus $\alpha = -d \frac{\ln(\varepsilon_2/\varepsilon_1)}{\ln(N_2/N_1)}$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 4: 2D cylinder flow at Reynolds 3900. Computed time step as a function of the physical time for space-time complexities 5M, 10M, 20M, 40M (from top to bottom). The convergence of the computed temporal discretization for the successive iterations of the fixed point is depicted.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: 2D cylinder flow at Reynolds 3900. Continuous spatial and temporal errors for complexities 5M, 10M, 20M, 40M (from top to bottom) at each fixed point iteration. The corresponding total error levels at convergence are respectively 0.034, 0.022, 0.014, 0.01.</figcaption></figure>
## 7.2 2D flow past a cylinder at Reynolds 3900 with multi-mesh adaptation

This section presents the computation of the previous test case when several different meshes on several time sub-intervals are used to adapt the space-time mesh, that is:

$$n_{adap} > 1.$$

The simulation parameters are then exactly the same as Section 7.1 (2D flow around a cylinder, Reynolds 3900, Spalart-Allmaras model, same initialisation, time interval is [0, 10](seconds), 15 adaptation cycles) except $n_{adap}$, which is set to $n_{adap} = 2$ when the 5 M complexity is run, $n_{adap} = 4$ when the 10 M complexity is run, $n_{adap} = 8$ when the 20 M complexity is run, and $n_{adap} = 16$ when the 40 M complexity is run.

<figure><img src="image_6.png" /><figcaption>Figure 6: 2D cylinder flow at Reynolds 3900 with multi-mesh adaptation. Final adapted meshes for a space-time complexity of 40M with 16 subintervals. View of the subinterval meshes 1, 4, 7, 10, 13, 16 (from left to right and from top to bottom) after the 15 fixed point iterations.</figcaption></figure>

A subset of the 16 final adapted meshes is shown in Figure 6. As $n_{adap} = 16$ and the 10 seconds simulation corresponds to two vortex shedding cycles, 8 different adapted meshes are used to simulate one vortex shedding cycle. As a result, we clearly see the vortex shedding phenomena inside the adapted meshes even if, in this case, the mesh size is on average 30K vertices while it was 60K vertices with $n_{adap} = 1$. Moreover, we observe also in Figure 6 the progressive refinement of the wake through the simulation of the last fixed point iteration.

In contrast to the $n_{adap} = 1$ calculation, a temporal discretization of 1 626 time steps has computed which is more or less twice the number (825) of time steps of the $n_{adap} = 1$ case. If we analyze the computed temporal discretization in Figure 7, the four periods are also visible but the curve is less smooth showing that some specific details in the solution are more accurately captured. When we examine the time dependence of the optimum time step length and compare to the one-mesh calculation, the multi-mesh optimum time step length seems importantly perturbed by the 15 interfaces between subintervals.
<figure><img src="image_7.png" /><figcaption>Figure 7: 2D cylinder flow at Reynolds 3900 with multi-mesh adaptation. Computed time step as a function of the physical time for space-time complexity 40M. The convergence of the computed temporal discretization for the successive iterations of the fixed point is depicted.</figcaption></figure>
<table><caption>Table 1 compares the numbers of both simulations.</caption><thead><tr><th>Method</th><th>Specified<br>space-time<br>complexity</th><th>Optimal<br>space<br>complexity</th><th>Optimal<br>number of<br>vertices</th><th>Optimal<br>number of<br>time steps</th></tr></thead><tbody><tr><td>One mesh</td><td>40M</td><td>48K</td><td>60K</td><td>825</td></tr><tr><td>Multi-mesh</td><td>40M</td><td>25K</td><td>30K</td><td>1626</td></tr></tbody></table>

<caption>Table 1: 2D cylinder flow at Reynolds 3900. Statistics of the one-mesh ($n_{adap} = 1$) and multi-mesh ($n_{adap} = 16$) simulations.</caption>

Figure 8 shows a plot of the total error $\mathcal{E}_{st}$ for our computations described above, in respect of the number of sub-intervals 2, 4, 8 and 16 for the respective complexity 5M, 10M, 20M and 40M. Table 2 compares the total error obtained when choosing a single mesh or several meshes for the experiment we have described. The total error is improved by the multimesh option but the gain is small due to difficulties with the present method to master the time error.

<table><caption>Table 2: 2D cylinder flow at Reynolds 3900. Total error for the one-mesh and multi-mesh simulations.</caption><thead><tr><th>Method</th><th>5M</th><th>10M</th><th>20M</th><th>40M</th></tr></thead><tbody><tr><td>One mesh</td><td>0.0344</td><td>0.0225</td><td>0.0149</td><td>0.00982</td></tr><tr><td>Multi-mesh</td><td>0.0294</td><td>0.0192</td><td>0.0131</td><td>0.00938</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 8: 2D cylinder flow at Reynolds 3900 with multi-mesh. Continuous total errors for complexities 5M and 2 sub-intervals, 10M and 4 sub-intervals, 20M and 8 sub-intervals, 40M and 16 sub-intervals for the last fixed point iteration (of each computation).</figcaption></figure>
## 7.3 2D flow past a cylinder at Reynolds 1M

The second test case is the 2D computation of a flow around a cylinder at Reynolds number 1M. The same circular computational domain as above is used. The solution is initialized by running the solution for 145 seconds in physical time on an initial radial mesh composed of 40K vertices. Then, the simulation consists in running a period of 9.3 seconds which corresponds to two vortex shedding cycles. The space-time mesh adaptation is carried out for this 9.3 seconds time frame by using a single adapted spatial mesh, that is $n_{adap} = 1$. Again, four space-time complexity values were considered, namely $N_{st}$ equal to 12.5M, 25M, 50M and 100M. Similar outputs are presented for the analysis.

Figure 9 shows the obtained adapted spatial mesh and the associated final solution for a space-time complexity of 12.5M and 100M. We note that the mesh is highly refined in the very thin boundary layer region and in the turbulent wake of the cylinder. The adapted mesh is highly anisotropic in the boundary layer (with aspect ratio $\simeq 10^2$) while it is almost isotropic in the wake. We note that the wake region is narrower near the cylinder compared to the previous case.

The computed adapted temporal discretization for the successive iterations of the fixed point for each complexity is depicted in Figure 11. As previously, two fixed point iterations are sufficient to converge the temporal adapted mesh. Whatever the complexity, the four temporal meshes are similar thus showing the same physics with five periods. The larger the complexity, the better the accuracy with smaller time steps. Note that for the finer mesh, a new detail in the solution is captured and it appears in the temporal mesh.

Figure 12 shows the evolution of the spatial $\mathcal{E}_{\text{space}}$ and the time $\mathcal{E}_{\text{time}}$ error functionals for the successive adaptation fixed point iterations. Again, the quick convergence of the process is clear as only two fixed point iterations are sufficient to reach the optimal errors values. As stated previously, in the

<figure><img src="image_9.png" /><figcaption>Figure 9: 2D cylinder flow at Reynolds 1M. Adapted mesh (left) and the associated Mach number (right) for a space-time complexity of 12.5M (top) and 100M (bottom).</figcaption></figure>
proposed method, the optimization is performed in a continuous context, producing an optimality condition, which we in turn discretize. We may have a difference between spatial and temporal errors due to the discretization (in space and time). The difference between the temporal and spatial errors is of the order of 20%, which is acceptable. Note that the difference between the two errors in absolute value decreases when the space-time complexity increases: $\delta\mathcal{E} \sim 0.0012$ for $N_{\text{st}} = 5M$ and $\delta\mathcal{E} \sim 0.0004$ for $N_{\text{st}} = 40M$. This points out the mesh convergence of the GSTTFP process.

For this case, we obtain a convergence order of 1.81 of the space-time error between the space-time complexity of 12.5M and 100M.

Finally, Figure 10 plots the simulation CFL (min, max and average) and the minimal triangle's height of the mesh for the four complexities. We note that the process automatically sets an average CFL close 6500 while preserving the solution accuracy. At these high CFL values, implicit time integration schemes are a lot more efficient in CPU time than explicit schemes. Note that, when the mesh complexity increases, the smallest height of the mesh decreases (as we increase the mesh accuracy) and the computed time step decreases. But, the CFL stays almost the same whatever the complexity meaning that, depending on the physics, the GSTTFP algorithm automatically finds the optimal CFL to run the simulation. This demonstrates the powerfulness of the proposed method.

<figure><img src="image_10.png" /><figcaption>Figure 10: 2D cylinder flow at Reynolds 1M. Evolution of CFL and the minima mesh height for the all the space-time complexities.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: 2D cylinder flow at Reynolds 1M. Computed time steps as a function of the physical time for space-time complexities 12.5M, 25M, 50M, 100M (from top to bottom). The convergence of the computed temporal discretization for the successive iterations of the fixed point is depicted.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 12: 2D cylinder flow at Reynolds 1M. Continuous spatial and temporal errors for complexities 12.5M, 25M, 50M, 100M (from top to bottom) at each fixed point iteration. The corresponding total error levels at convergence are respectively 0.014, 0.009, 0.006, 0.004.</figcaption></figure>
7.4 3D flow past a cylinder at Reynolds 3900

We consider again the first test case, but this time in 3D. The computational domain is cylindrical of
radius 40 diameter of cylinder and a span of 3.14 diameters. The solution is initialized by running the
solution for 140 seconds in physical time on an initial radial mesh composed of 840K vertices. The
simulation consists in running a time frame of 20 seconds from the initialization. The space-time mesh
adaptation is performed for this 20 seconds time frame by using a single adapted spatial mesh, that is
$n_{adap} = 1$. The space-time complexity values is set to $N_{st} = 750M$.

The final spatial adapted mesh is rather coarse, it is composed of 1M vertices and the final temporal
mesh has 1 700 time steps.

The spatial mesh accuracy is sufficient for the apparition of 3D features, as witnesses the examination
of the Q factor, see Figure 13. Vortices are propagated in the wake without dissipation.

Figure 14 shows the computed time steps as functions of the physical time, i.e., the computed adapted temporal mesh, for the successive iterations of the fixed point. As in the 2D case, the temporal mesh is quickly converged. Two fixed point iterations are sufficient. The periodicity of the flow is again captured in the temporal mesh.

Figure 15 shows the evolution of the spatial $\mathcal{E}_{\text{space}}$ and the time $\mathcal{E}_{\text{time}}$ error functionals for the successive adaptation fixed point iterations. Similarly to the 2D case, we note a fast convergence of the process as only two fixed point iterations are sufficient to reach the optimal errors values. The temporal and spatial error are not perfectly balanced, there is a 20% difference. This is due to the fact that the optimization is performed in a continuous context, producing an optimality condition, which we in turn discretize. The difference between both errors will decrease when the space-time complexity increases.
<figure><img src="image_13.png" /><figcaption>Figure 13: 3D cylinder flow at Reynolds 3900. Different views of the Q-criterion isosurface colored with the velocity magnitude for the final solution.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: 3D cylinder flow at Reynolds 3900. Computed time step as a function of the physical time for the space-time complexity 750M. The convergence of the computed temporal discretization for the successive iterations of the fixed point is depicted.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: 3D cylinder flow at Reynolds 3900. Continuous spatial and temporal errors for space-time complexity 750M at each fixed point iteration. The corresponding total error level at convergence is 0.011.</figcaption></figure>
## 7.5 3D flow past a cylinder at Reynolds 1M

Final, we consider the second test case in 3D. The computational domain is cylindrical of radius 40 diameter of cylinder and a span of 3.14 diameters. The solution is initialized by running the solution for 160 seconds in physical time on an initial radial mesh composed of 920K vertices. The simulation consists in running a time frame of 20 seconds from the initialization. The space-time mesh adaptation is performed for this 20 seconds time frame by using a single adapted spatial mesh, that is $n_{adap} = 1$. The space-time complexity values is set to $N_{st} = 725M$.

The final spatial adapted mesh is rather coarse, it is composed of 1.2M vertices and the final temporal mesh has 1300 time steps.

The spatial mesh accuracy is sufficient to capture the truly 3D flow. Figure 19 shows these 3D structures by plotting the Q-criterion. Vortices and other structures are propagated in the wake without dissipation.

Figure 16 shows the computed time steps as functions of the physical time, i.e., the computed adapted temporal mesh, for the successive iterations of the fixed point. We note that this 3D case is harder to converge, seven fixed point iteration are necessary to converge the temporal adapted mesh. The periodicity of the flow is again captured in the temporal mesh, we clearly see nine periods.

Figure 18 shows the evolution of the spatial $\mathcal{E}_{space}$ and the time $\mathcal{E}_{time}$ error functionals for the successive adaptation fixed point iterations. The convergence of the process is as only four fixed point iterations are sufficient to reach the optimal errors values. The temporal and spatial error are not perfectly balanced, there is a 30% difference again due to the discretization of the optimality conditions. But, the difference between both errors will decrease when the space-time complexity increases.

Figure 17 shows the evolution of the horizontal force applied to the cylinder and demonstrates the impact of the space-time adaptation. As a result, the drag coefficient, starting at a value of 0.58 after the first computation, takes a value of 0.52 after 10 space-time adaptation cycles.

<figure><img src="image_16.png" /><figcaption>Figure 16: 3D cylinder flow at Reynolds 1M. Computed time step as a function of the physical time for the space-time complexity 725M. The convergence of the computed temporal discretization for the successive iterations of the fixed point is depicted.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 17: 3D cylinder flow at Reynolds 1M. Drag coefficient evolution during the simulation time frame for the space-time complexity 725M. The convergence of the drag coefficient for the successive iterations of the fixed point is shown.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 18: 3D cylinder flow at Reynolds 1M. Continuous spatial and temporal errors for space-time complexity 725M at each fixed point iteration. The corresponding total error level at convergence is 0.012.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 19: 3D cylinder flow at Reynolds 1M. Different views of the Q-criterion isosurface colored with the velocity magnitude for the final solution.</figcaption></figure>
# 8 Concluding remarks

During the computation of an unsteady RANS flow with an implicit time advancing, the choice of the time step is delicate and much influences both efficiency and accuracy. Manual prescription of the time-step requires expertise and may lead to inaccurate or expensive simulations when unknown configurations are considered. An automatic prescription of the time steps is therefore of main interest in an industrial context.

In this work, we have presented a Global Space-Time Transient Fixed Point (GSTTFP) method which automatically provides the optimal space-time adapted mesh, i.e., it provides the optimal spatial adapted mesh and the optimal time steps, to minimize the considered space-time error model. Two error models have been considered: the feature-based error model based on a control or the space-time interpolation error in $L^p$ norm of a given sensor and the goal-oriented error model based on the control of the space-time approximation error of a given scalar output functional. This work extends the previous Global Transient Fixed Point methods (GTFP) proposed in [5] for the feature-based approach and [7] for the goal-oriented one. It should be noted that the new GSTTFP method can be directly used in place of the existing GTFP approach resulting in a direct CPU improvement induced by the optimal choice of time step.

The central principle of all these methods is to take into account the different approximation components in the combination of an unique error functional with a unique complexity constraint functional. This determines mathematically the optimal weighting between the different errors thanks to the derivation of optimality conditions. The coupling between the solver and the adaptation is then a discretization of the continuous optimality conditions. That discretization is in turn solved by a fixed-point iteration.

The proposed method, in the feature error model context, has been validated on 2D and 3D simulations of turbulent flow past a cylinder at different regimes. First, we have pointed out the fast convergence of the GSTTFP algorithm as just a few fixed point iterations are required to obtain the optimal spatial and temporal adapted meshes. In all cases, the spatial and temporal errors have been almost balanced by the algorithm. The gap between both error components decreases when the space-time complexity increases. We have also observed that the temporal adapted mesh captures the dynamic of the flow, in particular, the periodicity in time of the flow is clearly visible in the temporal adapted mesh. As regards the accuracy of the method, a space-time convergence order of ~ 1.8 has been observed for the 2D simulations using $n_{adap} = 1$. A higher order of convergence can be obtained by increasing progressively $n_{adap}$, see the analysis in [5].

The perspectives of this work are numerous.

First, the case where $n_{adap}$ higher than unity requires a detailed numerical study for other types of applications involving, for instance, a rapid dynamics, although some results have been presented in this paper. In that context, the variable $n_{adap}$ changes during the mesh convergence study (i.e., when the space-time complexity increases) in order to optimize the space-time mesh and increase the order of convergence according to the theoretical study described in [5, 7].

Second, the proposed error model is presented in a general formulation applying to either feature-based criteria or goal-oriented criteria. The important goal-oriented case needs further study and development as it requires the solution of the backward in time adjoint problem. This will be experimented in a future work.

Lastly, the proposed approach for quasi-steady flow using only one adapted mesh for the entire simulation time frame is very suitable for Large Eddy Simulation (LES) applications. The extension of this work to LES criteria will be considered in the future.
# 9 Acknowledgements

This work was supported by the ANR NORMA project, grant ANR-19-CE40-0020-01 of the French National Research Agency. The authors gratefully acknowledge GENCI for granted access to HPC resources through CINES (grant 2023-A0142A11523) and IDRIS (grant 2022-A0132A05067). The authors also acknowledge Adrien Loseille for providing the local remesher feflo.a.

# References

[1] F. Alauzet and L. Frazza. Feature-based and goal-oriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace. *J. Comp. Phys.*, 439:110340, 2021.

[2] F. Alauzet, P.J. Frey, P.L. George, and B. Mohammadi. 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations. *J. Comp. Phys.*, 222:592-623, 2007.

[3] F. Alauzet, P.J. Frey, and B. Mohammadi. Adaptation de maillages non structurés pour des problèmes instationnaires. *C.R. Acad. Sci. Paris Ser. I*, 303:773-778, 2002. in French.

[4] F. Alauzet, P.L. George, B. Mohammadi, P.J. Frey, and H. Borouchaki. Transient fixed point based unstructured mesh adaptation. *Int. J. Numer. Meth. Fluids*, 43(6-7):729-745, 2003.

[5] F. Alauzet, A. Loseille, and G. Olivier. Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD. *J. Comp. Phys.*, 373:28-63, 2018.

[6] A. Belme, F. Alauzet, and A. Dervieux. An a priori anisotropic goal-oriented estimate for viscous compressible flow and application to mesh adaptation. *J. Comp. Phys.*, 376:1051-1088, 2019.

[7] A. Belme, A. Dervieux, and F. Alauzet. Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows. *J. Comp. Phys.*, 231(19):6323-6348, 2012.

[8] T. Coupez, G. Jannoun, N. Nassif, H.C. Nguyen, H. Digonnet, and E. Hachem. Adaptive time-step with anisotropic meshing for incompressible flows. *J. Comp. Phys.*, 241:195-211, 2013.

[9] A. Dervieux, F. Alauzet, A. Loseille, and B. Koobus. *Mesh adaptation for Computational Fluid Dynamics, t.1 and t.2.* ISTE Ltd and John Wiley & Sons, 1st edition, 2023.

[10] F. Gatti, M. Fois, C. de Falco, S. Perotto, and L. Formaggia. Parallel simulations for fast-moving landslides: space-time mesh adaptation and sharp tracking of the wetting front. *Int. J. for Numer. Meth. Fluids*, 95:1286-1309, 2023.

[11] P.L. George, F. Hecht, and M.G. Vallet. Creation of internal points in Voronoi's type method. Control and adaptation. *Adv. Eng. Software*, 13(5-6):303-312, 1991.

[12] D. Guégan, O. Allain, A. Dervieux, and F. Alauzet. An $L^\infty-L^p$ mesh adaptive method for computing unsteady bi-fluid flows. *Int. J. Numer. Meth. Engng*, 84(11):1376-1406, 2010.

[13] B. Khara, K. Saurabh, R. Dyja, A. Sharma, and B. Ganapathysubramanian. Space-time finite element analysis of the advection-diffusion equation using Galerkin/least-square stabilization. *arXiv:2307.00822 [math.NA]*, https://arxiv.org/pdf/2307.00822.pdf, 2023.

[14] A. Loseille and F. Alauzet. Continuous mesh framework. Part I: well-posed continuous interpolation error. *SIAM Num. Anal.*, 49(1):38-60, 2011.
[15] A. Loseille and F. Alauzet. Continuous mesh framework. Part II: validations and applications. *SIAM Num. Anal.*, 49(1):61–86, 2011.

[16] A. Loseille, A. Dervieux, and F. Alauzet. Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations. *J. Comp. Phys.*, 229:2866–2897, 2010.

[17] A. Loseille and R. Löhner. Cavity-based operators for mesh adaptation. *51st AIAA Aerospace Sciences Meeting, Jan 2013*.

[18] Y. Luo and K. J. Fidkowski. Output-based space-time mesh adaptation for unsteady aerodynamics. *49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition 4 - 7 January 2011, Orlando, Florida, AIAA 2011-491, 2011.*

[19] K. Mani and D. J. Mavriplis. Discrete adjoint based time-step adaptation and error reduction in unsteady flow problems. *18th AIAA Computational Fluid Dynamics Conference 25 - 28 June 2007, Miami, FL, AIAA 2007-3944, 2007.*

[20] S. Micheletti and Simona Perotto. Space-time adaptation for purely diffusion problems in an anisotropic framework. *International Journal of Numerical Analysis and Modeling, 7(1):125–155, 2010.*

[21] S. Micheletti and Simona Perotto. Time adaptivity and anisotropic mesh adaptation for CFD applications. In J. P. Moitinho de Almeida, P. Diéz, C. Tiago, and N. Parés, editors, *VI International Conference on Adaptive Modeling and Simulation ADMOS 2013*, 2013.

[22] P.R. Spalart and S.R. Allmaras. A one-equation turbulence model for aerodynamic flows. In *30th AIAA Aerospace Sciences Meeting and Exhibit, AIAA-92-0439, Reno, NV, USA, Jan 1992.*