J Wackers

R Pellegrini

M Diez

A Serani

M Visonneau

Improving Active Learning in Multi-fidelity Hydrodynamic Optimization

Despite recent advances in machine learning, simulationdriven design optimization using high-fidelity simulations may still be prohibitively expensive for practical applications. This paper investigates improvements in multi-fidelity surrogate-based hydrodynamic optimization, which are intended to make the process faster and more efficient. Specific innovations are: a) the use of a reduced initial dataset with only one data point for all fidelity levels except the lowest, to reduce the computational cost of surrogate-model initialization; b) accounting for noise variance in the selection of the fidelity level to sample, to avoid oversampling well-resolved but noisy fidelity levels, c) improving the automatic mesh adaptation protocol for the CFD simulations, to further optimize mesh adaptation; and d) restarting highfidelity simulations from converged low-fidelity results, to improve the overall efficiency of the design optimization process. These methodological advancements are demonstrated for an analytical test problem, as well as the shape optimization of a NACA 4-digit airfoil and the DTMB 5415 for calm-water resistance. The results show that the reduced dataset drastically reduces the computational cost of initialization and favors efficient low-fidelity exploration of the design space. The noise-corrected fidelity selection encourages the selection of higher fidelities, to effectively determine the true optimum. Finally, the CFD solver advancements make high-fidelity simulations faster, up to eight times. Compared with previous work, the solution of the DTMB 5415 problem exhibits a more robust training process, providing a slightly improved optimum, in less than half the computational time.

INTRODUCTION

Simulation-driven design optimization (SDDO) has emerged also in the naval architecture and marine engi-neering arena as an efficient and effective approach to explore the potentiality of innovative designs, assessing and optimizing the design performance under a variety of operative and environmental conditions. SDDO leverages a tight integration of mathematical models to assess the design performance (from analytical models to numerical simulations) with accurate shape deformation techniques (from morphing methods to parametric CAD models), and with effective optimization algorithms (global/local and derivative-based/derivative free algorithms).

Achieving accurate performance predictions for innovative configurations or off-design conditions requires high-fidelity physics-based models, such as CFD simulations based (at least) on Reynolds-Averaged Navier Stokes equations (RANS) with large grids. Furthermore, the number of iterations required by an optimization algorithm to converge, may imply performing a significant number of simulations, especially if a global optimum is desired. As a consequence, the computational requirements (such as hardware and/or computational time) may easily become unaffordable for most designers/users and/or for highly complex applications. To give an example, recent research showed how an accurate Unsteady RANSE simulation to evaluate statistically significant performance of ship maneuvering in irregular waves may require up to 1M CPU hours on HPC systems, see Serani et al., (2021a).

Supervised machine learning via surrogate modeling methods may be applied to reduce the computational cost associated with SDDO for complex engineering problems in naval architecture and marine engineering, see e.g. Serani et al., (2021b). Among these methods, multi-fidelity approaches are gaining attention, due to their capability to combine the accuracy of high-fidelity solvers with the computational cost of lowfidelity solvers. See for instance the NATO Science and Technology Organization, Applied Vehicle Technology task group AVT-331, which is collaboratively assessing "Goal-Driven, Multi-Fidelity Approaches for Military Vehicle System-Level Design", including applications from air, space, and sea domains [START_REF] Beran | Comparison of Multi-Fidelity Approaches for Military Vehicle Design[END_REF]. Multi-fidelity methods leverage a fidelity spectrum of computational models (from low-to high-fidelity), with the objective of maximizing the model accuracy while minimizing the associated computational cost [START_REF] Fernández-Godino | Issues in deciding whether to use multifidelity surrogates[END_REF]. Different fidelities may be defined based on the physical model [START_REF] Anselma | Multidisciplinary design optimization for hybrid electric vehicles: component sizing and multi-fidelity frontal crashworthiness[END_REF], spatial and/or temporal discretization [START_REF] Song | Multifidelity local surrogate model for computationally efficient microwave component design optimization[END_REF]) (e.g, grid size and/or time step), multidisciplinary coupling (e.g., one-or two-way, tight or loose coupling, etc.) [START_REF] Volpi | Multidisciplinary design optimization of a 3D composite hydrofoil via variable accuracy architecture[END_REF], degree of solution convergence [START_REF] Palar | Global sensitivity analysis via multi-fidelity polynomial chaos expansion[END_REF], accuracy of geometry description [START_REF] Peart | Multi-fidelity surrogate models for VPP aerodynamic input data[END_REF]. Experimental data may also be used as the highest fidelity level [START_REF] Choi | Multi-Fidelity Surrogate Models for Predicting Averaged Heat Transfer Coefficients on Endwall of Turbine Blades[END_REF].

Surrogate-based multi-fidelity methods fuse information from multi-fidelity computations to construct efficiently a surrogate model approximating the highest available fidelity. Examples include polynomial chaos [START_REF] Rumpfkeil | Multi-Fidelity, Gradient-enhanced, and Locally Optimized Sparse Polynomial Chaos and Kriging Surrogate Models Applied to Test Problems[END_REF], Gaussian processes [START_REF] Coppedè | Hydrodynamic shape optimization by high fidelity CFD solver and Gaussian process based response surface method[END_REF], and radial basis functions [START_REF] Piazzola | Comparing Multi-Index Stochastic Collocation and Multi-Fidelity Stochastic Radial Basis Functions for Forward Uncertainty Quantification of Ship Resistance[END_REF] approaches. Multi-fidelity data are usually fused by correcting low-fidelity surrogate models with suitable bridge functions [START_REF] Han | Improving variable-fidelity surrogate modeling via gradientenhanced kriging and a generalized hybrid bridge function[END_REF], representing the difference between fidelity levels, ordered following a given hierarchy. The resulting multi-fidelity surrogate is then explored by the optimization algorithm. To enhance the overall design optimization process, optimization and surrogate-model training may be integrated using active learning approaches [START_REF] Liu | A survey of adaptive sampling for global metamodeling in support of simulation-based complex engineering design[END_REF][START_REF] Di Fiore | Multifidelity domain-aware learning for the design of re-entry vehicles[END_REF]. Furthermore, model auto-tuning capabilities may be used to increase the surrogate-model accuracy, see e.g. [START_REF] Volpi | Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification[END_REF]; [START_REF] Rumpfkeil | Multi-Fidelity, Gradient-enhanced, and Locally Optimized Sparse Polynomial Chaos and Kriging Surrogate Models Applied to Test Problems[END_REF]. Earlier work by the authors on multi-fidelity active learning and model auto tuning may be found in [START_REF] Serani | Adaptive multi-fidelity sampling for CFD-based optimisation via radial basis function metamodels[END_REF][START_REF] Pellegrini | Towards Automatic Parameter Selection for Multifidelity Metamodels[END_REF] and (Wackers et al., 2020;[START_REF] Pellegrini | A Multi-Fidelity Active Learning Method for Global Design Optimization Problems with Noisy Evaluations[END_REF], respectively.

Despite the effectiveness of multi-fidelity approaches, SDDO using high-fidelity simulations still represents a significant investment in computation time. Therefore, it is vital to continuously improve the efficiency of the overall design optimization procedure. The present work aims to do this in two ways: i) increase the efficacy of the MF method presented by [START_REF] Pellegrini | A Multi-Fidelity Active Learning Method for Global Design Optimization Problems with Noisy Evaluations[END_REF] and ii) improve the computational efficiency of the numerical solver.

To achieve the first objective, two methodological innovations are proposed: a) reducing the computational cost associated with the initial training set and therefore increasing the computational budget dedicated to the active learning process by introducing a new initialization approach; the new approach uses only one initial point for all the fidelities except for the lowest; b) improving the efficiency of the active learning approach, by defining new metrics to avoid sampling lower fidelities when they are no longer informative. This is achieved by penalizing the choice of the lower fidelities based on their associated noise variance, while selecting the fidelity to sample. For the second objective, the computational cost of achieving highly accurate simulations is reduced by fine-tuning the adaptive grid refinement procedure in the RANS solver ISIS-CFD [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF]. Furthermore, a restart procedure is introduced, initializing high-fidelity simulations from already existing lowfidelity solutions.

The proposed methodological and software advancements are demonstrated for an analytical test case and two SDDO problems, namely i) the drag-coefficient minimization of a NACA hydrofoil and ii) the calmwater resistance reduction of the DTMB 5415 model. For each problem, a comparison between the authors' previous initialization approach, the new initialization approach without noise penalization, and the new initialization approach with noise penalization is provided.

MULTI-FIDELITY ACTIVE LEARNING METHOD

Consider x ∈ R D as a design variables vector of dimension D. Let the true function f (x) be assessed by N fidelity levels: the highest-fidelity level is f 1 (x), the lowest-fidelity is f N (x), and the intermediate fidelity levels are {f i } N -1 i=2 (x). Using • to denote surrogate model prediction, the multi-fidelity (MF) approximation fi (x) of f i (x) (i = 1, . . . , N -1) is the sum of the lowest-fidelity surrogate and surrogates of the errors (inter-level errors or bridge-functions, ε) between subsequent levels

fi (x) = fN (x) + N -1 k=i εk (x).
(1)

For each i-th fidelity level the training set is T i = {y j , f i (y j )} Ji j=1 , with J i the training set size. The resulting inter-level error training set is defined as

E i = {y j , ε i (y j )} Ji j=1
, where

ε i (y j ) = f i (y j ) -fi+1 (y j ). ( 2
)
Consider the surrogate model providing the prediction and the associated uncertainty. The uncertainty U fN of the lowest-fidelity prediction is considered as uncorrelated with the uncertainty U εk of the inter-level error predictions. Therefore, the uncertainty U fi of the MF prediction can be evaluated as (i = 1, . . . , N -1) [START_REF] Piazzola | Comparing Multi-Index Stochastic Collocation and Multi-Fidelity Stochastic Radial Basis Functions for Forward Uncertainty Quantification of Ship Resistance[END_REF])

U fi (x) = U 2 fN (x) + N -1 k=i U 2 εk (x).
(3)

Stochastic Radial Basis Functions with Least Squares Approximation

Given a (single-fidelity) training set T = {y i , f (y i )} J i=1 , the RBF surrogate model prediction f (x) is computed as the expected value (EV) over a stochastic tuning parameter of the surrogate model [START_REF] Volpi | Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification[END_REF]

, τ ∼ unif[1, 3] f (x) = EV [g (x, τ )] τ , g (x, τ ) = EV [f ] + M j=1 w j ||x -c j || τ , (4)
where w j are unknown coefficients, || • || is the Euclidean norm and c j are the RBF centers, with j = 1, . . . , M and M ≤ J. If the training set is not affected by numerical noise then exact interpolation of the training set is imposed and the coefficients w j are computed by solving

Aw = (f -EV [f ]), with c j = y j (yielding M = J) and f = {f (y i )} J i=1 .
If numerical noise affects the training set then noise reduction is achieved by choosing a number of RBF centers M smaller than the number of training points J, and c j coordinates are defined via kmeans clustering [START_REF] Lloyd | Least squares quantization in PCM[END_REF] of the training points coordinates. Hence, w j are determined with least squares regression by solving w

= (A T A) -1 A T (f -EV [f ]).
The optimal number of stochastic RBF centers (M ⋆ ) is defined by minimizing a leave-one-out cross-validation (LOOCV) metrics (Wackers et al., 2020).

To avoid abrupt changes in the surrogate model prediction from one iteration to the next one, during the active learning procedure the search for M ⋆ can be constrained. In the present work,

M ⋆ k-1 -2 < M ⋆ k < M ⋆ k-1 + 2, with k the active learning iteration. Further- more, a minimum M ⋆ k > 2D + 1 is always imposed. The uncertainty U f (x)
of the SRBF prediction is quantified by the 95%-confidence band using the cumulative density function of g(x, τ ) [START_REF] Volpi | Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification[END_REF].

Initial Training Set and Bounded Surrogate Model

Here the first methodological advancement is described, it is a new approach for defining the initial training set [START_REF] Pellegrini | Towards Automatic Parameter Selection for Multifidelity Metamodels[END_REF]. This new method uses a reduced training set (RS) with only one sample point for all fidelity levels except the lowest, as opposed to the authors' previous work where a full training set (FS) was used. Both the RS and FS approaches can be used with single-or multi-fidelity methods. Table 1 summarizes the RS and FS approaches for the single-and multifidelity cases, respectively.

The challenge for the RS approach is to create a SRBF surrogate model which can handle extrapolation. Using the RS approach, at the first iteration of the active learning the surrogate model prediction is an extrapolation based on the single training point available. Since the SRBF with power kernel has a low accuracy when extrapolating, a bounded surrogate model prediction and the associated uncertainty (both identified with the B subscript) are defined as described in Algorithm 1. When one training point is available, the surrogate model prediction and the associated uncertainty are set equal to the function value in the training point f (x ′ ). This approach is consistent with Eq. 4, where the expected value of the training set is added to the radial basis functions thus providing non-zero prediction when only one training point is available. When more training points are available, the surrogate model prediction and the associated uncertainty are bounded only in regions of the domain far from these training points. Specifically, the definition of U εBi (x), as defined in Algorithm 1, stems from the consideration that the error surrogates represent errors in the multi-fidelity approximation f . Therefore the average error can be used as reference for the surrogate model prediction uncertainty when an extrapolation is performed.

In Algorithm 1 a sigmoid-like function s(r) is used to provide a smooth transition between the SRBF prediction and the bounded prediction

s(r) = 1 1 + e α(r-γ) , (5)
(x) = f (x ′ ) ; U fB (x) = f (x ′ ) ; else if J > 1 then // J training points available fB (x) = f (x) [1 -s(r)] + EV[f ]s(r) ; U fB (x) = min(U f , EV[f ]) ; end else if N > 1 then // Multi-fidelity case if J i = 1, i = 1, . . . , N -1 then // One training point available εBi (x) = ε(x ′ i ) ; U εBi (x) = ε(x ′ i ) ; else if J i > 1, i = 1, . . . , N -1 then // J i training points available εBi (x) = εi (x) [1 -s i (r)] + EV[ε i ]s i (r) ; U εBi (x) = min(U εi , EVε i ]) ; end end
where, for the present work, α = -75 and γ = 0.2. To define r, the smallest hyperrectangle (whose edges are parallel to the Cartesian coordinated axis) containing the training points is defined and r is the Euclidean distance of x from the hyperrectangle boundaries.

Active Learning Method

The multi-fidelity surrogate model is dynamically updated by adding new training points. First, a new training point x ⋆ is identified based on the aggregate-criteria active learning (ACAS, see Fig. 1) presented in [START_REF] Serani | Adaptive multi-fidelity sampling for CFD-based optimisation via radial basis function metamodels[END_REF]. It aims to find points with large prediction uncertainty and small objective function value. Accordingly, ACAS identifies a new training point by minimizing the acquisition function ψ(x)

x ⋆ = argmin x [ψ(x)] , ψ(x) = f (x) -U f (x) + P x (x) P x (x) = 1 ϵ d0-d(x) d0 if d(x) < d 0 , 0 if d(x) ≥ d 0 , (6)
where P x (x) is a penalization factor based on the distance from the existing training sets (all fidelities) to prevent the sampling of already sampled points and having matrix A ill-conditioned. ϵ is a coefficient here set equal to 1E -1, d(x) is the distance of the point x to the closest point identified among all the training sets, and d 0 = 5.0E -3 is the minimum acceptable distance to an existing training point. It may be noted that the ACAS method is a special case of the lower-confidence bounding [START_REF] Cox | A Statistical Method for Global Optimization[END_REF], with equally weighted contributions of f and U f . Once x ⋆ is identified, the fidelity to be evaluated needs to be selected. In this step of the MF-SRBF method the second methodological advancement is introduced. In the authors' previous work the fidelity selection vector is defined as

ϕ ≡ {U ε1 /β 1 , ..., U εN-1 /β N -1 , U fN /β N }, (7)
where β i = c i /c 1 with c i the computational cost associated to the i-th level and c 1 the computational cost of the highest-fidelity. Then, the fidelity level to sample is k = maxloc(ϕ), and the new training point is added to the k-th training set T k and to the lower-fidelity sets from k + 1 up to N . Here, the elements of ϕ for i = 1, . . . , N -1 are redefined as

ϕ i =    U 2 εi -MSEi βi , if MSE i < U 2 εi U εi βi , if MSE i ≥ U 2 εi ( 8
)
and for i = N as

ϕ N =    U 2 fN -MSE N β N if MSE N < U 2 fN U fN β N if MSE N ≥ U 2 fN (9)
where MSE i is the mean-squared error computed for i = 1, . . . , N -1 as and for i = N as

MSE i = 1 J i Ji j=1 [ε i (y j ) -εi (y j )] 2 , (10)
MSE N = 1 J N J N j=1 f N (y j ) -fN (y j ) 2 . ( 11
)
It should be noted that with an optimal least-squares regression of the SRBF, the MSE converges towards the average variance of the noise in the training data (under the hypothesis of zero-mean noise). This approach aims to better distribute the available budget of function evaluations among the fidelity levels. Indeed, when the least-squares regression is accurate for a fidelity and the surrogate prediction uncertainty decreases towards the average noise variance in the training set, it is not interesting to continue sampling that fidelity level, since the surrogate prediction is already accurate compared with the noise affecting the training set. Thus, adding more (noisy) training points will not improve the metamodel accuracy.

FLOW SOLVER AND IMPROVED ADAPTIVE MESH REFINEMENT

When performing multi-fidelity optimization, it is particularly beneficial to use flow solution methods which have a natural capability to produce multi-fidelity data. In this paper, adaptive mesh refinement is used to provide multifidelity data efficiently. The following subsections describe the flow solver ISIS-CFD and its mesh adaptation method. Then a series of tests is used to determine the most efficient mesh adaptation approach for evaluating ship resistance. Finally, it is shown that mesh adaptation allows a solution initialization procedure which produces multi-fidelity data faster than HF-only data.

ISIS-CFD hydrodynamics solver

The CFD simulations are performed with the Navier-Stokes solver ISIS-CFD developed at ECN -CNRS, available in the FINE™/Marine computing suite from Cadence Design Systems. ISIS-CFD is an incompressible unstructured finite-volume solver for multifluid flow. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint transformed into a pressure equation. These equations are similar to the Rhie and Chow SIMPLE method [START_REF] Rhie | A numerical Study of the Turbulent Flow Past an Isolated Airfoil With Trailing Edge Separation[END_REF], but have been adapted for flows with discontinuous density fields. Free-surface flow is simulated with a conservation equation for the volume fraction of water, discretized with specific compressive discretization schemes. The method features sophisticated turbulence models, such as an anisotropic EASM model and hybrid RANS/LES models of the DES family.

The unstructured discretization is face-based. While all unknown state variables are cell-centered, the systems of equations used in the implicit time stepping procedure are constructed face by face. Therefore, cells with an arbitrary number of arbitrarily-shaped constitutive faces are accepted, which enables for example adaptive mesh refinement. The code is fully parallel using the message passing interface (MPI) protocol. A detailed description of the solver is given by [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF]. Information on the interface-capturing scheme can also be found in [START_REF] Wackers | Free-surface viscous flow solution methods for ship hydrodynamics[END_REF].

Mesh Adaptation for Shape Optimization

The computational grids are created through adaptive grid refinement [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF][START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF]. Mesh adaptation is a versatile tool which can be used for example to increase the resolution of details in a simulation, to reduce the simulation costs, and to simplify the automatic creation of suitable grids. The latter two points are essential for automatic shape optimization.

Adaptive grid refinement adjusts the computational grid locally, during the computation, by dividing the cells of an original coarse grid. The decision where to refine comes from a refinement criterion, a tensor field C(x, y, z) computed from the flow. The tensor is based on the water surface position and on second derivatives of the flow variables, which give a crude indication of the local truncation errors. The grid is refined until the dimensions d p,j (j = 1, 2, 3) of each hexahedral cell p satisfy

∥C p d p,j ∥ = T r . ( 12
)
The refinement criterion based on the second derivatives of the flow is not very sensitive to grid refinement [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF], so the cell sizes everywhere are proportional to the constant threshold T r .

For the MF optimization, grid adaptation is used to take into account the need for several fidelities. The interest of this procedure is that different fidelity results can be obtained by running the same simulations and simply changing the threshold T r (Figure 2). Thus, it is straightforward to automate the MF simulations. The grids for the simulation of different geometries are obtained through grid deformation. Each simulation starts from the same original grid (see Figure 3a). This grid is then divided in layers around the hull. For each geometry g(ξ, x), the displacement of the hull faces with respect to g 0 (ξ) is propagated through these layers [START_REF] Durand | Light and flexible Fluid/Structure Interaction, application to sailing boats[END_REF]. The displacements are multiplied with a weighting factor which goes from 1 on the hull to 0 on the outer boundaries, so that the latter are not deformed (see Figure 3b). The original grid is coarse, since deforming coarse cells is easier and safer than the small cells of fine grids. The final grid, including all the refinement at the free surface, is then created using adaptive refinement (see Figure 3c).

Mesh Refinement Protocol for Resistance Evaluation

Recently [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] performed a search for the optimal parameter settings to perform adaptive mesh refinement specifically for ship resistance evaluation. The conclusions of this study are summarized here.

Refinement criterion

To create the fine meshes entirely using adaptive refinement, the refinement criterion must react to all the flow features which are relevant for ship resistance. This implies a combination of a freesurface capturing criterion and the Hessian of the pressure, which is a suitable indicator of the orbital flow fields in the waves. Also the accurate resolution of the boundary layer and wake is needed to correctly predict both viscous and pressure forces, which suggests refine-ment based on the velocity Hessians. Finally, any freesurface dynamics, such as waves but also sprays and jets, requires a good mesh resolution at the surface. The combined free-surface and flux-component Hessian criterion of [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF] is therefore a suitable choice.

Thresholds To formulate generally valid guidelines for the threshold T r in equation ( 12), the refinement criterion is non-dimensionalized using the reference length and velocity. As a result, the choice for T r no longer depends on the ship length and velocity; only a (weak) dependence on F r, Re, and the hull shape remains, which can be ignored for a specific class of ship hulls

For the free-surface refinement criterion, the threshold equals the target cell size at the free surface [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF], L/1000 according to standard ISIS-CFD guidelines. The choice of T r for the Hessian criteria has been studied by running several test cases, two of which are shown here. The first is the KRISO Container Ship (KCS) [START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF] in model-scale towed condition with free trim and sinkage at Re = 1.257 • 10 7 and F r = 0.260. The second is the DTMB 5415 model in towed condition, fixed at the experimental attitude, at F r = 0.28 and Re = 1.19 • 10 7 following the experiments of [START_REF] Olivieri | Towing tank, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[END_REF]. The kω SST turbulence model with a wall law is used for both cases. See [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] for the full test case data.

Figure 4 shows the convergence of the resistance with the Hessian threshold T rH for the two cases, as well as the numerical uncertainty according to Ec ¸a and Hoekstra (2014) (where T rH is used as the measure for the grid size, instead of an averaged cell size) and the evolution of the number of cells. While these two graphs [START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF], the blue line gives the 5415 experimental result [START_REF] Olivieri | Towing tank, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[END_REF].

are not identical, the orders of magnitude of all variables are the same for both cases, so the same range of T rH is appropriate for both. This indicates that generally valid guidelines for T rH can be formulated. Figure 5 shows that the free surface is nearly converged at T rH = 0.05, since the difference with the finest threshold is negligible. Thus, for MF ship optimization T rH = 0.2L, 0.1L, and 0.05L is proposed for coarse to fine grids, where the finest threshold gives around 1% numerical uncertainty for the resistance.

Refinement cutoff filter Adaptive refinement can create very small cells to capture flow details. However, for the accurate evaluation of ship resistance, especially for bare hulls, the smaller flow details may not be required.

The level of resolution that is required can be controlled by specifying a minimum cell size, below which cells are no longer refined. A large minimum cell size works like a cutoff filter: while the adaptive refinement still captures the main flow features, very fine refinement to resolve small details is not applied.

Figure 6 shows the dependence of the KCS grid convergence on the minimum cell size. This parameter (chosen relatively small in earlier tests, around 10 -4 L) has a much more pronounced influence than expected. The figure shows that for a larger minimum cell size, the number of cells is strongly reduced, for only a modest loss of precision. Also, the solutions for the different minimum cell sizes converge to the same value when Tr is reduced. This means that for resistance simulations, the largest minimum cell size tested, i.e. L/1000 like the free-surface threshold, is preferable. Wake field Mesh adaptation has the capacity to accurately resolve the ship wake and the wave field over a long distance behind the ship. Although this is a powerful capability, it may not be needed for resistance computations. Therefore, limiting boxes were tested which forbid horizontal refinement from a given distance behind the stern. Preventing refinement aft of 0.3L behind the stern removes the entire far wake field, but does not significantly alter the forces while reducing the number of cells by up to 40% (figure 7). An even shorter box (0.1L behind the stern) changes the converged resistance, so it is not adopted. The box at 0.3L is the best compromise. Comparison of meshing methods To conclude, the results with adaptive refinement are compared with two other techniques for generating series of meshes. The first approach is to create meshes with Hexpress using the C-Wizard automatic setup tool, which provides uniform cell sizes over the hull. The second series uses a combined Hexpress and systematic grid refinement (SGR) approach, where the three coarsest grids are generated by Hexpress, using manual refinement at the bow and stern to capture the pressure fields. The three finest grids are created from these coarser grids by refining all cells once, using the adaptive refinement algorithm.

Figure 8 compares the three results. The figures show that for both cases, the adapted series converges much faster to its final value than the other two. This means that the same numerical accuracy is obtained on coarser grids: the same distance from the converged value is typically obtained with 1.5 to 2 times fewer cells than the C-Wizard series. The figure also confirms that the MS-FCH and SGR series converge to roughly the same values for the resistance, which underlines the reliability of both series. The C-Wizard series may converge to somewhat lower values however, because the meshes are not similar (in Hexpress, the boundary layer meshes become thinner and thinner as the outer mesh is refined).

The MS-FCH series also performs better for the uncertainty estimation. A variable power-law fit (a sign of good grid convergence) is obtained for all meshes and the behaviour of the different power law fits is similar. Also, the estimated uncertainty is low (around 0.4% on the finest grids). The SGR series is perturbed by the switch from pure Hexpress meshes to SGR, but it also produces power-law convergence with low uncertainty (around 1% on the finest grids). The convergence of the C-Wizard series is perturbed by the change in thickness of the viscous layer: because the convergence is not asymptotic, the uncertainty is high in both cases. Thus, the Hessian-based refinement is the best meshing strategy, both for coarse-mesh and fine-mesh simulations.

Solution Initialization and Convergence Criterion

When computing nested multi-fidelity data, i.e. coarsegrid and fine-grid solutions for the same design parameters, the coarse-grid results can be used to initialize the fine-grid computations to accelerate the computation. If done right, this may reduce the computation time for all the samples to less than the fine-grid simulation time without initialization, i.e. we get the coarse-grid data for (less than) free! However, for general unstructured grids, it is not trivial to initialize a fine-grid solution from a coarsegrid one, since this implies an interpolation between two unrelated grids. Especially in the boundary layer, any imprecision in the interpolation leads to physically incorrect initial solutions (such as oscillating turbulence kinetic energy), which may actually take more time to converge than starting from zero. Adaptive grid refinement solves this problem. In our mesh adaptation algorithm, moving to a finer mesh is a simple matter of restarting a simulation with a smaller threshold T r ; the solution is interpolated naturally when the cells are divided. Thus, the initial solutions make sense and the computations on the finer adapted grids usually converge well.

To benefit from the better initial solutions, a convergence criterion is required to stop the simulations when the forces have stabilized (in earlier optimizations (Wackers et al., 2020), the number of time steps was kept fixed). The convergence criterion in FINE/Marine stops the computations when the force in a given window does not differ more than a specified tolerance from the average solution (computed over a longer window).

An example is provided for the DTMB 5415 at F r = 0.28. The flow is considered converged when in the last 150 time steps, it differs no more than 0.25% from the average over the last 250 time steps. The test compares full 2000 time-step simulations, singlelevel convergence-checked simulations where each grid is started from zero, and multi-level ones where each simulation starts from the next coarsest solution. Results are provided up to a mesh that is twice finer than the optimization HF meshes. Table 2 gives the computation times on a 28core Intel Xeon node, for the computations on each individual grid as well as the cumulative time to obtain nested MF data. The results show that the convergence checker has a dramatic influence: the cumulative time for the single level series is about 2.5 times less than for the 2000 time-step simulations. The multi-level simulations reduces the cumulative time by another factor two. On the two finest grids, the cumulative time is even lower than the single-grid times for the single level approach: with multi-level, evaluating the full MF data set is indeed less costly than computing the HF data only. The multi-level series shows slightly bigger differences, although the global convergence behaviour is similar; the convergence check may not be strict enough for the multi-level series. To be on the safe side, a somewhat stricter convergence test is used for the DTMB 5415 optimization exercise shown below: the computation stops if no point in the last 500 time steps differs by more than 0.2% from the average over that interval.

OPTIMIZATION PROBLEMS

Three implementations of the MF SRBF method are assessed here: the MF SRBF with FS initialization and the selection of the fidelity level to sample based on Eq. 7; the MF SRBF with the RS initialization and the selection of the fidelity level to sample based on Eq. 7; and the MF SRBF with the RS initialization and the selection of the fidelity level to sample based on Eqs. 10 and (11). The three implementations are referred as to FS, RS, and RS-MSE, respectively.

The assessment of the FS, RS, and RS-MSE implementations is based on an analytical test and two CFD-based design optimization problems, with design space dimensions D = 2. Although shape optimization problems usually require a larger number of design variables (in the order of tens), here one-and twodimensional problems are selected for their ease of representation and discussion of the results. Problems are solved with a number of fidelity levels N ranging from 2 to 3. These are reasonable numbers in multi-fidelity shape optimization, where several grid resolution levels and/or physical models may be considered.

In the present work, the term uncertainty always refers to the surrogate prediction uncertainty (see Eq. 3), whereas the noise is associated with the objective function evaluation and intrinsically related to the fidelity level: higher fidelities tend to be less noisy.

A deterministic single-objective formulation of the particle swarm optimization algorithm (Serani et al., 2016), is used for the surrogate-based optimizations, as well as for the solution of the minimization sampling problem of Eq. 6. The optimization is performed with a fixed budget of function evaluations: considering a normalized computational cost of a highest-fidelity evaluation (equal to 1), the overall computational cost CC is proportional to the training set sizes J l and is defined as:

CC = J 1 + N l=2 β l J l . (13)

Analytical Test Problems

This is an analytical test problem affected by artificial numerical noise with D = 2, defined as

minimize f (x) subject to l ≤ x ≤ u, (14)
where l i = 0.3 and u i = 1 (for i = 1, . . . , D) are the lower and upper bound for x, respectively, and f (x) is approximated by N = 2 fidelity levels (f 1 (x) and s 2 (x)) as

f 1 (x) = sin 1 x 1 x 2 f 2 (x) =f 1 (x) -9A 2 cos 1 x 1 x 2 s 2 (x) =f 2 (x) + N (x) (15)
with A 2 = 0.5, N ∼ unif[-0.1R 1 ; 0.1R 1 ] the noise associated to the 2-nd fidelity, and R 1 = 2 the function range of the highest fidelity level. f 1 (x) has two loci with the same lowest value for x 1 x 2 = 2/(3π) and x 1 x 2 = 2/(7π), see Fig. 10a.

The computational cost of the analytical test problem is negligible, therefore an artificial computational cost is defined as β 1 = 1 and β 2 = 0.2. The performance of the method is assessed using N = 2 fidelity levels. A computational budget equal to 100D is used. Since the noise is synthetically added to the analytical functions by a numerical generator of random numbers, a statistical analysis [START_REF] Ficini | Assessing the Performance of an Adaptive Multi-Fidelity Gaussian Process with Noisy Training Data: A Statistical Analysis[END_REF] is performed varying the seed of the random number generator for 25 repetitions.

NACA Hydrofoil

This problem addresses the drag coefficient minimization of a NACA four-digit airfoil. The following minimization problem is solved

minimize f (x) = C D (x) subject to C L (x) = 0.6 and to l ≤ x ≤ u, ( 16
)
where x is the design variable vector, C D and C L are respectively the drag and lift coefficient. The equality constraint on the lift coefficient is necessary in order to compare different geometries at the same lift force (equal to the weight of the object), since the drag depends strongly on the lift. The hydrofoil shape (see Figure 11) is defined by the general equation for four-digit NACA foils [START_REF] Moran | An introduction to theoretical and computational aerodynamics[END_REF]. The upper (y u ) and lower (y l ) hydrofoil surfaces are computed as

       ξ u = ξ -y t sin θ ξ l = ξ + y t sin θ y u = y c + y t cos θ y l = y c -y t cos θ (17) with yc =        m p 2 2p ξ c - ξ c 2 , 0 ≤ ξ < pc m (1 -p) 2 (1 -2p) + 2p ξ c - ξ c 2 , pc ≤ ξ ≤ c ( 18
)
where ξ is the position along the chord, c the chord length, y c the mean camber line, p the location of the maximum camber, m the maximum camber value, t the maximum thickness, and y t the half thickness:

y t = 5t 0.2969 ξ -0.1260ξ -0.3516ξ 2 +0.2843ξ 3 -0.1015ξ 4 . (19)
In this work, the D = 2 design variables vector is defined as x = {t, m} with t ∈ [0.030, 0.120] and m ∈ [0.025, 0.065]. The maximum camber position is fixed at p = 0.4.

The simulation conditions are: velocity U = 10 m/s, chord c = 1 m, fluid density ρ = 1, 026 kg/m 3 , with a chord based Reynolds number Re = 8.41 • 10 6 . The computational domain runs from 11c in front of the leading edge to 16c behind the hydrofoil and from -10c to 10c vertically. Dirichlet conditions on the velocity are imposed, except on the outflow side which has an imposed pressure. The hydrofoil surface is treated with a wall law and y + = 60 for the first layer. Turbulence is modeled with the standard kω SST model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. To obtain the same lift for all geometries (see Eq. 16), the angle of incidence α for the hydrofoil is adjusted dynamically during the simulations.

Tests are run with one, two, and three fidelity levels (N = 1, 2, 3). The optimization budget is fixed at CC = 55. The initial computational grid has 2,654 cells, the refinement threshold value T r is set equal to 0.1, 0.2, and 0.4 from highest-to lowest-fidelity. This results in a cell size ratio of 4 : 1 between the refined fine and coarse grids. The final grids (G) have about 12.8k, 5.7k, and 3.7k cells, respectively (see Figure 2). For consistency with our earlier work (Wackers et al., 2020), computational cost ratios of β 2 = 0.5477 and β 3 = 0.3 are used, although the actual costs have changed thanks to the solution initialization procedure. Statistical information from 9 repetitions is obtained by altering the initial meshes for the adaptive refinement, which changes the noise for each simulation.

DTMB 5415 Model

The shape of the DTMB 5415 destroyer is optimized for minimal resistance R T . The optimization problem reads

minimize f (x) = R T (x) subject to L pp (x) = L pp,0 and to l ≤ x ≤ u, (20)
where L pp,0 = 5.72 m (model scale) is the original length between perpendiculars. The ship is at even keel, with Froude number Fr = 0.30 and Re = 1.18 • 10 7 . The L pp constraint is automatically satisfied by the shape modification method. The modified geometries (g) are produced by the linear superposition of D orthonormal basis functions (ψ) on the original geometry (g 0 ), as follows

g(ξ, x) = g 0 (ξ) + δ(ξ, x), (21) with δ(ξ, x) = D k=1 x k ψ k (ξ), ( 22
)
where ξ are the geometry Cartesian coordinates, whereas -1.25 ≤ {x k } D k=1 ≤ 1.25 and {ψ k } D k=1 are the reduced design variables and the eigenfunctions, respectively, provided by the design-space augmented dimesionality reduction (ADR) procedure described in [START_REF] Serani | Shape Optimization under Stochastic Conditions by Design-space Augmented Dimensionality Reduction[END_REF]. Details about the original design space definition can be found in Serani et al. (2016). In this work, two design variables are used.

The optimization is performed with N = 3 fidelity levels. For the initial sample plane (only for this problem), CFD simulations for all fidelities were run in the center of the domain, and with each design variable at either +1 or -1.

Simulations of the DTMB 5415 are performed on half geometries. The domain runs from 1.5L pp in front of the bow to 3L pp behind the stern, up to 2L pp laterally, and from -1.5L pp to 0.5L pp vertically. Dirichlet conditions on the velocity are imposed on the inflow and side faces, pressure is imposed on the top, bottom, and outflow side. The hull is treated with a wall law and y + = 60 for the first layer. Turbulence is modelled with kω SST. The initial grid has 130k cells. The Hessian thresholds T rH for the simulations with different fidelities are T r = 0.2L, 0.1L, and 0.05L from coarse to fine. This implies a 4 : 1 cell size ratio between the coarsest and finest grids. The free surface mesh is kept the same (L/1000) for all meshes. Again the cost ratios of (Wackers et al., 2020), β 2 = 0.21 and β 3 = 0.06; actual computation times are discussed in the following section.

NUMERICAL RESULTS

The optimization results are assessed by three error metrics (Serani et al., 2016). Knowing the position of the global optimum x, these metrics characterize the normalized error in the design space, the objective function space, and Euclidean distance in the normalized x-f (x) hyperspace, respectively:

E x ≡ ∥x ⋆ -x∥ √ D , (23)
E f ≡ f (x ⋆ ) -f (x) R 1 , ( 24
)
E t ≡ E 2 x + E 2 f 2 , ( 25
)
where x is the array of the design variables (normalized in unit hypercube), x ⋆ is the location of the optimum of the approximation to f , and R 1 is the range of the highest-fidelity level computed considering the initial FS training set. The error metrics E x and E f evaluate design and goal accuracy, whereas the aggregate metric E t evaluates accuracy balancing the performance quantification of the method when optima are in very flat or very peaky portions of the design space. Eq. 24 uses an evaluation of the objective with the highest fidelity level at the point x ⋆ identified by the surrogate as the global optimum.

For the DTMB 5415 model problem the reference optimum is not available, therefore a different set of design-sensitive metrics are employed. These metrics quantify design point location and objective function, respectively:

∆ x ≡ ∥x ⋆ -x 0 ∥ √ D , (26)
∆ f ≡ f (x ⋆ ) -f (x 0 ) f (x 0 ) , ( 27
)
where x 0 is the original objective function value, meaning that ∆ x evaluate the distance of the global optimum position from the original design in the design variable space, whereas ∆ f provides the objective function variation with respect to the parent design. Additionally, the prediction error is used to quantify the error of the surrogate model in predicting the minimum value:

E p = f (x ⋆ ) -f (x ⋆ ) R 1 . ( 28
)

Analytical Test Problem

The challenge of this optimization test is twofold (Fig. 10): first, it has multiple optima located on the loci x 1 x 2 = 2/(3π) and x 1 x 2 = 2/(7π), and second, the LF and HF optima are in different locations, which is a danger for multi-fidelity methods which rely too heavily on LF data.

Considering the sampling approach of our three methods, in Fig. 14a-b the training set sizes are discussed using box plots. The box plot shows the q 1 , q 2 (median), and q 3 quartiles, while the lower and upper whiskers are given by q 1 -1.5IQR and q 3 +1.5IQR, respectively, with IQR = q 3q 1 the interquantile range. Since nested training sets are used, the number of low-fidelity samples is equal to the number of active learning iterations minus five (for the initialization). Thus, the RS approach without the forced HF sampling in the startset leads to a higher number of iterations with respect to FS, allowing a potential better exploration of the design space. Differently, the RS-MSE approach performs a lower number of iterations but uses a higher number of high-fidelity samples in comparison with FS and RS, since LF sampling is discouraged once the LF uncertainty approaches the noise level. The evolution of the sample sizes (Fig. 15c) shows the higher number of LF samples at the start (below CC = 10) for both RS and RS-MSE, indicating that the early exploration is performed with low-fidelity samples only. This is followed by a larger increase of HF samples for RS-MSE. Still, RS and RS-MSE have larger whiskers than FS, showing a significant variability and less consistency in the results.

Figure 13 shows the box plots of the three error metrics. The RS and RS-MSE approaches achieve lower median values than the FS approach, while RS-MSE performs the best; its median result is an almost exact optimum. Finally, the FS approach achieves the smallest IQR. Thus, FS is consistent, but it is consistently wrong: the information in the large initial sample set forces the optimization into a fixed, but suboptimal direction. RS with its greater freedom performs better, but its reliance on LF data also leads to some bad results, since the LF optimum does not correspond to the HF one; hence the large IQR for this approach. Finally, RS-MSE with its initial LF sampling and final emphasis on HF samples provides more consistently good results.

In Fig. 15 the placement of HF samples with respect to x 1 x 2 is studied: the quantity c i is the number of elements in each bin of width 0.017. The figure shows that FS requests most of its high-fidelity training point close to the x 1 x 2 = 2/(3π) locus and only a small quantity in the neighborhood of the x 1 x 2 = 2/(7π) locus, negatively affecting the final performance. The RS approach request almost the same quantity of highfidelity training points between the two loci, in the position of the LF optimum. The RS-MSE approach requests for the highest number of high-fidelity training points in the x 1 x 2 = 2/(3π) locus and almost the same number of high-fidelity samples in the second locus as the RS approach. This explains why the RS-MSE approach achieves the best performance overall. cases (see the lower whisker). It is, as yet, not clear why this is the case.

NACA Hydrofoil SDDO Problem

The small computational cost of the NACA hydrofoil simulations allows to perform a statistical analysis of the problem with a few repetitions of each optimization. Specifically, 7 repetitions are performed, randomness in the repetitions is introduced by slightly translating the hydrofoil within the computational domain to perturb the adaptive grid refinement process.

Figure 16 shows the box plots of the E x , E f , and E t metrics. The RS approach achieves the lowest position error and RS-MSE achieves a lower error than FS. RS achieves the smallest IQR and RS-MSE the largest. Nevertheless, the q 3 of the RS-MSE is lower than the median achieved for the FS approach.

Figure 17 shows the response surfaces and the training sets of the optimizations which provide the medians of E t , for the FS, RS, and RS-MSE approaches.

All the approaches are effective in identifying the region of the minimum, but the RS and RS-MSE provide a better representation of the region than FS. Both RS and RS-MSE focused the medium-and high-fidelity samples in the neighborhood of the minimum. Figure 18 shows the training sets size versus computational cost. Like the analytical test case, the RS and RS-MSE approaches use a higher number of lowfidelity samples in the beginning, which suggest lowfidelity exploration of the design space. The RS approach then increases the number of medium-fidelity samples for CC>10. Both the RS and RS-MSE approaches use a very limited number of high-fidelity samples. However, contrary to the analytical test case, RS-MSE has less HF samples than RS.

Thus, for the NACA SDDO problem, the RS-MSE approach performs worse than RS. This could be caused by the limited number of low-fidelity samples in comparison with the analytical test problem. Another explanation is that the MSE becomes larger than the uncertainty due to the significant noise for this test case, which deactivates the procedure (Eqs. (8)-( 9)). As a preliminary check, Eqs. (8) and (9) were changed to have zero uncertainty (e.g.,

U εi = 0 or U fN = 0) if M SE i > U 2
i ; the results show that the median numbers of high-, medium-, and low-fidelity samples change from 2, 13, 134 to 6, 15, 96, respectively. Thus, this modification seems effective in forcing the active learning towards higher fidelities in presence of noise. Nevertheless, when applied to the analytical test with two fidelities it completely prevented the sampling of the lower fidelity. The optimal formulation of RS-MSE is therefore an open problem.

DTMB 5415 SDDO problem

Since this optimization of 3D free-surface flows is characterized by high computational costs, the problem is solved only with the RS-MSE approach. The comparison with the FS approach is provided taking the results presented in (Wackers et al., 2020). Since both the adaptive metamodelling strategy and the CFD simulations have changed, this does not provide a detailed assessment of one topic; rather, the comparison globally shows the progress that has been achieved in the last two years.

Figure 19a-b present the multi-fidelity surrogate models at the last iteration of the active learning approach, the MF datasets, and the predicted optima. The sampling strategies for the two approaches are radically different: RS-MSE performed an exploration of the domain using only low-fidelity samples, correctly identifying the minimum region. The precision in this region is then increased using mainly medium-fidelity evaluations; only two high-fidelity points are sampled, one of which is almost in the optimum location. Near the end of the sampling, most points are added around the optimum. FS on the other hand, uses more HF points spread around the parameter space. Not all these points are useful; note for example the set of points in the top left corner, where a second minimum was suspected in the initial stages of the sampling. And while the data points are clustered, none are placed directly around the optimum.

These differences are also reflected in the x ⋆ convergence (Fig. 19c). The CC of evaluating the startset is 7.35 for FS and only 1.24 for RS-MSE which, combined with the efficient initial LF exploration, means that RS-MSE has globally identified the optimum before FS finishes half its startset. The subsequent RS-MSE convergence is fast and without oscillations, as medium-and high-fidelity points are added around the optimum. The optimization has converged around CC = 15. The FS convergence is much more irregular, as it identifies two incorrect optima before finally settling on the correct one around CC = 24. This explains the lack of sampling in the optimum position (Fig. 19b).

Table 3 summarizes the performance of the FS and RS-MSE approaches. Although the CC is lower, the RS-MSE approach correctly identified the region of the minimum, using more low-and medium-fidelity data than the FS approach. The prediction error is twice lower, which confirms that the metamodel is accurate arround the optimum. The ∆ x % value is larger for RS-MSE than for FS, meaning that the exploration for the identification of the minimum moved further. Finally, the RS-MSE approach achieves a lower resistance than the FS approach.

Furthermore, the computational costs per simulation have decreased significantly. For the RS-MSE data points, Fig. 20 evaluates the number of cells, time steps, and simulation times on a 20-core Intel Xeon workstation. The latter, with a median value of 92 (mean 92), 107 (mean 113) and 129 (mean 127) minutes for coarse, medium, and high fidelity, can be compared with the mean time of (Wackers et al., 2020) which is about 70, 240, and 1140 minutes. There are two reasons for this acceleration:

• The number of cells in the meshes has decreased, for the same precision, thanks to the new simu- lation protocol and notably the higher minimum cell size. For the HF simulations, it has gone from 3.4M to 1.7M and for the medium fidelity, from 860k to about 700k. Only for LF, there is an increase from 280k to 400k because the free surface is meshed finer than before. Especially for HF, the gain is significant.

• The convergence checker and the solution initialization reduce the number of time steps required.

While this has little influence on the LF simulations, which generally require the full 2000 time steps, the median for the (highly variable) medium fidelity is 1152 time steps, while the HF requires only 600 time steps for convergence.

Of course, the different innovations interact. For example, the finer low-fidelity free-surface mesh may reduce the noise in the LF data and thus speed up the convergence of the optimization. Also, both the faster convergence and the lower times per simulation, as well as code optimization of the active learning algorithm, contribute to an increased efficiency. Altogether, the RS-MSE optimization produces a similar optimum as the old FS result, in a more robust manner, for a total wall clock time that has reduced from 25 to about 11 days.

CONCLUSIONS AND FUTURE WORK

This paper focuses on improving the efficiency of a multi-fidelity (MF) method for simulation-driven design optimization. The efficiency and accuracy of the MF surrogate approach is increased by: a) Reducing the computational cost of the initial training set by a new reduced startset (RS) with one single point for all the fidelities except the lowest one (and less expensive to evaluate).

b) Improving the fidelity sampling, by subtracting the mean-squared error (MSE) between the training set and the least-squares regression, from the prediction uncertainty on each fidelity level while selecting the fidelity to sample. Indeed, when the least-squares regression is accurate and the surrogate prediction uncertainty decreases towards the average noise variance in the training set it is less interesting continuing sampling that fidelity level.

Furthermore, the computational efficiency of the RANS solver ISIS-CFD is increased by: c) Optimizing the adaptive mesh refinement protocol, such that the same simulation accuracy can be obtained with up to 2 times less cells than before.

d) Implementing a restart procedure for evaluating high-fidelity solutions from already existing lowfidelity simulations and a convergence checker to reduce the number of time steps required for the convergence of the simulation.

The proposed methodological and software advancements are assessed with an analytical test problem and two SDDO problems: the drag minimization of a NACA hydrofoil and the DTMB 5415 model. Numerical results show that the RS approach is effective in improving the identification of the minima and better distributing the high-fidelity samples in interesting regions of the domain. The reasons for this are twofold. First, the reduced startset encourages early exploration with low-fidelity samples, unconstrained by high-fidelity data which may wrongly indicate sub-optimal regions. This early exploration leads to a reliable identification of the minimum region. Later on, thanks to the limits imposed on the extrapolated highfidelity uncertainty, sampling of high-fidelity data in the observed minimum region is encouraged, which efficiently increases the precision in the minimum region. This behavior is observed for all test cases.

For the analytical test case, the RS-MSE approach achieves even better results than the RS approach, since its formulation is effective in forcing the fidelity selection towards the high-fidelity after a large number of low-fidelity evaluations have been performed. This is necessary for this case, since the low-and high-fidelity optima are in different locations. For the NACA SDDO problem, the RS-MSE approach performs worse than RS, either due the limited number of low-fidelity samples or due to the MSE becoming larger than the uncertainty due to the significant noise for this test case. This implies that the RS-MSE method, despite its potential, requires modifications to be suitable for all test cases. Finally, the DTMB 5415 SDDO problem shows that the RS-MSE method is highly effective in identifying the region of the minimum, performing high-and medium-fidelity evaluations almost exclusively in this region. Although few high-fidelity points are added, it is likely that the sampling with (accurate yet cost-effective) medium-fidelity points is encouraged by the RS-MSE.

The new simulation protocol, convergence checker and multilevel simulation are demonstrated as very effective: for the DTMB 5415, they allow a reduction of the wall-clock time for a medium-and highfidelity simulation of 50% and 88% respectively. This confirms the conclusion in [START_REF] Pellegrini | A Multi-Fidelity Active Learning Method for Global Design Optimization Problems with Noisy Evaluations[END_REF], that it is advantageous to use MF optimization with CFD simulation methods which naturally produce multi-fidelity data. The total wall-clock time for the DTMB 5415 optimization was reduced from 25 to about 11 days.

In conclusion, the test results are promising, showing that the reduction of initial training set is beneficial and that an effective approach to limit the sampling of the lower fidelities contributes to improve the MF method performance. Nevertheless, three main issues need to be addressed: i) the robustness of RS-MSE can be improved, so that it becomes beneficial for any test case; ii) the active learning approach shows significant clustering of the samples, because the uncertainty of the noise-filtered LS-RBF may not be reduced when a sampling point is added to an already densely sampled area; and iii) the selection of the number of kernels for the RBF needs to be improved, since the analytical test problem showed that the MF method overestimates the noise variance in the training set. While this is not an issue for the RS-MSE approach, since it would mean that the method is "more" effective when enough data are available, the accuracy of the data regression is obviously important.

Thus, future work will focus among others on improving the RBF kernel estimation to improve the regression of noisy data and a new acquisition function, to improve the active learning process and reduce data clustering.

DISCUSSION

Pr. Markus P. Rumpfkeil, Hans von Ohain Endowed Chair, Department of Mechanical and Aerospace Engineering, University of Dayton, Dayton, Ohio. Markus.Rumpfkeil@udayton.edu

The paper is well written and demonstrates that by using various acceleration strategies for a high-fidelity CFD solver in combination with improvements to the underlying machine learning strategy significant computational savings can be achieved for performing optimization. It would be great to see the performance of the framework for more input design variables in future work. Also, the multimodal analytical test function example showed that the suggested improvements are not always beneficial and it would be nice to see this algorithm applied to a wider range of test cases to be able to better assess the overall utility. Lastly, does the ISIS-CFD solver have adjoint capabilities or are they being developed? If so, adjoint-based mesh adaptation can replace a lot of the heuristics currently used for the adaptation. Additionally, adjoint-gradients would enable more powerful gradient-based optimization methods (which scale much better with the number of design variables) or can possibly enhance the construction of the surrogate model itself.

AUTHORS' REPLY

We thank the discusser for his kind words about our work and for his interesting questions.

The discusser's suggestion to assess the proposed approaches for a larger number of design variables and for a wider range of test cases is of great interest and will be completely addressed in future work.

Obviously, the performance of the multi-fidelity metamodeling methods are problem dependent and not always better than the single high-fidelity counterpart. Our approach to overcome this issue is to work on increasing the adaptivity of the proposed multi-fidelity method, in order to make it more and more robust to different problems. Since the RS approach is developed to give as much freedom as possible to the adaptive sampling procedure, it is expected that its performance will remain promising also on a wider range of test cases.

When increasing the number of variables, it is expected that the advantage of using the RS approach will become more evident. Indeed, the bounded values of the metamodel uncertainty away from the training set do not change with the design space dimensions, therefore the likelihood of sampling the corners will not increase. Differently, with the FS approach the metamodel uncertainty will increase in the corners as the design space dimension increases. This is because the distance of the corners from the FS initial training set increases with the dimensions. Therefore, with a larger uncertainty in the corners it is more likely that the active learning method samples in the corners, eventually "wasting" objective function evaluations. Differently, the RS approach is shown to place the high-fidelity samples almost exclusively in the region of the optimum. To validate this point a test will be performed and included in the final version of the paper, testing the Paciorek function also with 5 and 10 variables.

Goal-oriented refinement criteria are ideal for our mesh adaptation applications. However, the difficulty with adjoint solvers for naval applications is the presence of the free surface, which should be taken into account for the adjoint solution. While this is possible, it is a difficult matter1 . We have performed initial investigations using ISIS-CFD but have not been successful so far. The current refinement criteria for ISIS-CFD are not goal-oriented, but they are not fully heuristic either: the Hessian tensors provide an indication of the truncation errors. Comparisons have shown that such feature-based criteria form a reasonable alternative when adjoint-based refinement criteria are not available2 . Furthermore, the current method is easier to apply to time-dependent problems than adjoint approaches (although this paper does not exploit this capacity).

Regarding the optimization algorithm, for this work we aimed at a method which does not require any specific input from the CFD solver, other than force values. Indeed, the objective is to produce a metamodelling and optimization approach which automatically adjusts itself to the behaviour of these data, independently of their origin, to remain efficient for higher dimensions through metamodel adaptivity. We think that this approach is interesting for the many practical situations where adjoint solutions are not available. Still, as the discusser mentions, adjoint solutions could make a positive contribution to all aspects of the method we presented.

DISCUSSION

Dr. Thomas Scholcz, MARIN, Wageningen, The Netherlands t.p.scholcz@marin.nl A very novel computational method is described that aims to improve multi-fidelity active learning applied to hydrodynamic optimisation problems. Especially practical aspects such as how to treat noise in objectives and the use of a budget driven multi-fidelity strategy makes this research highly relevant to the hydrodynamic research community.

Q1: In Equation (2), the training set is denoted by y j instead of x j . From the context it is clear that the vectors y j are points in R D that form the training set for each fidelity level. However, the vectors y are not introduced and as f i is a function of x, one would expect that these points are denoted by x j . This may cause confusion while reading the section. Is there a specific reason to use y j instead of x j ? Q2: On page 3, an extrapolation procedure is described when only a single training point is available. The prediction and uncertainty are kept constant in this case. The quality of this approach will deteriorate with the distance in D w.r.t to this point. This distance will grow with increasing number of dimensions. Will the advantage of the RS approach w.r.t the FS maintain for a high number of dimensions?

Q3: Using one-leave-out cross-validation to find hyper parameters introduces a dependency on the initial DoE. For example, a DoE can be too regular for cross-validation to work well. Could it be that the speedup of the RS approach w.r.t the FS approach is partly caused by the fact that points are added in an irregular fashion, thereby facilitating the hyper parameter procedure in a very early stage of the optimisation?

AUTHORS' REPLY

We thank the discusser for his thoughtful comments and for several interesting questions.

Regarding Q1, two different symbols are used to denote coordinates because y j indicates sample locations, which are fixed data once each sample is per-formed. On the contrary, x are coordinates which vary during operations on the metamodel, such as minimum searches. However, we have used x j in other papers to denote sample locations, and following the discusser's suggestion we will probably do so again in the future.

Q2 concerns an important issue. The constant RS extrapolation is based on the idea that the error metamodels represent numerical errors, which are more or less constant for similar designs. Since the distance from the domain center to its corners increases for higher dimensions, the accuracy of this assumption diminishes indeed. However, with respect to RS, the FS approach only adds data in the face centers, whose distance to the corners increases in the same way. Therefore, its accuracy may deteriorate as much as for RS.

And furthermore, the RS extrapolation does not have to be accurate. The main purpose of RS (besides reducing the cost of the startset) is to encourage a thorough LF exploration of the domain before adding any HF points. Contrary to the (unlimited) FS uncertainty, whose values in the domain corners grow for higher dimensions, the limiting of the RS uncertainty ensures that no HF sampling will take place in the corners. This idea to concentrate HF points only around the optimum region instead of covering the entire design space, becomes more important in higher dimensions, where the design space is large. Therefore, we expect that the advantage of RS will actually increase for higher dimensions. We will investigate this in the near future.

Finally, regarding Q3 the discusser is right, the initial DoE affects the one-leave-out cross-validation as well as the active learning procedure. Nevertheless, it is worth mentioning that the initial number of samples is small (2D + 1) for the FS approach and only these initial points are placed in regular fashion using a central composite design. Therefore, as other points are added by the active learning method, the samples distribution will rapidly become mostly irregular for the FS approach as well as for the RS. It is our opinion, that the speed-up of the RS approach is mainly due to the better handling of the high-fidelity samples, as discussed in the answer for Q2. Furthermore, with a fixed budget the use of the RS approach allows more active learning iterations, thus increasing the chance to improve the results.

Figure 1 :

 1 Figure 1: Example of the active learning method using one fidelity: (a) shows the initial surrogate model with the associated prediction uncertainty and training set; (b) shows the position of the new training point and the new surrogate model prediction and its uncertainty.

Figure 2 :Figure 3 :

 23 Figure 2: NACA hydrofoil computational grids (G) for ISIS-CFD: (a) G1, 12.8k cells, (b) G2, 5.7k cells, and (c) G3, 3.6k cells.

Figure 5 :

 5 Figure 5: DTMB 5415, free surface at F r = 0.28. Top to bottom: TrH = 0.2L, 0.1L, 0.05L, and 0.025L.

  (a) KCS (b) DTMB 5415

Figure 4 :

 4 Figure 4: Resistance, grid convergence with TrH and estimated numerical uncertainty. The box shows the uncertainty interval of the KCS experimental results[START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF], the blue line gives the 5415 experimental result[START_REF] Olivieri | Towing tank, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[END_REF].

Figure 6 :

 6 Figure 6: KCS, influence of the minimum cell size on the convergence of the resistance.

Figure 7 :

 7 Figure 7: KCS, convergence as a function of TrH for three different limiting box sizes.

Figure 8 :

 8 Figure 8: Convergence in terms of the number of cells for three meshing strategies and estimated uncertainties. The x-axis uses a log scale.

Figure 9 :

 9 Figure9: DTMB 5415 at F r = 0.28, grid convergence for three simulation approaches.

Figure 9

 9 Figure9compares the resistance for the three series. The 2000-step simulations and the single level series are very close, even the uncertainty estimations are similar. The multi-level series shows slightly bigger differences, although the global convergence behaviour is similar; the convergence check may not be strict enough for the multi-level series. To be on the safe side, a somewhat stricter convergence test is used for the DTMB 5415 optimization exercise shown below: the computation stops if no point in the last 500 time steps differs by more than 0.2% from the average over that interval.

  Figure 10: Analytical test problem without artificial numerical noise. The dashed lines show the two loci with the lowest f1(x) value (the x1x2 = 2/(7π) locus is in the neighborhood of the bottom-right corner).

Figure 11 :

 11 Figure 11: NACA 4-digit hydrofoil.

Figure 13 :

 13 Figure 13: Analytical test, box plots of the Ex, E f , and Et metrics.

Figure 14 :

 14 Figure 14: Analytical test, box plot of the training sets size and training sets sizes versus computational cost.

Figure 15 :

 15 Figure 15: Analytical test, histogram of the x1x2 quantity for the high-fidelity training set (with ci the number of elements in each bin).

Figure 12 :

 12 Figure 12: Analytical test, box plot of the MSE value of the low-fidelity training set.Finally, Fig. 12 shows the value of the predicted MSE for the low-fidelity training set with respect to the

Figure 16 :

 16 Figure 16: NACA problem, box plots of the Ex, E f , and Et metrics.

Figure 17 :

 17 Figure 17: NACA problem, response surface and training sets for the three approaches. The optimization corresponding to the median Et value is selected for each approach.

Figure 18 :

 18 Figure 18: NACA problem, median of the training sets size versus CC.

  Figure 19: DTMB 5415 SDDO problem, multi-fidelity surrogate model prediction and x ⋆ convergence for the FS and RS-MSE approaches.

( a )

 a Number of cells (b) Time steps (c) Simulation time

Figure 20 :

 20 Figure 20: DTMB 5415 SDDO problem, box plots for the number of cells, time steps, and simulation time.

Table 1 :

 1 Comparison between the new reduced training set (RS) and the full training set (FS).

	Approach	N 1	Fidelity level i = 1	|J i | 1	Training points placement Center of domain
	RS	> 1	i = 1, . . . , N -1 i = N	1 2D+1 Center of domain and of the boundaries Center of domain
		1	i = 1	2D+1 Center of domain and of the boundaries
	FS	> 1	i = 1, . . . , N -1 2D+1 Center of domain and of the boundaries i = N 2D+1 Center of domain and of the boundaries

Table 2 :

 2 Computation time (in minutes) for three simulation approaches. The number of time steps is given between parentheses. Top table: time to compute only on the given mesh. Bottom table: all the meshes up to the given mesh.

	Tr,H	Full simulation	Single level	Multi-level
	Current grid only
	0.2L	47.5 (2000)	20.4 (915)	20.4 (915)
	0.1L	94.8 (2000)	35.1 (857)	15.7 (356)
	0.05L	296.4 (2000) 144.6 (1107)	54.3 (358)
	0.025L	767.6 (2000)	279.9 (910)	183.9 (466)
	Cumulative
	0.2L	47.5 (2000)	20.4 (915)	20.4 (915)
	0.1L	142.3 (4000)	55.5 (1772)	36.1 (1271)
	0.05L	438.7 (6000) 200.1 (2879)	90.4 (1629)
	0.025L	1206.3 (8000) 480.0 (3789) 274.3 (2095)

Table 3 :

 3 DTMB 5415 SDDO problem, summary of the results.

	Approach CC	x1	x2	∆x% ∆ f % |Ep|% J1 J2	J3
	FS	24.0 0.5506 0.1330	26.2	-4.5	1.73	16 18	72
	RS-MSE	18.4 0.5043 0.1525	37.2	-4.9	0.87	3	44 103

N. Kühl et al., Adjoint complement to the volume-of-fluid method for immiscible flows, J. Comput.Phys. 440, pp. 110411 (2021).

F. Alauzet and L. Frazza, Feature-based and goal-oriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace, J.Comput. Phys. 439, pp. 110340 (2021).

ACKNOWLEDGMENTS CNR-INM is partially supported by the Office of Naval Research through NICOP grant N62909-21-1-2042, administered by Dr. Elena McCarthy and by Dr. Woei-Min Lin of the Office of Naval Research Gobal and the Office of Naval Research, respectively.

The work is conducted in collaboration with NATO STO AVT task group on "Goal-driven, multi-fidelity approaches for military vehicle system-level design" (AVT-331). ECN/CNRS received access to HPC resources of CINES under the allocation A0092A01308 made by GENCI (Grand Equipement National de Calcul Intensif), which is gratefully acknowledged.