A stochastic domain decomposition method for time dependent mesh generation
Alexander Bihlo and Ronald D. Haynes 00footnotetext: Memorial University of Newfoundland, St. John’s, NL, Canada firstname.lastname@example.org Memorial University of Newfoundland, St. John’s, NL, Canada email@example.com
We are interested in PDE mesh generation where the mesh is computed by solving a mesh PDE which is coupled to the physical PDE of interest. In [bihl13a] we proposed a stochastic domain decomposition (DD) method to find adaptive meshes by solving a linear elliptic mesh generator. The stochastic DD approach, as originally formulated in [aceb05a], relies on an accurate numerical solution using the probabilistic form of the exact solution of the linear elliptic boundary value problem. Monte–Carlo simulations are used to evaluate this probabilistic form of the solution only at the sub-domain interfaces. These interface approximations can be computed independently and are then used as Dirichlet boundary conditions for the deterministic sub-domain solves. Within the framework of grid adaptation it is generally not necessary to solve the mesh PDEs with high accuracy. The reason is that the mesh equations are only a means to an end. Only a good quality mesh, one that allows an accurate representation of the physical PDE, is required. This relatively low accuracy requirement makes the proposed stochastic DD method computationally more attractive, reducing the number of Monte–Carlo simulations required.
Grid adaptation by a stochastic DD approach does generate interesting issues in its own right. Grid quality should be monitored during the interface solves to give a suitable stopping criteria for the stochastic portion of the algorithm. In [bihl13a] only the steady grid generation problem was considered. Of course, in practice, the problem of grid generation is coupled with the process of solving the system of physical, usually time dependent, PDEs. It is this latter issue that we begin to explore in this paper.
We are interested in time dependent PDEs whose solutions evolve on disparate space and time scales. The solution behaviour lends itself to the use of time dependent meshes which automatically adapt and evolve to efficiently resolve the solution features. The generation of these time dependent grids can be done either by statically applying an elliptic mesh generator using the physical solution obtained at the previous time step or by employing a time relaxation of the static mesh PDE resulting in a parabolic mesh equation, as in [huan10a].
In [bihl13a] we applied the stochastic DD algorithm to a linear elliptic mesh generator. Here we consider the extension to time dependent mesh PDEs. This extension to (linear) parabolic mesh generators is possible due to the existence of a stochastic representation of the exact solution of such linear parabolic problems. For the sake of illustration, we will work with the time-relaxed form of the Winslow–Crowley variable diffusion mesh generation method, first described in [wins66a].
2 Winslow’s method
The equipotential method of mesh generation in 2D, [crowley62], found the mesh lines in the physical co-ordinates and as the level curves of the potentials and satisfying Laplace’s equations
and appropriate boundary conditions which ensure grid lines lie along the boundary of the domain. Here derivatives are with respect to the physical co–ordinates. The physical mesh transformation, and , in the physical domain , can be found by (inverse) interpolation of the solution of (1) onto a (say) uniform grid. In practice, the inversion to the physical co–ordinates is not necessary. Instead one could transform the physical PDE of interest into the computational co–ordinate system.
Winslow [wins81a] generalized (1) by adding a diffusion coefficient depending on the gradient or other aspects of the solution. This gives the linear elliptic mesh generator
The function , known as a mesh density function, characterizes regions where additional mesh resolution is needed and in general depends on the solution of the physical PDE.
Here we assume the solution of the physical PDE is time dependent and hence the mesh density function is changing with time, . One could still use (2) to solve the mesh transformation at each time . For time dependent PDEs this would result a system of differential–algebraic equations for the physical solution and the mesh. Instead, we choose to relax (2) to obtain a parabolic linear mesh generator of the form
This gives a mesh PDE which depends explicitly on the mesh speed and provides a degree of temporal smoothing for the mesh, cf. [huan94a].
3 Linear parabolic differential equations and stochastic domain decomposition
The system of mesh PDEs (3) is of the form
where and are the computational coordinates defined over , where is the spatial domain in physical coordinates. In system (4), is a linear elliptic operator of the form
with continuous coefficient matrix , , and drift vector . Here we employ the summation convention over repeated indices.
System (4) is accompanied by boundary and initial conditions
The solution of such linear parabolic problems can be described using the tools of stochastic calculus [aceb10a, kara91a]. The point–wise solution of system (4) at is given probabilistically as
where the process satisfies, in the Îto sense, the stochastic differential equation (SDE)
The relation between and is given through
for all . The solution for is completely analogous.
In (4), the denotes the expected value, is the time when the stochastic path starting at first hits the boundary of the physical domain , is two-dimensional Brownian motion and is the indicator function.
where is the identity matrix.
For our two dimensional mesh generator we choose the initial conditions and , corresponding to an initial uniform mesh, and the static boundary conditions , , and . This ensures we use the standard computational domain and the rectangular physical domain . The remaining boundary conditions for and are determined by solving the 1D version of (2) along the boundaries. Collectively, we use and to denote these boundary conditions for and .
Hence we have to solve the SDE
|for the single path . The stochastic form of the exact solution of Eq. (3) for is then obtained by evaluating|
Thepoint–wise solution for is obtained in an analogous fashion.
In principle, the probabilistic solution (7) allows one to determine the computational coordinates and at each point in the space–time domain . However, this is prohibitively expensive. A more efficient approach is to evaluate the solution (7) only at certain points in space and time which then serve as boundary points for a DD implementation. This stochastic DD approach for parabolic problems has been studied by Acebrón et. al. [aceb10a].
In the mesh generation context it is not possible to obtain the solution of (5) at all times, as the solution of the mesh PDE is coupled to the physical solution. That is, rather than solving (5) for a time , it is generally only be possible to use this stochastic solution to advance the solution of (4) over one single time step from to . In this case, and should be interpreted as the values of and at time and the monitor function, , is given at either or and remains constant over the time step.
4 The numerical method
Stochastic solver and domain decomposition. The use of the stochastic solution (5) for the time-relaxed Winslow mesh generator with parameters (6) is straightforward. We solve (7a) using the classical Euler–Maruyama scheme, i.e. we employ linear time-stepping. An alternative would be to use exponential time-stepping as advocated e.g. in [aceb05a, bihl13a, jans03a]. In our results linear time-stepping gives sufficient accuracy. The components of the Brownian motion are computed as , where is a normally distributed random number with mean zero and variance one [kara91a].
The time dependent weight only becomes available only at each time step (due to a possible coupling with a physical PDE). Hence we are only able to employ formula (7b) to integrate over a single time step, i.e. from to . Over this time step, the weight function is evaluated at and held constant, i.e. we have in (7a). Accordingly, in Eq. (7b) is to be interpreted as , i.e. the values of the computational coordinates at the current time . Moreoever, the boundary functions and are updated at each time to reflect changes in the physical solution.
We then numerically integrate the SDE (7a) from to . The drift vector is required everywhere along the path of the stochastic process but is only available at the grid points of the domain. Bilinear interpolation is used to obtain the values of in between these grid points. The derivatives in are approximated using finite differences.
In the DD context, the stochastic solution is only required at a selection of points, , which live on the interfaces between sub-domains. One time step is split into several smaller sub-time steps in order to numerically integrate the SDE (7a) from to . We found this splitting of into sub-time steps necessary to determine sufficiently accurate whether the stochastic processes started at an interface point has left the domain during . This is not unlike the approach for mesh generation discussed in [huan10a]. At each sub-time step, a boundary test is performed to determine whether the stochastic process has left the domain . If this is the case, the process contributes via the second term in Eq. (7b) to the approximation of . If the stochastic process did not leave the domain until is reached, it contributes to the first term in the approximation of in Eq. (7b). to the approximation of . The computation of is handled in the analogous way. The expected values are then replaced by arithmetic means and approximated using the Monte-Carlo method. Note, it is not desirable to make itself smaller, as this would degrade the efficiency of the (deterministic) implicit single domain solver, which is described below.
Deterministic single-domain solver. The values of and along the subdomain interfaces serve as boundary conditions for the single-domain solver. The single-domain solver we employ is an implicit finite-difference discretization of Eq. (3). The matrix system is solved using LU-factorization.
Parallelization and further speed-up. It is well-known that Monte-Carlo techniques converge rather slowly [pres07a] and are usually most competitive for problems in high dimensions. The use of the stochastic solution to obtain the interface values of a DD problem only, however, is considerably more efficient and provides a fully parallel grid generation algorithm. In particular, it is not required to pass information from one sub-domain to another. Moreover, the stochastic solutions on the interfaces can be determined at each point separately and each Monte-Carlo simulation is independent. Additionally, each sub-domain solution could potentially be assigned to a single processor once the interface solutions are obtained, yielding excellent scalability. Due to the fully parallel nature of the algorithm, the method is also fault tolerant. This renders the method suitable for an implementation on massively parallel computing architectures, cf. [aceb05a, aceb10a, bihl13a].
A further source of improvement stems from the fact that not all values of and on the interfaces have to be computed using the stochastic solution (7). As proposed in [aceb05a] it may be sufficient to use the stochastic solution only at few points on the interface and recover the solution at the remaining interface points using interpolation. In [bihl13a] we have used a relatively simple optimal placement strategy to determine the most important locations on the interface where the stochastic solution should be computed. We use the same strategy in the present algorithm, i.e. the stochastic solution is computed near the maxima and minima of and along the horizontal interfaces and and along the vertical interfaces.
5 Numerical Results
We present an example our combined deterministic-stochastic DD method to generate an adaptive (moving) mesh for the weight function , where
We choose the parameters and used in [huan10a]. Both the physical and computational domain are the unit square. The grid we generate has nodes and is divided into four sub-domains. On the interfaces we determine the stochastic solution at the key points using the optimal placement strategy mentioned in the previous section. Piecewise cubic Hermite interpolation is used to determine the remaining interface points. We integrate (3) up to using . Each time step is split into sub-time steps while solving the SDE (7a) and Monte-Carlo simulations are used to estimate the expected value in (7b). The resulting meshes at , , and are depicted in Fig. 5.