Dimension reduction in heterogeneous parametric spaces with application to naval engineering shape design problems
Abstract
We present the results of the first application in the naval architecture field of a methodology based on active subspaces properties for parameters space reduction. The physical problem considered is the one of the simulation of the hydrodynamic flow past the hull of a ship advancing in calm water. Such problem is extremely relevant at the preliminary stages of the ship design, when several flow simulations are tyipically carried out by the engineers to assess the dependence of the hull total resistance on the geometrical parameters of the hull, and others related with flows and hull properties. Given the high number of geometric and physical parameters which might affect the total ship drag, the main idea of this work is to employ the active subspaces properties to identify possible lower dimensional structures in the parameter space. Thus, a fully automated procedure has been implemented to produce several small shape perturbations of an original hull CAD geometry, in order to exploit the resulting shapes to run high fidelity flow simulations with different structural and physical parameters as well, and then collect data for the active subspaces analysis. The free form deformation procedure used to morph the hull shapes, the high fidelity solver based on potential flow theory with fully nonlinear free surface treatment, and the active subspaces analysis tool employed in this work have all been developed and integrated within SISSA mathLab as open source tools. The contribution will also discuss several details of the implementation of such tools, as well as the results of their application to the selected target engineering problem.
Contents
1 Introduction
Nowadays engineering simulations present a wide range of different parameters. When it comes to find an optimal solution with respect to the physical constraints it is easy to be affected by the curse of dimensionality, when the number of parameters makes the simulation unfeasible. This problem arises quite easily even with a small parameters space dimension (depending on the simulation, even ten parameters could take months to be optimized). In this framework reducing the dimension of this space becomes crucial and a priority. Some techniques arose in the last years to tackle it, such as the methods grouped under the name of manifold learning algorithms and the active subspaces approach.
Manifold learning is an approach to nonlinear dimensionality reduction. They try to identify a lower dimension structure in the behaviour of the output of interest with respect to the parameters involved [47]. The most interesting algorithms are isometric mapping (isomap) and local linear embedding (LLE [39]). Engineering applications include manufactoring processes (see [27], where they used a manifold walking algorithm, and [38]), mechanical tests [30] and structural optimization problems [37].
We will focus instead on the active subspaces properties (see [11]) to carry out a technique applied on a naval engineering problem, that is the computation of the total wave resistance of a hull advancing in calm water. In the framework of simulationbased design and shape optimization we cite among others [17, 46, 48, 15]. The computational pipeline we are going to present is composed first by a geometrical parametrization and deformation of the hull through free form deformation (see [44]). Then the use of a high fidelity solver based on Boundary Elements Method (BEM) to get the wave resistance with respect to the geometrical parameters. We consider also a structural parameter — the initial displacement of the hull — and a physical one — the velocity of the hull —. Subsequently active subspaces are identified thanks to the data collected from the high fidelity solver, and finally a proper reduced response surface is constructed. The result allows the final user to have an estimate of the wave resistance under a certain threshold and within a time of one second with respect to the hours needed for a single classic full simulation. Moreover during the process is possible to identify the most important parameters and have insights on how they influence the output of interest. Figure 1 summarizes the proposed computational pipeline.
In section 2 we present the active subspaces properties and its features, with a numerical recipe to identify them. In section 3 we recall the free form deformation technique and we show the main features of the developed tool to manage parametric shapes. Section 4 aimes the purpose of introducing what high fidelity solver we used and how. Then section 5 shows the numerical results obtained by coupling the three methods in sequence. Finally conclusions and perspectives are drawn in section 6. Section 7 is an appendix that clarifies the pipeline in a more simple context.
2 Active subspaces properties
The active subspaces (AS) approach represents one of the emerging ideas for dimension reduction in the parameter studies and it is based on the homonymous properties. The concept was introduced by Paul Constantine in [11], for example, and employed in different real world problems. We mention, among others, aerodynamic shape optimization [29], the parameter reduction for the HyShot II scramjet model [13] and active subspaces for integrated hydrologic model [24].
A characteristic of the active subspaces is that they identify a set of important directions in the space of all inputs, instead of identifying a subset of the inputs as important. If the simulation output does not change as the inputs move along a particular direction, then we can safely ignore that direction in the parameter study. In Figure 2 it is possible to capture the main idea behind the active subspaces approach: we try to rotate the inputs domain in such a way lower dimension behavior of the output function is revealed. When an active subspace is identified for the problem of interest, then it is possible to perform different parameter studies such as response surfaces [6], integration, optimization and statistical inversion [25].
There are some main ingredients in order to employ active subspaces. The first is a scalar function smooth enough depending on the inputs x that represents the quantity of interest. Moreover we need the gradients of this map with respect to the inputs or an approximation of them, a sampling density , and a gap between eigenvalues of the covariance matrix associated to the gradients. That is a symmetric, positive semidefinite matrix whose elements are the average products of partial derivatives of the simulations’ input/output map, that reads
(1) 
where is the expected value, a probability density function — usually a uniform one —, and the socalled uncentered covariance matrix of the gradients of (for a more deep understanding of these operators see for example [16]). Usually a Monte Carlo method (see [31]) is used in order to approximate the eigenpairs of this matrix (see [9]) as follows
(2) 
where we draw samples independently from the measure and where . This matrix is symmetric and positive semidefinite, so it admits a real eigenvalue decomposition
(3) 
where is a column matrix of eigenvectors, and is a diagonal matrix of eigenvalues. Then we order the eigenpairs in descending order. We will select the first eigenvectors to form a reducedorder basis. This is where we reduce the dimensionality of the design problem. We will attempt to describe the behaviour of the objective function by projecting the fullspace design variables onto this active subspace. On average, perturbations in the first set of coordinates change more than perturbations in the second set of coordinates. Low eigenvalues suggest that the corresponding vector is in the nullspace of the covariance matrix, and to form an approximation we can discard these vectors. We select the basis as follow
where with , and contains the first eigenvectors. The active subspace is the span of the first few eigenvectors of the covariance matrix. We define the active variables to be the linear combinations of the input parameters with weights from the important eigenvectors. In particular we define the active subspace to be the range of . The inactive subspace is the range of the remaining eigenvectors in . With the basis identified, we can map forward to the active subspace. So is the active variable and the inactive one, respectively:
(4) 
In particular any point in the parameter space can be expressed in terms of and :
Having this decomposition in mind we can rewrite
and construct a surrogate model discarding the inactive variables
We underline that the size of the eigenvalue problem is the limiting factor. We need to compute eigenvalue decompositions with matrices, where is the dimension of the simulation, that is the number of inputs.
Active subspaces can be seen in the more general context of ridge approximation (see [36, 26]). In particular it can be proved that, under certain conditions, the active subspace is nearly stationary and it is a good starting point in optimal ridge approximation as shown in [12, 23].
Now we are going to briefly present the free form deformation for shape parametrization and deformation, and why we have chosen this particular technique. Moreover we underline the improvements and contributions in this research direction.
3 Free form deformation techniques
As already mentioned, we are interested in problems characterized by both physical and geometrical parameters. Concerning the geometrical ones, we adopt a free form deformation (FFD) approach in order to deform the hull of our test case.
A very detailed description of FFD is beyond the scope of the present work. In the following we will give only a brief overview. For a further insight see [44] for the original formulation and [28, 40, 42, 19, 43] for more recent works.
We decided to adopt free form deformation among other possibilities (including, for instance, Radial Basis Functions or Inverse Distance Weighting) because it allows to have big and global deformations with a few parameters. For the complexity of the problem at hand, by trying to reduce the number of parameters starting from hundreds of them can be infeasible for the number of Monte Carlo simulations required. One of the possible drawbacks of FFD is generally that the parameters do not have a specific geometric meaning, like, for instance, a prescribed length or distance. In the case of application to AS this is not a problem since AS itself identifies new parameters, obtained by combination of the original ones, meaningless from the geometric and physical point of view.
FFD consists basically in three different step, as shown in Figure 3:

Mapping the physical domain to the reference one with the map .

Moving some control points to deform the lattice with . The movement of the control points is given by the weights of FFD, and represent our parameters .

Mapping back to the physical domain with the map .
So FFD map is the composition of the three maps, i.e.
(5) 
In this way, in the three dimensional case, every unperturbed point inside the FFD box changes its position according to
(6) 
where are Bernstein polynomials of degree , , in each direction, respectively. In Figure 4, we show the application of the FFD morphing on a very simple geometry, a sphere.
We implemented this algorithm in a python package called PyGeM [1] in order to deal with the major industrial CAD file formats. It handles iges, stl, step, vtk, unv, keyword, and openfoam files. It extracts the coordinates of the control points, deforms them according to the inputs given by the user and writes a new file with the morphed CAD. We improve the traditional version of the algorithm by allowing a rotation of the FFD lattice in order to give more flexibility to the tool. In general with our package it is possible to have a generic bounding box (not only a cube) as long as the map is affine.
We will see now how to integrate this parametrized geometrical setting with our problem of interest, used to introduce a complete computational embedded pipeline.
4 Three dimensional FluidStructure Interaction simulations based on rigid hull dynamics and potential flow theory
In this work, the free form deformation tool described in section 3 has been employed for the generation of a very large number of hull geometries obtained from the morphing of the DTMB5415 naval combatant hull. Each geometry generated has been used to set up a high fidelity hydrodynamic simulation also accounting for rigid hull motion. Ultimately, the output of such simulations is represented by the total resistance of the considered hull advancing in calm water at the specified speed. In section 5 we will provide a detailed report of the geometric parameters involved in the DTMB5415 shape modifications, along with the numerical experiments carried out in this work.
In the present section, we will instead provide a brief description of the unsteady fully nonlinear potential fluid dynamic model used to carry out the high fidelity simulations. In addition, we will describe the rigid body equations based on hull quaternions used to compute the hull linear and angular displacements. A few details on the spatial and time discretization, along with the fluid and structure problem coupling strategies will also be reported. We refer the interested reader to [33], [34] and [32] for further information on the fully nonlinear potential free surface model, on its application to complex hull geometries, and on the treatment of hull rigid motions respectively.
4.1 Fully nonlinear potential model
In this work, we are only considering the motion of a ship advancing at constant speed in calm water. For such reason, we solve the problem in a global, translating reference frame , which is moving with the constant horizontal velocity of the boat . Thus, the axis of the reference frame is aligned with , the axis is directed vertically (positive upwards), while the axis is directed laterally (positive port side).
Under the assumptions of irrotational flow, non viscous fluid, and simply connected domain , the velocity field admits a representation through a scalar potential function , namely
(7) 
in which is perturbation potential. Under the present assumptions, the equations of motion simplify to the unsteady Bernoulli equation and to the Laplace equation for the perturbation potential:
(8a)  
(8b) 
where is an arbitrary function of time, and , is the gravity acceleration vector, directed along the axis. The unknowns of such mathematical problem and are uncoupled. This means that the solutuon of Poisson problem (8b) can be obtained independently of the pressure field. Once such solution is obtained, the pressure can be obtained through a postprocessing step based on Bernoulli’s equation (8a). Thus, the Laplace equation is the governing equation of our model. Such equation is complemented by non penetration boundary conditions on the hull surface and water basin bottom boundary , and by homogeneous Neumann boundary conditions on the truncation boundaries of the numerical domain. On the water free surface , we employ the kinematic and dynamic semiLagrangian fully nonlinear boundary conditions, which respectively read
(9)  
(10) 
The former equation expresses the fact that a material point moving on the free surface will stay on the free surface —here assumed to be a single valued function of the horizontal coordinates and . The latter condition is derived from Bernoulli’s equation (8a), under the assumption of constant atmospheric pressure on the water surface. This peculiar form of the fully nonlinear boundary conditions was proposed by Beck [5]. Eq. 9 allows for the computation of the vertical velocity of markers which move on the water free surface with a prescribed horizontal speed . Eq. 10 is used to obtain the velocity potential values in correspondence with such markers. The resulting vector is the time derivative of the position of the free surface markers. In this work, such free surface markers are chosen as the free surface nodes of the computational grid. To avoid an undesirable mesh nodes drift along the water stream, the markers arbitrary horizontal velocity is set to 0 along the direction. The component of the water nodes in contact with the ship —which is moved according with the computed linear and angular displacements— is chosen so as to keep such nodes on the hull surface. As for the remaining water nodes, the lateral velocity value is set to preserve mesh quality.
4.2 Three dimensional hull rigid motions
The ship hull is assumed to be a rigid body. A second, hullattached reference frame , which follows the hull in its translations and rotations is employed to study the ship motions. The center of such reference frame is located at the ship center of gravity, which in the global reference frame reads , where , , are the unit vectors along the global system axes.
The rotation matrix is used to convert the coordinates of a point written in the hullattached reference frame, to those in the global frame , namely
(11) 
The global frame velocity of a point having coordinates in the hullattached frame is obtained as
(12) 
in which is the angular velocity vector.
Eqs. 11 and 12 imply that once , , and are known at time , the position and velocity of each point of the hull can be obtained. For this reason, writing the time evolution equations for , , and is sufficient to fully determine the hull dynamics. In the next sections, we will present such evolution equations written in the global reference frame, as presented in [2] and [18].
4.2.1 Linear momentum conservation.
The evolution equation for is obtained via the linear momentum conservation equation, which in the case of our hydrodynamics simulation framework reads
(13) 
In Eq. 13, is the mass of the ship, while the hydrodynamic force vector is in principle obtained as the sum of the pressure and viscous forces on the hull, propeller and appendages.
4.2.2 Angular momentum conservation.
The evolution equation for is obtained writing the angular momentum conservation, namely
(14) 
where is the matrix of inertia of the ship in the hullattached reference frame, and hydrodynamic moment vector is the sum of the moment about the ship center of gravity of the pressure and viscous forces on hull, propeller and appendages.
4.2.3 Rotation matrix and hull quaternions.
To write an evolution equation for , we first the introduce the angular velocity tensor associated to , namely
(15) 
Note that tensor will act on a vector as if the operator were applied to :
(16) 
Making use of such tensor, an evolution equation for the rotation matrix reads
(17) 
which can be advanced in time to obtain the components of and close the equations of motions of a rigid body in three dimension. Yet, in common the practice of rigid body simulations, direct numerical integration of Eq. 17 is avoided. The most important reason for this is related to numerical drift. If we in fact keep track of the orientation of a rigid body integrating Eq. 17, numerical error will build up in the entries of , so that it will no longer be a rotation matrix, losing its properties of orthogonality and of having determinant equal to 1. Physically, the effect would be that applying to a body would cause a skewing effect.
A better way to represent the orientation of a rigid body in three dimensions (even with large rotations) is represented by the use of unit quaternions (see the work of [45] for details). For our purposes, quaternions can be considered as a particular type of four element vector, normalized to unit length. If we indicate the quaternion as , the internal product of two quaternions and is defined as
(18) 
The norm of a quaternion is defined as . Unit quaternions can be used to represent rotations in a three dimensional space. In fact, given a quaternion , we can obtain the corresponding rotation matrix as
(19) 
Finally, the equation needed to describe the time evolution for the hull quaternion is
(20) 
where is the quaternion associated with the angular velocity vector . As quaternions only have four entries, there only is one extra variable used to describe the three degrees freedoms of a three dimensional rotation. A rotation matrix instead employs nine parameters for the same three degrees of freedom; thus, the quaternions present far less redundancy than rotation matrices. Consequently, quaternions experience far less numerical drift than rotation matrices. The only possible source of drift in a quaternion occurs when the quaternion has lost its unit magnitude. This can be easily corrected by periodically renormalizing the quaternion to unit length.
4.3 Discretization and numerical solution
The boundary value problem described is governed by the linear Laplace operator. Yet, it is nonlinear due to the presence of the boundary conditions in equations 9 and 10. Further sources of nonlinearity are given by continuous change of the domain shape over time and by the arbitrary shape of the ship hull. Thus, for each time instant, we will look for the correct values of the unknown potential and node displacement fields by solving a specific nonlinear problem resulting from the spatial and time discretization of the original boundary value problem. The spatial discretization of the Laplace problem is based upon a boundary integral formulation described in [20]. In this framework, the domain boundary is subdivided into quadrilateral cells, on which bilinear shape functions are used to approximate the surface, the flow potential values, and the normal component of its surface gradient. The resulting Boundary Element Method (BEM, see [7]) consists in collocating a Boundary Intergal Equation (BIE) in correpondence with each node of the numerical grid, and computing the integrals appearing in such equation by means of the described isoparametric formulation. The linear algebraic equations obtained from such discretization method are combined with the Ordinary Differential Equations (ODE) derived from the Finite Element spatial discretization of the fully nonlinear free surface boundary conditions 9 and 10. The spatial discretization described is carried out making use of the deal.II open source library for C++ implementation of finite element discretizations ([4, 3]). To enforce a strong coupling between the fluid and structural problem, the aforementioned system of Differential Algrbraic Equations (DAE) is complemented by the equations of the rigid hull dynamics (Eqs.13, 14 and 20). The DAE system solution is time integrated by means of an arbitrary order and arbitrary times step implicit Backward Difference Formula (BDF) scheme implemented in the IDA package of the SUNDIALS open source C++ library SUNDIALS ([22]). The potential flow model described has been implemented in a stand alone C++ software, the main features of which are described in [33].
The solver is complemented with a mesh module directly interfaced with CAD data structures [34]. Such feature allows for fully automated mesh generation once a hull shape is assigned at the start of each simulation by means of a — possibly non watertight — CAD geometry. Figure 5, on the left, displays the mesh generated on the surface of a DTMB5415 Navy Combatant hull. At each time step of the simulation, the wave resistance is computed as
(21) 
where the is the pressure value obtained introducing the computed potential in equation (8a). The inviscid fluid dynamic model drag prediction is then corrected by adding a viscous drag contribution obtained by means of the ITTC57 formula. Figure 5, on the right, depicts a comparison between the computed total resistance curve and the corresponding one measured by Olivieri et al. [35]. As can be appreciated in the plot, for all the Froude numbers tested the computed total drag difference with respect to the measurements is less then 6%.
5 Numerical results
In this section we are going to introduce the particular naval problem and how we merged together the free form deformation tool, the high fidelity solver, the active subspaces properties approach and the creation of the response surface.
For a simpler test case with only geometrical parameters see section 7. There the reader can clear any doubts and questions about the pipeline proposed seeing a much familiar problem of interest involving only geometrical deformation.
5.1 The hull
Given the relatively accurate numerical results obtained for the DTMB5415 benchmark hull, we selected such geometry as the reference hull for our simulations campaign. The DTMB5415 was originally conceived for preliminary design of a US Navy Combatant in 1980. The hull geometry includes both a sonar dome and a transom stern (see [35]), as can be appreciated in Figure 6 which represents a side view of the hull.
5.2 The parameters
As geometrical parameters we select 6 components of 4 different control points of a FFD lattice over one side wall of the hull. Then we apply the same deformation to the other side. This because one of our hypothesis is the symmetry of the deformed hull. In particular table 1 summarizes the set of design variables, the associated FFDnode coordinate modified ( is the span of the hull, its length and its depth) and the lower and upper bound of the modification. The structural parameter is the displacement of the hull and the physical one is the velocity. We underline that the original shape coincides with the domain center. From now on we denote with the column vector of the parameters, where are defined in table 1. Figure 7 shows the FFD lattice and the semihull.
Parameter  Nature  Lower bound  Upper bound 

FFD Point 1  0.2  0.3  
FFD Point 2  0.2  0.3  
FFD Point 3  0.2  0.3  
FFD Point 4  0.2  0.3  
FFD Point 3  0.2  0.5  
FFD Point 4  0.2  0.5  
weight (kg)  500  800  
velocity (m/s)  1.87  2.70 
To create a dataset with 130 simulations we take the original iges file of the hull and deforme it with an in house developed python package called PyGeM (Python Geometrical Morphing, see [1]) distributed via GitHub. The deformations are generated randomly with an uniform distribution. All the inputs are bounded from above and from below in order to satisfy physical constraints.
5.3 The solver
Then the high fidelity solver reads the iges and simulates the behaviour of the hull for 30 seconds. As solver we use the inhouse developed Boundary Element Method solver WaveBEM (as described in section 4). To further speedup the computations, the total resistance computed as in equation 21 is extrapolated to obtain the total resistance at the final, steady state regime, with an error in the order of 0.1%. In Figure 8 we can see the original simulation of the total resistance for the first 30 seconds and then the extrapolation we have done after a proper filtering of the data. We fit the maximums using the following function: . For the minimums we use . Then we set the approximated wave resistance to the average of the two at infinity.
At this point we have the scalar output function computed for all the 130 samples in the parameter space. We split the dataset in two, creating a train dataset with 80% of the data, and a test dataset with the remaining 20%.
5.4 Eigenvalues analysis
In order to construct the uncentered covariance matrix , defined in Eq. (1), we use a Monte Carlo method as shown in Eq. (2), employing the software in [10]. Since we have only pairs of input/output data we need to approximate the gradients of the total wave resistance with respect to the parameters, that is . To do so we use a local linear model that approximates the gradients with the best linear approximation using 14 nearest neighbors. After constructing the matrix we calculate its real eigenvalue decomposition. Recalling section 2 here the number of parameters is equal to 8, so .
In Figure (a)a we see the eigenvalues of the matrix and the bootstrap intervals. Bootstrapping is the practice of estimating properties of a quantity (such as its variance) by measuring those properties when sampling from an approximating distribution. It relies on random sampling with replacement. The bootstrap intervals in grey are computed using 1000 bootstrap replicates randomly extracted from the original gradient samples. We underline a major gap between the first and the second eigenvalue and a minor one between the second and the third. This suggests the existence of a one dimensional subspace and possibly the presence of a two dimensional one. To better investigate the first case, in Figure (b)b we present the components of the eigenvector with index 1 that correspondes to the greatest eigenvalue, that is the matrix — in this case it is a vector — of Eq. (3). Since they are the weights of the linear combination used to construct the active direction they provide useful informations. We can see that the major contributions are due to the velocity, the weight, and the depth of the hull. We underline that a weight that is almost zero means that the output function is almost flat along the direction identified by the corresponding parameter. This is a very useful information for a designer because in such a way he can deform the shape along prescribed directions without affecting the total wave resistance, allowing for example to store more goods inside the hull preserving the performances.
In Figure 10 we explore the training dataset, plotting the sufficient summary plot (see [14]) for one and two active variables. Sufficient summary plots are powerful visualization tools to identifying lowdimensional structure in a quantity that depends on several input variables. A scatter plot that contains all available regression information is called a sufficient summary plot. Recalling Eq. (4), Figure 10 shows against , where contains the first one and two eigenvectors respectively. In particular each point represents the value of the output function for a particular choice of the parameters, mapped in the active subspace. The two plots confirm the presence of an active subspace of dimension one and two. The latter seems to capture the output function in a much finer way, but as we are going to show the gain in terms of error committed is not so big.
5.5 Error analysis
We can compare the decay of the eigenvalues with the decay of surrogate model error on the test dataset shown in Figure (a)a. We project the data onto active subspaces of varying dimension, and construct a surrogate model with a leastsquaresfit, global, multivariate polynomial approximation of different orders. Then we calculate the rootmeansquare error of the test data against the surrogate. This error is scaled with respect to the range of the function evaluations, making it a relative error. We repeat this procedure 20 times constructing every time the uncentered covariance matrix of Eq. (2), since a Monte Carlo approximation is involved. Finally we take the average of the errors computed. Because we have a large amount of training data, we can expect the surrogate model constructed in a low dimension to be accurate if the data collapses into a manifold. Thus the test error is an indication of how well the active subspace has collapsed the data. In Figure (a)a are depicted the errors with respect to the active subspace dimension and the order of the response surface. The subspace dimension varies from 1 to 3, while the order of the response surface from 1 to 4. The minimum error is achieved using a two dimensional active variable and a response surface of degree 4 and it is . Further investigations show that increasing the dimension of the active variable does not decrease significantly the error committed while the time to construct the corresponding response surface increases. This is confirmed by the marginal gains in the decay of the eigenvalues for active variables of dimension greater then three as shown in Figure (a)a. We can affirm that the active subspace of dimension one is sufficient to model the wave resistance of the DTMB 5415 if we can afford an error of approximately 4.5%. In particular in Figure (b)b we can see the predictions made with the surrogate model of dimension one and the actual observations. Otherwise, we can achieve a error if we take advantage of two active dimensions and a response surface of order four, preserving a fast evaluation of the surrogate model.
We want to stress the fact that the result is remarkable if we consider the heterogeneous nature of all the parameters involved. In the case of only geometrical parameters one can easily expect such a behaviour but considering also physical and structural ones make the result not straightforward at all. Moreover the evaluation of the response surface takes less than one second compared to the 12 hours of a full simulation per single set of parameters on the same computing machine. This opens new potential approaches to optimization problems.
6 Conclusions
In this work we presented a numerical framework for the reduction of the parameters space and the construction of an optimized response surface to calculate the total wave resistance of the DTMB 5415 advancing in calm water. We integrate heterogeneous parameters in order to have insights on the more important parameters. The reduction both in terms of cost and time, remaining below the 4.2% of error on new unprocessed data, is very remarkable and promising as a new design interpreted tool. The methodological and computational pipeline is also easily extensible to different hulls and/or different parameters. This allows a fast preprocessing and a very good starting point to minimize the quantities of interest in the field of optimal shape design.
This work is directed in the development of reduced order methods (ROMs) and efficient parametric studies. Among others we would like to cite [8, 21, 41, 43] for a comprehensive overview on ROM and geometrical deformation. Future developments involve the application of the POD after the reduction of the parameters space through the active subspaces approach.
Acknowledgement
This work was partially supported by the project OpenViewSHIP, “Sviluppo di un ecosistema computazionale per la progettazione idrodinamica del sistema elicacarena” and “Underwater Blue Efficiency”, supported by Regione FVG  PAR FSC 20072013, Fondo per lo Sviluppo e la Coesione, by the project “TRIM â Tecnologia e Ricerca Industriale per la Mobilità Marina”, CTN0100176163601, supported by MIUR, the italian Ministry of Education, University and Research, by the INDAMGNCS 2016 projects “Tecniche di riduzione della complessità computazionale per le scienze applicate” and “Numerical methods for model order reduction of PDEs” and by European Union Funding for Research and Innovation — Horizon 2020 Program — in the framework of European Research Council Executive Agency: H2020 ERC CoG 2015 AROMACFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” P.I. Gianluigi Rozza.
7 Appendix: the geometrical test case
Here we are going to present the pipeline proposed using only geometrical parameters. The hull, the solver and the free form lattice are the same presented in section 5. The output function is slightly different since it accounts also the immersed volume of the hull. The immersed volume is obtained using the vertical component of the hydrostatic force acting on a submerged surface that reads
The output function is the total wave resistance divided by the immersed volume .
As parameter inputs we select 8 components of 4 different control points of a FFD lattice over one side wall of the hull as shown in table 2. Then we apply the same deformation to the other side. The displacement is set to 550 kg and the Froude number is fixed at 0.28.
The train dataset is composed by the 80% of the data with the other 20% is used as test to perform error analysis. To calcolate the gradients of the output function with respect to the parameters we use a local linear model that approximates them with the best linear approximation using 19 nearest neighbors.
Parameter  Nature  Lower bound  Upper bound 

FFD Point 2  0.2  0.3  
FFD Point 4  0.2  0.3  
FFD Point 1  0.2  0.3  
FFD Point 2  0.2  0.3  
FFD Point 3  0.2  0.3  
FFD Point 4  0.2  0.3  
FFD Point 3  0.2  0.5  
FFD Point 4  0.2  0.5 
In Figure 12 we can see the eigenvalues estimates and the eigenvector corresponding to the major eigenvalue. For the details we refer to section 5.4. The main eigenvalues gap is between the first and the second eigenvalue while a minor one is present between the second and the third one. We underline that the first eigenvector suggests that the fourth and the sixth parameter variation does not contribute significantly to the output function.
The sufficient summary plots corresponding to the active subspace of dimension 1 and 2 are depicted in figure 13. The dimension of the active variable is dictated by the gaps between the eigenvalues.
In Figure 14 are depicted the errors with respect to the active subspace dimension and the order of the response surface. The subspace dimension varies from 1 to 3, while the order of the response surface from 1 to 4. The minimum error, , is achieved using a two dimensional active variable and a response surface of degree 3. For more details on how the error is computed see section 5.5.
References
 [1] PyGeM: Python Geometrical Morphing. Available at: https://github.com/mathLab/PyGeM.
 [2] R. Azcueta. Computation of turbulent freesurface flows around ships and floating bodies. Schriftenreihe Schiffbau, Bericht Nr. 612, June 2001. ISBN: 3892206120.
 [3] W. Bangerth, D. Davydov, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier, B. Turcksin, and D. Wells. The deal.II library, Version 8.4. Journal of Numerical Mathematics, 24(3):135–141, 2016.
 [4] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II – a general purpose object oriented finite element library. ACM Transactions on Mathematical Software, 33(4):24/1–24/27, 2007.
 [5] R. F. Beck. Timedomain computations for floating bodies. Applied Ocean Research, 16:267–282, 1994.
 [6] G. E. Box and N. R. Draper. Empirical modelbuilding and response surfaces, volume 424. Wiley, New York, 1987.
 [7] C. A. Brebbia. The Boundary Element Method for Engineers. Pentech Press, 1978.
 [8] F. Chinesta, A. Huerta, G. Rozza, and K. Willcox. Model order reduction. Encyclopedia of Computational Mechanics. Second Edition. John Wiley & Sons, Ltd., New York, 2017.
 [9] P. Constantine and D. Gleich. Computing active subspaces with Monte Carlo. arXiv preprint arXiv:1408.0545, 2015.
 [10] P. Constantine, R. Howard, A. Glaws, Z. Grey, P. Diaz, and L. Fletcher. Python activesubspaces utility library. The Journal of Open Source Software, 1(5), September 2016.
 [11] P. G. Constantine. Active subspaces: Emerging ideas for dimension reduction in parameter studies, volume 2. SIAM, Philadelphia, 2015.
 [12] P. G. Constantine, A. Eftekhari, and R. Ward. A nearstationary subspace for ridge approximation. arXiv preprint arXiv:1606.01929, 2016.
 [13] P. G. Constantine, M. Emory, J. Larsson, and G. Iaccarino. Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet. Journal of Computational Physics, 302:1–20, 2015.
 [14] R. D. Cook. Regression graphics: ideas for studying regressions through graphics, volume 482. John Wiley & Sons, New York, 2009.
 [15] J. Dambrine, M. Pierre, and G. Rousseaux. A theoretical and numerical determination of optimal ship forms based on michellÂs wave resistance. ESAIM: Control, Optimisation and Calculus of Variations, 22(1):88–111, 2016.
 [16] J. L. Devore. Probability and Statistics for Engineering and the Sciences. Cengage Learning, 2015.
 [17] M. Diez, E. F. Campana, and F. Stern. Designspace dimensionality reduction in shape optimization by Karhunen–Loève expansion. Computer Methods in Applied Mechanics and Engineering, 283:1525–1544, 2015.
 [18] L. Formaggia, A. Mola, N. Parolini, and M. Pischiutta. A threedimensional model for the dynamics and hydrodynamics of rowing boats. Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology, 224(1):51–61, 2010.
 [19] D. Forti and G. Rozza. Efficient geometrical parametrisation techniques of interfaces for reducedorder modelling: application to fluid–structure interaction coupling problems. International Journal of Computational Fluid Dynamics, 28(34):158–169, 2014.
 [20] N. Giuliani, A. Mola, L. Heltai, and L. Formaggia. FEM SUPG stabilisation of mixed isoparametric BEMs: application to linearised free surface flows. Engineering Analysis with Boundary Elements, 8–22:59, 2015.
 [21] J. S. Hesthaven, G. Rozza, and B. Stamm. Certified reduced basis methods for parametrized partial differential equations. SpringerBriefs in Mathematics. Springer, 2016.
 [22] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward. Sundials: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software, 31(3):363–396, 2005.
 [23] J. M. Hokanson and P. G. Constantine. Datadriven polynomial ridge approximation using variable projection. arXiv preprint arXiv:1702.05859, 2017.
 [24] J. L. Jefferson, J. M. Gilbert, P. G. Constantine, and R. M. Maxwell. Reprint of: Active subspaces for sensitivity analysis and dimension reduction of an integrated hydrologic model. Computers & Geosciences, 90:78–89, 2016.
 [25] J. Kaipio and E. Somersalo. Statistical and computational inverse problems, volume 160. Springer Science & Business Media, 2006.
 [26] S. Keiper. Analysis of generalized ridge functions in high dimensions. In Sampling Theory and Applications (SampTA), 2015 International Conference on, pages 259–263. IEEE, 2015.
 [27] G. Le Quilliec, B. Raghavan, and P. Breitkopf. A manifold learningbased reduced order model for springback shape characterization and optimization in sheet metal forming. Computer Methods in Applied Mechanics and Engineering, 285:621–638, 2015.
 [28] M. Lombardi, N. Parolini, A. Quarteroni, and G. Rozza. Numerical simulation of sailing boats: Dynamics, FSI, and shape optimization. In Variational Analysis and Aerospace Engineering: Mathematical Challenges for Aerospace Design, pages 339–377. Springer, 2012.
 [29] T. W. Lukaczyk, P. Constantine, F. Palacios, and J. J. Alonso. Active subspaces for shape optimization. In 10th AIAA Multidisciplinary Design Optimization Conference, page 1171, 2014.
 [30] L. Meng, P. Breitkopf, B. Raghavan, G. Mauvoisin, O. Bartier, and X. Hernot. Identification of material properties using indentation test and shape manifold learning approach. Computer Methods in Applied Mechanics and Engineering, 297:239–257, 2015.
 [31] N. Metropolis and S. Ulam. The monte carlo method. Journal of the American statistical association, 44(247):335–341, 1949.
 [32] A. Mola, L. Heltai, and A. De Simone. Wet and Dry Transom Stern Treatment for Unsteady and Nonlinear Potential Flow Model for Naval Hydrodynamics Simulations. Journal of Ship Research, 61(1):1–14, 2017.
 [33] A. Mola, L. Heltai, and A. DeSimone. A stable and adaptive semilagrangian potential model for unsteady and nonlinear shipwave interactions. Engineering Analysis with Boundary Elements, 128–143:37, 2013.
 [34] A. Mola, L. Heltai, A. DeSimone, et al. Ship Sinkage and Trim Predictions Based on a CAD Interfaced Fully Nonlinear Potential Model. In The 26th International Ocean and Polar Engineering Conference, volume 3, pages 511–518. International Society of Offshore and Polar Engineers, 2016.
 [35] A. Olivieri, F. Pistani, A. Avanzini, F. Stern, and R. Penna. Towing tank experiments of resistance, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant insean 2340 model. Technical report, DTIC Document, 2001.
 [36] A. Pinkus. Ridge functions, volume 205. Cambridge University Press, 2015.
 [37] B. Raghavan, P. Breitkopf, Y. Tourbier, and P. Villon. Towards a space reduction approach for efficient structural shape optimization. Structural and Multidisciplinary Optimization, 48(5):987–1000, 2013.
 [38] B. Raghavan, G. Le Quilliec, P. Breitkopf, A. Rassineux, J.M. Roelandt, and P. Villon. Numerical assessment of springback for the deep drawing process by level set interpolation using shape manifolds. International Journal of Material Forming, 7(4):487–501, 2014.
 [39] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323–2326, 2000.
 [40] G. Rozza, A. Koshakji, and A. Quarteroni. Free Form Deformation techniques applied to 3D shape optimization problems. Communications in Applied and Industrial Mathematics, 4, 2013.
 [41] F. Salmoiraghi, F. Ballarin, G. Corsi, A. Mola, M. Tezzele, and G. Rozza. Advances in geometrical parametrization and reduced order models and methods for computational fluid dynamics problems in applied sciences and engineering: overview and perspectives. In M. Papadrakakis, V. Papadopoulos, G. Stefanou, and V. Plevris, editors, Proceedings of the VII European Congress on Computational Methods in Applied Sciences and Engineering, Crete, Greece, volume 1, pages 1013–1031, 510 June 2016.
 [42] F. Salmoiraghi, F. Ballarin, L. Heltai, and G. Rozza. Isogeometric analysisbased reduced order modelling for incompressible linear viscous flows in parametrized shapes. Advanced Modeling and Simulation in Engineering Sciences, 3(1):21, 2016.
 [43] F. Salmoiraghi, A. Scardigli, H. Telib, and G. Rozza. Free form deformation, mesh morphing and reduced order methods: enablers for efficient aerodynamic shape optimization. Submitted, 2017.
 [44] T. W. Sederberg and S. R. Parry. Freeform deformation of solid geometric models. In ACM SIGGRAPH computer graphics, volume 20, pages 151–160. ACM, 1986.
 [45] K. Shoemake. Animating rotation with quaternion curves. In ACM Computer Graphics (Proc. SIGGRAPH), pages 245–254, 1985.
 [46] Y. Tahara, H. Kobayashi, M. Kandasamy, W. He, D. Peri, M. Diez, E. Campana, and F. Stern. CFDbased multiobjective stochastic optimization of a waterjet propelled high speed ship. In 29th Symposium on Naval Hydrodynamics, Gothenburg, Sweden, 2012.
 [47] J. B. Tenenbaum, V. De Silva, and J. C. Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319–2323, 2000.
 [48] S. Volpi, M. Diez, N. J. Gaul, H. Song, U. Iemma, K. Choi, E. F. Campana, and F. Stern. Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification. Structural and Multidisciplinary Optimization, 51(2):347–368, 2014.