In the next subsections an approach is explained, which reduces this factor extremely, and is labeled ’EOS-quadtree’ in Table 3.1.
Table 3.1:Comparison of evaluation time (µs) of perfect gas and realistic EOS
Perfect Gas 0.01712
Vapor Two-Phase Liquid direct EOS 232.015 2392.149 260.016 EOS-quadtree 0.32402 0.36802 0.30002
poly-nomial representation of the thermodynamic variables in the unit element is built by
qhN(ξ1, ξ2) =
N
X
i,j=0
ˆ
qijψij(ξ1, ξ2), ψij=`Ni (ξ1)`Nj (ξ2). (3.1) This is again a nodal interpolation ansatz (see Section 2.3) with a tensor prod-uct of 1-D Lagrange polynomials `N of degreeN as basis functions. The nodal valueqˆij is evaluated with the help of the EOS library and two input values
ˆ
qij=q(aelemroot +ξ1i∆la, belemroot +ξ2j∆lb). (3.2) Nis again the polynomial degree. The root of each element is the lower left corner which is denoted by the minimum values foraandbin this element.
The used nodal Lagrange basis functionsψ(ξ1, ξ2)are build with Chebyshev-Gauß nodes. These are defined in the interval[−1,1]and a mapping to the desired unit space is given by
ξ1= ξ1CG+ 1
2 and ξ2= ξ2CG + 1
2 . (3.3)
The polynomial degreeNcan vary for each element. The maximum polyno-mial degree of an element can be choosen from 1 to 10 at the beginning of the building process. This chosen degree is only the upper limit. For each polynomial degree from 1 to the chosen upper limit theL∞-error is calculated for every element
L∞= max
qh(ξ1, ξ2)−q(ξ1, ξ2) q(ξ1, ξ2)
< t, ξ1, ξ2∈[0,1]. (3.4) The lowest polynomial degree, which satisfies theL∞-error, is used in the corresponding element. By adapting the polynomial degree for each ele-ment the quadtree storage size is reduced. The infinity norm error-criterion (L∞) is checked on 20×20 equidistant points in each element, by using the interpolation (qh) (see Equation (3.1)) and evaluating the EOS directly (q) at these points. The equidistant points also include the element bound-aries, which ensures that the jump of the solution over element interfaces is
lower thant. In case Equation (3.4) is not satisfied by any polynomial de-gree for an element, this element is going to be refined by dividing it in2×2 new quadratic elements. All four new elements are in the next level (l+ 1) compared to the refined element, which is in level l. This procedure is re-peated for each element until the error is lower thantor the maximum level is reached. With this technique not every level contains 2l ×2l elements, instead the amount depends on t and the chosen upper limit for the poly-nomial degree. In this implementationLq = 32is the maximum level be-cause the element localization, e.g., finding the root of an element, is done with a 64-bit binary number. Figure 3.1 shows the 1st and 2nd level of a
00 01
10 11
00 00 00 01
00 10 00 11 01 00 01 01
01 10 01 11
10 00 10 01
10 10 10 11 11 00 11 01
11 10 11 11
Figure 3.1:Level 1 and 2 with identification numbers
table with the corresponding identification (ID) binary numberB. This fig-ure illustrates an example where all elements in level 1 needed refinement.
The numberB is read from left to right. The first two digits correspond to level one. The second pair of numbers belongs to the second level and so on.
This approach is similar to the Morton numbering [40]. Algorithm 2 shows a
way to calculate the root values for each element by using the ID numberB.
aelemroot =amin; belemroot =bmin; forl←1toledo
∆la= amax2−al min;
∆lb= bmax2−blmin;
aelemroot =aroot+B[2l−1]∆la; belemroot =broot+B[2l]∆lb; end
Algorithm 2:Finding the root for an identification numberN In this algorithmle is the level of the element. B[i] denotes the bit at the ithposition and the value ofB[i]is either 0 or 1. A schematic example of a quadtree is shown in Figure 3.2 with the connections to the previous levels.
The top plane represents a levellwith four equidistant elements. The middle
Top view Top view Top view
Figure 3.2:Schematic quadtree with level connections
plane shows levell+ 1. For this schematic example only three elements need to be refined from levell to level l+ 1and one satisfies the error criterion stated by Eq. (3.4). The bottom plane shows levell + 2and displays the refined elements from the previous level. The top view is visible in the lower
right of the picture. From level l+ 1to level l+ 2four elements do need refinement in this example.
Figure 3.3 shows a quadtree for water build with the above described method for(a, b)→q = (ρ, e)→T, withρ∈ [ρmin = 1, ρmax = 1330]kg/m3and e ∈ [emin = −9.9, emax = 4.056·103]kJ/kg.The building error was set to t = 1·10−7 andLq = 17for this quadtree. The black dot in Figure 3.3 is situated in the two-phase region, which is enveloped by the saturation lines and shows one root of an element given by following ID numberB
B = 00
|{z}
l=1
11
|{z}
2
10
|{z}
3
11
|{z}
4
01
|{z}
5
→ρelemroot = 582.44 kg
m3 , eelemroot = 1387.76 kJ kg. (3.5) Because of very high gradients at the liquid and vapor saturation line the algo-rithm has to refine along these lines to reach the desiredL∞-error. To reduce the amount of elements due to refinement at these lines, a cut-cell approach is implemented and explained in Section 3.4, which was also used to build Figure 3.3. Every element of such a quadtree is saved during the build process even if it does not satisfy Equation (3.4). This is done to be able to load the quadtree with a higher error than for the build process to save memory during runtime. The memory usage is lower for a higher error since a lower amount of elements has to be loaded.
To evaluate the quadtree during computation the first step is to load the saved quadtree into memory. For the loading process the error can be set higher than the build error l ≥ t. Starting at level one each element in the first level is checked if the desired error is reached by the polynomial ap-proximation inside the element. For each element where the desired error is not reached the next level is examined and the connection to the upper level is built. This is repeated until the maximum level of the table is reached or all elements satisfy the errorl. After the table is loaded into memory a fast quadtree bisection method shown in Algorithm 3 is used to find the corre-sponding element for a given input set (ag, bg). In this algorithm INT() symbolizes a integer devision, (aquadtreemid , bquadtreemid )is the middle point of the quadtree. In the first level the element is set to the corresponding i, j ele-ment. The row ’elem→elem%connectionij’ changes the current element to the element one level below the current element by using thei, jconnection.
The point(aupper elemmid , bupper elemmid )is the middle point of the element one level
ρ e
1 333.25 665.5 997.75 1330
499.38 582.44 623.97
−9.9 4056
2023.05
1006.58
1514.82 1260.7 1387.76
3039.53
Figure 3.3:A visualization of the quadtree for water (Lq=17).
ρ∈[1,1330]kg/m3ande∈[−9.9,4.056×103]kJ/kg
up, which is connected to the current element. An example for Algorithm 3 is given in Table 3.2 with the given valuesag = ρg = 582.5kg/m3and bg = eg = 1338kJ/kg. The element ID corresponding to the given values isB = 00 11 10 11 01. This element can have two states, the first is that it contains a approximated solution to the desired variableqand the second state is that this element does not contain a approximated solution. If the first state is true than the approximated solution is evaluated at(ξ1g, ξ2g)with Equation (3.1), whereξ1g andξ2g are given by the transformation from physical space to reference space
ξ1g = ag−aelemroot
∆la , ξ2g = bg−belemroot
∆lb . (3.6)
i=INT a
g
aquadtreemid
; j=INT b
g
bquadtreemid
; elem→elemij;
whileelem has connectionsdo i=INT a
g
aupper elemmid
; j=INT b
g
bupper elemmid
;
elem→elem%connectionij; end
Algorithm 3:Finding the element for a given input set(ag, bg)
To accelerate the evaluation of Equation (3.1) the following property is used ψ(ξ) =V−Tφ(ξ) or ψi=
N
X
n=0
Vin−Tφn(ξ) (3.7) whereφ(x) = (1, x, x2, . . . , xN)T is the monomial basis andV the so called Vandermonde matrix
Vmn=φn(ξm). (3.8)
Rewriting Equation (3.1) the result is:
qhN ξ1g, ξ2g
=
N
X
i,j=0
ˆ qij
N
X
n=0
Vin−Tφn(ξ1g)
N
X
m=0
Vmj−Tφm(ξ2g). (3.9) The termsφn(ξ1g)andφm(ξ2g)can be precomputed. The evaluation of Equa-tion (3.9) is by a factor of 1.3 faster compared to the evaluaEqua-tion of EquaEqua-tion (3.1) with the Lagrange basis function. This is due the higher computational cost of evaluating the Lagrange basis (see Equation (2.32)). The time to cal-culate the needed variables temperature, pressure, speed of sound, viscosity and heat conductivity from density and specific internal energy with Equation (3.9) for the quadtree approach can be seen in Table 3.1 compared to the per-fect gas equation and the EOS library. If the second state of the element is
Table 3.2:Example for finding the corresponding element for ag=ρg= 582.5kg/m3andbg=eg= 1338kJ/kg
Level Middle value integer division Element ID
1 ρquadtreemid = 665.5kg/m3 i= 0
equadtreemid = 2023.05kJ/kg j= 0 00 2 ρ00mid= 333.25kg/m3 i= 1 e00mid = 1006.58kJ/kg j= 1 11
3 ρ00 11mid = 499.38kg/m3 i= 1
e00 11mid = 1514.82kJ/kg j= 0 10
4 ρ00 11 10mid = 582.44kg/m3 i= 1
e00 11 10mid = 1260.7kJ/kg j= 1 11 5 ρ00 11 10 11
mid = 623.97kg/m3 i= 0 e00 11 10 11 01
mid = 1387.76kJ/kg j= 1
true, which means that it does not contain a approximate solution, then the calculation stops and a better resolved table has to be used. This second state occurs if the maximum level during loading the quadtree is reached and the approximated solution of the element does not satisfy the loading-errorl.