Field Projections

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 Nas 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 Mthe arrays [a_k], [b_k], [c_k] of coefficients that minimize the distance between the projected field \psi field and the original field \varphiThus, 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 Nas 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>1implies 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_kthe 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

Influence radius R_q in the N field

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_kwe 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_wof 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.

Execution times for three minimisation algorithms (from (Klinyam 2019))