The role of data transfer on the selection of a single vs. multiple mesh architecture for tightly coupled multiphysics applications

Applied Mathematics and Computation 217 (2011) 8943–8962

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The role of data transfer on the selection of a single vs. multiple mesh architecture for tightly coupled multiphysics applications Richard W. Johnson a,⇑, Glen Hansen a, Chris Newman b a b

Nuclear Science and Technology, Idaho National Laboratory, Idaho Falls, ID 83415-3840, United States T-3 Fluid Dynamics and Solid Mechanics, Los Alamos National Laboratory, Los Alamos, NM 87545, United States

a r t i c l e

i n f o

Keywords: Data transfer Multiphysics spatial integration Multiple mesh simulation

a b s t r a c t Data transfer from one mesh to another may be necessary in a number of situations including spatial adaptation, remeshing, arbitrary Lagrangian–Eulerian (ALE), and multiphysics simulation. Data transfer has the potential to introduce error into a simulation; the magnitude and impact of which depends on the application, transfer scenario, and the algorithm used to perform the data transfer. During a transient simulation, data transfer may occur many times, with the potential of error accumulation at each transfer. This paper examines data transfer error and its impact on a set of simple multiphysics problems. Data transfer error is examined using analytical functions to compare schemes based on interpolation, area-weighted averaging, and L2 minimization. An example error analysis is performed to illustrate data transfer error and behavior for a simple problem. Data transfer error is also investigated for a one-dimensional time-dependent system of partial differential equations. This study concludes that data transfer error can be significant in coupled multiphysics systems. These numerical experiments suggest that error is a function of data transfer scheme, and characteristics of the field data and mesh. If there are significant differences in the meshes in a multiple mesh simulation, this study suggests that data transfer may lead to error and instability if care is not taken. Further, this work motivates that data transfer error should be included in the estimation of numerical error when data transfer is employed in a simulation. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Computational physics simulation often involves the solution of discretized partial differential equations (PDEs) describing some physical process. The PDE is discretized spatially on a mesh with sufficiently small elements so that the solution lies within the asymptotic range of convergence of the discrete approximation. During the course of computing the solution to the PDE one might employ mesh adaptation to keep the quality of the mesh high as features develop. Such an evolving mesh would require data transfer from the previous to the revised mesh, incurring data transfer error each time the transfer is applied, which will impact the overall accuracy of the solution. As numerical simulation has become more sophisticated, so has the ability to analyze more than one phenomenon simultaneously. Models for each phenomenon might be contained in separate simulation codes, each implementing a subset of the desired physics, or the models might be within the same code. In cases where the phenomena are interdependent, either the meshes used by each model are identical or data transfer operations are needed to communicate the data from the mesh of one equation set to the mesh of the other. This transfer of data, typically represented by point values at a finite ⇑ Corresponding author. E-mail addresses: [email protected] (R.W. Johnson), [email protected] (G. Hansen), [email protected] (C. Newman). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.03.101

8944

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

number of points, to a different set of locations on a different mesh is also referred to as remapping, rezoning, or interpolation. There are advantages and disadvantages to using a single mesh vs. multiple meshes. Using a common shared mesh for coupled physics avoids data transfer altogether and is attractive for both simplicity and accuracy. However, if the spatial discretization requirements of the physics are significantly different, for example if the different solutions require refinement in different locations of the problem, one must make a compromise when a single mesh is used. Unless the single mesh is generated using a composite error metric that equidistributes the global solution error [1], a single mesh potentially overresolves some phenomena in selected regions and underresolves them other areas. In the first case, a loss of efficiency will result as more degrees-of-freedom (DOF) are being used in the calculation than needed. In the second case, underresolution of the physics will result in loss of accuracy and perhaps stability. One solution to this compromise is the use of a separate mesh to solve each problem, along with a data transfer algorithm to couple the solutions between meshes. In this case, if the solution is adapted on each mesh to support the requirements of the solution calculated there, the minimum number of elements will be used for the solution. Unfortunately, the data transfer operation between meshes is not error free, and the multiple mesh method introduces complexity in that multiple meshes must be generated and maintained, and load balanced in the case of a parallel multiphysics application. A third approach is to use a single mesh as a common basis of refinement [2]; using mesh refinement (e.g. h-FEM), order elevation (e.g. p-FEM), or both (e.g. hp-FEM) to adapt this base mesh to form multiple mesh instantiations, each tailored to the requirements of the phenomena that it hosts. In this vein, Li [3] proposes an h-FEM multiple mesh integration procedure, based on multiple hierarchically related meshes, and demonstrates it on a weakly coupled elliptic system and on two optimal control problems; one in two-dimensions and the second in three dimensions. Huang et al. [1] suggest an a posteriori error estimator where the meshes are refined such that the error estimators are equally distributed over the mesh. In other words, one needs to define a global error indicator that includes error contributions from all meshes considered as a group, and then refine the meshes for each physical model such that the global error indicator is equidistributed. Unfortunately, refinement on each mesh may affect the refinement needed on the others at particular spatial locations, which may lead to iteration of the refinement process. Li [3] presents a semi-regularization algorithm to produce a mesh on which a finite element space may be built, as well as an approach of coarsening the mesh over time. Solin et al. [4,5] present a general hp-FEM based multiple mesh integration strategy that is based on a master mesh and a set of mutually independent refinements for each physics (that all are derived using hp refinement from the master). In this method, data is not transferred from one mesh (or element) to another; it is integrated in place and the integral is assembled on the master mesh. Note that this method uses adaptation to support both the solution and physics integration; sufficient order is employed to satisfy the error criteria of both. The adaptation metric used by Solin et al. [4,5] selects elements to be refined based on the maximum error across all models and across all elements. Dubcova et al. [6] extend this work by comparing the hp multiple mesh assembly strategy to an exact approach, linear interpolation, and linear and quadratic common refinement on a multigroup neutron diffusion problem. Wang et al. [7] discuss an h-FEM multiple mesh procedure that is similar to [3], on a multigroup neutron diffusion problem. In [8], they present an hp-FEM approach similar to [4,5], and present multigroup neutron diffusion results that employ an equidistribution error metric [1]. The study by Wang et al. [7] is based on an hp-FEM code developed by Demkowicz and Kim [9]. While the use of a single common mesh and some of the multiple mesh approaches above do not require data transfer, a number of other applications are based on intermesh data transfer operations. Data transfer is used with the arbitrary Lagrangian–Eulerian (ALE) method where the mesh moves with the transient simulation but does not necessarily follow the material motion directly, requiring continual rezoning [10]. Other applications include data transfer between overlaid meshes that are used to capture different features in the simulation. Jaiman et al. [11] provide a comparative study of several conservative data transfer schemes for a problem involving fluid–structure interaction. A conservative method is one that approximates the underlying conservation law faithfully. For example, a mass conservative approximation is one that neither loses or gains mass when it is employed. In a finite element method, element conservation means that the weighting function can be set exactly to one in the element and zero elsewhere [12]. In the case of the finite volume method, a similar statement holds for the volumes. A detailed discussion of conservation in the context of the finite element method is presented by Hughes et al. [12]. Both discretization and data transfer methods are typically classified as conservative or nonconservative. It is argued by Jiao and Heath [13] that data transfer methods need to be both numerically accurate and physically conservative. Conservation and accuracy are not necessarily congruent at a given level of discretization; the higher-order (but nonconservative) cubic spline interpolation is shown to have greater accuracy than many of the other schemes applied. The present study takes the view that conservation is generally desirable, but results will be presented that show that accuracy may be more important than the selection of a conservative method when it comes to data transfer between meshes. Our objective is to illustrate potential issues that may arise from the application of data transfer operations through numerical experiments in simplified scenarios that involve disparate length scales. Our experiments suggest that in the application of data transfer operations, error estimates should capture data transfer error as well as discretization error. Section 2 reviews and classifies several data transfer schemes; Section 3 applies these schemes to the data transfer of static functions in one and two dimensions; Section 4 provides an example to illustrate the nature of data transfer error; Section 5 examines the effects of data transfer in a coupled, unsteady system of two partial differential equations that represent a multiphysics reaction diffusion system; and Section 6 presents conclusions from the work.

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

8945

2. Data transfer methods Jiao and Heath [13] present three schemes based on L2 minimization and compare to consistent (linear), thin plate spline and cubic spline interpolation schemes using multiple data transfer cycles of analytical functions. The schemes are used to transfer data from a source mesh to a target mesh and back multiple times. Our study adds first order area-weighted averaging to the methods they presented and examines the effects of failing to capture extrema present in the source solution. The approach employed by Jiao and Heath [13] demonstrates cumulative data transfer error over the course of multiple data transfer cycles which provides an appreciation for the error that may be encountered for a transient problem where data transfer occurs at each time step. We classify these schemes into three categories; interpolation, area-weighted averaging and L2 minimization. 2.1. Interpolation Interpolation methods for data transfer are based on the construction of interpolants using nodal values from the source mesh. If the interpolation curve is of the same order as the basis function used to compute the source mesh data, the method is called consistent interpolation [13]. If the interpolant constructed is of a different order than used on the source grid, the method is called alternative interpolation. A common data transfer method used in applications is simple linear interpolation, which is a consistent interpolation method if linear basis functions are used to compute the values on the source mesh. Examples of alternative interpolation are quadratic and cubic spline interpolation using nodal values on the source grid. Another example is a B-spline curve of higher order fitted to the source data. As a class, interpolation schemes are not considered to be conservative. 2.2. Area-weighted averaging Area-weighted averaging schemes are typically associated with the finite volume method where conservation is required for each cell in the mesh. Values on the target mesh, usually at cell centers, are determined by integration over a common refinement of the two meshes, computing values for new cells by integration over the subcells of each new cell. The integrated value of the mesh data is preserved on the domain. First and second order area-weighted averaging schemes for application to arbitrary Lagrangian Eulerian methods are discussed in [10,14–18]. These schemes use the divergence theorem to convert area integration into line integration on edges of old and new cells, which are summed on each new cell to obtain new values. First order rezoning methods assume a constant value for the conserved quantity of the cell domain and are discussed in [14,15,17]; second order methods, which assume a linearly varying conserved quantity are discussed in [10,16,18]. 2.3. L2 minimization methods L2 minimization for data transfer is typically associated with the use of the finite element method [13]. Both the source and target functions are assumed to be expanded in a finite element basis; the source function is given by P Pn f ¼ m i¼1 fi /i and the target function by g ¼ i¼1 g i wi , where /i and wi are basis functions, fi and gi are the interpolation values, and m and n are the number of nodes in the source and target meshes, respectively. The L2 norm, kg  fk2, is minimized when

@ @g i

Z

ðg  f Þ2 dx ¼ 0;

i ¼ 1; . . . ; n;

ð1Þ

X

where X is the function domain. Eq. (1) is satisfied when n Z X j¼1

X

wi wj dx g i ¼

Z

wi f dx;

i ¼ 1; . . . ; n:

ð2Þ

X

The integral on the left hand side of (2) is computed over elements of the target mesh in all cases, as the coefficients gi are the unknown physics values on the target mesh. The integral on the right hand side involves function values from the source mesh and basis functions from the target mesh. This integral can be integrated over either elements of the source or target mesh. Integration of the right hand side term of (2) over the source mesh results in the source scheme and integration of the right hand side term of (2) over the target mesh results in the target scheme [13]. One objective in [13] is to present results for a scheme where the term on the right hand side of (2) is integrated over elements of a common refinement of the two meshes. This is accomplished by constructing a commonly refined mesh, defining the integration points on the commonly refined mesh, then performing the numerical quadrature by evaluating the source mesh function and target mesh basis functions at the integration points. The order of quadrature used is high enough to yield exact integration of the integrands of (2). Fig. 1 illustrates the finite element basis functions on the source and target R meshes and provides a schematic of the computation of X wi f dx on the commonly refined mesh. If higher order basis functions are present, higher order quadrature is typically used. Arguments are presented in [13] that show that source and common refinement approaches are conservative while target-based integration is not conservative. They further show that the common refinement scheme delivers more accurate

8946

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

Source basis functions

Target basis functions

Load vector on commonly refined grid Gauss pts.

+

+

+

+ ++

+

+

+

+ +

+

+

+

++ +

+

+

+

Fig. 1. Linear basis functions for the source mesh (top), the target mesh (middle), and a schematic of the computation of the right hand side of (2) using values of f from the source mesh, the basis functions of the target mesh and the integration points of the commonly refined mesh (indicated by ‘‘+’’ signs along the bottom).

results over many data transfer cycles than either source or target methods. Table 1 illustrates the classes of schemes investigated in the present study. 3. Data transfer error: static functions The mesh discretization philosophy used in this study was designed to mimic common practice for large multidimensional problems where coarse discretization near the upper limit of the asymptotic range of convergence is typically employed. For such problems, the use of excessively fine discretization is avoided due to both run time and memory considerations. Neither of these is of practical concern for the simple problems considered here. However, the use of excessively fine discretization here might well render the study useless for the target application space. Further, this study employed both coarse and fine mesh resolutions. The fine mesh was chosen to contain from 5% to 50% more nodes than the coarse, to prevent the fine mesh from getting too fine and violating the above considerations. One approach for error quantification, proposed in [13], employs analytical functions for one- and two-dimensional meshes that are transferred from a coarse to a fine mesh and back, defining one transfer cycle. Multiple cycles are used to emphasize the relative accuracy of the several schemes employed. The technique of applying many cycles to a static function is also performed to examine possible cumulative effects of data transfer applied over a simulation spanning many time steps. As in [13], we examine the Runge function,

f ðxÞ ¼ 1=ð1 þ 25x2 Þ;

1 6 x 6 1;

ð3Þ

as a one-dimensional example. Fig. 2 illustrates the results of 64 transfer cycles for source, target, linear interpolation, areaweighted averaging, common refinement using both linear and quadratic basis functions and cubic spline schemes. Results for the source scheme are limited to 8 cycles, as the results using this method deteriorate quickly. Note that the mesh elements for the quadratic common refinement scheme use quadratic basis functions (Com ref 2) with mesh spacing exactly twice the length of the elements using linear basis functions. The fine mesh consists of 45 uniformly spaced nodes, and the coarse meshes consist of 32 and 31 uniformly spaced nodes. The notation [45, 32] indicates the numbers of nodes for the fine and coarse meshes, respectively. Fig. 2 shows that both the conservative source scheme and the nonconservative target scheme produce erratic results between the meshes consisting of 32 and 31 nodes. For nonconservative linear interpolation between the 32 and 31 node coarse meshes, the coarse mesh of 31 nodes locates a node exactly at the peak of the Runge function (3), and shows very different results from the mesh where there is no node at the peak. This result shows the importance of capturing extrema

Table 1 Categorization of data transfer methods. Category

Interpolation

Conservative Nonconservative

Linear Cubic spline

Area-weighted avg.

L2 minimization

First order Second order

Source Common refinement Target

8947

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 -1

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

0.6

0.8

0 -1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 -1

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

0.6

0.8

0 -1

1

Runge fctn. Com ref [45,31] Com ref [45,32] Com ref2 [45,31]

1

0.6

0.6

y

0.8

0.4

0.4

0.2

0.2

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

0.6

0.8

1

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

0 -1

0.6

0.8

1

Runge fctn. Area-wt [45,31] Area-wt [45,32]

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

1

0.8

0 -1

Runge fctn. Target [45,31] Target [45,32]

1

Runge fctn. Linear [45,31] Linear [45,32]

y

y

1

y

1

Runge fctn. Source [45,31] Source [45,32]

y

y

1

0.6

0.8

1

Runge fctn. Cubic [45,31] Cubic [45,32]

-0.8 -0.6 -0.4 -0.2

0

x

0.2

0.4

0.6

0.8

1

Fig. 2. Data transfer for (3) using source, target, linear interpolation, area-weighted averaging, common refinement using both linear (Com ref) and quadratic basis functions (Com ref 2), and cubic spline schemes for 64 cycles. Number of nodes used in the fine and coarse meshes are indicated in brackets.

in the data being transferred by locating nodes very close to the peak. Note that conservative area-weighted averaging diffuses the peak when a sampling location coincides with the peak (sampling locations for the area-weighted averaging are halfway between nodes.) The data transfer error for linear interpolation is examined in more detail in Section 4. To study the effect of the function shape on these observations, a second set of results are generated using the sine function,

f ðxÞ ¼ 0:6 þ 0:1 sinðpxÞ:

ð4Þ

8948

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

Data transfer operations are performed on uniform meshes [23, 16], which have similar node spacing as for the Runge function study, but on the domain 0 6 x 6 1. Nodal values are used to compute normalized error for the Runge and sine functions using kfnumerical  fexactk2/kfexactk2. Fig. 3 illustrates the normalized error as a function of data transfer cycles for the previous data transfer schemes using both the Runge (for [45, 32] grids) and sine functions as input. Comparing linear versus quadratic common refinement schemes and linear versus cubic interpolation, it is clear that using a higher-order basis function or interpolation method decreases the data transfer errors incurred. Further, the magnitude of the normalized error depends on the function being transferred. The conservative linear and quadratic common refinement schemes, which are based on L2 minimization, perform very well over the 64 data transfer cycles. Finally, the nonconservative cubic spline scheme shows the best accuracy of all the methods considered. While it is an interpolation method, it performs much better than linear interpolation in terms of accuracy (exhibiting the least error for a given number of data transfer operations). Next, a cyclical data transfer process between fine and coarse grids for the two-dimensional peaks function, which exhibits multiple peaks and valleys on a square domain is performed. The approach taken here is identical to [13], but with the addition of first and second order area-weighted averaging schemes. The peaks function is 2 ðyþ1Þ2

f ðx; yÞ ¼ 3ð1  xÞ2 ex

x  2 2 1 2 2  10  x3  y5 ex y  eðxþ1Þ y 5 3

ð5Þ

for 3 6 x, y 6 3. The peaks function is illustrated in Fig. 4. The 2-D coarse and fine meshes are defined by 32  32 and 45  45 uniformly spaced nodes. The peaks function is transferred from the coarse to the fine mesh and back to form one cycle. A measure of the accuracy of the 2-D data transfer schemes is provided by contour plots of the original function and after data transfer. Fig. 5 shows the results for linear and cubic interpolation as well as first and second order area-weighted averaging schemes after 32 transfer cycles. Fig. 6 illustrates the normalized error defined as above for seven 2-D data transfer schemes. It is interesting to compare the 1st and 2nd order area-weighted averaging schemes. The latter is significantly superior to the former and compares well with common refinement for less than 8 cycles. Area-weighted averaging schemes might be considered ‘‘finite-volumebased common refinement methods;’’ the second-order area-weighted scheme assumes linear variation of the cell value, as does the linear common refinement approach. Experiments show that cubic spline interpolation is the most accurate data transfer method; common-refinement approaches do not perform as well. The two other conservative approaches, source and first-order area weighted averaging, did not perform as well as the nonconservative cubic spline method. It is noted that both higher order as well as L2 minimization schemes may produce overshoot and undershoot for source functions that exhibit large gradients and/or curvatures [13]. The investigation in the present study is limited to the behavior of data transfer for smooth functions. 4. An example of data transfer error In this section, a simple Taylor series error analysis using a one-dimensional problem is performed to motivate the potential for the introduction of error when data is transferred from one mesh to another. Assume that two functions f and g are defined on the identical interval [xi, xi+1], as shown in Fig. 7. Here, there is no distinction between whether these are approximate functions obtained from a numerical solution or analytical functions. For this example, f is increasing monotonically over the interval and g has a local maximum. The function values are known at the ends of the interval. These locations and values could be nodal points on a finite element mesh or face values on a finite volume mesh. The goal is to

Normalized error

10

0

10

-1

10

-2

10

-3

Source, Runge Linear, Runge Target, Runge Com ref, Runge Com ref 2, Runge Cubic, Runge Source, sine Linear, sine Target, sine Com ref, sine Com ref 2, sine Cubic, sine

10-4 10-5 10-6 10

0

1 10 Data transfer cycles

10

2

Fig. 3. Normalized error as a function of data transfer cycles for (3) and (4) using common refinement, cubic spline, common refinement 2, source, target, and linear interpolation.

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

8949

8 6 4 2 0 -2 -4 -6 -8

3

2

1

0

y

-1

-2

-3

-3

-2

-1

0

1

2

3

x

Fig. 4. The peaks function (5) used for 2-D data transfer.

calculate the approximate value for a function at the spatial location xr for transfer to another mesh which has a node located at xr. First, define the distances h1 = xr  xi and h2 = xi+1  xr, where Dx = xi+1  xi. Further, expand f(xr) in Taylor series about xi, 2

f ðxr Þ ¼ f ðxi þ h1 Þ ¼ f ðxi Þ þ h1 f 0 ðxi Þ þ

h1 00 3 f ðxi Þ þ Oðh1 Þ; 2

ð6Þ

or more compactly, 2

fr ¼ fi þ h1 fi0 þ

h1 00 3 f þ Oðh1 Þ: 2 i

ð7Þ

Similarly, an expansion about xi+1 results in 2

0 fr ¼ fiþ1  h2 fiþ1 þ

h2 00 3 f þ Oðh2 Þ: 2 iþ1

ð8Þ

Next, the value of function f at location xr is estimated using a given data transfer scheme (e.g. one in Fig. 3). The form of this approximation depends on the data transfer approach employed. For linear interpolation the approximation is

fDT ðxr Þ ¼ fi

h2 h1 þ fiþ1 : Dx Dx

ð9Þ

The error in estimating the function using the selected data transfer scheme is

Eðxr Þ ¼ f ðxr Þ  fDT :

ð10Þ

Substitution of (7) and (9) into (10) yields 2

Eðxr Þ ¼

h1 h 3 ðfi  fiþ1 Þ þ h1 fi0 þ 1 fi00 þ Oðh1 Þ: Dx 2

ð11Þ

Note that the leading term is linear in Dx. This demonstrates that data transfer error for linear interpolation is first order. Subtracting (8) from (7) and substituting into (11) yields 2

2

h1 h h 00 0 Eðxr Þ ¼  1 fi00 þ 2 fiþ1 h1 f10  h2 fiþ1 Dx 2 2

!

2

þ h1 fi0 þ

h1 00 f þ OðDx3 Þ: 2 i

ð12Þ

Simplifying,

Eðxr Þ ¼

   h1 h2 0 1 0 00 þ h1 fi00 þ h2 fiþ1 fi  fiþ1 þ OðDx3 Þ: 2 Dx

ð13Þ

This result represents the data transfer error at xr for the functional form f using linear interpolation as the data transfer mechanism. Examining this result in detail, E(xr) ? 0 as one approaches the nodes xi or xi+1 (as xr ? xi, h1 ? 0). Further, it is intuitive that if f is well approximated by the order of the interpolation scheme over the interval, the error over the interval is min0 00 imized. This is clear for linear interpolation as if f is linear, fi0 ¼ fiþ1 and fi00 ¼ fiþ1 ¼ 0, and thus E(xr) = 0.

8950

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

2

2

1

1

0

0

-1

-1

-2

-2

-3 -3

-2

-1

0

x

3

1

2

-3 -3

3

2

2

1

1

0

0

-1

-1

-2

-2

-3 -3

-2

-1

0

x

1

2

3

Cubic 32 c.

-2

-1

3

Area-wt 1st ord. 32c.

y

y

3

Linear 32 c.

y

y

3

-3 -3

0

x

1

2

3

Area-wt 2nd ord. 32 c.

-2

-1

0

x

1

2

3

Fig. 5. Contour maps showing the application of 32 data transfer cycles using linear interpolation, cubic spline, and 1st and 2nd order area weighted schemes. The solid lines show contours of the peaks function, the dashed lines show the locations of these lines after 32 transfer cycles. The amount of error inherent in each method is evident by how close the dashed contours are to the original data.

10

Normalized error

10

0

-1

10

-2

10

-3

10

-4

10

Linear Cubic Area-wt1st Area-wt2nd Source Target Comref

0

1 10 Data transfer cycles

10

2

Fig. 6. Normalized errors for the 2-D data transfer schemes for the peaks (5) function.

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

fi xi

h1

xr

h2

f i +1

gi

x i +1

xi

8951

g i +1 h1

xr

h2

x i +1

Fig. 7. Functions f and g on the interval [xi, xi+1]. The point xr is the location at which the value of g will be transferred to the equations of f, and vice versa. 0 Taking this analysis further, if an extremum exists on the interval, fi0 and fiþ1 will be of different sign. Thus, if one ignores the second derivative terms, the error on the interval depends on h1 and h2,

Eðxr Þ ¼

h1 h2 h 1 ð Dx  h 1 Þ C¼ C; Dx Dx

ð14Þ

where C is some constant. This exhibits extrema at

E0 ðxr Þ ¼ 1 

2h1 ¼ 0; Dx

ð15Þ

or

h1 ¼ h2 ¼

Dx ; 2

ð16Þ

i.e. the point of maximum error occurs at the midpoint of the interval. This example suggests the importance of effectively capturing the behavior of the solution being transferred by use of an appropriate data transfer function and discretization. Important considerations include:  Use of mesh refinement as a strategy to support accurate data transfer reduces Dx and thus the error in (13). This may be very important in areas where the data transfer algorithm does not match the order of f.  In intervals with significant extrema, placing an additional node point at the peak of f serves to decrease the magnitude of 0 fi0  fiþ1 significantly. In the results that follow the importance of locating nodes near extreme values is demonstrated numerically. 5. Multiphysics numerical experiments This section explores the effects of data transfer applied to the Brusselator problem [19], a coupled, nonlinear parabolic partial differential equation system. To most efficiently solve this problem, one would likely use a single mesh implementation. Here, each equation is computed on a different mesh to serve as a simple data transfer study. In this example, the solution for each equation is transferred from one mesh to the other each time step, using cubic spline, common refinement, source, target and linear interpolation data transfer schemes. Exchanging data in this manner, between the coupled problems at the beginning of each timestep, is usually termed explicit coupling of the multiphysics problem. Explicit coupling of multiphysics problems is often used for efficiency reasons, thus it is employed here. The accuracy of the approach can be very good, especially when small timesteps are used. 5.1. Brusselator system The Brusselator is a system used to model a chemical reaction diffusion system, coupled through the nonlinear source terms [20],

@T @2T ¼ D1 2 þ a  ðb þ 1ÞT þ T 2 C; @t @x @C @2C ¼ D2 2 þ bT  T 2 C; @t @x

ð17Þ

where X is the one dimensional domain a 6 x 6 b, T(a, t) = T(b, t) = a and C(a, t) = C(b, t) = b/a. The system (17) has steady state, oscillatory and chaotic solutions. With the choice of a = 0.6, b = 2.0 and D1 = D2 = 0.025, T and C are oscillatory [19]. The initial conditions, T(x, 0) = f(x), C(x, 0) = g(x), can be specified as any smooth functions, which allows the study of the effects of the initial conditions and evolving solutions on aggregate numerical error, including data transfer error. This study applies three sets of initial conditions to investigate their effects on data transfer error:

f ðxÞ ¼ 0:6 þ 0:1 sinðpxÞ; 10 þ 0:1 sinðpxÞ; gðxÞ ¼ 3

ð18Þ

8952

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

1 ; 1 þ 25x2 10 1 ; þ gðxÞ ¼ 3 1 þ 25x2

ð19Þ

expðaxÞ  1 ; expðaÞ  1 10 1 þ ; gðxÞ ¼ 3 1 þ 25x2

ð20Þ

f ðxÞ ¼ 0:6 þ

f ðxÞ ¼ 0:6 þ

where a = 15 in (20). 5.2. Numerical solution and convergence Eq. (17) are spatially discretized using either linear or quadratic finite elements [21] and temporally using the explicit Euler method [22]. The resulting linear equations are solved using the GMRES linear iterative solver in Trilinos [23], using a relative convergence tolerance of 1  104. To verify the temporal and spatial orders of convergence and determine the asymptotic ranges of convergence, numerical tests are performed by varying spatial and temporal discretizations. Uniform common meshes for T and C are used for the convergence tests. The three pairs of initial conditions, given by (18)–(20), are used in separate convergence tests. The domain intervals used are 0 6 x 6 1 for the first set of initial conditions (18) and 1 6 x 6 1 for (19) and (20). The solutions are compared to reference solutions Tref and Cref that employ fine discretizations of Dx = 5.0  104 and Dt = 5.0  105. Figs. 8 and 9 show the error curves for the temporal and spatial discretizations. The normalized L2 norms of the error are computed using

R e ¼ 100 

X ½Tðx;

sÞ  T ref ðx; sÞ2 dx

R

X

T 2ref ðx; sÞdx

1=2

1=2

ð21Þ

;

where s is the time at which the solutions are compared; s = 1 for (18) and (19), and s = 3 for (20). Convergence rates are first order in time; spatial convergence is second order for linear and third order for quadratic finite elements. All subsequent calculations use Dt = 5.0  105, which is well within the temporal asymptotic convergence range. 5.3. Sine and Runge function initial conditions

L2 norm relative error

In this section, results for computations based on using similar functions for both initial conditions are presented for two sets of similar functions, sine and Runge functions; results are then compared and discussed. Calculations are performed for

10

-2

10

-3

10

-4

T 1-sine C 2-sine T 1-Runge C 2-Runge T Expon. C Runge 1st order 2nd order

10-5 10-6 10

-7

10

-8

10

-3

Δt

10

-2

Fig. 8. Relative error (21) vs. timestep size for the Brusselator (17), where the source terms are (18) (sine, sine), (19) (Runge, Runge) and (20) (Expon., Runge). Results shown use Dx = 5.0  104 and are relative to a reference solution using a small timestep (Dt = 5.0  105). Uniform common meshes are used for these results, thus there is no error from data transfer.

8953

10

-1

10

-1

10

-2

10

-2

10

-3

10

-3

10-4 T sine C sine T Runge C Runge Second order Third order

10-5 10

-6

10

-7

10

-2

10

Δx

L2 norm relative error

L2 norm relative error

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

-1

10-4 T Exp. C Runge Second order Third order T Exp. quad. C Runge quad.

10-5 10

-6

10

-7

10

-2

Δx

10

-1

Fig. 9. Relative error (21) vs. mesh size for the Brusselator (17), where the source terms are (18) (sine, sine), (19) (Runge, Runge) and (20) (Expon., Runge) for both linear and quadratic finite elements. Results shown use Dt = 5.0  105 and are relative to a reference solution using small mesh elements (Dx = 5.0  104). Uniform common meshes are used for these results, thus there is no error from data transfer.

(17) using sine-function based initial conditions given by (18), on intervals 0 6 x 6 1 and 0 6 t 6 1 for uniform meshes of 16 (Dx  0.067) and 23 (Dx  0.045) nodes for T and C, respectively. Also, calculations are made using the Runge function-based initial conditions, (19), on intervals 1 6 x 6 1 and 0 6 t 6 1. These meshes consist of 32 (Dx  0.065) and 45 (Dx  0.045) nodes for T and C, respectively. Dt = 5.0  105 is used. All spatial and temporal discretizations are within the asymptotic ranges shown earlier. Mesh node locations are selected to produce similar mesh spacings as used for the data transfer of similar functions in Section 3. Data transfer is performed at each timestep using cubic spline, common refinement, source, target, and linear interpolation. Error is then estimated for the several data transfer schemes using the L2 error (21) relative to the fine mesh reference solution. In a coupled, multiphysics solution such as this one, the L2 errors do not reflect error solely due to data transfer. They also include the usual discretization error from the discrete solution in addition to propagation of errors that may occur from the interaction of the discretization and data transfer errors and the equations being solved. It is not straightforward to isolate the contribution to the aggregate L2 error solely due to data transfer, but one can compare the aggregate error with the discretization error. In order to compare aggregate error to the discretization error, two reference solutions are computed. The reference solutions all employ common meshes for the two equation variables (no data transfer is used). The first reference solution uses a very fine mesh of 4001 nodes for both sets of initial conditions. The second reference solutions using 23 and 45 nodes for the sine and Runge function, respectively. The second (coarse) reference solutions are used to estimate the discretization error compared to the fine mesh reference solutions. The total aggregate error, including discretization and data transfer error, are

0.9

T

0.8

Cubic Com ref Source Target Linear fine mesh ref sol

0.7

0.6

0

0.2

0.4

x

0.6

0.8

1

Fig. 10. T at s = 1 for sine initial conditions (18) for cubic spline, common refinement, source, target and linear interpolation data transfer schemes compared to the fine mesh reference solution.

8954

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

Table 2 Relative error (21) for the sine and Runge initial conditions explored in Figs. 10 and 11.

a

Com ref

Cubic

Linear

Source

Target

Coarse mesha

Sine fctn. Net error

0.17 0.08

0.18 0.10

0.11 0.02

1.48 1.40

0.17 0.08

0.09

Runge fctn. Net error

0.31 0.17

0.63 0.48

1.36 1.22

2.27 2.13

0.32 0.17

0.14

No data transfer employed.

Fig. 11. T at s = 1 for Runge initial condition (19) using cubic spline, common refinement, source, target, and linear interpolation schemes. The figure in the rectangle is a magnification of the region 0 6 x 6 0.3. The initial condition is also shown.

computed from the solutions involving data transfer compared to the fine reference solution. The net error including data transfer error can then be estimated from the differences between these two errors. Fig. 10 illustrates the results for T at s = 1 for cubic spline, common refinement, source, target, and linear interpolation data transfer schemes compared to the fine reference solution for the sine function-based initial conditions (18). Note that the symbols for the fine mesh solution are spaced 100 nodes apart. The results for the several schemes are very similar, except for the source scheme, which falls significantly below the other solutions. The L2 errors for the common mesh problem and for those involving data transfer are tabulated in Table 2. Also shown is the net error including data transfer error obtained by subtracting the estimated discretization error from the aggregate error for each data transfer scheme. Fig. 11 shows T at s = 1 using the data transfer schemes for use of the Runge function initial condition (19). Also shown are the initial condition and the fine mesh reference solution (with symbols spaced 100 nodes apart). There is a clear variation of accuracy for the data transfer schemes considered. Common refinement and target schemes show similar results that match the reference solution closely. Both source and linear interpolation schemes show poor accuracy. Cubic spline behaves slightly poorer than common refinement and target. Table 2 shows a summary of the L2 error for the above. Table 2 shows that transfer errors are larger using the Runge function initial conditions than the sine function. Additionally, the various data transfer schemes perform differently. For instance, linear interpolation performs extremely well for the sine function, but second poorest for the Runge function. However, the source scheme is seen to perform the poorest for both. These observations lead to the conclusion that the nature of the function has an important influence on the data transfer error, as well as the particular scheme used. Eq. (13) suggests that the second derivative of the solution may be an important consideration in the data transfer error. Towards this end, we chose to investigate this numerically by considering the curvature [24] of the solution, 2

j¼h

jd f ðxÞ=dx2 j 1 þ ðdf ðxÞ=dxÞ2

i3=2 ;

where f(x) is any twice-differentiable function. The curvatures of (18)–(20) over 0 6 x 6 1 are shown in Fig. 12. Note that x 2 [0, 1] in the figure whereas (19) and (20) are defined on x 2 [1, 1]. The curvature for (20) can be adjusted by modifying the parameter a. The maximum curvature of the sine function is the lowest, while that for Runge function is the greatest. It can be seen that the curvature of the Runge function is almost always greater than that for the sine function. This is one

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

8955

reason that all of the data transfer schemes show greater data transfer error for the Runge function than for the sine function. Linear interpolation will be exact for a linear function and less accurate as the curvature increases. This leads one to conclude that the curvature of a function will have an important influence on the accuracy of a particular data transfer scheme used for a multiphysics problem requiring data transfer. 5.4. Combined exponential and Runge function initial conditions One challenge of multiphysics simulation is that solutions that are strongly coupled to each other may exhibit different length scale phenomena at different spatial locations in the mesh. To explore the potential for data transfer error in such situations, in this study we use diverse functions for the initial conditions. Towards this goal, an exponential function along with a Runge function (20) over 1 6 x 6 1 are chosen. These functions are shown in Fig. 13; f(x) exhibits a ramp near x = 1 while g(x) has a pronounced ‘‘hump’’ at x = 0. The objectives are to examine the effects of data transfer scheme, mesh coarseness, and mesh disparity on numerical error given these initial conditions. Mesh disparity refers to how different the node spacings are between the two meshes where data transfer is taking place.

Curvature

5.4.1. Uniform mesh results The first study examines data transfer between uniform meshes of 40 (Dx  0.051) and 39 (Dx  0.053) nodes (indicated by [40, 39]), both in the asymptotic range, for T and C, respectively, given in (17). Because the source scheme showed poor results and because the target scheme showed accuracy similar to common refinement in Section 5.3, these studies examine only cubic spline, common refinement, and linear interpolation schemes. A fine mesh reference solution was calculated using Dx = 5.0  104 for both meshes, and using a constant Dt = 5.0  105. Fig. 14 details the results for the data transfer calculations for T at s = 3 s, compared to the reference solution and the initial condition. The L2 errors, given by (21) using the fine mesh solution as the reference, are 0.71, 0.73 and 0.94 for common refinement, cubic spline and linear interpolation, respectively. Computing a solution at s = 3 for common meshes [40, 40], results in an L2 error of 0.73, compared to the fine reference solution. This is of the same magnitude as the above aggregate errors. Given this result, one might conclude that for this problem the data transfer error using all three methods is inconsequential. To investigate the effects of mesh coarsening on the data transfer error, new solutions are computed using uniform meshes [30, 29] (Dx  0.069, 0.071), which are both in the asymptotic range. Fig. 15 illustrates the results for T and C for the above data transfer methods in comparison to the reference and initial solutions. The L2 errors (21) based on the fine mesh solution, are 1.28, 1.31 and 1.73 for common refinement, cubic spline and linear interpolation. The L2 error for the nontransferred solution using [30, 30] meshes is 1.31. Again, the aggregate error is similar to the error due to the discretization error. Further, the aggregate error has increased due to mesh coarsening as would be expected. The next study examines results obtained by refining the meshes and changing their spatial disparity to [45, 33] (Dx  0.045, 0.063). Fig. 16 shows results using common refinement, cubic, common refinement 2 and linear interpolation schemes for both T and C at s = 3. Note that quadratic finite elements are used in conjunction with the common refinement 2 scheme; the quadratic finite elements are twice as large as the linear elements. In this result, the aggregate error has increased for linear interpolation for T at the sharp peak near x = 0.9 as compared to Figs. 14 and 15. The L2 norm of the error in T is 0.55 for common refinement, 0.64 for cubic spline, 0.14 for common refinement 2 and 3.17 for linear interpolation compared to the fine mesh reference solution. While the results for linear interpolation for the refined T mesh are better for x 6 0.6, the overestimation of the solution for T near x = 0.9 is significantly worse than for the [30, 29] mesh (cf Fig. 15). Both

10

2

10

1

10

0

10

-1

10-2

10-3

Sine fctn Runge fctn Exponential fctn

0

0.2

0.4

x

0.6

0.8

Fig. 12. Curvature of functions in (18)–(20) (for x 2 [0, 1]).

1

8956

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

1.2

T, initial cond., Exponential fctn C, initial cond., Runge fctn

1

T, C

0.8 0.6 0.4 0.2 0 -0.2 -1

-0.5

0

0.5

1

x Fig. 13. Eq. (20) used for initial conditions.

Fig. 14. T at s = 3 for common refinement, cubic spline and linear interpolation schemes for uniform meshes [40, 39] compared to the fine mesh reference solution and the initial condition. The right figure is a magnification of the area near the peak.

meshes are finer than the [30, 29] meshes. Results for C also indicate that the error for linear interpolation near x = 0.8 on this finer mesh is noticeably greater than for the [30, 29] mesh. This suggests that the disparity between the 45 and 33 node meshes has a significant effect on the aggregate numerical error, more so than mesh coarseness. To summarize the uniform mesh results, aggregate error was of the same magnitude as the discretization error for all calculations except linear interpolation on the [45, 33] problem. One can conclude from this that, given the data transfer choices considered, data transfer error significantly degraded solution quality only in the [45, 33] mesh. Further, given these limited results, even within the spatial asymptotic range of the physics approximations one needs to verify that the data transfer scheme chosen for the problem supports the disparity in meshes present (or expected) during that simulation. 5.4.2. Nonuniform mesh results The next phase of this study examines effects of using nonuniform meshes optimized for the initial conditions on the aggregate error. Optimized meshes consisting of 67 and 45 nodes for the exponential and modified Runge functions (20) are used. The mesh optimization to the initial source data was accomplished by specifying a minimum cell size located at the point of greatest curvature [25]. A smoothing algorithm was used to increase element size as curvature decreases; this smooth expansion of cells from the minimum was limited to a factor of 1.15. The minimum cell size at the location of highest curvature in the initial conditions (20), the expansion factor and forcing terms (method of manufactured solutions) are selected such that solution norms of the ordinary differential equation system

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

8957

Fig. 15. T and C at s = 3 for common refinement, cubic spline and linear interpolation schemes for uniform meshes [30, 29] compared to a fine mesh reference solution and the initial condition.

Fig. 16. T and C at s = 3 for common refinement, cubic spline, common refinement 2, and linear interpolation for uniform meshes [45, 33] compared to a fine mesh reference solution and the initial condition.

2

d T 2

dx

2

¼ ST ;

d C 2

dx

¼ SC ;

ð22Þ

where

ST ¼ f 00 ðxÞ; 00

SC ¼ g ðxÞ;

f ðxÞ ¼ ½expðaxÞ  1=½expðaÞ  1; 2

gðxÞ ¼ 1=ð1  25x Þ;

a ¼ 15;

ð23Þ

on 1 6 x 6 1, T(±1) = f(±1), C(±1) = g(±1), are within 0.32% of the exact solution norms. The maximum element sizes (Dx  0.096, 0.091) are within the spatial asymptotic convergence range computed above for uniform meshes. Note that this curvature adaptation metric adapts each mesh to the curvature of the source function, f or g, that it hosts. While it is important to adapt each instance to the spatial resolution requirements of the physics it supports, the following results show that this strategy does not guarantee accuracy of the coupled problem. Solutions for T and C using common refinement, cubic spline, common refinement 2 and linear interpolation compared to the fine mesh and initial solutions for s = 3 are shown in Fig. 17. Note that for the quadratic finite elements used with the common refinement 2 scheme, the elements consist of three nodes with the intermediate node located at the center of the element; the nonuniform meshes used for quadratic finite elements are constructed from the meshes used for linear finite elements with intermediate nodes adjusted to lie at element centers. The solution for T at the location of the peak near x = 0.9 significantly overshoots the fine mesh solution for linear interpolation and common refinement schemes; common refine-

8958

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

ment 2 results are slightly high and cubic spline results are slightly low at the peak. The L2 norm of the error for T is 2.51 for common refinement, 0.57 for cubic spline, 0.51 for common refinement 2 and 6.59 for linear interpolation compared to the fine mesh reference solution. Comparing the results for C in the vicinity of the minimum near x = 0.8, it is obvious that the error is worse than for the previous uniform [45, 33] mesh. The disparity in the two meshes is greater in this case because the T mesh has been refined while the C mesh has been coarsened in this vicinity. Again, this supports the conclusion that mesh disparity is the reason for the large errors for the less accurate data transfer schemes. While the results for cubic spline and common refinement 2 approaches are of the same magnitude as discretization error, both the linear interpolation and common refinement solutions depart significantly from the peak near x = 0.9, likely by an unacceptable amount. Even though the maximum element sizes on both meshes fall within the asymptotic ranges of the solutions, the disparity of meshes near the peak combine with the tight coupling between the equations and the data transfer algorithm to result in significant spatial error. This leads to a hypothesis that even though the separate physics equations are spatially resolved for their initial conditions, the coupled multiphysics problem (which includes the chosen data transfer algorithm), is not resolved at the location of the T peak. Indeed, even though the meshes are adaptively refined to their respective source functions, the C mesh is probably too coarse near the T peak to support the data transfer operation there. To provide additional evidence in support of this hypothesis, a similar nonuniform problem was considered where the C mesh was refined to provide results for meshes [67, 67]. The refinement in the C mesh, although still targeting only the source term for that equation, now provides enough sampling near the peak in the T solution to support data transfer operations performed at that location. The L2 norm of the error is 0.35 for common refinement, 0.42 for cubic spline, 0.05 for common refinement 2 and 0.91 for linear interpolation compared to the fine mesh reference solution. 5.4.3. Coarse nonuniform mesh results While a numerical solution that lies outside the asymptotic convergence range can be suspect, the proper element size for a new problem is not always clear, initially. Further, large three-dimensional multiphysics analysis may be initiated with a mesh that is understood to be too coarse to verify the construction of the input data set without spending an appreciable amount of resources on a flawed specification. Poor input data may lead to instability or the simulation failing; it is routine to check the validity of input data using an abbreviated problem. The previous result appears to indicate that, for a nonuniform coupled problem with data transfer, one must exercise caution when the mesh is coarsened. If one coarsens the nonuniform grids such that the largest elements are just larger than the element size needed to fall within the spatial range of asymptotic convergence, one obtains undesirable results. Consider an example of nonuniform grids [55, 37] with maximum element sizes of (Dx  0.158, 0.130). The largest elements in the T mesh (adapted to the exponential function of (23)) occurs for the element adjacent to x = 1, whereas the largest elements in the C mesh (adapted to the Runge function of (23)) lie at both ends of the domain. Fig. 18 illustrates the results for T at s = 1 and 3 for common refinement, cubic spline, common refinement 2 and linear interpolation schemes (results for linear interpolation are not shown for s = 3 as the method has become unstable). Results for common refinement, common refinement 2 and linear interpolation exhibit considerable error near the peak at x = 0.9. The accuracy shown by the schemes is in the same relative order as the order (degree) of the scheme with the linear common refinement more accurate than linear interpolation. Given previous results, the C mesh in this area is now clearly too coarse. These results reinforce an earlier point stressed in Sections 3 and 4, that it is important to capture the extrema in the source data with the nodes of the target mesh, but also reinforces the dangers inherent in under-resolved calculations. While

Fig. 17. T and C at s = 3 for common refinement, cubic spline, common refinement 2 and linear interpolation for nonuniform meshes of [67, 45] nodes compared to a fine mesh reference solution and the initial condition.

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

3

T

2.5 2

8959

T, Initial condition T (55,37) 1 sec, nonunif., Com ref T (55,37) 1 sec, nonunif., Cubic T (55,37) 1 sec, nonunif., Com ref 2 T (55,37) 1 sec, nonunif., Linear T, fine mesh, 1 sec, no data transfer T (55,37) 3 sec, nonunif., Com ref T (55,37) 3 sec, nonunif., Cubic T (55,37) 3 sec, nonunif., Com ref 2 T, fine mesh, 3 sec, no data transfer

1.5 1 0.5 -1

-0.5

0

0.5

1

x Fig. 18. T at s = 1 and s = 3 using common refinement, cubic spline, common refinement 2 and linear interpolation on nonuniform meshes of [55, 37] nodes compared to the fine mesh reference solution and the initial condition (linear interpolation results at s = 3 are not shown due to instability).

C is under-resolved globally, it is greatly under-resolved near the peak in T. This example shows that coupled multiphysics problems with data transfer, whether on uniform or nonuniform meshes, may produce substantial error near the coarse asymptotic limit, depending on the data transfer strategy employed. The nonuniform mesh problems [55, 37] and [67, 45] are both under-resolved multiphysics problems, where uniform mesh [40, 39] and [30, 29] results show aggregate error results indistinguishable from discretization error (showing negligible data transfer error). These results suggest that, when the coupled multiphysics problem is considered, more mesh nodes may sometimes be necessary to spatially resolve the problem using a nonuniform mesh than are required for a uniform mesh. These results run contrary to one’s intuition that fewer elements would be needed if a nonuniform mesh is employed. In actuality, these results indicate the importance of selecting a proper strategy for mesh refinement. Recall that the C mesh was refined solely using the curvature of the g initial function. Likewise, the T mesh was refined only using the curvature of the f initial function. This refinement strategy completely ignores the need to refine the T mesh where it is necessary to host the coupling (and data transfer) with C, and vice versa. Indeed, one should consider use of a global error indicator that insures that the composite error across the simulation is factored into the refinement applied to each individual element on each instance of the mesh used. While a global error equidistribution strategy such as proposed in [1] is a good first step, there is considerable research opportunity in this area. Figs. 19 and 20 are designed to illustrate a rudimentary global refinement strategy. Here the meshes from the nonuniform [67, 45] problem, Fig. 17, were modified by coarsening the T grid (exponential function, 67 nodes) near x = 0.9 and refining it near x = 0, such that the total number of nodes remains the same. The C grid (Runge function, 45 nodes) was refined near x = 0.9 and coarsened near x = 0, again keeping the total number of nodes the same. The L2 errors for T obtained previously using the old meshes were 2.51 and 6.59 for common refinement and linear interpolation; for the newly adjusted meshes, the L2 errors are 0.87 and 2.18, respectively. The L2 errors using the exponential and Runge function initial conditions (20) are summarized in Table 3. Note that, even though the spatial refinement was reduced at the highly curved areas of the respective functions, the accuracy of the solutions more than double. These results clearly indicate that a significant increase of accuracy is achieved by performing adaptation considering the global multiphysics solution, even though less refinement is used in the highly curved areas of the initial function on its respective grid. Secondly, the use of a global indicator will tend to refine each mesh instance toward the others. If the accuracy needs of the simulation are high, this equidistribution of mesh refinement across the meshes may effectively nullify any reduction in overall DOF achieved by using the multiple mesh approach. In closing, these results indicate that a careful assessment should be performed when comparing a multiphysics approach using a single mesh vs. a multiple mesh approach. While the multiple mesh approach may have a theoretical advantage in the minimization of the DOF required for the simulation, this advantage may not be significant in practice due to the need to perform global refinement to accurately couple the problem. Further, the complexity of implementation of the multiple mesh approach does not work in its favor, either. 6. Summary The present study has involved the application of numerical experimentation to investigate potential problems that may arise in the course of data transfer operations in a multiple mesh multiphysics calculation. While such investigations are valuable to provide information and insight into numerical processes, the extension to other problems will not necessarily

8960

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

1.2 1

Runge fctn, old grid Runge fctn, new grid Exp. fctn, old grid Exp. fctn, new grid

T, C

0.8 0.6 0.4 0.2 0 -0.2 -1

-0.5

0

0.5

1

x Fig. 19. Original and adapted [67, 45] T and C meshes supporting the calculation in Fig. 20. Note that each new grid possesses refinement in highly curved areas of the other.

Fig. 20. Results for T and C using the new adapted grids compared to results for the original [67, 45] grids.

produce the same results. The study discussed a selection of data transfer algorithms and the impact of these algorithms on static functions and on a coupled multiphysics example, the latter employing both uniform and nonuniform meshes. The nonuniform meshes were adapted to the initial conditions on each mesh. Examining the selected data transfer algorithms on a cyclical data transfer problem led to the conclusions that: (a) the selection of a conservative data transfer scheme is not sufficient, by itself, to guarantee acceptable accuracy for a particular problem, (b) it is usually important that the nodes

Table 3 Relative error (21) for exponential and Runge function initial conditions (20). Mesh [40, 39] [30, 29] [45, 33] [67, 45] [67, 45] [67, 67] a

uniform uniform uniform nonuniform old nonuniform new nonuniform

No data transfer employed.

Com ref

Cubic

Linear

0.71 1.28 0.55 2.51 0.87 0.35

0.73 1.31 0.64 0.57

0.94 1.73 3.17 6.59 2.18 0.91

0.42

Com ref 2

Coarse mesha [40, 40] 0.73 [30, 30] 1.31

0.14 0.51 0.05

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

8961

on the target mesh capture the extrema in the source data (by locating nodes at the extrema), and (c) higher-order data transfer methods may be of benefit when considering the accuracy of the data transfer procedure alone. Indeed, the nonconservative cubic-spline interpolation data transfer procedure consistently performed very well across this study; the conservative quadratic common refinement scheme also performed very well. When the coupled multiphysics problem was considered, several additional observations were made. First, the fact that each mesh lies within the range of asymptotic convergence for the equations it is hosting does not necessarily mean that the mesh is fine enough for accurate data transfer operations. Depending on the data transfer algorithm employed, both of these meshes may need to be finer than one would expect. Indeed, if the mesh is too coarse, discretization and data transfer error can combine in a coupled multiphysics problem given the selection of a suboptimal data transfer method to deliver a very inaccurate solution. Secondly, one needs to exercise caution if mesh adaptation is used; application of the chosen mesh adaptation algorithm must not result in excessive coarsening of one mesh locally that leads to insufficient resolution of the physics on the opposite mesh at that location. Further, the use of adaptation does not guarantee that fewer elements will be needed to obtain a resolved simulation in a multiphysics application. One solution to the problem of inadequate resolution of physical features important to a simulation involving data transfer is to employ simultaneous, dynamic adaptive mesh refinement on all of the meshes involved. Further, the error metric used to drive this adaptation should be derived from a composition of the errors of the discretizations on their respective meshes as well as errors introduced by the data transfer process. Such a dynamic adaptation method should reduce the aggregate error by varying the number and spacing of the nodes on each mesh as the calculation proceeds. Several authors are developing advanced approaches for multiphysics spatial integration along these lines. In closing, the present results indicate that data transfer error can be quite significant in multiphysics simulations, and that combining a poorly chosen data transfer strategy with meshes that do not capture the spatial features in a simulation can lead to undesirable results. This study concludes that: 1. Multiphysics application developers and users should carefully contemplate the selection of a multiple mesh architecture. One should not be too quick to abandon the simplicity and potential accuracy of a single mesh implementation unless the advantages in doing so are clear on the class of problems to be solved. 2. Source and target mesh resolutions, capturing extrema in the solution with mesh refinement, solution element order, and data transfer method (particularly the order of the method), all influence the accuracy of the final solution in a multiple mesh application that employs data transfer. There is considerable research opportunity in developing error indicators that quantify these sources of error to drive spatial adaptation approaches that reduce error in both the solution process and in the data transfer between the solutions. Higher-order methods might be particularly promising given that the higher order data transfer approaches (cubic spline and common refinement 2) performed very well in this study. 3. Issues with both common mesh approaches and multiple mesh approaches using low-order data transfer operations for coupled problems leave open the question as to whether there is a better approach. Monolithic hp-FEM multiple mesh methods that do not employ data transfer such as those advanced by Solin et al. [5], Dubcova et al. [6] look quite promising for multiphysics problems with considerable spatial diversity.

Acknowledgments The authors wish to express our gratitude to the three reviewers of our manuscript. We appreciate the reviewers’ suggestions towards improving this work. Further, we wish to explicitly thank the second reviewer for suggestions in improving Section 4 and the example contained therein. The submitted manuscript has been authored by a contractor of the US Government under Contract No. DE-AC0705ID14517 (INL/JOU-08-14009). Accordingly, the US Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for US Government purposes. References [1] Y.-Q. Huang, R. Li, W.-B. Liu, N.-N. Yan, Efficient discretization for finite element approximation of constrained optimal control problems SIAM J. Control Optimization, Submitted for Publication. PDF at: . [2] G. Hansen, S. Owen, Mesh generation technology for nuclear reactor simulation; barriers and opportunities, Nucl. Eng. Design 238 (10) (2008) 2590– 2605. [3] R. Li, On multi-mesh H-adaptive methods, J. Sci. Comput. 24 (3) (2005) 321–341. [4] P. Solin, J. Cerveny, L. Dubcova, I. Dolezel, Multi-mesh hp-FEM for thermally conductive incompressible flow, in: M. Papadrakakis, E. Onate, B. Schrefler (Eds.), ECCOMAS Conference on Coupled Problems, CIMNE, Barcelona, 2007, pp. 547–550. [5] P. Solin, J. Cerveny, L. Dubcova, D. Andrs, Monolithic discretization of linear thermoelasticity problems via adaptive multimesh hp-FEM, J. Comput. Appl. Math. 234 (2010) 2350–2357. [6] L. Dubcova, P. Solin, G. Hansen, H. Park, Comparison of multimesh hp-FEM to interpolation and projection methods for spatial coupling of thermal and neutron diffusion calculations, J. Comput. Phys. 230 (4) (2011) 1182–1197. . [7] Y. Wang, W. Bangerth, J. Ragusa, Three-dimensional h-adaptivity for the multigroup neutron diffusion equations, Prog. Nucl. Energy 51 (2009) 543– 555. [8] Y. Wang, J. Ragusa, Application of hp adaptivity to the multigroup diffusion equations, Nucl. Sci. Eng. 161 (2009) 22–48.

8962

R.W. Johnson et al. / Applied Mathematics and Computation 217 (2011) 8943–8962

[9] L. Demkowicz, C.-W. Kim, 1D hp-adaptive finite element package, Fortran-90 implementation (1Dhp90), Tech. Rep. TICAM Report 99-38, Texas Institute of Computational and Applied Mathematics, University of Texas, Austin, 1999. [10] J.K. Dukowicz, J.W. Kodis, Accurate conservative remapping (rezoning) for arbitrary Lagrangian–Eulerian computations, SIAM J. Sci. Statist. Comput. 8 (3) (1987) 305–321. [11] R.K. Jaiman, X. Jiao, P.H. Geubelle, E. Loth, Assessment of conservative load transfer for fluid-solid interface with non-matching meshes, Int. J. Numer. Methods Eng. 65 (15) (2005) 2014–2038. [12] T.J.R. Hughes, G. Engel, L. Mazzei, M.G. Larson, The continuous Galerkin method is locally conservative, J. Comput. Phys. 163 (2) (2000) 467–488, doi:10.1006/jcph.2000.6577. [13] X. Jiao, M.T. Heath, Common-refinement-based data transfer between non-matching meshes in multiphysics simulations, Int. J. Numer. Methods Eng. 61 (14) (2004) 2402–2427. [14] J.K. Dukowicz, Conservative rezoning (remapping) for general quadrilateral meshes, J. Comput. Phys. 54 (1984) 411–424. [15] J. Grandy, Conservative remapping and region overlays by intersecting arbitrary polyhedra, J. Comput. Phys. 148 (1999) 433–466. [16] D.S. Miller, D.E. Burton, J.S. Oliviera, Efficient second order remapping on arbitrary two dimensional meshes, Tech. Rep. UCRL-ID-123530, Lawrence Livermore National Laboratory, 1996. [17] J.D. Ramshaw, Conservative rezoning algorithm for generalized two-dimensional meshes, J. Comput. Phys. 59 (1985) 193–199. [18] J.D. Ramshaw, Simplified second-order rezoning algorithm for generalized two-dimensional meshes, J. Comput. Phys. 67 (1986) 214–222. [19] D.L. Ropp, J.N. Shadid, C.C. Ober, Studies of the accuracy of time integration methods for reaction-diffusion equations, J. Comput. Phys. 194 (2) (2004) 544–574. [20] I. Prigogine, R. Lefever, Symmetry breaking instabilities in dissipative systems. II, J. Chem. Phys. 48 (4) (1967) 1695–1700. [21] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fourth ed., Butterworth-Heinemann, Oxford, 2000. [22] U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations and Differential–Algebraic Equations, SIAM, Philadelphia, PA, 1998. [23] M. Heroux et al., Trilinos: an object-oriented software framework for the solution of large-scale, complex multi-physics engineering and scientific problems. Available from: , 2008. [24] D. Zwillinger, CRC Standard Mathematical Tables and Formulae, Chapman & Hall/CRC, 2003. [25] R.W. Johnson, Progress on the development of B-spline collocation for the solution of differential model equations: a novel algorithm for adaptive knot insertion, in: C.A. Brebbia, G.M. Carlomagno, P. Anagnostopoulos (Eds.), Computational Methods and Experimental Measurements XI, WIT Press, Southhampton, UK, 2003, pp. 547–556.