ON THE SPACE SEPARATED REPRESENTATION WHEN ADDRESSING THE SOLUTION OF PDE IN COMPLEX DOMAINS

. Separated representations allow impressive computational CPU time savings when applied in diﬀerent ﬁelds of computational mechanics. They have been extensively used for solving models deﬁned in multidimensional spaces coming from (i) its proper physics, (ii) model parameters that were introduced as extra-coordinates and (iii) 3D models when the solution can be separated as a ﬁnite sum of functional products involving lower dimensional spaces. The last route is especially suitable when models are deﬁned in hexa-hedral domains. When it is not the case, diﬀerent possibilities exist and were considered in our former works. In the present work, we are analyzing two alternative routes. The ﬁrst one consists of immersing the real non-separable domain into a fully separable hexahedral domain. The second procedure consists in applying a geometrical transformation able to transform the real domain into a hexahedra in which the model is solved by using a fully separated representation of the unknown ﬁeld.


1.
Introduction. Today many problems in science and engineering remain intractable, in spite of the impressive progresses attained in modeling, numerical analysis, discretization techniques and computer science during the last decade, because their numerical complexity, or the restrictions imposed by different requirements (real-time on deployed platforms, for instance) make them unaffordable for today's technologies. Among them we can consider: • Models that are defined in high dimensional spaces, usually encountered in quantum chemistry describing the structure and mechanics of materials [15], the kinetic theory description of complex materials [10], social dynamics and economic systems, vehicular traffic flow phenomena, complex biological systems involving mutation and immune competition, crowds and swarms encountered in congested and panic flows, among many other unimaginable possibilities (see [8] and the references therein); the chemical modeling in too diluted systems where the concept of concentration cannot be used, that results in the so-called chemical master equation governing for example cell signaling and other phenomena in molecular biology. • Many problems in parametric modeling, inverse identification, and process or shape optimization, usually require, when approached with standard techniques, the direct computation of a very large number of solutions of the concerned model for particular values of the problem parameters. When the number of parameters increases such a procedure becomes inapplicable. On the other hand a new paradigm in the field of Applied Sciences and Engineering has emerged in the last decade. Dynamic Data-Driven Application Systems (DDDAS) constitute nowadays one of the most challenging applications of simulation-based Engineering Sciences. By DDDAS we mean a set of techniques that allow the linkage of simulation tools with measurement devices for real-time control of simulations. • Finally many models require a fine (3D) solution event when they are defined in degenerated domains like beams, plates or shells like domains. In those cases, when 3D solutions are compulsory, standard mesh-based discretizations are challenging in order to prevent the existence of too distorted elements.
While the previous list is by no means exhaustive, it includes a set of problems with no apparent relationship between them that can however be treated in a unified manner.
The solution of these three kind of models, assumed, without loss of generality, scalar, can be expressed respectively from: • Solution involving physical and conformational coordinates: where x and t are the standard space and time coordinates, x ∈ Ω ⊂ R d (d = 1, 2, 3) and t ∈ I ⊂ R. Coordinates c 1 to c D denotes the conformational coordinates depending on the involved physics. When c j ∈ C j , j = 1, · · · , D, then one can approximate the unknown field in the separated form: • Solution involving the parameters (related to the model, to the geometry or to the initial or boundary conditions) as model extra-coordinates: where again x and t are the standard space and time coordinates and the extracoordinates p 1 to p D denotes the parameters considered now as coordinates. When p j ∈ P j , j = 1, · · · , D, then one can approximate the unknown field in the separated form: • Solution involving a partial or fully separated representation of the space coordinates: u(x, t) = u(x, y, z, t) (5) When the model is defined in a hexahedral domain, the solution can we searched as In the case of extruded profiles, plates or shells the most appropriate representation reads: These separated representations were considered extensively in our former works: [16] [6] and [31] in the case of multidimensional physics; [5] [23] [24] [27], [39] and [18]in the case of parametric models and in [11] [32] [19] [12] and [25] in the case of geometrically degenerated domains. Separated representations in computational mechanics come from the pioneering works of P. Ladeveze that proposed using space-time separated representations as one of the main ingredients of the LATIN nonlinear/multiscale solver [29] [30] and the references therein. Later stochastic and deterministic coordinates were separated in the works of A. Nouy [37] and the references therein. The invested reader can refer to the recent review [17] that contains a large list of valuable references.
There exist in the context of model order reduction alternatives to the separated representation just considered. Proper Orthogonal Decomposition and Reduced Bases are the two most valuable alternatives.
Proper orthogonal decompositions (POD) allows extracting the most significant characteristic of the solution, that can be then applied for solving models slightly different to the ones that served to defined the reduced approximation bases. There is an extensive literature. The interested readers can reefer to [38] [43] and the numerous references therein. The extraction of the reduced basis is the tricky point when using POD-based model order reduction, as well its adaptivity when addressing scenarios far from the ones considered when constructing the reduced basis [41] [42].
On the other hand, reduced bases (RB) are constructed by combining a greedy algorithm and a priori error indicator. After some off-line work, the resulting reduced basis can be used on-line for solving different models with a perfect control of the solution accuracy because the availability of error bounds. When the error is inadmissible, the reduced basis can be enriched by invoking again the same greedy algorithm. The interested readers can refer to [33] [34] [44] [40] and the references therein.
In the context of the three computational challenges described above, POD and RB can only address efficient solutions to the second one, that is, the efficient solution of parametric models and real time simulation of many model encountered in computational physics and engineering, however, both procedures become inefficient for addressing multidimensional physics involving conformational coordinates for representing fine physics or for addressing models defined in degenerated domains where a space decomposition seems compulsory.
In this paper we focuses in the efficient treatment of models involving the usual space coordinates (x, y) where for the sake of simplicity and without loss of generality the time dependence is ignored and only 2D models are addressed. As just discussed when the model is defined in a rectangular domain the decomposition becomes the most natural choice. Its construction is revisited in section 2.
In the case of non rectangular domains (hexahedral in the 3D case) different possibilities exist, among them: • Immerse the real domain ω into a rectangular domain Ω whilst enforcing the boundary condition on γ := ∂ω. • Consider an appropriate geometrical transformations between the real and a reference domain, the last fully separable (rectangular in 2D and hexahedral in 3D). The first approach was considered in [26] whereas a particular case of the second one was successfully applied in [7] for addressing shape optimization.
In this work we propose two new approaches, one belonging to the first class of techniques that proceeds by immersing the non-separable domain into a separable one and the second one associated to a appropriate mapping able to transform some classes of non-separable geometries into separable domains. Both will be described and validated in section 3 and 4 respectively.

2.
On the space separated representation in hexahedral domains. It is now time to detail the inner workings of the PGD. We begin with a simple but illustrative case study. Consider the solution of the Poisson equation in a two-dimensional rectangular domain Ω = Ω x × Ω y = (0, L) × (0, H). We specify homogeneous Dirichlet boundary conditions for the unknown field u(x, y), i.e. u(x, y) vanishes at the domain boundary Γ. Furthermore, we assume that the source term f is constant over the domain Ω.
For all suitable test functions u * , the weighted residual form of (9) reads or more explicitly Our goal is to obtain a PGD approximate solution to (9) in the separated form We shall do so by computing each term of the expansion one at a time, thus enriching the PGD approximation until a suitable convergence criterion is satisfied.
2.1. Progressive construction of the separated representation. At each enrichment step n (n ≥ 1), we have already computed the n − 1 first terms of the PGD approximation (12): We now wish to compute the next term X n (x) · Y n (y) to obtain the enriched PGD solution Both functions X n (x) and Y n (y) are unknown at the current enrichment step n, and they appear in the form of a product. The resulting problem is thus non-linear and a suitable iterative scheme is required. We shall use the index p to denote a particular iteration.
At enrichment step n, the PGD approximation u n,p obtained at iteration p thus reads u n,p (x, y) = u n−1 (x, y) + X p n (x) · Y p n (y). (15) The simplest iterative scheme is an alternating direction strategy that computes X p n (x) from Y p−1 n (y), and then Y p n (x) from X p n (x). An arbitrary initial guess Y 0 n (y) is specified to start the iterative process. The non-linear iterations proceed until reaching a fixed point within a user-specified tolerance , i.e.
where · is a suitable norm. The enrichment step n thus ends with the assignments X n (x) ← X p n (x) and Y n (y) ← Y p n (y). The enrichment process itself stops when an appropriate measure of error E(n) (based on the norm of the residual or on a quantity of interest among many other possible choices) becomes small enough, i.e E(n) <˜ . Several stopping criteria are available [20]. We now describe in more detail one particular alternating direction iteration at a given enrichment step.

2.2.
Alternating direction strategy. Each iteration of the alternating direction scheme consists in the following two steps: • Calculating X p n (x) from Y p−1 n (y) In this case, the approximation reads where all functions are known except X p n (x).
The simplest choice for the weight function u * in the weighted residual formulation (11) is which amounts to select the Galerkin weighted residual form of the Poisson equation.
Injecting (17) and (18) into (11), we obtain Here comes a crucial point: since all functions of y are known in the above expression, we can compute the following one-dimensional integrals over Ω y : Equation (19) becomes We have thus obtained the weighted residual form of a one-dimensional problem defined over Ω x that can be solved (e.g. by the finite element method) to obtain the function X p n we are looking for. Alternatively, we can return to the corresponding strong formulation and solve it numerically by means of any suitable numerical method (e.g. finite differences, pseudo-spectral techniques, . . . ). The strong form (22) is a second-order ordinary differential equation for X p n . This is due to the fact that the original Poisson equation involves a second-order x-derivative of the unknown field u.
With either the weighted residual or strong formulations, the homogeneous Dirichlet boundary conditions X p n (x = 0) = X p n (x = L) = 0 are readily specified.
Having thus computed X p n (x), we are now ready to proceed with the second step of iteration p.
• Calculating Y p n (y) from the just-computed X p n (x) The procedure exactly mirrors what we have done above. Indeed, we simply exchange the roles played by all relevant functions of x and y. The interested reader can refer to [20] for a detailed description and analysis. We have thus completed iteration p at enrichment step n. It is important to realize that the original two-dimensional Poisson equation defined over Ω = Ω x × Ω y has been transformed within the PGD framework into a series of decoupled one-dimensional problems formulated in Ω x and Ω y .
3. Immersed formulation for non-separable domains. We consider a model defined in a complex domain ω ⊂ R 2 where its space decomposition (8) does not apply straightforward, that is, ω is not rectangular. One possibility consists in immersing ω into a rectangular domain Ω, ω ⊂ Ω ⊂ R 2 , where Ω = Ω x × Ω y . The boundaries of domains ω and Ω, ∂ω and ∂Ω, are noted respectively by γ and Γ.
Consider the problem defined by: with n the unit outwards vector defined on the domain boundary. γ D and γ N is a partition of γ according to the specified boundary conditions, Dirichlet and Neuman respectively, veryfying γ D ∪ γ N = γ and γ D ∩ γ N = ∅. Moreover the given functions f (x) and g(x) are assumed regular enough. The solution of problem (23), u(x), is supposed to be the restriction to ω of a function U (x) defined in Ω that verifies ∆U = 0. Note that this consideration exclude solutions of problem (23) involving singularities in the neighborhood of γ.
First we assume that γ D = γ, that is, the problem only contains Dirichlet boundary conditions (more complex situation will be considered later). If we consider H a Hilbert space defined from H(Ω) = {U : Ω → R/∆U = 0 in Ω}, u = U | ω := U (x ∈ ω) is solution of problem (23) if and only if: where · γ is a norm on H(γ), associated with the scalar product (·|·) γ . Two natural choices for H(γ) are: Thus (24) defines a well posed problem. Its solution determines the values of the solution on γ as well as the values of its normal derivative. Thus, the calculation of U (x) in Ω − ω defines a Cauchy's problem that in principle is ill posed. For solving it we consider the evanescent regularization procedure introduced in [21] [22] that proceed by iterating until reaching convergence according to: with k ≥ 1, c a positive constant small enough and U 0 = 0. In order to proceed with the just described algorithm, and because the linearity of the considered model, we start by considering M n nodes on Γ, with coordinates X r , r = 1, · · · M n , for approximating U | Γ := U (x ∈ Γ). Now, we compute the solution U r (x), solution of ∆U = 0 with U r (X s ) = δ rs , ∀r, s = 1, · · · , M n , with δ rs the Kroenecker's delta function. Being Ω = Ω x × Ω y the space decomposition (8) can be easily defined within the PGD framework: being its extension to higher dimensional spaces straightforward. Now, the general solution of ∆U = 0 in Ω can by expressed from: The unknowns coefficients α i , i = 1, · · · , M n , α i = U (X i ), are computed in order to fulfill Eq. (27). For this purpose we consider m n nodes on γ, with coordinates x s , s = 1, · · · , m n , able to describe accurately f (x ∈ γ), and we define the residuals r s 1 and r t 2 at iteration k, related to the tentative solution U k : The cost function C(α 1 , · · · , α M n ) to be minimized reads: and the minimization is performed by using the Levenberg -Marquardt's strategy. The different components of the Jacobian J have simple expressions: for s = 1, · · · , m n , and ∂r t for t = 1, · · · , M n . Finally, the Levenberg-Marquardt updating reads: where U k is the vector containing the tentative alpha coefficients at iteration k, i.e. U k = (α k 1 , · · · α k Mn ), λ is a parameter which improves convergence and r is the vector of residuals.
3.1. Numerical verification. In order to check the proposed strategy we compute the solution of problem (23) in ω, where the nodes used for representing the solution on γ ≡ ∂ω are represented in Fig. 1. Fig. 2 compares f (x ∈ γ) with U | γ := U (x ∈ γ), whose gap measured using the L 2 norm is lower than 1% after 3 iterations. A more complex domain ω is shown in Fig. 3, with again f (x ∈ γ) and U | γ compared in Fig. 4. In this case again the error was about 1% after 3 iterations.  Figure 2. f (x ∈ γ) (in red) versus U | γ := U (x ∈ γ) (blue) 3.2. Extension to more complex boundary conditions. We applied the same strategy to a similar problem ∆u(x) = 0, x ∈ ω, involving both essential (Dirichlet) and Robin boundary conditions, on γ D and γ R respectively, with γ = γ D ∪ γ R .  Figure 4. A more complex geometry: f (x ∈ γ) (in red) versus U | γ := U (x ∈ γ) (blue) reduced number of nodes on Γ. As expected the accuracy improves when increasing the number of nodes on Γ.
Robin's boundary condition reads: being g(x) a given constant value, i.e. g(x) ≡ g. On the other hand we prescribe u(x) on γ D from u(x ∈ γ D ) = f (x). In this case, the evanescent regularization writes: Fig. 6, computed after three iterations, contains an error using the L 2 norm lower than 1%.
3.3. Application to shape optimization. Shape optimization requires in general solving a model in a given geometry, compute the cost function related to such solution and if it is not considered good enough, the domain geometry is modified and the model is solved again. The iteration procedure continues until reaching the convergence, that is, an acceptable value of the cost function. The main issue Figure 7. Parametrized geometry of domain ω related to such an approach is the necessity of constructing a mesh in each tentative domain before solving the model by using a standard mesh-based discretization technique.
In this section the immersed domain technique previously proposed leads to an efficient solution procedure. For illustrating the technique we consider the problem just analyzed in section 3.2 but now the domain is parametrized by two geometrical parameters L 1 and L 2 indicated in figure 7.
The cost function related to the shape optimization aims minimizing the part area A = L 1 × L 2 whilst maximizing the thermal flux exchanges Q through γ R Both competing criteria are expressed respectively from the two cost functions and Because the problem is multi-objective we decided to proceed by determining the Pareto's front. We consider a tentative domain ω(L 1 , L 2 ), from which we can extract its boundary γ(L 1 , L 2 ). Now, as each tentative domain is immersed into Ω, i.e. ω(L 1 , L 2 ) ⊂ Ω and solutions U i (x ∈ Ω) were recalculated in a separated form and offline for i = 1, · · · , M n , we only need few iterations for minimizing the boundary gap according to (36). Thus in few seconds and without remeshing needs we can evaluate many choices of the parameters defining the domain geometry. Each tentative geometry consists in a point in the (J 1 , J 2 ) representation. The Pareto's front is easily determined from the envelop of the resulting cloud of points. Other strategy for determining it exists, but they do not offer additional computing time savings. Figure 8 shows the cloud of points related to different evaluations as well as the extracted Pareto's front. However, an appropriate geometrical transformation allows its extension.
Consider for example the domainω, with boundaryγ ≡ ∂ω, shown in Fig. 9 (in which ∆u = 0 with u(x ∈γ) = f (x) is defined) that using polar coordinates results in the domain ω (with boundary γ) depicted in Fig. 10. Now, again its immersion in a rectangular one Ω (Ω = Ω r × Ω θ with boundary Γ), as depicted in Fig. 11, allows using the strategy previously described.
For addressing the solution of the model in the new coordinate system it suffices of expressing the Laplace's differential operator in polar coordinates and proceed as previously. The only difference concerns the boundary conditions. Now, on the boundary θ = 0 and θ = 2π where γ and Γ superpose, we must enforce ∀r ∈ Ω r the solution periodicity that writes: Figure 11 shows the rectangular domain Ω, with the boundary's partition Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 . In this case boundary conditions are enforced on Γ 1 and Γ 3 , while periodicity applies on Γ 2 and Γ 4 . A solution related to a choice of coefficients α 1 , · · · α Mn is depicted in Fig. 12. Fig. 13 illustrates the nodal distribution considered for approximating the solution in γ. Fig. 14 shows the computed solution in polar coordinates, being its Cartesian counterpart the one depicted in Fig. 15. 4. Geometrical transformation for non-separable domains. In this section we consider a non separable domain ω ⊂ R 2 with boundary γ ≡ ∂ω, where we with Dirichlet boundary conditions on γ u(x ∈ γ) = f (x) (43) where s(x) represents the source term. Figure 11. Immersion of ω into the larger separable (rectangular) domain Ω Figure 12. Solution U (r, θ) in Ω related to a particular choice of parameters α 1 , · · · , α Mn If the domain boundary is expressed from (θ, R(θ)), with θ ∈ [0, 2π) as illustrated in Fig. 16, we can define the mapping: that transform the complex domain ω into a separable one Ω = Ω r × Ω θ defined in the coordinate system (r, θ), with r ∈ Ω r = (0, 1) and θ ∈ Ω θ = [0, 2π). Again, on the boundaries related to θ = 0 and θ = 2π periodicity boundary conditions apply. On the boundary r = 1 we must apply the boundary conditions specified on γ and finally at the pole r = 0, and in absence of localized sources at this position, the conservation balance implies [13]: lim r→0 2π 0 ∇u · n rdθ = 0 (45) Figure 13. Nodal distribution used for approximating the solution on γ: blue stars are nodes where boundary conditions must be enforced while red lines are related to the nodes in which solution periodicity applies Figure 14.
The partial differential equation related to the initial domain ω can be easily transformed into the one applying in Ω by considering the coordinates transformation (44). If we define From Eq. (48) we can calculate second order derivatives leading to the problem strong formulation in Ω: that introducing the notation can be rewritten as It can be noticed that as soon as R(θ) = 1, Eq. (51) reduces to the Laplace equation expressed in polar coordinates (r, θ). where J denotes the mapping Jacobian and ∆ (r,θ) the Laplace differential operator when considering the (r, θ) coordinate system. The 2D resulting problem can be easily solved by integrating by parts and discretizing the resulting week form by using any standard and well experienced discretization technique like the finite element method.
Being Ω separable, i.e. Ω = Ω r ×Ω θ and because the differential operator involved in Eq. (51) is also separable, one could easily compute the separated form of the problem solution u(r, θ) from The construction is carried out by enriching until convergence the finite sum representation, assumed known at iteration n, u n (r, θ): looking for the solution u n+1 (r, θ) according to: The test function u * (r, θ) in Eq. (54) is taken as By introducing the trial and test functions, (57) and (58) respectively, into the problem weak form (54) and following the rationale deeply described in section 2, we obtain two boundary value problems for determining functions R n+1 (r) and Θ n+1 (θ).
The error was defined from the L 2 norm of u exact −u N , i.e. E(N ) = u exact (r, θ)− u N (r, θ) , where u N (r, θ) refers to the PGD solution consisting of N functional products.
The evolution of the error with the number of functional products involved in the finite sum (N ) is depicted if Fig. 17 where we can notice that as expected the solution can be represented with a single mode because the extant solution consists of a single product of functions, one depending on r, r 3 , and the other on θ, cos θ. The solutions in both systems of coordinates u(x, y) and u(r, θ) are depicted in Fig. 18, being the error distribution E(x, y) = u exact (x, y) − u P GD (x, y) and E(r, θ) = u exact (r, θ) − u P GD (r, θ) shown in Fig. 19. Here as just indicated u P GD ≡ u N =1 .

4.2.2.
A more complex domain. Now, we are solving again ∆u(x, y) = s(x, y) in the domain Ω defined from ω depicted in Fig. 20. The domain boundary γ is composed of an internal and external part, γ i and γ e respectively, γ = γ i ∪ γ e , where homogeneous Dirichlet boundary conditions are enforced, i.e. u(x ∈ γ) = 0.
The source term s(x, y) writes in the present case In this case points on γ i and γ e can be expressed from and x = r e R(θ) cos θ y = r e R(θ) sin θ , θ ∈ [0, 2π) respectively, with r i = 0.1, r e = 1 and R(θ) = 0.6 + 0.25sin(5θ) Thus the domain ω can be transformed into a separable domain Ω by using the (r, θ) coordinates, i.e. Ω = Ω r × Ω θ , with Ω r = (0.1, 1) and Ω θ = [0, 2π). The x = rR(θ) cos θ y = rR(θ) sin θ , r ∈ Ω r , θ ∈ Ω θ (63) This problem was first solved by using the FEM and a fine enough mesh of Ω, whose solution was considered as reference solution for quantifying the errors associated with the PGD solution when different number of terms N were involved in the finite sum decomposition. 200 nodes where uniformly distributed in the domains Ω r and Ω θ where the PGD solution was computed.
The error was defined from the L 2 norm of u F EM −u N , i.e. E(N ) = u F EM (r, θ)− u N (r, θ) , where u N (r, θ) refers to the PGD solution consisting of N functional products.
The evolution of the error with the number of functional products involved in the finite sum (N ) is depicted if Fig. 21 where we can notice that as expected the solution needs more terms for approximating it because the complexity of both the geometry and the source term.
The solutions in both systems of coordinates u(x, y) and u(r, θ) are depicted in Fig. 22, being the error distribution E(r, θ) = u F EM (r, θ) − u P GD (r, θ) shown in Fig. 23. Finally Fig. 24 illustrates the evolution of the CPU time when increasing the number of nodes N r involved in the discretization of domain Ω r for a given precision of the PGD solution.

5.
Conclusions. In this paper we analyzed two different approaches for performing space separated representations when considering models defined in non-separable domains. The analysis carried out proved their ability to efficiently address this type of scenarios.  The first approach consists of immersing the complex domain into a hexahedral one that allows the use of fully separated space representations. The algorithm proceeds by identifying the boundary conditions in that domain in order to ensure the verification of the real boundary conditions on the boundary of the non-separable one.
The second technique consists of performing an adequate mapping able to transform the non-separable domain into a separable one. The price to pay is the transformation of all the differential operators to the new system of coordinates, but in many cases this is not a major issue.
We noticed that the solution procedure when considering the immersed formulation for non-separable domains described in Section 3 involved the solution of many problems, whose solution does not represented a big challenge due to the use of Figure 23. E(r, θ) = u F EM (r, θ) − u P GD (r, θ) Figure 24. FEM versus PGD CPU time when increasing the number of nodes N r involved in the discretization of Ω r separated representations within the PGD framework. However, after being proved that the solution in a circular (spherical) domain can be easily addressed thanks to the geometrical transformation for non-separable domains described in Section 4, one could imagine solving a single problem in the circular geometry, being the solutions related to the different boundary conditions considered in Section 3 obtained from the first one by applying an adequate rotation.