Considerable effort has been devoted in recent years to develop robust and efficient computational strategies for the simulation of physical phenomena, that take into account shape uncertainty. Often, the problem is formulated as an elliptic partial differential equation (PDE) whose domain boundaries are uncertain. Such problems arise due to imperfections in manufacturing processes, e.g., in nano-optics where the production of nano particles is, often, inaccurate relatively to nano-scale electromagnetic wave lengths bejan2000shape . Other examples arise in the context of inverse problems, such as tomography where the visual representation of some hidden object is constructed by partial, and possibly noisy, measurements sethian1999level .
The common practice for quantifying uncertainty in computational models, is to employ the parametric method. The theoretical basis for the method was laid down by Wiener gPC-Wiener-HomoChaos . The method itself was initially developed by Ghanem and Spanos gPC-GS-StochasticFE , and later generalized by Xiu and Karniadakis gPC-WXK-FlowSimulations ; gPC-XK-SteadyStateDiffusion ; gPC-XK-Wiener-Askey
. In this approach the uncertain parameters are replaced by random quantities, and the problem is recast as a system with random input. The solution is estimated via global expansion of the random variables into a basis of uncorrelated functions. Thus, the stochastic problem is transformed into a deterministic system in higher dimension. The most popular expansions employed are theKarhunen-Loève expansion, and the generalized Polynomial Chaos (gPC) expansion.
The parametric approach, often, demonstrates superior performance in terms of computational effort over other traditional methods, see gPC-Xiu-Book for a detailed review. However, when the physical domain of the problem is uncertain the quantification by the parametric approach becomes much more challenging. The main difficulty stems from the fact that the problem is not characterized by smooth coefficients whose dependence on the random parameters is known. Thus, an accurate discretization which captures desired features of the solution for any realization of the random shape is not readily available.
The stochastic collocation method and Monte Carlo sampling, which rely on samplings of the random parameters and the solution of each realization deterministically are well established. However, for random domain problems each realization is, essentially, characterized by a different geometry and requires a custom discretization scheme. Generally, when the variations of the random domain are large and undergo complicated changes as a function of the random parameters, these methods become extremely expensive to employ with prohibitive computational costs.
To overcome the difficulties associated with the quantification of a random shape or domain, various techniques have been proposed. Typically, these are classified as one of the following: perturbation, fictitious domain, level-set or random domain mapping. Perturbation techniquesShape-Perturbation are straightforward and simple to apply, however, their applicability is limited to small shape deformations. The fictitious domain Shape-Fictitious and level-set methods Shape-LevelSet-Nouy20084663 ; Shape-LevelSet-XFEM are based on embedding the random domain in a larger, deterministic domain containing all possible realizations. These methods are capable of handling very irregular non-smooth geometries. However, the embedding introduces non-smoothness in the spatial region, that intersects with the random boundary. Thus, high-order convergence is only partially ensured in the entire computational domain.
The random domain mapping method Shape-TX2006b ; Shape-TX2006a is the most common tool used for solving PDEs on uncertain domains. The method is based on a realization-dependent coordinate transformation uniformly mapping all the realizations of the domain to a fixed, reference configuration. The variational formulation of the PDE on the random domain can then be posed on the reference domain, reducing the problem to a PDE on a fixed domain with stochastic coefficients. The transformed PDE whose domain is fixed is solved using standard techniques. However, the method is highly sensitive to the non-linear dependence of the problem on the random boundary. In case of complex evolution of the shape, the random coefficients are difficult to obtain and typically exhibit highly varying behavior. The common practice to overcome this difficulty is to impose a highly accurate discretization grid, often combined with dimensionality reduction techniques, e.g., sparse grids, to ensure reasonable computational effort. See Shape-DomainMapping-Regularity ; Shape-Harbrecht2016 ; Shape-Hptmr2015 for further details.
In this work an alternative method that attempts to mitigate the difficulties associated with the more standard techniques for PDEs on uncertain domains is proposed. The method is of boundary integral type BIE-CK-InverseScattering-Book and, thus, can handle large shape deformations. The analysis is based on two observations. First, that as a function of the random boundary of the domain the solution is piece-wise smooth in the spatial domain. Second, that in practice we seek to approximate the outcome of a predetermined set of linear output functionals operating on the random boundary. The key idea is of this work is to construct an integration grid whose associated weight function encompasses the irregularities and non-smoothness imposed by the random boundary. Thus, the outcome of the functionals can be evaluated accurately with relatively low number of integration gridpoints. This idea is similar to certain classic numerical techniques for estimating integrals of highly oscillatory functions, which rely on oscillatory weighted Gaussian integration formulae.
The proposed method constructs a discretization grid of the random surface for all possible realizations in two stages. In the first stage a spatial low-dimensional embedding of the family of random surfaces is constructed via the Coarea formula Federer . The embedding, essentially, captures any irregular behavior of the random surface and a discretization is applied only on a compact region in the spatial domain. In the second stage a parametric grid corresponding to the low-dimensional spatial grid is imposed. A sparse or hierarchical parametric grid can be applied for dealing with high dimensionality, while the spatial grid effectively ensures that the bulk variation of output functionals defined on the boundary is captured. In general, the method allows the handling of non-trivial geometries without the loss of accuracy in the region intersecting with the random interface.
Since this is a first case study and for the ease of presentation, the discussion has been limited to time-harmonic wave scattering by star-shaped obstacles. In the analysis and numerical study a scattering object and low-dimensional parametric space are assumed. More complicated examples in higher spatial and parametric dimensions are discussed. However, in-depth study of this topic is deferred to future work. For its simplicity, acoustic fluid-structure interaction has been chosen as the physical application. In that case, the solution represents small oscillations of pressure in a compressible ideal fluid. The method and ideas presented in this work can also be applied to electrodynamics and elastodynamics.
This work employs the null-field approach NFM-Martin-Book ; NFM-Waterman-MatrixFormulation ; NFM-Waterman-T-Matrix ; NFM-Wriedt-Review , which in contrary to the better known boundary element method (BEM) BIE-sauter2010boundary and the Nyström method kress2013linear , does not involve singular integrals. Null-field methods are fast and much easier to implement compared to BEM and the Nyström method. Their applicability range is, however, more limited. The null-field reconstruction technique NFM-EM-Random-Shape-1D ; NFM-Optimal-Reconstruction is inherently stable, admits a-priori error evaluation, and facilitates the extraction of features of interest without prior estimation of the entire solution. The method enables us to perform analysis from a purely geometric point of view, which avoids the additional complications associated with integration of weakly singular kernels. Combining low-dimensional surface embedding with BEM and Nyström method can be foreseen in a future study.
The paper is organized as follows. The fundamentals of the null-field reconstruction method for the time-harmonic wave scattering problem is presented in Section 2. Section 3 reviews the procedure of optimal reconstruction from a numerical linear algebra point of view. Section 4 consists of the main theoretical results of this work and includes the formulation of the problem. In Section 5 the proposed method is applied to a class of randomly shaped polygonal cylinders, as a proof of concept that the suggested method can, indeed, handle complex non-smooth shapes. Summary of the results, conclusions, and suggestions for applying the method in more complicated scenarios are given in Section 6.
2 Null-Field Reconstruction for Time Harmonic Wave Scattering
In this section a brief review on the null-field reconstruction method for the time-harmonic wave scattering problem is given. The time-harmonic acoustic scattering problem is presented, followed by a review of the fundamental theory of null-field methods. The main idea of the null-field reconstruction technique for time-harmonic wave scattering is presented in the concluding subsection.
2.1 Acoustic Scattering by Impenetrable Obstacles
Let denote a bounded domain in representing an impenetrable obstacle with boundary . We denote by the closure of . Let be the unbounded exterior region occupied by a uniform medium. Let denote a general spatial point,
For an incident time-harmonic field ’illuminating’ the obstacle, the scattered field satisfies the following exterior boundary value problem:
where . Equation (1) is known as the Helmholtz equation, where is the Laplacian and is the wavenumber. Equation (2) specifies the boundary condition, depending on the physical problem: adopting the acoustic terminology, it is sound-soft for Dirichlet problems and sound-hard for Neumann problems. Here is the unit outward normal to and is the normal derivative of . The last condition (3), known as the Sommerfeld radiation condition, ensures that the scattered field propagates from the obstacle to infinity. The solution of the exterior scattering problem is unique. A solution to the Helmholtz equation is called a wavefunction. A wavefunction satisfying the Sommerfeld condition (3) is called an outgoing wavefunction.
2.2 Null-Field Theory Fundamentals
Null-field methods for the acoustic scattering problem (1) are based on Green’s second theorem
which holds for any bounded domain with a Lipschitz piecewise smooth boundary , where and are scalar fields, and and denote corresponding normal derivatives.
Let be an outgoing wavefunction and let denote the total field, . Assuming is analytic in , it can be shown by (4) that
Thus, for a sound-soft obstacle ( on )
while for a sound-hard obstacle ( on )
Using (5) or (6), an infinite set of equations can be produced from which or on are approximated. In practice, one chooses a finite subset of equations of the form of (5) and (6) which are employed to optimally reconstruct the scattered field without an explicit estimation of or on . The core idea of reconstruction by functionals is presented in the next subsection, while the numerical procedure for its practical implementation is covered in Section 3.
2.3 Reconstruction of Surface Functionals
Typically we are interested in estimating features of interest which are expressed by the unknown surface density, or on . Often such features are the outcomes of functionals in an appropriate Hilbert space. Indeed, let denote the complex conjugate of or on . Then for a general non-smooth surface , the surface density belongs to the complex Hilbert space whose inner-product and norm are defined by
where denotes the complex conjugate of and is the induced volume form on the surface. Recall that by Riesz representation theorem any bounded linear functional operating on surface densities, is of the form . We call such functionals surface functionals.
In this work we focus on the estimation of the scattering coefficients of the expansion of to cylinder harmonics in the case. These coefficients, denoted by , satisfy
where denotes the th-order Hankel function of the first kind. See BIE-CK-InverseScattering-Book for further details. The scattering coefficients are very useful features of the surface density, since they can easily express other important quantities such as the far-field pattern and the radar cross section NFM-Thesis .
Consider the sound-soft case (5). Using the Hilbert space notation, it follows that the scattering coefficients satisfy
where denotes the th-order Bessel function of the first kind, and denotes the complex conjugate of the surface density . In practice only satisfying
are required for an accurate description of , see NFM-Optimal-Reconstruction for more details. Similar expressions can be derived for the sound-hard case.
The core idea of the reconstruction procedure is to approximate the outcome of target functionals (7) and without producing an explicit approximation of the surface density on . This, generally, allows us to handle complex geometries as well as irregular or singular surface densities much more accurately . Explicitly, we approximate the each outcome (7) by a linear combination of the following form
where are predetermined sets of functionals, whose outputs are either known or can be calculated directly. We call such functionals the information functionals.
As shown in (5) the outcome of the information functionals are readily available if are outgoing wavefunctions whose singularities are located in . Hence, given the information
the outcome of the target functional can be approximate by
Obtaining the coefficients while ensuring measurable error bounds of the estimations of the scattering coefficients can be achieved by reconstruction kernel approximation which is the main topic of Section (3).
3 Optimal Reconstruction with A-priori Error Estimate
In this section we review the procedure for the recovery of target functionals by information functionals in general Hilbert space. Reconstruction problems often involve regularization parameters which govern the stability and accuracy of the procedure. The result of the optimization of the reconstruction with respect to the regularization parameters is referred to as optimal reconstruction. Optimal reconstruction can be traced back to the notion of optimal recovery OR-GW-OptimalApproximation ; OR-MR-OptimalEstimation . A more modern analysis from an inverse problem point of view can be found in Moment-Louis-FeatureReconstruction ; Moment-Louis-UnifiedApproach ; Moment-LM-MollifierMethod .
We begin with the definition of the reconstruction problem and the notion of reconstruction kernel. This is followed by a brief description of the numerical procedure including error analysis. The final part elaborates on proper numerical integration rules, that are needed for the error estimates. The method and error analysis presented here, as well as further technical details have been initially introduced in NFM-Optimal-Reconstruction .
3.1 The Reconstruction Problem and Reconstruction Kernels
Let be a complex Hilbert space, whose inner-product is denoted by . The reconstruction problem is to approximate a finite set of target functionals
by a given finite set of information functionals,
where the element is unknown.
Let denote the norm induced by in , and let denote a closed convex subset of . A linear combination whose coefficients satisfy the minimality condition
is called an optimal reconstruction kernel of the target functional by the information functionals over .
Clearly, (11) is a projection on a convex set. The key point which is addressed later, is how to determine the convex set . Note that almost no prior knowledge on the element is assumed.
We will show in the next subsection, that obtaining (11) vastly exceeds our needs. In practice, it is sufficient to obtain an approximation satisfying
with respect to some predetermined threshold, . In that case the linear combination is simply called a reconstruction kernel (i.e., not optimal).
3.2 The Discrete Reconstruction Procedure with Error Analysis
For the evaluation of the reconstruction kernel, we assume a finite dimensional discretization satisfying the following definition.
Let be a bounded subset of a Hilbert space whose inner-product is denoted by . A mapping,
is called an inner-product preserving discretization of of accuracy if
. The vectorsand are called the corresponding inner-product preserving discretizations of and on .
Let denote inner-product preserving discretizations of some , respectively. Let denote the orthogonal projection of on the subspace spanned by . By definition (2) we obtain
for all and . Hence, given the information (10) and an approximation of ,
we can reconstruct the unknown target coefficients (9) via
To evaluate the error of the reconstruction (15) we denote for each and the discretization errors
and obtain the following estimate
where our assumption ensure that . Note that is the projection error which can not be reduced if the set of information functionals, , is predetermined.
To control the error we impose the following regularization constraint
where is a chosen or given evaluation error bound. Thus, we obtain
Often, the summation of evaluation errors is not cumulative. Thus, the overall error is typically assuming . This facilitates an efficient approximation technique which is based performing successive singular value decompositions on subsets of information functionals. The technique was presented in NFM-Optimal-Reconstruction and demonstrated high stability and good convergence properties. Further details including different variants of the technique can be found in NFM-Thesis .
3.3 Inner-Product Preserving Discretization and Numerical Integration
Obtaining inner-product preserving discretizations of surface functionals is a fundamental issue. Let us focus on the case, as in this work, where and are smooth functions in some space with an inner-product,
where is compact and Jordan measurable, is the complex conjugate of and is a proper weight function.
To numerically compute the integrals 18, we observe that it is sufficient to employ an integration rule which is accurate on the finite dimensional subspace of smooth functions spanned by . Thus, we assume the availability of a standard rule of the following form
with integration nodes contained in and real positive weights . Discretizing an element as a weighted gridfunction
essentially, satisfies the inner-product preserving assumption (13) if is sufficiently large. The number of elements required for an effective inner-product preserving discretization depends on the convergence rate of the numerical integration formula and, typically, under some smoothness assumption of the integrands. Indeed, if the weight function encompasses all the singularities while and are analytic, a Gaussian numerical integration rule with respect to ensures exponential convergence. Note that the weights of Gaussian rules are always positive and uniformly bounded. See davis2007methods for more details.
4 Surface Embedding of Random Star-Shaped Obstacles
In this section the main theoretical contribution of this paper is presented. The first two subsections cover the setting of the problem, where Subsection (4.1) defines the random shape properties, and Subsection (4.2) covers relevant components of the generalized Polynomial Chaos (gPC) expansion theory. The chosen framework leads to a reconstruction problem in a Hilbert space. A concise discussion on the disadvantages of naive discretization of the reconstruction problem concludes Subsection (4.2).
In Subsection (4.3) we present an analytic approach for overcoming the difficulties associated with the naive discretization approach. Using the Coarea formula we construct a low-dimensional spatial embedding within the family of random surfaces, which facilitates a natural choice for setting a cubature rule in a compact region of
. The chosen integration weight function is a strictly positive minimal variance quantity encompassing the irregularities of the family of random surfaces.
In Subsections (4.4) and (4.5) we focus on the case of a single random variable describing the randomness of the object. Using the implicit function theorem we obtain explicit formulas including full characterization of the singular behavior of the integration weight function. The usage of the single random variable formulation as a building block for the more general case of multiple random variables is considered and discussed in Section (6).
Subsections (4.6) and (4.7) are devoted to the demonstration of the preceding theoretical parts on a model problem of a randomly oriented elliptic cylinder. The random orientation problem is a very simple ’toy’ problem. However, it allows us us to demonstrate in an affable fashion the implementation of the theory.
4.1 The Random Shape Setting
For brevity, we focus on the sound-soft case and assume that represents a star-shaped obstacle in whose boundary,
, depends smoothly on a real valued vector of mutually independent and continuous random variables
The boundary is, however, not assumed to be uniformly smooth in the spatial domain. We assume that each random variable
has finite even moments
where is the support of and is the probability density function of . Property (20) effectively ensures the existence of surface functionals suitable for the reconstruction of the scattering coefficients.
Our assumption that the obstacle is star-shaped for any realization of the random vector , ensures that its boundary possesses a polar representation,
and the existence of two positive radial bounds, and , satisfying
Thus, as illustrated in Figure 1, is confined to the transition region,
4.2 Random Shape and Generalized Polynomial Chaos Expansion
Given our assumptions we observe that the scattering coefficients (7) are finite dimensional random fields,
A common method to approximate these fields is to obtain their generalized Polynomial Chaos (gPC) expansions,
where is a multi-index and is an orthogonal basis of the inner-product space induced by the probability density function of ,
whose support is . Hence, the expansion coefficients are readily available by the orthogonality via
For a randomly shaped obstacle each coefficient in the gPC expansion (24) is a target functional of the following form
where the normalized metric tensor in polar coordinates is given by
In principle, we need to devise a discretization scheme for (26) and apply the optimal reconstruction procedure of Section (3). However, inherits any irregularity of family of surfaces; e.g., lack of smoothness and oscillatory behaviour, which often necessitates specialized high-order discretization of the surface . Additionally, discretizing the random surface integral with a grid of numerical integration nodes has to be realized for every grid point in the parameters domain, . Hence, in general, the practical implementation of an inner-product preserving discretization satisfying (13) is a difficult task. An analytic approach for overcoming this fundamental difficulty is presented in the next subsection.
4.3 Random Surface Embedding and the Coarea Formula
In this subsection we present an analytic approach for producing inner-product preserving discretizations of functionals of the form of (26). The key idea is to apply a change of variables transforming (26) to the following equivalent representation
where the weight function , is proportional to the conditional expectation of given the information . Thus, has minimal variance while, essentially, encompassing the irregularities of the family of random surfaces, . The term is a linear functional uniformly bounded in operating on . A high-order numerical integration rule with respect to would serve as a discretization satisfying (13). The transformed representation (27) is obtained by the so-called Coarea Formula Federer which allows us to express the surface integral in terms of the integral of the level sets of another function.
(The Coarea Formula)
Let be an open Jordan measurable subset of where and is a non-negative integer. Let be a piecewise smooth Lipschitz function, such that the level set,
is a piecewise smooth -dimensional manifold in . Then for any integrable function, , we have
where is the Jacobian of , and denotes surface measure of .
The coarea formula expresses the integral of a function over in terms of the level sets of the function . The level sets, , are called fibers of the domain . The formula is a kind of ”curvilinear” version of Fubini’s theorem.
Let us consider the function
By direct calculations we obtain
contains at least one smooth -dimensional manifold in . Explicitly, these points satisfy where is non-empty. Thus, we obtain the following representation,
where the domain of integration of the inner integral in (29) is given by
and the domain of integration of (30) is defined by the subset of irregular points of the set ,
Note that or (for certain values of ) can be empty sets, and in that case the associated integral is taken to be zero.
The advantage of the representation (29,30) is that the wave function, , in (29) is no longer composed with the boundary and does not inherit its irregular properties. The subset of irregular points, (32), defines portions of the random surface, (21), which are independet on ; i.e., non-random, thus reduces to an integral of the following general form,
whose discretization is straitforward.
A major challenge is to efficiently evaluate the inner integral in (29),
This integral can become infinite since and are, essentially, singular. Accordingly, using the following weight function
we have that
where denotes the irregular component (30) and
To show that is a bounded linear functional, let us assume for simplicity that (32) is an empty set. In that case, we we observe that (34) is, in fact, proportional to the conditional expectation of given ,
where is the probability density function of . It is well known that conditional expectation is a minimum variance predictor as a function of the given information. Hence, the choice (34) implies minimization of oscillatory behaviour of as a function of . Employing a similar argument we obtain that the linear functional (36) is, in fact,
Hence, the Cauchy-Schwarz inequality implies that
which shows that the linear functional, , is, indeed, bounded. A similar argument can be applied to show that the functional remains bounded when is not an empty set.
Setting the integration grid for (35) can be done in two stages. First we obtain a cubature rule with nodes and corresponding cubature weights with respect to the weight function , regardless of . In the second stage we identify the integration grid in which corresponds for the evaluation of .
4.4 Explicit Representation for a Single Random Variable
Let us consider a simplified case of one-dimensional random vector,
We will show that in this case the functional (36) can be explicitly represented. This approach can serve as a building block for the case of a general random vector, which is discussed in Section (6). The assumption does not imply a loss of generality, since the discretization of the irregular part (30) is carried out directly without applying the coarea formula.
Our assumptions imply that for any the fiber set (31) is either an empty set or composed of a finite set of discrete points; i.e., a zero-dimensional sub-surface. Thus, we obtain the following explicit representation of (29) for a single random variable,
where is the reduction of to a zero dimensional subset of ,