Dimension reduction in heterogeneous parametric spaces with application to naval engineering shape design problems

Dimension reduction in heterogeneous parametric spaces with application to naval engineering shape design problems

Marco Tezzele, Filippo Salmoiraghi, Andrea Mola and
Gianluigi Rozza

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.

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

Figure 1: Scheme of the structure of the pipeline proposed. The geometrical deformation is performed via free from deformation (FFD), then a BEM solver computes the wave resistance, the active subspaces (AS) are detected and finally a response surface method (RSM) is employed.

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

(a) Original output function
(b) First rotation
(c) Second rotation
(d) Output function with respect to the active variable
Figure 2: Example of a bivariate output function LABEL:sub@subfig:output_funct, intermediate rotations of the domain LABEL:sub@subfig:output_funct_rot1 and LABEL:sub@subfig:output_funct_rot2, and the final state LABEL:sub@subfig:output_funct_active, where we can see the variation of the function along the active variable.

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


where is the expected value, a probability density function — usually a uniform one —, and the so-called 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


where we draw samples independently from the measure and where . This matrix is symmetric and positive semidefinite, so it admits a real eigenvalue decomposition


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 reduced-order 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 full-space 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:


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.


In this way, in the three dimensional case, every unperturbed point inside the FFD box changes its position according to


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.

Figure 3: Sketch of the FFD map construction.
Figure 4: FFD morphing on a simple geometry: the sphere case. Here we move only one FFD control point in the lattice.

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 Fluid-Structure 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 DTMB-5415 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 DTMB-5415 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


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:


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 semi-Lagrangian fully nonlinear boundary conditions, which respectively read


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, hull-attached 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 hull-attached reference frame, to those in the global frame , namely


The global frame velocity of a point having coordinates in the hull-attached frame is obtained as


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


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


where is the matrix of inertia of the ship in the hull-attached 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


Note that tensor will act on a vector as if the operator were applied to :


Making use of such tensor, an evolution equation for the rotation matrix reads


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


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


Finally, the equation needed to describe the time evolution for the hull quaternion is


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 bi-linear 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 iso-parametric 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].

Figure 5: On the left, the mesh automatically generated on the surface of the DTMB-5415 Navy Combatant Hull. On the right, the total resistance of the DTMB-5415 hull as a function of the Froude number associated with the surge velocity imposed in the simulations. The blue continuous line represents the experimental values presented in Olivieri et al. [35]. The values obtained in this work are represented by the dashed magenta line.

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 water-tight — CAD geometry. Figure 5, on the left, displays the mesh generated on the surface of a DTMB-5415 Navy Combatant hull. At each time step of the simulation, the wave resistance is computed as


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 ITTC-57 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 DTMB-5415 benchmark hull, we selected such geometry as the reference hull for our simulations campaign. The DTMB-5415 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.

Figure 6: Original DTMB 5415 semi-hull on which we performed the FFD.

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 FFD-node 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 semi-hull.

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
Table 1: Design space (FFD lattice ) with 8 design parameters. 6 geometrical parameters chosen among the FFD control points, 1 structural parameter that is the initial displacement of the hull and finally 1 physical parameter given by the velocity.

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.

(a) Front view
(b) Back view
Figure 7: Plot LABEL:sub@subfig:ffd_hull_1 shows the FFD lattice over one side wall of the hull from the front, while plot LABEL:sub@subfig:ffd_hull_2 depicts the hull and the lattice from the back together with the numbers that identify the FFD points.

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

(a) Original wave resistance
(b) Fitted wave resistance
Figure 8: Plot LABEL:sub@subfig:wave_res shows the original wave resistance simulated for 30 seconds. Then plot LABEL:sub@subfig:wave_res_fitted depicts, after a filter has been applied, the exponential fitting of the maximums and minimums and the average at regime for 60 seconds.

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

Now the data are ready to be analyzed exploiting the active subspaces properties [11] we have seen in section 2.

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.

(a) Eigenvalues estimates.
(b) First eigenvector estimate.
Figure 9: Plot LABEL:sub@subfig:eigenvalues_hull shows the eigenvalue estimates in block circles with the bootstrap intervals (grey region). The order-of-magnitude gaps between the eigenvalues suggest confidence in the dominance of the active subspace. Plot LABEL:sub@subfig:eigenvector_hull shows the components of the eigenvector correspondent to the greatest eigenvalue.

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

(a) One active variable.
(b) Two active variables.
Figure 10: Sufficient summary plots for LABEL:sub@subfig:ssp1_hull one and LABEL:sub@subfig:ssp2_hull two active variables using the training dataset.

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 least-squares-fit, global, multivariate polynomial approximation of different orders. Then we calculate the root-mean-square 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.

(a) Surrogate model error with respect to the active subspace dimension and the order of the response surface.
(b) Observations and the corresponding predictions using a polynomial response surface of order three and dimension one.
Figure 11: Plot LABEL:sub@subfig:heat_map shows the relative error on the test dataset with respect to the active subspace dimension and the order of the response surface; plot LABEL:sub@subfig:predictions_deg3_dim1 shows the observations of the test dataset and the corresponding predictions using a polynomial response surface of order three.

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.


This work was partially supported by the project OpenViewSHIP, “Sviluppo di un ecosistema computazionale per la progettazione idrodinamica del sistema elica-carena” and “Underwater Blue Efficiency”, supported by Regione FVG - PAR FSC 2007-2013, Fondo per lo Sviluppo e la Coesione, by the project “TRIM – Tecnologia e Ricerca Industriale per la Mobilità Marina”, CTN01-00176-163601, supported by MIUR, the italian Ministry of Education, University and Research, by the INDAM-GNCS 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 AROMA-CFD 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
Table 2: Design space (FFD lattice ) with 8 design parameters. is the span of the hull, its length and its depth. Refer to Figure (b)b for the exact collocation of the control points.

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.

(a) Eigenvalues estimates.
(b) First eigenvector estimate.
Figure 12: Plot LABEL:sub@subfig:eigenvalues_hull shows the eigenvalue estimates in block circles with the bootstrap intervals (grey region). The order-of-magnitude gaps between the eigenvalues suggest confidence in the dominance of the active subspace. Plot LABEL:sub@subfig:eigenvector_hull shows the components of the eigenvector correspondent to the greatest eigenvalue.

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.

(a) One active variable.
(b) Two active variables.
Figure 13: Sufficient summary plots for LABEL:sub@subfig:ssp1_hull_geo one and LABEL:sub@subfig:ssp2_hull_geo two active variables using the training dataset.

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.

Figure 14: Surrogate model error with respect to the active subspace dimension and the order of the response surface.


  • [1] PyGeM: Python Geometrical Morphing. Available at: https://github.com/mathLab/PyGeM.
  • [2] R. Azcueta. Computation of turbulent free-surface flows around ships and floating bodies. Schriftenreihe Schiffbau, Bericht Nr. 612, June 2001. ISBN: 3-89220-612-0.
  • [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. Time-domain computations for floating bodies. Applied Ocean Research, 16:267–282, 1994.
  • [6] G. E. Box and N. R. Draper. Empirical model-building 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 active-subspaces 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 near-stationary 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. Design-space 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 three-dimensional 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 reduced-order modelling: application to fluid–structure interaction coupling problems. International Journal of Computational Fluid Dynamics, 28(3-4):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. Data-driven 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 learning-based 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 semi-lagrangian potential model for unsteady and nonlinear ship-wave 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, 5-10 June 2016.
  • [42] F. Salmoiraghi, F. Ballarin, L. Heltai, and G. Rozza. Isogeometric analysis-based 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. Free-form 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. CFD-based 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description