Projection de champs - L'hypercube

1. Presentation of the problem

A common pitfall when solving problems with multiple physics is that the use of multiple tools often results in the creation of several meshes (see figure below). It is then a question of projecting the fields on a common grid in order to be able to carry out calculations at the same point in the space for the two calculated quantities.

Difference in mesh size between the tool for calculating solar fluxes and air velocities over a geometry.

In the following section, a mathematical tool is presented to project one field onto another, regular to no.

2. Projection methods

There are several algorithms of varying complexity to perform a field projection. Some variants are presented in this section.

2.1. Nearest neighbors

One technique for projecting a field onto another is to solve a minimization problem. The brief method is explained here given in code Aster in dimension two for the algorithm of the closest neighbors.

The aim is to project a cloud of points $N$ bearing the field $\varphi$ on a number of points $M$ different of $N$ . Let us call $\psi$ la projection de $\varphi$ on $M$ .

Two clouds of points N and M, respectively initial and projection

For each point $k$ belonging $M$ with coordinates $(x_k,y_k)$ , we're looking for its $n$ nearest neighbors belonging to N $N$ as shown in following Figure.

Illustrations of the n=5 closest neighbours of point k

The value of the field on the projected point $\psi_k$ is computed with coefficients $a_k,b_k,c_k$ such that:

\psi_k = a_k + b_k x_k + c_k y_k

It is therefore a question of finding, for all points $k$ of the cloud $M$ the arrays $[a_k], [b_k], [c_k]$ of coefficients that minimize the distance between the projected field $\psi$ field and the original field $\varphi$ Thus, for all point $k$ of the cloud $M$ we try to minimize the least squares $\varepsilon_k$ between the $\psi_k$ field and the original field $\varphi_i$ on the $n$ nearest neighbours

\varepsilon_k = \sum_{i=1}^{n} w(d_i) \big(\psi_k-\varphi_i \big)^2

Where the weighting function $w(d_i)$ depends on distance $d_i$ between point $(x_i,y_i)$ of the initial cloud and point $k$ of the projection grid $(x_k,y_k)$ :

w(d_i) = e^{-(\frac{d_i}{d_{r}})^\beta}

Where
– $d_i$ is the distance between $(i-k)$ , i. e. $d_i=\sqrt{(x_i-x_k)^2+(y_i-y_k)^2}$ ,
– $d_r$ is the reference distance of the mesh. This is defined as the radius of the smallest circle centred at a point of $M$ and containing three points of the cloud $N$ as shown in the figure below (this radius can be determined with a kdtree from the scipy Python library for example). As a first approximation, if the distance $\Delta_x$ between 2 points is known (e.g. in a regular mesh) we may consider $d_r \sim \frac{\Delta x}{2}$ . Otherwise, on can take the distance $d_r \sim \frac{D^2}{N} = \frac{2D}{\sqrt{N}}$ where $D$ is the greatest distance between any two points of the cloud.

Illustration of the reference distance d_r

The choice of the parameter set $(n,M,\beta,d_{r})$ will influence the global algorithm performance: given the cloud point distances one have to set a number $M$ big enough, while keeping in mind that the calculation time is proportional to $M$ . Furthermore, chosing $\beta>1$ implies that the influence of distant points is quickly soften when their distance from the reference point is greater than $d_{r}$ .

Tests conducted on analytical fields with a homogeneous point density have shown however that increasing the number of neighbors does not affect the quality of the projection. For such fields, the use of $n \in [4,8]$ is recommended, with a weighting function coefficient of $\beta=1.5$ .

Example of a field projected on a regular grid using the nearest neighbour method.

2.2 Shepard's method

(adapted from T. Kliyam internship at AREP)

This method is a form of inverse distance weighting which is a spatial interpolation method. The main effort here will be to build an interpolation function $Q(x,y)$ for the considered field, instead of using directly the value of the field

In the analytical fields studied, this method is $\sim 2 \times$ more performant in terms of precision, for an execution time increasing by an order of magnitude. If the gain may seem small compared to the computation time, the robustness of Shepard's method for the projection of real fields has been proven, in particular because it induces less smoothing on them.

Construction of an interpolation function on the initial field $N$

The simplest form of the method is called "Shepard's method" (Shepard 1968), with the following equation:
$f\left(x,y\right)=\sum_{k=1}^{N}{w_k\left(x,y\right) F_k}$

Where $w_k$ is a weighting function generally dependent on distance and $F_k$ is a function that calculates the value at the point $k$ . In the simplest cases, such as the nearest neighbor method, we define $F_k=\varphi_k$ the value of the initial field at point $k$ .

In Shepard's method, the weighting is defined as:

w_k=\frac{\left[\frac{\left(R_q-d_k\right)_+}{R_q\cdot d_k}\right]^2}{\sum_{k=1}^{N}\left[\frac{\left(R_q-d_k\right)_+}{R_q\cdot d_k}\right]^2}

Where

– $d_k$ is the distance between point $(x,y)$ and point $(x_k,y_k)$ , i. e. ${ d_k=\sqrt{\left(x-x_k\right)^2+\left(y-y_k\right)^2}}$ , with $\left(x,y\right)$ the coordiantes of the interpolated point and $(x_k,y_k)$ the coordinates of the cloud of points of the initial field.

– $R_q$ is the influence radius of a circle centered on the point $(x,y)$ (represented on the figure below). In practice we choose $N_q \sim 40 - 50$ . The radius of influence is calculated so that

$(R_q-d_k)_+= R_q - d_k \ \ \text{ pour }\ \ d_k < R_q$
$(R_q-d_k)_+=0 \text{ pour } d_k \ge R_q$

Let $D$ be the distance between the 2 most distant points of $N$ . We define $R_q$ through the nodal function parameter $N_q$ such that
$R_q=\frac{D}{2}\sqrt{\frac{N_q}{N}}$

As for the value of the function $F_k$ we then construct an interpolation law $Q_k(x,y)$ defined as follows:

Q_k(x,y)=\varphi_k+a_{k2}(x-x_k)+a_{k3}(y-y_k)+a_{k4}(x-x_k)^2+a_{k5}(x-x_k)(y-y_k)\\+a_{k6}(y-y_k)^2

Where $\varphi_k$ is the point value at $(x_k,y_k)$ and the coefficients $(a_{k2}, ...,\ a_{k6})$ are specific to each point $(x_k,y_k)$ . They are determined by the following least squares minimization, for $k \in [1,N^*]$ , where $N^*$ is the number of points in the radius of influence $R_q$ (in order not to overload the notations we will write indifferently $N^*$ and $N$ in the following, indeed, the weighting function cancels the contribution of the other points of the cloud):

\varepsilon_k = \sum_{ i=1\ et\ i\neq k}^{N^*}\left [ \frac{Q_k(x_i,y_i)-\varphi_i}{\rho_k(x_i,y_i)} \right ]^2

Explicitly:

\varepsilon_k = \sum_{ i=1\ et\ i\neq k}^{N^*} \bigg [ \frac{(R_q-d_i)_+}{R_qd_i} \bigg( \varphi_k+a_{k2}(x_i-x_k)+a_{k3}(y_i-y_k)+a_{k4}(x_i-x_k)^2 \\ +a_{k5}(x_i-x_k)(y_i-y_k)+a_{k6}(y_i-y_k)^2 - \varphi_i \bigg) \bigg ]^2

With $d_i$ the distance between points $i$ and $k$ : $d_i= \sqrt{(x_k-x_i)^2+(y_k-y_i)^2}$

The interpolation equation then becomes :

f(x,y)=\sum_{k=1}^{N^*}{w_k(x,y) \cdot Q_k(x,y)}

Note that the computational effort is located in the identification step of the parameters $a_{2k},..., a_{6k}$ : the greater the number $N_q$ is, the larger the number of points involved, which increases the number of coefficients to be determined for the minimization function.

Calculation of the projection on the field $M$

We then define the radius of influence $R_w$ of which an illustration is given in the Figure below, such that:

R_w=\frac{D}{2}\cdot \sqrt{\frac{N_w}{N}}

In practice, an effective choice is to take $N_w=N_q/2$ .

Radius of influence Rw and distance dk (from (Klinyam 2019))

The projection on the point $(x_m,y_m)$ of the cloud $M$ is calculated by summing the interpolation functions for each point in the radius of influence $R_w$ (see figure below):

\psi_m=\sum_{k=1}^{N}{w_k(x_m,y_m)}\cdot Q_k(x_m,y_m)

Projected field calculation with Q_k interpolation functions of points within the radius of influence (from (Klinyam 2019)).

It is then necessary to calculate the projection $\psi_m$ with the weighting function $w_k$ defined analogously to the equation \eqref{eq:shep1_wk2} whose radius of influence $R_w$ is defined according to the equation \eqref{eq:shep_rq} using the influence parameter $N_w$ .

An applied example of Shepard's method is shown below.

Applied example with urban air flows: amplification factor map - Initial field (left) and projected (right) (from (Klinyam 2019)).

Results of the various minimisation algorithms on the execution time

Lastly, we try out some of the internal functions from the Optimize library from Scipy for an analytical 2D field projection problem. Interestingly, the L-BFGS-B method is the fastest. Above a threshold of about 2000 points on which we compute the interpollation function $Q(x,y)$ , one can notice a linear increase of the execution time.