Abstract
An efficient parallel multiscale numerical algorithm is proposed for a parabolic equation with rapidly oscillating coefficients representing heat conduction in composite material with periodic configuration. Instead of following the classical multiscale asymptotic expansion method, the Fourier transform in time is first applied to obtain a set of complex-valued elliptic problems in frequency domain. The multiscale asymptotic analysis is presented and multiscale asymptotic solutions are obtained in frequency domain which can be solved in parallel essentially without data communications. The inverse Fourier transform will then recover the approximate solution in time domain. Convergence result is established. Finally, a novel parallel multiscale FEM algorithm is proposed and some numerical examples are reported.
1. Introduction
Consider the parabolic equation with rapidly oscillating coefficients as follows: where is a bounded smooth domain or a bounded polygonal convex domain with a periodic microstructure and denotes the time domain. The unknown is the temperature increment . and are some given functions which denote heat source and initial temperature, respectively. is a small parameter and it represents the relative size of a periodic cell.
Let and the heat conductivity coefficient is assumed to satisfy the following.) is -periodic in , , where is the reference cell.() is symmetric and satisfies the uniform elliptic condition; that is, where is an arbitrary vector in , and are some constants independent of .()Let and . Suppose that the boundaries for some and , for some ; see [1, 2].
In (1), the rapid spatial oscillations in the coefficients translate the periodic structure of the body which comes from assembling -scaled versions of the reference cell . This body has to be thought of as being made of a composite material.
An efficient analysis approach to such a problem induces a need for homogenization. The mathematical theory has been presented in a lot of works; see, for example, [3–6]. Several efficient numerical methods have been proposed and analyzed; for instance, see [7–9]. From these literatures, we found that the homogenized equation happens to be a system of the same type with constant coefficients. So the normal step-by-step time-marching algorithms such as the backward-Euler, Crank-Nicolson methods, or any higher order methods need to be applied to solve the homogenized equation step by step in time. Although these algorithms are effective in solving many practical problems, one of major drawbacks is that they are not easily parallelizable along the time axis, since they require the information of the solutions at previous time steps in order to advance to the next time step.
In this work, in order to propose a parallelizable multiscale numerical algorithm, we apply the Fourier transform and inverse Fourier transform to present the multiscale asymptotic analysis for (1). Numerical methods for time-dependent problems based on the use of Fourier transform have been considered in [10, 11], providing us with the starting point for using Fourier transform in multiscale asymptotic analysis. In [10, 11], Douglas et al. introduced an efficient parallel method for solving wave equations in the space-time domain after taking the Fourier transforms. The parallel method had been further developed and extended to viscoelasticity, parabolic problems, and Navier-Stokes equations; see [12, 13]. Later Sheen et al. also applied Laplace transform to solve parabolic problems; see [14, 15]. In addition, the early study of homogenization by employing means of integral transformation is [16], where Fourier and Laplace transforms were applied to prove the weak convergence theorems of the homogenization procedure for the time-dependent equations.
The main procedures in this paper are now briefly described. First, we take the Fourier transform with respect to the time variable in the parabolic equation and then present the multiscale asymptotic expansions of the solutions for the resulting equations in space-frequency domain at specific frequency points of interest, which are independent and may therefore be done in parallel. Finally, the approximate solution of parabolic equation in space-time domain is retrieved by the discrete inverse Fourier transform. It should be stated that our purpose of using Fourier transform is not just to prove the homogenization convergence theorems as done in [16], but to propose a parallelizable multiscale numerical algorithm. The new contributions presented in this study are the following. The strong convergence theorem of the second-order multiscale asymptotic solution in space-frequency domain is derived where we carefully estimate the bounds associated with frequency variable as it has a singularity for both zero and infinity in inverse Fourier transform. The convergence result for the approximate solution in time domain is obtained. A novel parallel multiscale FEM algorithm and numerical examples in 3D composite materials are given.
The remainder of this paper is organized as follows. In Section 2, the original parabolic equation is transformed to a set of complex-valued elliptic problems. The multiscale asymptotic expansions in space-frequency domain are then discussed. In Section 3, the multiscale truncation error estimates in space-frequency domain are derived under some assumptions. The approximate solution in space-time domain is inverted back by discrete inverse Fourier transform and the convergence result is proved in Section 4. Finally, a multiscale algorithm is proposed and numerical simulations are carried out to validate the presented method.
Throughout the paper, let denote a positive constant independent of and . All complex-valued functions are assumed to have values in the complex field and standard notations for function spaces and their associated norms will be used in this paper.
2. Fourier Transform and Multiscale Asymptotic Expansions
Denote the Fourier transform of a function by , where
In order to take the Fourier transform of (1), we first extend and by zero when and then we transfer (1) into the following complex-valued elliptic problem in space-frequency domain, where is some fixed parameter in frequency domain.
It is not difficult to check the existence and uniqueness of the problem (4) with an application of the Lax-Milgram lemma; see [17].
Proposition 1. Let , under the assumptions ()-(); then, for any , , (4) has a unique weak solution .
The method requires solving a finite set of problem (4) at specific frequency point . It is natural to treat them simultaneously because they are independent. Moreover, observe that (4) is an equation with rapidly oscillating coefficients in spatial domain. The direct accurate numerical computation of the solution is difficult especially in 3D problems because it would require huge computational cost. We will apply the method of multiscale asymptotic expansion. Let the derivative with respect to be where denotes the macroscopic variable and describes the microscopic one. Two scales describe the model: gives the position of a point in domain while gives its position in the reference cell . On expanding in asymptotic expansion of the form where is -periodic in , the space-frequency (4) becomes
On equating orders at , with the well-known theorem of the existence and uniqueness of the solutions in quotient space, we deduce that
Equating terms at , we have with the method of separation of variables and we can obtain where are the solutions of the following problems:
At , we have On averaging (12) over , we get the homogenized equation where the homogenized coefficients are defined by
Remark 2. Under the assumptions ()-() of the coefficients for the parabolic equation, one can prove that the homogenized coefficients are also symmetric and uniform elliptic for some constants ; see, for example, [3, 5]. Therefore the existence and uniqueness of the solution for the homogenized (13) can be established.
Subtracting (13) from (12), we have and using the method of separation of variables, we can obtain where are the solutions of the following problems:
Remark 3. Note that , defined in (11) and (17) can be obtained by resolving the following cell problems in : and then extending their solutions by periodicity (see, e.g., [5]), so they are called cell functions as usual.
Therefore, we can obtain the first-order and second-order multiscale asymptotic solution for the space-frequency (4) as
Remark 4. It is worthwhile to note that all the cell problems are real valued, while the homogenized equation is complex valued. Therefore, when we do numerical computation, we need to repeatedly solve the homogenized problem at a set of specific frequency points but solve all the cell problems only once.
3. Multiscale Truncation Error Estimate in Frequency Domain
In order to derive the truncation error estimate for multiscale asymptotic solutions, we first discuss the regularity for the homogenized solution . The key point is to derive precise estimates for the regularity bounds associated with since the behavior of these bounds as approaches zero and infinity is of critical importance in obtaining the multiscale truncation error estimate and when inverting the multiscale asymptotic solution back to time domain by the inverse Fourier transform.
Proposition 5. Let be the weak solution of the homogenized equation (13) whose coefficients are constants defined in (14). Then,(i) for any , , it holds where the constant depending only on and the coefficients;(ii) for any , and any subdomain , it holds where the constant is depending only on , and the coefficients. Moreover, if is , then we have the global estimate (22) on .
Proof. Define the conjugate bilinear form :
where denotes the complex conjugate of . Then the variational formulation of the homogenized equation (13) is to find such that
Set in (24) and take the imaginary part; one obtains
which yields
Let in (24) and consider the real part; with the coercivity of and the Poincare’s inequality, we have
and thus
together with (26); we have
Next we investigate the interior regularity of the homogenized solution. Fix any open set , and choose an open set such that . Then select a cutoff function satisfying
Rewrite (24) as
and set , where denotes the difference quotient,
By a standard regularity technique (see, e.g., [18], pages 310–313), we have
Cauchy’s inequality with yields
and use the properties of the difference quotient for ,
and choose that yields the bound
Combining (33) and (36) implied
In view of the properties of the difference quotient, with (29) and (28), we deduce , with the estimate
For higher regularity of , insert into (24), , where , and perform integration by parts; we have
In view of (38), applying the interior regularity result to , then
which yields
Finally we can obtain (22) by induction on . Given , following the standard boundary regularity technique we get the global estimate
Next we can derive the truncation error estimate for multiscale asymptotic solutions (20) in frequency domain.
Theorem 6. Let be the weak solution of (4), and let be the multiscale asymptotic solutions given in (20). Under the assumptions ()(), with , , , there exists a constant independent of and such that
Proof. For , set ; one can check that satisfy
where
Under the assumptions , it follows from Theorem 1.1 of [2] (also see [1]) that
Together with the regularity estimate of , one can obtain
and therefore, we get and
In order to get the boundary estimate, introducing the following cutoff function such that
and setting , , with the regularity estimate of , we can obtain
Using the following estimate on boundary layer , (see [6])
we have
By the trace theorem, we know
Similar to Proposition 5 for (44) with nonhomogeneous boundary condition, we can drive the following estimate:
together with (48) and (53) and we have
where is independent of and .
For , by means of Theorem 1.2 in [6], we can also derive the estimate (43). Therefore the proof is complete.
Remark 7. It is worth noting that, for any , it is well known that homogenized solution if , and . If does not have the local regularity near the boundary , follow the way in [19]; we can use the asymptotic expansion in the interior subdomain and construct the boundary layer near .
4. Inverse Fourier Transform
In this section, we will construct the approximate solution of the parabolic equation (1) in space-time domain by the inverse Fourier transform. Recall the formula of the inverse Fourier transform as provided that is a real function.
Fix a sufficiently large such that is negligible for . Set , where is a selected segment number for the interval and the space-time solution of (1) can be approximated by where , .
In order to get an estimate for , let and we get
For , from Proposition 5, we can check
For the estimate of , the following lemma was proved in [17].
Lemma 8. Let and suppose that for , for all . Then
For , by Theorem 6, we have
Combining (60), (62), and (63), we have the following theorem for the full error estimate.
Theorem 9. Let be the solution of (1), and let be the approximate solution given in (57). Under the assumptions ()(), , , , and suppose that for , defined in (61) belong to for all ; there exists a constant independent of and such that
5. Parallel Multiscale FEM Algorithm and Numerical Tests
In simulation, we will apply the following first-order and second-order difference quotients of some functions given by where , the set of elements with the node ; , is the number of elements of ; is the number of nodes on the element ; are the nodes of the element ; are the Lagrangian shape functions; denotes the FEM solution of ; is the value of at the node associated with the element .
The following parallel multiscale FEM algorithm are presented to solve (1).
Step 1. Solve the first-order cell problems of defined in (18) for and the second-order cell problems of defined in (19) for in local fine grid of the periodic cell by FEM.
Step 2. Compute the homogenized coefficients defined in (14) in the local fine grid of the periodic cell by numerical integration.
Step 3. Fix a sufficiently large and set , and simultaneously solve the complex-valued homogenized equation (13) in the whole domain in a coarse mesh at , . Then compute the first-order and second-order partial derivatives of the homogenized solutions by high order difference quotients in the same coarse mesh.
Step 4. Extend all the cell functions by periodicity from the local fine grid of the periodic cell to the gird of the whole domain . Then compute the multiscale asymptotic solutions given by (20) at frequency points , .
Step 5. Calculate the approximation solutions for any fixed time for the original equation (1) according to the formula (57).
Remark 10. Notice that cell problems are independent of the time variable or the frequency variable , so they only need to be solved once.
Remark 11. It should be stated that Step 1 to Step 5 in the algorithm are serial. However, for the time-dependent problems, the main computation of the existing multiscale algorithms is in computing the time-dependent homogenized equation step by step. In our presented algorithm, we can simultaneously solve the homogenized equations at different frequent points, so the main computation of the algorithm is parallel. Also, notice that each step can be done simultaneously. In this respect, different from other multiscale algorithms for the time-dependent equation, we called it parallel algorithm.
Remark 12. If some conditions of the geometric symmetry are satisfied (see [19]), then we can employ the homogeneous Dirichlet boundary conditions on the boundary for the cell problems (18) and (19). Following the way in [19], one can also obtain a similar convergence result to that of Theorem 6.
Now we discuss the finite element computation for the complex-valued homogenized Equation (13) in . Let and , where the subscript denotes the real part of a function and the subscript denotes the imaginary part.
The variational form of (13) is to find that satisfied for any , and is given by (23). Let be the finite element space, the discrete variational form of (66) is to find that satisfied for any . Let be a set of basis of , and set and therefore, we solve the linear system in dimension as follows: where
To validate the proposed method and to confirm the theoretical analysis reported in this paper, we present numerical simulations for the following cases.
Example. Consider the parabolic equation (1) in 3D composite materials, where whole domain and the periodic cell are illustrated in Figures 1(a) and 1(b). In the configuration, we take and the blue cube denotes the inclusion and other parts are the matrix. The parameters of the composite materials for two cases will be considered, where . Let and .
(a)
(b)
To assess the validity of the presented method, the exact solution of the original parabolic equation (1) must be available. It is extremely difficult and even impossible to find the exact solution; we replace it by the numerical solution fully resolved by FEM in space and FDM in time using fine meshes (FD-FE method). It should be noted that more considerable computing time will be required to solve (1) using fully resolved FD-FE method when is small, so we only present numerical results for . Here, the linear tetrahedral elements are employed for the semidiscrete problem of (1) using fine meshes, and the solution is computed by the Euler midpoint scheme.
Without confusion, let denote the numerical solution for the original problem (1) using the fully resolved FD-FE method. Let , denote the numerical solutions inverted back by the first-order and second-order multiscale asymptotic solutions in accordance with (57), respectively. Also, let denote the numerical solution inverted back by the homogenized solutions according to the following formula:
All the codes are implemented on parallel hierarchical grid (PHG), which is a toolbox for developing parallel adaptive finite element programs. PHG is currently under active development at State Key Laboratory of Scientific and Engineering Computing of CAS. We use linear tetrahedron elements to solve the related problems numerically.
The numerical results in the cross-section and at time for Case 1 and Case 2 are displayed in Figures 2 and 3, respectively. Although the first-order and second-order multiscale asymptotic solutions have the same convergence rate as stated in Theorem 9, numerical results displayed in Figure 2 confirm that the second-order corrector terms are crucial in some cases.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The relative errors vary with the segment number in the numerical integration formula (57) for homogenized method (HM), first-order multiscale method (1st MsM), and second-order multiscale method (2nd MsM) corresponding to , , and that are listed in Figure 4. We observe that the relative errors become stable after . Also, the relative errors evolution over time is listed in Figure 5. The results show that the accuracy of the proposed method especially the 2nd MsM is in good agreement with the fully resolved case. It is worthwhile to note that the relative errors do not grow rapidly as time increases.
(a)
(b)
(a)
(b)
The costs required by the developed multiscale finite element method are not sensitive to the time variable but lie on the number of space-frequency equations which is determined by the segment number , and each space-frequency equation is solved by multiscale finite element method with the computational cost listed in Table 1; for the fully resolved case, we need to compute steps and the cost for each time step is also listed in Table 1. This clearly demonstrates that the present algorithm is very effective and tremendous saving in computing time is achieved.
Conflict of Interests
The authors declare that they have no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by National Natural Science Foundation of China (Grant no. 11301329) and the Grant of the First-class Discipline of Universities in Shanghai.