1 Introduction

The alternating direction solver [1, 2] has been recently applied for numerical simulations of non-linear flow in heterogeneous media using the explicit dynamics [3, 4].

The problem of extraction of liquid fossil fuels with hydraulic fracturing technique has been considered there. During the simulation two (contradictory) goals i.e., the maximization of the fuel extraction and the minimization of the ground water contamination have been considered [4, 14]. The numerical simulations considered there are performed using the explicit dynamics with B-spline basis functions from isogeometric analysis [5] for approximation of the solution [6, 7]. The resulting computational cost of a single time step is linear, however the number of time steps is large due to the Courant-Fredrichs-Lewy (CFL) condition [8]. In other words, the number of time steps grows along with the mesh dimensions.

Our ultimate goal is to extend our simulator for implicit dynamics case, following the idea of the implicit dynamics isogeometric solver proposed in [9]. The problem is that the extension is possible only if the permeability coefficients of the elliptic operator are expressed as the tensor product structure. Thus, we focus on the algorithm approximating the permeability coefficients with tensor products iteratively.

The algorithm is designed to be a preconditioner for the implicit dynamics solver. With such the preconditioner the number of time steps of the non-stationary problem can be reduced, while the numerical accuracy preserved.

Our method presented in this paper is an alternative for other methods available for approximating coefficients of the model, e.g., adaptive cross approximation [15].

2 Explicit and Implicit Dynamics Simulations

Following the model of the non-linear flow in heterogeneous media presented in [1] we start with our explicit dynamics formulation of the problem of non-linear flow in heterogeneous media where we seek for the pressure scalar fieldu:

$$\begin{aligned} \begin{aligned} \bigg (\frac{\partial u(x,y,z)}{\partial t},\upsilon (x,y,z)\bigg ) = \bigg ( \Big ( K(x,y,z)e^{\mu u(x,y,z)} \Big ) \nabla u(x,y,z), \nabla \upsilon (x,y,z) \Big ) \\ +\, \big ( f(x,y,z), \upsilon (x,y,z) \big ) \quad \forall \upsilon \in V \end{aligned} \end{aligned}$$
(1)

Here \(\mu \) stands for the dynamic permeability constant, K(x, y, z) is a given permeability map, and f(x, y, z) represents sinks and sources of the pressure, modeling pumps and sinks during the exploration process.

The model of non-linear flow in heterogeneous media is called exponential model [12] and is taken from [10, 11].

In the model, the permeability consists of two parts, i.e., the static one depending on the terrain properties, and the dynamic one reflecting the influence of the actual pressure.

The broad range of the variable known as the saturated hydraulic conductivity along with the functional forms presented above, confirm the nonlinear behavior of the process.

The number of time steps of the resulting explicit dynamics simulations are bounded by the CFL condition [8], requesting to reduce the time step size when increasing the mesh size. This is important limitation of the method, and can be overcome by deriving the implicit dynamics solver.

Following the idea of the implicit dynamics solvers presented in [9], we move the operator to the left-hand side:

$$\begin{aligned} \begin{aligned} \bigg (\frac{\partial u}{\partial t},\upsilon \bigg ) - \bigg ( \Big ( K(x,y,z)e^{\mu u(x,y,z)} \Big ) \nabla u, \nabla \upsilon \Big ) = \big ( f, \upsilon \big ) \quad \forall \upsilon \in V, \end{aligned} \end{aligned}$$
(2)

where we skip all arguments but the permeability operator.

In order to proceed with the alternating directions solver, the operator on the left-hand-side needs to be expressed as a tensor product:

$$\begin{aligned} \begin{aligned} \bigg (\frac{\partial u}{\partial t},\upsilon \bigg ) - \bigg ( \Big ( K(x)e^{\mu \overline{u}(x)} K(y)e^{\mu \overline{u}(y)} K(z)e^{\mu \overline{u}(z)} \Big ) \nabla u, \nabla \upsilon \Big ) = \\ \big ( f, \upsilon \big ) + \Big ( K(x)K(y)K(z) e^{\mu \overline{u}(x)} e^{\mu \overline{u}(y)} e^{\mu \overline{u}(z)} - K(x,y,z) e^{\mu u(x,y,z)} \nabla u , \nabla \upsilon \Big ) \quad \forall \upsilon \in V \end{aligned} \end{aligned}$$
(3)

It is possible if we express the static permeability in a tensor product form:

$$\begin{aligned} K(x,y,z) = K(x)K(y)K(z) \end{aligned}$$
(4)

using our tensor product approximation algorithm described in Sect. 3.

Additionally, we need to replace the dynamic permeability with an arbitrary selected tensor product representation:

$$\begin{aligned} u(x,y,z) = \overline{u}(x)\overline{u}(y)\overline{u}(z) \end{aligned}$$
(5)

It can be done by adding and subtracting from the left and the right hand sides the selected tensor product representation.

One simple way to do that is to compute the average values of u along particular cross-sections, namely using:

$$\begin{aligned} u(x,y,z) = \sum _{i=1}^{N_x} \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( d_{ijk}B_{i,p}(x)B_{j,p}(y)B_{k,p}(z) \Big ) \Big ) \Big ) \end{aligned}$$
(6)

so we define:

$$\begin{aligned} \overline{u}(x) = \sum _{i=1}^{N_x}\overline{u}_iB_{i,p}(x) \end{aligned}$$
(7)
$$\begin{aligned} \overline{u}(y) = \sum _{j=1}^{N_y}\overline{u}_jB_{j,p}(y) \end{aligned}$$
(8)
$$\begin{aligned} \overline{u}(z) = \sum _{k=1}^{N_z}\overline{u}_kB_{k,p}(z) \end{aligned}$$
(9)

and

$$\begin{aligned} \overline{u}_i = \frac{ \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z}(d_{ijk}) \Big ) }{N_yN_z}; \quad \overline{u}_j = \frac{ \sum _{i=1}^{N_x} \Big ( \sum _{k=1}^{N_x}(d_{ijk}) \Big ) }{N_xN_z} ; \quad \overline{u}_k = \frac{ \sum _{i=1}^{N_x} \Big ( \sum _{j=1}^{N_y}(d_{ijk}) \Big ) }{N_xN_y} \end{aligned}$$
(10)

In other words, we approximate the static permeability and we replace the dynamic permeability.

Finally we introduce the time steps, so we deal with the dynamic permeability explicitly, and with the static permeability implicitly:

$$\begin{aligned} \begin{aligned} \bigg (u_{t+1},\upsilon \bigg ) - \bigg ( \Big ( K(x)e^{\mu \overline{u}_t(x)} K(y)e^{\mu \overline{u}_t(y)} K(z)e^{\mu \overline{u}_t(z)} \Big ) \nabla u_{t+1}, \nabla \upsilon \Big ) = \\ \big ( f, \upsilon \big ) + \Big ( K(x)K(y)K(z) e^{\mu \overline{u}(x)} e^{\mu \overline{u}(y)} e^{\mu \overline{u}(z)} - K(x,y,z) e^{\mu u_t(xyz)} \nabla u_t , \nabla \upsilon \Big ) \quad \forall \upsilon \in V \end{aligned} \end{aligned}$$
(11)

In the following part of the paper the algorithm for expression of an arbitrary material data function as the tensor product of one dimensional functions that can be utilized in the implicit dynamics simulator is presented.

3 Kronecker Product Approximation

As an input of our algorithm we take a scalar function defined over the cube shape three-dimensional domain. We call this function a bitmap, since often the material data is given in a form of a discrete 3D bitmap.

First, we approximate this bitmap with B-spline basis functions using fast, linear computational cost isogeometric L2 projections algorithm.

$$\begin{aligned} Bitmap(x,y,z) \approx \sum _{i=1}^{N_x} \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( d_{ijk}B_{i,p}(x)B_{j,p}(y)B_{k,p}(z) \Big ) \Big ) \Big ) \end{aligned}$$
(12)

Now, our computational problem can be stated as follows:

Problem 1

We seek coefficients \(a_1^x, \ldots , a_{N_x}^x\),\(b_1^y, \ldots , b_{N_y}^y\), \(c_1^z, \ldots , c_{N_z}^z\) to get the minimum of

$$\begin{aligned} \begin{aligned} F(a_1^x, \ldots ,a_{N_x}^x,b_1^y, \ldots , b_{N_y}^y, c_1^z, \ldots , c_{N_z}^z) \\ =\displaystyle \int _{\varOmega }\Big [ \Big ( \sum _{i=1}^{N_x}a_iB_{i,p}^x \Big ) \Big ( \sum _{j=1}^{N_y}b_jB_{j,p}^y \Big ) \Big ( \sum _{k=1}^{N_z}c_kB_{k,p}^z - \sum _{i=1}^{N_x} \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \big ( d_{ijk}B_{i,p}(x)B_{j,p}(y) B_{k,p}(z) \big ) \Big ) \Big ) \Big ) \Big ]^2 \\ =\displaystyle \int _{\varOmega } \big [ \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} \big ( a_ib_jc_k-d_{ijk}B_{i,p}(x)B_{j,p}(y) B_{k,p}(z) \big ) \big ) \big ) \big ]^2 \end{aligned} \end{aligned}$$
(13)

The minimum is realized when the partial derivatives are equal to zero:

$$\begin{aligned} \frac{\partial F}{\partial a_l^x}(a_1^x, \ldots , a_{N_x}^x,b_1^y, \ldots , b_{N_y}^y, c_1^z, \ldots , c_{N_z}^z)=0 \end{aligned}$$
(14)
$$\begin{aligned} \frac{\partial F}{\partial b_l^y}(a_1^x, \ldots , a_{N_x}^x,b_1^y, \ldots , b_{N_y}^y, c_1^z, \ldots , c_{N_z}^z)=0 \end{aligned}$$
(15)
$$\begin{aligned} \frac{\partial F}{\partial c_l^z}(a_1^x, \ldots , a_{N_x}^x,b_1^y, \ldots , b_{N_y}^y, c_1^z, \ldots , c_{N_z}^z)=0 \end{aligned}$$
(16)

We compute these partial derivatives:

$$\begin{aligned} \begin{aligned} \frac{\partial F}{\partial a_l^x}(a_1^x, \ldots , a_{N_x}^x,b_1^y, \ldots , b_{N_y}^y, c_1^z, \ldots , c_{N_z}^z)=0 \\ =\displaystyle \int _{\varOmega } \big [ \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} \big ( 2(a_lb_jc_k-d_{ljk} \big ) \Big ( \frac{\partial (a_ib_jc_k)}{\partial a_l^x} - \frac{\partial (d_{ijk})}{\partial a_l^x} \Big ) B_{l,p}^xB_{j,p}^y B_{k,p}^z) \big ) \big ]=0, \end{aligned} \end{aligned}$$
(17)

where the internal term:

$$\begin{aligned} \frac{\partial (a_ib_jc_k)}{\partial a_l^x} = \frac{\partial (a_i)b_jc_k}{\partial a_l^x} + a_i \frac{\partial (b_jc_k)}{\partial a_l^x} = b_jc_k\delta _{il}+0, \end{aligned}$$
(18)

thus

$$\begin{aligned} =\displaystyle \int _{\varOmega } \big [ \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} \big ( 2(a_lb_jc_k-d_{ljk} \big ) b_jc_k B_{l,p}^xB_{j,p}^y B_{k,p}^z \big ) \big ]=0, \quad l=1,\ldots ,N_x \end{aligned}$$
(19)

Similarly we proceed with the rest of partial derivatives to obtain:

$$\begin{aligned} =\displaystyle \int _{\varOmega } \big [ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} \big ( 2(a_ib_lc_k-d_{ilk} \big ) a_ic_k B_{i,p}^xB_{l,p}^y B_{k,p}^z \big ) \big ]=0, \quad l=1,\ldots ,N_y \end{aligned}$$
(20)
$$\begin{aligned} =\displaystyle \int _{\varOmega } \big [ \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} \big ( 2(a_ib_jc_l-d_{ijl} \big ) a_ib_j B_{i,p}^xB_{j,p}^y B_{l,p}^z \big ) \big ]=0, \quad l=1,\ldots ,N_z \end{aligned}$$
(21)

This is equivalent to the following system of equations:

$$\begin{aligned} \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} 2 \big ( a_lb_jc_k-d_{ljk} \big )b_jc_k \big ) =0 \end{aligned}$$
(22)
$$\begin{aligned} \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} 2 \big ( a_ib_lc_k-d_{ilk} \big )a_ic_k \big ) =0 \end{aligned}$$
(23)
$$\begin{aligned} \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} 2 \big ( a_ib_jc_c-d_{ijl} \big )a_ib_j \big ) =0 \end{aligned}$$
(24)

We have just got a non-linear system of \(N_x+N_y+N_z\) equations with \(N_x+N_y+N_z\) unknowns:

$$\begin{aligned} a_l\Big ( \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} \big ( b_jc_k \big ) b_jc_k \big ) \Big ) = \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \big ( d_{ljk} b_jc_k \big ) \Big ) \end{aligned}$$
(25)
$$\begin{aligned} b_l\Big ( \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} \big ( a_ic_k \big ) a_ic_k \big ) \Big ) = \sum _{i=1}^{N_x} \Big ( \sum _{k=1}^{N_z} \big ( d_{ilk} a_ic_k \big ) \Big ) \end{aligned}$$
(26)
$$\begin{aligned} c_l\Big ( \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} \big ( a_ib_j \big ) a_ib_j \big ) \Big ) = \sum _{i=1}^{N_x} \Big ( \sum _{j=1}^{N_y}\big ( d_{ijl} a_ib_j \big ) \Big ), \end{aligned}$$
(27)

what implies:

$$\begin{aligned} a_l = \frac{ \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} d_{ljk} b_jc_k \big )}{ \sum _{j=1}^{N_y} \big ( \sum _{k=1}^{N_z} \big ( b_jc_k \big )^2 \big ) } \end{aligned}$$
(28)
$$\begin{aligned} b_l = \frac{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} d_{ilk} a_ic_k \big )}{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} \big ( a_ic_k \big )^2 \big ) } \end{aligned}$$
(29)

We insert these coefficients into the third equation:

$$\begin{aligned} \begin{aligned} c_l\sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \Big ( \frac{ \sum _{m=1}^{N_y} \bigg ( \sum _{n=1}^{N_z}d_{imn}b_mc_n \bigg ) }{ \sum _{m=1}^{N_y} \bigg ( \sum _{n=1}^{N_z}(b_mc_n)^2 \bigg ) } \Big )^2 \Big ( \frac{ \sum _{m=1}^{N_x} \bigg ( \sum _{n=1}^{N_z}d_{mjn}a_mc_n \bigg ) }{ \sum _{m=1}^{N_x} \bigg ( \sum _{n=1}^{N_z}(a_mc_n)^2 \bigg ) } \Big )^2 \Bigg ) \\ = \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y}d_{ijl} \frac{ \sum _{m=1}^{N_y} \bigg ( \sum _{n=1}^{N_z}d_{imn}b_mc_n \bigg ) }{ \sum _{m=1}^{N_y} \bigg ( \sum _{n=1}^{N_z}(b_mc_n)^2 \bigg ) } \frac{ \sum _{m=1}^{N_x} \bigg ( \sum _{n=1}^{N_z}d_{mjn}a_mc_n \bigg ) }{ \sum _{m=1}^{N_x} \bigg ( \sum _{n=1}^{N_z}(a_mc_n)^2 \bigg ) } \Bigg ) \end{aligned} \end{aligned}$$
(30)
$$\begin{aligned} \begin{aligned} c_l\sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \bigg ( \sum _{m=1}^{N_y} \Big ( \sum _{n=1}^{N_z}d_{imn}b_mc_n \Big ) \bigg ) \bigg ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z}d_{mjn}a_mc_n \Big ) \bigg ) \Bigg ) \\ = \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y}d_{ijl} \Bigg ) \Bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y} \big ( b_mc_n \big )^2 \Big ) \Bigg ) \Bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \big ( a_mc_n \big )^2 \Big ) \Bigg ) \end{aligned} \end{aligned}$$
(31)
$$\begin{aligned} \begin{aligned} c_l\sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y}d_{imn}b_mc_n \Big ) \bigg ) \bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_x}d_{mjn}a_mc_n \Big ) \bigg ) \Bigg ) \\ = \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y}d_{ijl} \Bigg ) \Bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y} \big ( b_mc_n \big )^2 \Big ) \Bigg ) \Bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \big ( a_mc_n \big )^2 \Big ) \Bigg ) \end{aligned} \end{aligned}$$
(32)
Fig. 1.
figure 1

The original configuration of static permeability

Full size image
Fig. 2.
figure 2

The result obtained from the heuristic algorithm (a) and from the heuristic plus genetic algorithms (b).

Full size image
Fig. 3.
figure 3

The tensor product approximation after one (a) and five (b) iterations of Algorithm 1.

Full size image
$$\begin{aligned} \begin{aligned} c_l\sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y}d_{imn}b_mc_n \Big ) \bigg ) \bigg ( \sum _{m=1}^{N_x}d_{mjn}a_mc_n \bigg ) \Bigg ) \\ = \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y}d_{ijl} \Bigg ) \Bigg ( \sum _{m=1}^{N_y} \big ( b_mc_n \big )^2 \Bigg ) \Bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \big ( a_mc_n \big )^2 \Big ) \Bigg ) \end{aligned} \end{aligned}$$
(33)
$$\begin{aligned} \begin{aligned} \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y} \big ( \sum _{o=1}^{N_x} d_{ojn}a_oc_nd_{imn}b_mc_nc_l \big ) \Big ) \bigg ) \Bigg ) \\ = \sum _{i=1}^{N_x} \Bigg ( \sum _{j=1}^{N_y} \bigg ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_y} \big ( \sum _{o=1}^{N_x} (a_oc_nb_mc_n )^2d_{ijl} \big ) \Big ) \bigg ) \Bigg ) \end{aligned} \end{aligned}$$
(34)

The above is true when

$$\begin{aligned} \begin{aligned} d_{imn}b_mc_nc_ld_{ojn}a_oc_n= (a_oc_nb_mc_n )^2d_{ijl}, \end{aligned} \end{aligned}$$
(35)

so:

$$\begin{aligned} \begin{aligned} d_{imn}c_ld_{ojn}= a_oc_nb_mc_nd_{ijl} \end{aligned} \end{aligned}$$
(36)

thus:

$$\begin{aligned} \begin{aligned} \frac{d_{ojn}d_{imn}}{d_{ijl}}= \frac{a_oc_nb_mc_n}{c_l} \end{aligned} \end{aligned}$$
(37)
Fig. 4.
figure 4

The tensor product approximation after ten (a) and fifty (b) iterations of Algorithm 1.

Full size image
Fig. 5.
figure 5

The error of the tensor product approximation after one (a), and five (b) iterations of Algorithm 1.

Full size image
Fig. 6.
figure 6

The error of the tensor product approximation after ten (a), and fifty (b) iterations of Algorithm 1.

Full size image

We can setup now \(a_1\), \(b_1\), and \(c_1\) arbitrary and compute \(c_l\) using the derived proportions.

In a similar way we compute \(a_l\), namely we insert:

$$\begin{aligned} b_l = \frac{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} d_{ilk} a_ic_k \big )}{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} \big ( a_ic_k \big )^2 \big ) } \end{aligned}$$
(38)
$$\begin{aligned} c_l = \frac{ \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} d_{ijl} a_ib_j \big )}{ \sum _{i=1}^{N_x} \big ( \sum _{j=1}^{N_y} \big ( a_ib_j \big )^2 \big ) } \end{aligned}$$
(39)

into

$$\begin{aligned} \begin{aligned} a_l \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z} (d_{mjn}a_mc_n) \Big ) \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_y} (d_{mnk}a_mb_n) \Big ) \Big ) \Big ) \Big ) \Big ) \\ = \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z}d_{ljk} \Big ) \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z}(a_mc_n)^2 \Big ) \Big ) \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_y}(a_mb_n)^2 \Big ) \Big ) \Big ), \end{aligned} \end{aligned}$$
(40)

then:

$$\begin{aligned} \begin{aligned} a_l \Big ( \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z}d_{mjn}a_mc_n \Big ) \Big ( \sum _{o=1}^{N_y}d_{mok}a_mb_o \Big ) \Big ) \Big ) \Big ) \Big ) \\ = \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z}d_{ljk} \Big ) \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z}(a_mc_n)^2 \Big ) \Big ) \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{o=1}^{N_y}(a_mb_o)^2 \Big ) \Big ), \end{aligned} \end{aligned}$$
(41)

and finally:

$$\begin{aligned} \begin{aligned} \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z} \Big ( \sum _{o=1}^{N_y}a_ld_{mok}a_mb_od_{mjn}a_mc_n \Big ) \Big ) \Big ) \Big ) \Big ) \\ = \sum _{j=1}^{N_y} \Big ( \sum _{k=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z} \Big ( \sum _{m=1}^{N_x} \Big ( \sum _{n=1}^{N_z}(a_mb_oa_mc_n)^2d_{ljk} \Big ) \Big ) \Big ) \Big ) \Big ), \end{aligned} \end{aligned}$$
(42)

what results in:

$$\begin{aligned} \begin{aligned} a_ld_{mok}a_mb_od_{mjn}a_mc_n= (a_mb_oa_mc_n )^2d_{ljk}, \end{aligned} \end{aligned}$$
(43)

so:

$$\begin{aligned} \begin{aligned} a_ld_{mok}d_{mjn}= a_mb_oa_mc_n d_{ljk}, \end{aligned} \end{aligned}$$
(44)

thus:

$$\begin{aligned} \begin{aligned} \frac{d_{mok}d_{mjn}}{d_{ljk}}= \frac{a_mb_oa_mc_n}{a_l} \end{aligned} \end{aligned}$$
(45)

We compute \(b_l\) from (we already have \(a_i\) and \(c_k\)):

$$\begin{aligned} b_l = \frac{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} d_{ilk} a_ic_k \big )}{ \sum _{i=1}^{N_x} \big ( \sum _{k=1}^{N_z} \big ( a_ic_k \big )^2 \big ) } \end{aligned}$$
(46)

The just analyzed Problem 1 has multiple solutions, and the algorithm presented above finds one exemplary solution, for the assumed values of \(a_1\), \(b_1\), and \(c_1\).

This however may not be the optimal solution, in the sense of equation (13), and thus we may improve the quality of the solution executing simple genetic algorithm, with the individuals representing the parameters \(a_1^x,\ldots ,a_{N_x}^x,b_1^y,\ldots ,b_{N_y}^y,c_1^z,\ldots ,c_{N_z}^z\), and with the fitness function defined as (13).

4 Iterative Algorithm with Evolutionary Computations

The heuristic algorithm mixed with the genetic algorithm, as presented in Sect. 3, is not able to find the solution with 0 error, for non-tensor product structures, since we approximate \(N*N\) data with \(2*N\) unknowns. Thus, the iterative algorithm presented in 1 is proposed, with the assumed accuracy \(\epsilon \).

figure a

In the aforementioned algorithm we approximate the static permeability as a sequence of tensor product approximations:

$$\begin{aligned} K(x,y,z)=\sum _{m=1}^M K_m^x(x)K_m^y(y)K_m^z(z) \end{aligned}$$
(47)

Practically, it is realized according to the following equations:

$$\begin{aligned} \begin{aligned} \Big (u_{t+m},\upsilon \Big )- \Big ( K_m^x(x)e^{\mu \overline{u}_{t+m-1}(x)} \Big ) \Big ( K_m^y(x)e^{\mu \overline{u}_{t+m-1}(y)} \Big ) \Big ( K_m^z(x)e^{\mu \overline{u}_{t+m-1}(z)} \Big )\nabla u_{t+m},\nabla \upsilon \Big ) \\ = - \sum _{n=1,m\ne n} \Bigg ( K_n^x(x)e^{\mu \overline{u}_{t+n}(x)} K_n^y(y)e^{\mu \overline{u}_{t+n}(y)} K_n^z(z)e^{\mu \overline{u}_{t+n}(z)} \nabla u_{t+n}, \nabla \upsilon \Bigg ) \\ + \big (f,\upsilon \big ) +K_m^x(x)K_m^y(y)K_m^z(z) \\ \Big [\Big ( e^{\mu \overline{u}_{t+m}(x)} e^{\mu \overline{u}_{t+m-1}(y)} e^{\mu \overline{u}_{t+m-1}(z)} \Big ) - e^{\mu \overline{u}_{t+m-1}(x,y,z)} \Big ]\nabla u, \nabla \upsilon \Big ) \quad \forall \upsilon \in V \end{aligned} \end{aligned}$$
(48)

5 Numerical Results

We conclude the paper with the numerical results concerning the approximation of the static permeability map. The original static permeability map is presented in Fig. 1. The first approximation has been obtained from the heuristic algorithm described in Sect. 3. We used the formulas (25)–(27) with the suitable substitutions. In the first approach we first compute the values of a, next, the values of b and finally the values of c. As the initial values we picked \(\root 3 \of {d_{111}}\).

Deriving this method further we decided to compute particular points in the order of \(a_2\), \(b_2\), \(c_2\), \(a_3\), \(b_3\) and so on. This gave us the final result presented in Fig. 2a.

We have improved the approximation by post-processing with the generational genetic algorithm as implemented in jMetal package [13] with variables from [0,1] intervals. The fitness function was defined as:

$$\begin{aligned} \begin{aligned} f(a_1,\ldots ,a_{N_x},b_1,\ldots ,b_{N_y},c_1,\ldots ,c_{N_z})=\sum _{i=1}^{N_x}\sum _{l=1}^{N_y}\sum _{k=1}^{N_z}\Big (d_{ilk}-a_ib_lc_k\Big )^2 \end{aligned} \end{aligned}$$
(49)

The results are summarized in Fig. 2b.

To improve the numerical results we have employed the Algorithm 1. In Figs. 3 and 4 results obtained after 1, 5, 10 and 50 iterations of Algorithm 1 are presented.

In order to analyze the accuracy of the tensor product approximation, we also present in Figs. 5 and 6 the error after 1, 5, 10, 50 iterations. We can read from these Figures, how the error decreases when adding particular components.

6 Conclusions and the Future Work

In the paper the heuristic algorithm for tensor product approximation of material data for implicit dynamics simulations of non-linear flow in heterogeneous media is presented.

The algorithm can be used as a generator of initial configurations for a genetic algorithm, improving the quality of the approximation. The future work will involve the implementation of the implicit scheme and utilizing the proposed algorithms as a preconditioner for obtaining tensor product structure of the material data.

We have analyzed the convergence of our tensor product approximation method but assessing how the convergence influences the reduction of the iteration number of the explicit method will be the matter of our future experiments.

Our intuition is that 100 iterations (100 components of the tensor product approximation) should give a well approximation, and thus we can use the implicit method not bounded by the CFL condition, which will require 100 sub-steps in every time step.