CC BY
10
MSC 35J25
DOI: 10.14529/ mmp210407
HIGH ACCURACY NUMERICAL SOLUTION OF ELLIPTIC EQUATIONS WITH DISCONTINUOUS COEFFICIENTS
V.P. Shapeev1'2, V.A. Belyaev1, L.S. Bryndin1'2
1Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, Russian Federation
2Novosibirsk State University, Novosibirsk, Russian Federation
E-mails: vshapeev@gmail.com, va.belyaev@itam.nsc.ru, l.bryndin@g.nsu.ru
We develop an approach to constructing a new high-accuracy hp-version of the least-squares collocation (LSC) method for the numerical solution of boundary value problems for elliptic equations with a coefficient discontinuity on lines of different shapes in a problem solution domain. In order to approximate the equation and the conditions on the discontinuity of its coefficient, it is proposed to use the external parts and irregular cells (i-cells) of the computational grid which are cut off by the line of discontinuity from regular rectangular cells. The proposed approach allows to obtain solutions with a high order of convergence and high accuracy by grid refining and/or increasing the degree of the approximating polynomials both in the case of the Dirichlet conditions on the boundary of the domain and in the case of the presence of Neumann conditions on a large part of the boundary. Also, we consider the case of the problem with a discontinuity of the second derivatives of the desired solution in addition to the coefficient discontinuity at the corner points of the domain. We simulate the heat transfer process in the domain where particles of the medium move in a plane-parallel manner with a phase transition and heat release at the front of the discontinuity line. An effective combination of the LSC method with various methods of accelerating the iterative process is demonstrated: the acceleration algorithm based on Krylov subspaces; the operation of prolongation along the ascending branch of the V-cycle on a multigrid complex; parallelization. The results are compared with those of other authors on solving the considered problems.
Keywords: numerical method; elliptic equations; coefficient discontinuity; conservation law; high accuracy.
Introduction
Numerical simulation based on the solution of equations describing processes of various nature has achieved significant success. In particular, when simulating some of them taking place in physically homogeneous media, a high accuracy of solving the equations describing them is often achieved. However, in reality, many physical media are heterogeneous. They may contain discontinuities in physical parameters that affect the quantitative characteristics of the solution results. Simulation of heat transfer in structures made of dissimilar materials, processes with abrupt phase transitions are common applied problems with a discontinuity in the thermal conductivity coefficient [1,2]. Interest in such problems has been growing recently in connection with the creation of new materials and their research [2]. Thermophysical processes and their effect on the components of composite materials play an essential role in matters of their strength. Examples of solving elliptic equations with a discontinuous coefficient can be found in electrostatics, when simulating multiphase flows, describing solidification processes, and in other applied problems [3].
In the presence of discontinuities in the coefficients of equations in the computational domain, many numerical methods significantly lose their accuracy or do not allow at
all to obtain a convergent solution using the same computational resources as in the absence of such features. In order to achieve increased accuracy in solving problems with discontinuities in physical parameters, the numerical and other approximate methods must have certain properties that differ from those of methods suitable for constructing classical smooth solutions. It is known that one of them is the requirement that the applied method for solving the problem should approximate the conservation laws with sufficient accuracy, the consequences of which are the solvable differential equations [1,4]. For example, the approximation of the conservation law by the finite difference method (FDM) allows to obtain second-order numerical solutions of the heat conduction equation with discontinuous coefficients [1]. This is relatively easy to achieve in the case of a discontinuity front in the form of a straight line passing through the grid nodes. To achieve the same accuracy of numerical solutions to problems in the case of a curvilinear front, the numerical algorithms include additional techniques, for example, "of front straightening". Numerical simulation of processes in physical media in the presence of discontinuities in their parameters, the search for new algorithms in this direction, and the development of existing ones [3,5] remain an urgent task.
This paper is devoted to the development of a new conservative hp-version of the LSC method with the approximation of the integral conservation law, the possibility of refining the grid steps and/or increasing the degree of the approximating polynomials. To approximate the equation and conditions for the discontinuity of its coefficient, which in the case of the heat transfer process take place from the requirement that the heat conservation law is fulfilled, the proposed algorithm uses irregular parts of the cells of the computational grid, which are called i-cells here, cut off by the discontinuity line from regular rectangular cells. At that, we use the possibility available in the LSC method to require that neighboring pieces of an approximate piecewise polynomial solution satisfy simultaneously different matching conditions, which makes it easy to take into account the features of the problem. For simplicity, we consider only the case of a square domain and the solution of several problems when the line of discontinuity of the equation coefficient is straight or curvilinear. A similar idea of using irregular cells was successfully applied in the LSC method for solving boundary value problems in irregular domains with increased accuracy when constructing a solution in cells near the outer boundaries and, in the case of multi-connectivity, inner boundaries of the domain [6,7]. In many cases, this approach avoids the appearance of small and/or elongated cells, which is essential for the numerical solution of ill-conditioned problems [6].
1. Conservation Law for Elliptic Equation
First of all, consider the Dirichlet problem for the elliptic equation in the square domain Q with the boundary 5Q
where k(x1,x2) > 0, f(x^x2), and ub(x1,x2) are given, and u(x1, x2) is the required function. The integral conservation law corresponding to equation (1) has the form
(1) (2)
uL„ = Ub(x i,x2),
du f f
k(x\, x2) — \dS\ = // f(xi,x2)dxidx2, (3)
S V
where S is the boundary of the subdomain V C Q, n is the outer unit normal to S, |dS| is the length of an element of its arc. Suppose that Q contains a continuous line r dividing Q into two subdomains Q1 and Q2, on which coefficient k(x^x2) has a discontinuity with the values k1(x1,x2) in Q1 and k2(x1 , x2) in Q2. In generalized statement of problem (1), (2), it is required to find the function u(x1;x2) that satisfies boundary condition (2) and equation (3) for any subdomain V C Q.
2. Defining Equations of Approximate Problem
The algorithm for solving the problem by the LSC method and the corresponding computer program in the case of a rectangular domain, apart from additional notations, does not fundamentally differ from the case under consideration. We cover the domain Q by a regular grid with NxN square cells of the size 2hx2h with sequential numbering of all cells j = 1, 2,..., N2. In each j-th cell, first of all, for the convenience of implementing the algorithm, we introduce a local coordinate system ^1(x1) = (x1 — x1j)/h, £2(x2) = (x2 — x2j)/h, j = 1,..., N2, where (x1j, x2j) are the coordinates of its center. Here, in each j-th cell, the approximate solution to problem (3), (2) is sought in the form of a linear combination of basic monomials in the linear space of polynomials of the degree K
K K-ii
uij= £ E Ciii2j£11 £22 (4)
il =0 i2=0
with undefined coefficients Cili2,j. To find their values in the LSC method, a local overdetermined system of linear algebraic equations (SLAE) is written for each grid cell, which defines a piece of the approximate solution in the vicinity of the cell center. The system consists of the collocation equations, which represent the requirement for the solution to fulfill the conservation law (3) in the cell, the matching conditions, which ensure the continuity of the solution and equality of flows on the common sides of neighboring cells, and the boundary conditions if the cell side belongs to the boundary of the computational domain.
To derive the collocation equations, we divide each j-th cell into K2 equal square partial cells with the coordinates of the vertices (x1;1, x2,s), (x1;1+1, x2,s), (x1;1+1, x2>s+1), (x1>i ,x2>s+1), where l, s = 0, ...,K — 1. We require that (3) is fulfilled for each partial cell on the desired approximate solution (4). After substituting the representation of approximate solution (4) into (3) we have K2 integral collocation equations
xl,l+l X2,s + l
Fl+ltS+i - FltS+i + Fl+hs+l - Fl+hs = J J f(x1,x2)dx1dx2, (5)
Xl,l X2,s
?l(xi,i+l) ?2(x2,s+l)
Fl+iyS= J k{xux2)-T^r-dCi,FlyS+i= J k{xux2)-^-
?l(xi,l) ?2(x2,s )
in each j-th cell. As a result of calculations of the integrals along the contour on the left side and the integral over the area of a partial cell on the right side, we obtain linear algebraic equations for the unknowns Ci1i2 j with numerical coefficients. Here, it is often possible to take the integrals analytically, otherwise for calculating two-dimensional integrals on the right-hand side (5), we can use, for example, the direct product of one-dimensional Gauss quadrature formulas. At that, the coefficient k(x1, x2) is taken in (5) with the index of the domain in which the center of the square cell under consideration is located.
To formulate the boundary and matching conditions, we divide each side of each cell into K equal partial segments. In the j-th cell, consider the flow through an arbitrary arc of the line 7 given parametrically by x1(t),x2(t),t G [T1;T2], with the outer unit normal n to it (here, 7 is a part of the straight side of the cell or a part of r)
t2
duhij(xi(t),x2(t))
F=-J k(xux2y^K^^KV,yx?(t)+x^(t)\dtl (6)
ti
If the centers of neighboring cells are in the same subdomain Qj, then K matching conditions of solution pieces on the common sides of neighboring cells ensure the continuity of the linear combination of the required solution and the flow
uh j + F = Uh + F. (7)
Here, uhj and Uh are the limits of the approximate problem solution when approaching the side of the j-th cell from inside and outside, which are written for the middle of a partial segment; F and F are the flows through this partial segment calculated using uhj and Uh, respectively.
If the centers of neighboring cells belong to different subdomains Qj, then the part of the line r that belongs to the union of two neighboring cells is approximately evenly divided into K parts and K ratios are written similar to (7). The values of the coefficient k(x1, x2) under the matching conditions for F and F are chosen depending on the subdomain, which contains the centers of the cells in which the approximate solutions uhj and Uh are constructed, respectively.
If the side of the j-th cell is a segment of 5Q, then, in the midpoints of the K partial segments, we formulate the following relations in the local variables of this cell:
uh(Ci,C2)|5Q = ub(x1 (Ci),x2(C2)). (8)
Combining collocation equations (5), matching conditions (7), and boundary conditions (8) in the j-th cell, if it is a boundary one, we obtain a local overdetermined SLAE. The collection of all local SLAEs constitutes a global SLAE determining an approximate global solution to problem (1), (2). The global SLAE can be solved using direct methods or iterations over subdomains, similar to solving various equations using the LSC method presented in [6, 7]. Here, the latter method was applied using QR decomposition of local matrices when constructing a piecewise polynomial solution. At that, the iterations continue until the following inequality is fulfilled:
max |Cn+j - Cn j2j | < e, 0 < n < K, 0 < 12 < K - ¿1, j = 1,...,N2,
j1 ,j2
where n is the iteration number, e is a small value called pseudo-error and specified manually.
3. Algorithm for Constructing Solution in the Case of Curvilinear Front of Coefficient Discontinuity
Let us consider in more detail the case when the discontinuity line of the coefficient k(xi,x2) is a curve (Fig. 1 a). Crossing the square cells, the curve divides them into two parts, which are called i-cells above. Note that with respect to the subdomain Qi, the i-cells located on its other side in the subdomain Q2 are external ("outside-the-contour"). And vice versa, i-cells located in Q2 are outside-the-contour in relation to Qi.
The square cell that r crosses is said to be the parent. The i-cell containing the center of the parent cell is said to be independent (Fig. 1 b, cells 2-6), otherwise, non-independent. All regular square cells are also considered to be independent and parent (Fig. 1 b, cells 1,7). The independent cells and the parent cells containing them are formally considered to be parts of Qi,i = 1, 2, if the parent's center belongs to Qi. As before, we divide each square cell in Q into K2 equal partial cells and divide each cell side in Q into K equal partial segments. In all square cells, as above, the collocation equations are written (Fig. 1 b). Similarly, the equations take the coefficient k(xi ,x2) with the index of the domain in which the center of the parent cell under consideration is located. The collocation equations written for the parent cell are included in the local SLAE of the independent cell contained in it. The equations are valid only for the solution in the independent cell and for a subdomain in which its center is located. In other words, to construct the solution in the independent cell, formally and in fact, we use the entire parent cell and its part that belongs to another subdomain ("outside-the-contour") and is non-independent in it.
On the common sides of neighboring parent cells belonging to the same subdomain, as before (in the case of a straight line of the coefficient discontinuity), we formulate the matching conditions, which we attach to the local SLAEs of independent cells contained
—B-B— [ ] ) 1 [] ) 1—B-B— V V 1—B-B—| V V —B—|—B— c 1 [] C | [] V ' y
[] 2 V V 3 ; . o A A c n2 ; V V A A C [] < [] V V
[] ) [] ) V V c N o c 4 V v y 5 C [] < [] V V
A A [] ) [] ) —B-B— A A c Hi ; —B-B— o 1 : 6 I —b-a— A A [] 7 [] —B-B—
a b
Fig. 1. An example of the problem solution in a domain with a curve line of the coefficient discontinuity (a). The computational domain at K = 2 (b). Here, the dotted lines show the borders of the partial cells for formulating (5); o denotes the centers of parent cells; x, solid diamond, o, A. A, solid rectangle, and empty rectangle denote the centers of partial line segments for writing (7) between cells located in one subdomain Ql or Q2, 1 and 2, 2 and 3, 3 and 4, 4 and 5, 5 and 6, 6 and 7, respectively; □ are points for formulating (8)
in the parent cells under consideration. For matching in independent cells, the centers of which are located in different subdomains, the part of the line r located in the union of these two cells is used (Fig. 1 b). On the sides of each boundary cell located on the border of the domain, as above, we uniformly place K points to approximate (8) (Fig. 1 b).
Therefore, in all regular and independent cells with respect to the coefficients Ci1i2j, to construct a local solution, we write K2 + 4K equations. In the implemented version, an exception was made for the independent boundary i-cells: in the outside-the-contour parts of their parent cells, points on were not used to formulate the boundary conditions, because the solution does not converge with a high order and does not have a high accuracy, which corresponds to the physics and the essence of this approach.
To construct a solution in a non-independent i-cell, at each specific point, an independent cell is sought in the same subdomain, the center of which is closest to it. For the solution at this point, the one continued from the found independent cell is taken. Hence, a part of the i-cell joins one independent cell, and the other part joins another one.
4. Numerical Experiments
In all examples with a known exact solution for calculating the order of convergence R of the relative error ||Er of the solution in infinite norm, at least 100 control points were uniformly distributed in each cell. In the numerical experiments, to significantly accelerate the process of convergence of iterations and reduce the computation time, we used the combined application [8] of the operation of continuation along the ascending branch of the V-cycle on a multigrid complex in the Fedorenko method [9], Krylov subspaces on the intermediate grids of the complex [10], and parallelization using the OpenMP open standard with subdomain traversal based on the red-black ordering [11].
Example 1. Let the line of discontinuity of the coefficient fc(xi, x2) be the straight line x1 = 0,5. Consider the test problem (1), (2) in Q = [0,1]x[0,1], in which the solution u(x1,x2) and the coefficient k(x1 , x2) are given by the expressions
i 10ex2(ex1 - e0'5), xi < 0,5, f 1, xi < 0,5, (q)
u 1 ex2(exi - e0'5), xi > 0,5, k [10, xi > 0,5, (9)
with the inhomogeneous equation term f (xi, x2) = 20exiex2 — 10e0'5ex2. In Dirichlet condition (2), the values of the exact solution (9) on Ш were taken as ub(x1 ,x2).
The results of numerical experiments on a sequence of grids are given in Table 1. We can see a significant superiority in accuracy and order of convergence of the approximate solution in the cases of using polynomials of high degrees for its representation in comparison with the cases of polynomials of low degrees. Here we present the number of iterations Niter required for convergence for different e, and the time in seconds spent on the iteration process. The calculations were carried out on the computer Intel Core i5-8265U CPU 1.6 GHz, 4 Cores, DIMM DDR4-2400 1200 MHz 8 Gb. Also, we present the ratios of the number of iterations AFj and time AFt in the case of solving the problem without using their acceleration according to Krylov, the operation of continuation according to Fedorenko, and parallelization to the corresponding values in the case of their combined application. Note that the error values hardly differ in both cases (with and without acceleration), therefore, they are given here once. It is found that with an increase in the degree of the polynomials K and a grid refinement, the combined application of various
Table 1
Results of numerical experiments in Example 1
K = 2 (e= 10"10)
NxN 11 Er 11 oo R without acceleration with acceleration AFj AFt
Niter time, s Niter time, s
20x20 2,88E-5 2,91 390 1,828 83 0,172 4,70 10,62
40x40 3,78E-6 2,92 1491 24,39 161 1,5 9,26 16,26
80x80 5,35E-7 2,82 5597 360,047 319 10,453 17,55 34,41
160x160 l,18E-7 2,18 20834 5784,22 601 88,265 34,67 65,53
K = 4 (e= 10"14)
NxN 11 11 qo R without acceleration with acceleration AFj AFt
Niter time, s Niter time, s
10x10 2,53E-8 — 286 3,172 90 0,453 3,17 7,00
20x20 8,15E-10 4,96 1107 48,344 108 1,703 10,25 28,38
40x40 2,58E-11 4,98 4250 720,609 202 12,625 21,03 57,07
80x80 8,13E-13 4,99 16267 13086,7 414 124,297 39,29 105,28
K = 6 (e= 10"14)
NxN 11 Er 11 oo R without acceleration with acceleration AFj AFt
Niter time, s Niter time, s
2x2 l,16E-7 — 85 0,203 41 0,078 2,07 2,60
4x4 l,21E-9 6,58 135 1,297 42 0,125 3,21 10,37
8x8 1,05E-11 6,85 282 10,172 42 0,531 6,71 19,16
16x16 8,89E-14 6,88 1067 152,016 51 2,453 20,92 61,97
acceleration methods becomes more efficient. Table 1 shows that, at K = 4, the number of iterations can be reduced by almost 40 times and the calculation time can be reduced by more than 100 times.
Remark 1. It is known that if the Dirichlet conditions are replaced by the Neumann conditions on a part of the boundary in a differential problem, then the problem becomes worse conditioned. It is found that when solving the boundary value problem from Example 1 with the Neumann conditions on some parts of the solution accuracy obtained by the LSC method decreases. However, the use of polynomials of high degrees (K = 4, 6) allows to achieve an accuracy of the order of 10-9-10-12 on grids of a fairly moderate size (NxN < 40x40). The study of the values of the conditionality numbers v of local SLAEs in the spectral norm showed: 1) in the inner cells, depending on K, the values of v belong to the interval (10,102); 2) larger values of K correspond to larger values of v; 3) the conditionality of the SLAE of the boundary cells slightly differs from the conditionality of the SLAE of the cells inside the domain in the case of the problem with the Dirichlet conditions only; in the case of the Neumann conditions, the values of v in the boundary cells belong to the interval (102,104). At that, the larger the part of the boundary with the Neumann conditions, the slower the iterative process converges, and at large K it diverges.
Remark 2. When solving practical problems, their exact solutions are often not known. In addition, in such cases, additional difficulties associated with the presence of features often arise. Within the framework of this study, an example of solving problem (1), (2) was considered, in which, simultaneously with the coefficient discontinuity in the domain, there
was a discontinuity of the second derivative of the solution at the corner points of Q. The results of numerical experiments showed the presence of the second order of convergence in the entire domain with a good accuracy at K = 2, 4, 6. The restriction of the order of convergence by the order of the solution smoothness was also observed when solving the Poisson equation with a similar feature by the FDM [12] and the LSC method [13]. It was also established here that the solution accuracy and the order of its convergence (for K > 2), calculated outside the neighborhood of the corner points, were significantly higher than in the entire domain Q. For this, we observed the convergence of the norm of the solutions difference on two successive grids.
Example 2. Consider an example of a typical heat transfer process in a domain in which particles of the medium move plane-parallel with a constant velocity v along the axis x1 in the presence of a discontinuity in the thermal conductivity coefficient with a phase transition and heat release (absorption) at the front of the discontinuity line. A part of such a problem occurs in more general problems of heat and mass transfer, for example, when welding metals, melting ice, permafrost, etc. In many cases, the model of such a problem in the domain Q with the boundary 5Q and the discontinuity line of the thermal conductivity coefficient is determined by the equation
ki
d2 u d2 u dxf dx 2
+ VCi рг
du dx i
fi (xi,X2 ), (xi ,X2) G П,
(10)
with Dirichlet boundary conditions (2). Let the straight line x\ = 0,5 be the line of discontinuity of the coefficient k(xi,x2) in Q = [0,1]x[0,1], with different physical parameters, where the constants ki, ci, pi and the functions fi(x1 ,x2) be the thermal conductivity coefficient, specific heat, density of the medium, and distributed heat sources in the i-th subdomains, i = 1, 2. Suppose that the particles of the medium from the subdomain Q1 with the first phase of the medium move through the discontinuity line into the subdomain Q2 with the second phase of the medium. According to the heat conservation law in the simulated process on the discontinuity line, the following relation must be fulfilled:
k2
du dxi
- k
du dx i
pi KV,
where k is the latent specific heat of the phase transition, i k2
du 'ôxi
and ( k
du dxi
(11)
are
the heat fluxes calculated from different sides of the discontinuity line.
Slight differences of this example from the previous one are that when writing the collocation equations, the conservation law for equation (10), similar to (3), is approximated, and when writing the matching conditions, (11) is taken into account.
We check the capabilities of the LSC method by solving a test problem. Table 2 shows the results of a numerical experiment for solving the problem when the discontinuity line is the straight line x\ = 0, 5. The test solution has the form u = 2048(xi — 1)2x6+64(x2 — 1)3x2 and u = 64(x1 — 1)2x1(0,5 + x1 )3 + 64(x2 — 1)3x2 in Q1 and Q2, respectively. As boundary conditions and heat sources f in (10), we took the corresponding functions calculated from the test solution. The physical parameters were k1 = 10, c1 = 1, p1 = 2, k2 = 1, c2 = 1, p2 = 1, k = 632. As in the previous case, it can be seen from the analysis of the numerical results that the approximate problem solution converges with a higher order.
2
i
2
i
Table 2
Results of numerical experiments in Example 2
NxN K = 2 K = 4 K = 6
11 Fr 11 oo R 11 Fr 11 oo R 11 Fr 11 oo R
10x10 l,96E-2 — 3,19E-5 — 4,49E-8 —
20x20 5,02E-3 1,96 l,92E-6 4,06 7,17E-10 5,98
40x40 l,25E-3 2,01 l,18E-7 4,02 1,10E-11 6,03
80x80 3,14E-4 1,99 7,37E-9 4,00 2,21E-12 2,32
Example 3. Consider problem (1), (2) in Q = [0,1] x [0,1] with the test solution u(xi, x2) (Fig. 1 a) and the functions k(x1 ,x2) and f (x1 ,x2), which are given by the formulas
= f 10 (er - e0'8), r< 0,8, = f 1, r< 0,8, U = \ er - e0'8, r> 0,8, k [10, r> 0,8,
f(xi,x2) = 10er(l + r)/r, where r = ^J(X\ + 0,1)2 + {x2 + 0,1)2. In this example, un is divided on the curved line r, but the flow F (6) through r is continuous.
The results of numerical experiments are given in Table 3, from which it follows that when using polynomials of high degrees to approximate the problem, as in the previous examples, its solution achieves high accuracy already on coarse grids.
Table 3
Results of numerical experiments in Example 3
NxN K = 2 K = 4 K = 6
11 Fr 11 oo R 11 Fr 11 oo R 11 Fr 11 oo R
10x10 l,46E-3 — 4,92E-5 — 2,69E-6 —
20x20 2,48E-4 2,55 2,97E-6 4,05 6,18E-8 5,44
40x40 4,92E-5 2,33 l,35E-7 4,46 9,08E-10 6,09
80x80 l,20E-5 2,04 5,28E-9 4,68 9,95E-12 6,51
Example 4. Depending on the problem statement, jumps may occur not only in the functions un, but also in kun and u. In the LSC method, such features are approximated with good accuracy by specifying matching conditions for neighboring pieces of a piecewise polynomial solution to the problem. Note that from the perspective of software implementation, such modifications of the algorithm are very simple. Let us explore the capabilities of the LSC method when solving several problems considered in [3,5]. Here, Q1 = Q\Q2, where Q2 is the inner circle.
Example 4.1. Consider (1), (2) in Q = [-1,1]x[-1,1] with k1 = 2, k2 = 1 (Fig. 2 a) [3]:
u
r4 - 0,54 - 0,1 log(2r) 0,52 , , ^
k1
r2
k
2f (X1,X2) E Q2,
16r2, (x1,x2) E Q1, 4, (x1,X2) E Q2,
a b
Fig. 2. Solution plots in Example 4.1 (a) and Example 4.3 (b)
r = \J x\ + x2, r is the circle x2 + x\ = 0,52. Here, un and kun suffer discontinuity on T. Example 4.2. Consider (1), (2) in ft = [0,1]x[0,1] with ki = 1, k2 = 2 [5]:
f 0, (xi,x2) e fti, f i 0, (xi ,x2) e fti,
u = \ e-x?-x2, (xi, X2) e ft2, f I 8(xi + x2 - 1)e-x1-x2, (xi, X2) e ft2,
r is the circle (xi — 0,5)2 + (x2 — 0,5)2 = 0,252. Here, un, kun and u suffer discontinuity on r.
Example 4.3. Consider (1), (2) in ft = [0,1]x[0,1] with ki = 0,02, k2 = 1 (Fig. 2 b) [5]:
= \ e Xl x2, (xi, x2) e fti, f = \ 4(x2 + x2 — 1)e Xl x2, (xi, x2) e fti,
U [ ex2+x2, (xi, x2) e ft2, f [ 8(x2 + x2 + 1)ex?+x2, (xi, x2) e ft2,
r is the circle (xi — 0,5)2 + (x2 — 0,5)2 = 0,252. Here, un, kun and u suffer discontinuity on r. Note that k2/ki ^ 1, i.e. there is a small parameter at the highest derivative.
In [3,5], the authors managed to achieve the second order of the solution convergence using the FDM. On 320x320 grid, the accuracy of 2,98E-5 was achieved in Example 4.1, the accuracy was about 10 -6-10 -7 in Example 4.2, and finally, in Example 4.3, the accuracy was about 10 -5-10 -6. Table 4 shows the high-precision results of numerical experiments obtained using the LSC method at K = 4. It can be seen that R is no worse than the third. Moreover, in Example 4.3, more iterations were required for convergence compared to the other two. Therefore, the use of polynomials of high degrees in the LSC method often allows achieving high accuracy of the solution to the problems with singularities, while the scheme of its construction is quite simple. However, it is necessary to take into account the type of features and other factors that affect the convergence of solutions. Examples 4.1, 4.2, and 4.3 deal with more complex cases than Examples 1-3 in this paper. It is found that when polynomials of the degree K = 6 are used, divergence of the iterative process can be observed. In addition, calculations using high degree polynomials behave less robustly than those using lower degree polynomials. Nevertheless, they converge on coarse grids. In particular, in Example 4.2, the accuracy of 2,14E-9 was achieved on a 20x20 grid at K = 6.
Table 4
Results of numerical experiments in Example 4 at K = 4
NxN Example 4.1 Example 4.2 Example 4.3
11 Fr 11 oo R 11 Fr 11 oo R 11 Fr 11 oo R
10x10 3,87E-5 — 7,67E-6 — l,85E-4 —
20x20 l,16E-5 1,74 3,08E-7 4,64 2,29E-5 3,01
40x40 7,81E-7 3,89 2,13E-8 3,85 l,00E-6 4,52
80x80 2,22E-8 5,14 l,48E-9 3,85 l,23E-7 3.02
Conclusion
We proposed the high-precision algorithm for solving elliptic equations with a coefficient discontinuity on lines of various shapes in the computational domain. The algorithm allows obtaining the solution with a high convergence order. On moderately coarse grids, in many cases it is possible to construct the solution practically with accuracy of the rounding error. The efficiency of the proposed approach was tested by solving various problems. It can be seen that the scheme of the numerical algorithm is the same, which confirms the universality of the LSC method. In addition, the use of the LSC method with a combination of different methods of accelerating iterations makes it possible to reduce their number by almost 40 times and reduce the calculation time by more than 100 times.
Acknowledgements. The research was carried out within the state assignment of Ministry of Science and Higher Education of the Russian Federation (project no. 121030500137-5 and AAAA-A19-119051590004-5).
References
1. Samarskii A.A., Vabishchevich P.N. Vychislitel'naya teploperedacha [Computational Heat Transfer]. Moscow, Editorial URSS, 2003. (in Russian)
2. Isaev V.I., Cherepanov A.N., Shapeev V.P. Numerical Study of Heat Modes of Laser Welding of Dissimilar Metals with an Intermediate Insert. International Journal of Heat and Mass Transfer, 2016, vol. 99, pp. 711-720. DOI: 10.1016/j.ijheatmasstransfer.2016.04.019
3. Zhilin L. A Fast Iterative Algorithm for Elliptic Interface Problems. SIAM Journal on Numerical Analysis, 1998, vol. 35, no. 1, pp. 230-254. DOI: 10.1137/S0036142995291329
4. Godunov S.K., Zabrodin A.V., Ivanov M.I., Kraiko A.N., Prokopov G.P. Chislennoe reshenie mnogomernyh zadach gazovoy dinamiki [Numerical Solution of Multidimensional Problems of Gas Dynamics]. Moscow, Nauka, 1976. (in Russian)
5. Tzou C., Stechmann S.N. Simple Second-Order Finite Differences for Elliptic PDEs with Discontinuous Coefficients and Interfaces. Communications in Applied Mathematics and Computational Science, 2019, vol. 14, no. 2., pp. 121-147. DOI: 10.2140/camcos.2019.14.121
6. Belyaev V.A., Shapeev V.P. The Versions of Collocation and Least Residuals Method for Solving Problems of Mathematical Physics in the Trapezoidal Domains. Computational Technologies, 2017, vol. 22, no. 4, pp. 22-42. (in Russian)
7. Shapeev V.P., Golushko S.K., Bryndin L.S., Belyaev V.A. The Least Squares Collocation Method for the Biharmonic Equation in Irregular and Multiply-Connected Domains. Journal of Physics, 2019, vol. 1268, article ID: 012076. DOI: 10.1088/1742-6596/1268/1/012076
8. Vorozhtsov E.V., Shapeev V.P. On the Efficiency of Combining Different Methods for Acceleration of Iterations at the Solution of PDEs by the Method of Collocations and Least Residuals. Applied Mathematics and Computation, 2019, vol. 363, pp. 1-19. DOI: 10.1016/j.amc.2019.124644
9. Fedorenko R.P. Vvedenie v vychislitel'nuyu fiziku [Introduction to Computational Physics]. Moscow, Moscow Institute of Physics and Technology, 1994. (in Russian)
10. Saad Y. Numerical Methods for Large Eigenvalue Problems. Philadelphia, Society for Industrial and Applied Mathematics, 2011.
11. Degi D.V., Starchenko A.V. Numerical Solution of Navier-Stokes Equations on Computers with Parallel Architecture. Tomsk State University Journal of Mathematics and Mechanics, 2012, no. 2(18), pp. 88-98. (in Russian)
12. Shapeev V.P., Shapeev A.V. Solutions of the Elliptic Problems with Singularities Using Finite Difference Schemes with High Order of Approximation. Computational Technologies, 2006, vol. 11, no. 2, pp. 84-91. (in Russian)
13. Shapeev V.P. Bryndin L.S., Belyaev V.A. Solving Elliptic Equations in Polygonal Domains by the Least Squares Collocation Method. Bulletin of the South Ural State University. Series Mathematical Modelling, Programming and Computer Software, 2019, vol. 12, no. 3. pp. 140-152. (in Russian)
Received May 26, 2021
УДК 519.632.4 БЭТ: 10.14529/mmp210407
ЧИСЛЕННОЕ РЕШЕНИЕ С ПОВЫШЕННОЙ ТОЧНОСТЬЮ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ С РАЗРЫВНЫМИ КОЭФФИЦИЕНТАМИ
В.П. Шапеее1'2, В.А. Беляев1, Л.С. Брындин1'2
Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, г. Новосибирск, Российская Федерация
2Новосибирский государственный университет, г. Новосибирск, Российская Федерация
Разработан подход построения нового hp-варианта метода коллокации и наименьших квадратов (КНК) численного решения с повышенной точностью краевых задач для эллиптических уравнений с разрывом коэффициента на линиях различных форм в области решения задачи. Для аппроксимации уравнения и условий на разрыве его коэффициента в алгоритме предложено использовать законтурные части и нерегулярные ячейки (н-ячейки) расчетной сетки, отсеченные линией разрыва от регулярных
прямоугольных ячеек. Предложенный подход позволил получить решения с повышенным порядком сходимости и высокой точности при измельчении шага сетки и/или увеличении степени аппроксимирующих полиномов как в случае условий Дирихле на границе области, так и в случае наличия условий Неймана на значительной части границы. Рассмотрен также случай задачи, когда кроме разрыва коэффициента в угловых точках области имеется разрыв вторых производных искомого решения. Проведено моделирование процесса теплопереноса в области, в которой частицы среды перемещаются плоскопараллельно с фазовым переходом и выделением тепла на фронте линии разрыва. Продемонстрировано эффективное сочетание метода КНК с различными способами ускорения итерационного процесса: алгоритм ускорения, основанный на подпространствах Крылова; операция продолжения вдоль восходящей ветви V-цикла на многосеточном комплексе; распараллеливание. Проведено сравнение с результатами других авторов по решению рассмотренных задач.
Ключевые слова: численный метод; эллиптические уравнения; разрыв коэффициента; закон сохранения; повышенная точность.
Литература
1. Самарский, А.А. Вычислительная теплопередача / А.А. Самарский, П.Н. Вабищевич. -М.: Едиториал УРСС, 2003.
2. Isaev, V.I. Numerical Study of Heat Modes of Laser Welding of Dissimilar Metals with an Intermediate Insert / V.I. Isaev, A.N. Cherepanov, V.P. Shapeev // International Journal of Heat and Mass Transfer. - 2016. - V. 99. - P. 711-720.
3. Zhilin, Li. A Fast Iterative Algorithm for Elliptic Interface Problems / Li Zhilin // SIAM Journal on Numerical Analysis. - 1998. - V. 35, № 1. - P. 230-254.
4. Годунов, С.К. Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, М.Я. Иванов, А.Н. Крайко, Г.П. Прокопов. - М.: Наука, 1976.
5. Tzou, C. Simple Second-Order Finite Differences for Elliptic PDEs with Discontinuous Coefficients and Interfaces / C. Tzou, S.N. Stechmann // Communications in Applied Mathematics and Computational Science. - 2019. - V. 14, № 2. - P. 121-147.
6. Беляев, В.А. Варианты метода коллокации и наименьших невязок для решения задач математической физики в трапециевидных областях / В.А. Беляев, В.П. Шапеев // Вычислительные технологии. - 2017. - Т. 22, № 4. - С. 22-42.
7. Shapeev, V.P. The Least Squares Collocation Method for the Biharmonic Equation in Irregular and Multiply-Connected Domains / V.P. Shapeev, S.K. Golushko, L.S. Bryndin, V.A. Belyaev // Journal of Physics. - 2019. - V. 1268, article ID: 012076.
8. Vorozhtsov, E.V. On the Efficiency of Combining Different Methods for Acceleration of Iterations at the Solution of PDEs by the Method of Collocations and Least Residuals / E.V. Vorozhtsov, V.P. Shapeev // Applied Mathematics and Computation. - 2019. -V. 363, № 124. - P. 1-19.
9. Федоренко, Р.П. Введение в вычислительную физику / Р.П. Федоренко. - М.: Московский физико-технический институт, 1994.
10. Saad, Y. Numerical Methods for Large Eigenvalue Problems / Y. Saad. - Philadelphia: Society for Industrial and Applied Mathematics, 2011.
11. Деги, Д.В. Численное решение уравнений Навье-Стокса на компьютерах с параллельной архитектурой / Д.В. Деги, А.В. Старченко // Вестник Томского государственного университета. Математика и механика. - 2012. - № 2 (18). - С. 88-98.
12. Шапеев, В.П. Решение эллиптических задач с особенностями по схемам высокого порядка аппроксимации / В.П. Шапеев, А.В. Шапеев // Вычислительные технологии. -2006. - Т. 11, № 2. - С. 84-91.
13. Шапеев, В.П. Решение эллиптических уравнений в полигональных областях методом коллокации и наименьших квадратов / В.П. Шапеев, Л.С. Брындин, В.А. Беляев // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. -2019. - Т. 12, № 3. - С. 140-152.
Василий Павлович Шапеев, доктор физико-математических наук, главный научный сотрудник, Институт теоретической и прикладной механики им. С.А. Христиа-новича СО РАН (г. Новосибирск, Российская Федерация); профессор, кафедра математического моделирования, Новосибирский государственный университет (г. Новосибирск, Российская Федерация), vshapeev@gmail.com.
Василий Алексеевич Беляев, аспирант, младший научный сотрудник, Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН (г. Новосибирск, Российская Федерация), va.belyaev@itam.nsc.ru.
Лука Сергеевич Брындин, старший лаборант с высшим профессиональным образованием, Институт теоретической и прикладной механики им. С.А. Христи-ановича СО РАН (г. Новосибирск, Российская Федерация); аспирант механико-математического факультета, кафедра математического моделирования, Новосибирский государственный университет (г. Новосибирск, Российская Федерация), l.bryndin@g.nsu.ru.
Поступила в редакцию 26 мая 2021 г.
