Projection de champs

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 cloud $M$ the arrays $[a_k], [b_k], [c_k]$ of coefficients that minimize the distance between the projected field $\psi$ and the original field $\varphi$ Thus, for all point $k$ of cloud $M$ on cherche à minimiser les moindres carrés $\varepsilon_k$ entre la fonction $\psi_k$ 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 du package scipy par exemple). En première approximation, si la distance $\Delta_x$ entre deux points du nuage est connue (cas d’un maillage régulier par exemple), on peut prendre $d_r \sim \frac{\Delta x}{2}$ . Autrement, on pourra aussi utiliser la distance $d_r \sim \frac{D^2}{N} = \frac{2D}{\sqrt{N}}$ where $D$ est la distance entre les deux points les plus éloignés du nuage.

Illustration of the reference distance d_r

The choice of the parameter set $(n,M,\beta,d_{r})$ influe sur la performance globale de l’algorithme : selon la distance entre points du maillage il convient d’avoir un nombre de points de projection $M$ while keeping in mind that the calculation time is proportional to $M$ . What is more, chosing $\beta>1$ implies that the influence of distant points is quickly amortized when their distance from the reference point is greater than $d_{r}$ .

Les tests menés sur champs analytiques avec un densité de points homogène ont montré cependant que l’augmentation du nombre de voisins n’affecte pas la qualité de la projection. Pour de tels champs l’utilisation de $n \in [4,8]$ est conseillée, avec pour coefficient de fonction de pondération $\beta=1.5$ .

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

2.2 La méthode de Shepard

(élaboré à partir des travaux de T. Klinyam au sein d’AREP)

Cette méthode est une forme de pondération inverse à la distance qui est une méthode d’interpolation spatiale. L’effort principal ici va consister à construire une fonction d’interpolation $Q(x,y)$ pour le champ considéré, au lieu d’utiliser directement la valeur du champ.

Sur les champs analytiques étudiés, cette méthode est $\sim 2 \times$ plus performante en termes de précision, pour un temps d’exécution augmentant d’un ordre de grandeur. Si le gain peut paraître faible par rapport au temps de calcul, la robustesse de la méthode de Shepard pour la projection de champs réels a été éprouvée, notamment car elle induit moins de lissage sur ceux-ci.

Construction d’une fonction d’interpolation sur le champ initial $N$

La forme la plus simple de la méthode est appelée « Méthode de Shepard » (Shepard 1968) avec l’équation suivante :
$f\left(x,y\right)=\sum_{k=1}^{N}{w_k\left(x,y\right) F_k}$

Where $w_k$ est une fonction de pondération généralement dépendante de la distance et $F_k$ est une fonction qui permet de calculer la valeur au point $k$ . Dans les cas les plus simples, comme la méthode des plus proches voisins, on définit $F_k=\varphi_k$ , soit la valeur du champ initial au point $k$ .

Dans la méthode de Shepard, on définit la pondération telle que :

$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$ est la distance entre le 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}}$ , dont $\left(x,y\right)$ sont les coordonnées du point interpolé et $(x_k,y_k)$ sont des coordonnées du nuage des points du champ de départ

– $R_q$ est le rayon d’influence d’un cercle centré sur le point $(x,y)$ (représenté sur la figure ci-dessous), En pratique on choisit $N_q \sim 40 - 50$ . Le rayon d’influence est calculé de sorte que :

$(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$

Rayon d’influence Rq dans le champ N (tiré de (Klinyam 2019))

Soit $D$ la distance entre 2 points les plus éloignés du nuage $N$ . On définit $R_q$ via le paramètre de fonction nodale $N_q$ such that
$R_q=\frac{D}{2}\sqrt{\frac{N_q}{N}}$

Pour ce qui est de la valeur de la fonction $F_k$ , on construit ensuite une loi d’interpolation $Q_k(x,y)$ définie de la manière suivante :

$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$ est la valeur sur le point $(x_k,y_k)$ et les coefficients $(a_{k2}, ...,\ a_{k6})$ sont propres à chaque point $(x_k,y_k)$ . On les détermine par la minimisation des moindres carrés suivants, pour $k \in [1,N^*]$ , where $N^*$ est le nombre de points dans le rayon d’influence $R_q$ (afin de ne pas surcharger les notations on écrira indifféremment $N^*$ and $N$ dans la suite, en effet la fonction de pondération annule la contribution des autres points du nuage) :

$\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$

Soit, explicitement :

$\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$ la distance entre les points $i$ and $k$ : $d_i= \sqrt{(x_k-x_i)^2+(y_k-y_i)^2}$

L’équation d’interpolation devient alors :

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

On notera que l’effort calculatoire est situé dans l’étape d’identification des paramètres $a_{2k},..., a_{6k}$ : plus le nombre $N_q$ est grand, plus le nombre de points concerné est important, ce qui augmente la quantité de coefficients à déterminer pour la fonction de minimisation.

Calcul de la projection sur le champ $M$

On définit alors le rayon d’influence $R_w$ , dont une illustration est donnée sur la Figure ci-après, tel que :

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

En pratique, un choix efficace consiste à prendre $N_w=N_q/2$ .

Rayon d’influence Rw et distance dk (tiré de (Klinyam 2019))

La projection sur le point $(x_m,y_m)$ of cloud $M$ est calculée en faisant la somme des fonctions d’interpolations pour chaque point dans le rayon d’influence $R_w$ (voir figure ci-après) :

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

Calcul du champ projeté avec les fonctions d’interpolation Q_k des points dans le rayon d’influence (tiré de (Klinyam 2019)).

Il s’agit ensuite de calculer la projection $\psi_m$ avec la fonction de pondération $w_k$ définie de manière analogue à l’équation \eqref{eq:shep1_wk2} dont le rayon d’influence $R_w$ est défini selon l’équation \eqref{eq:shep_rq} en utilisant le paramètre d’influence $N_w$ .

Un exemple d’application de la projection selon la méthode de Shepard est donné ci-dessous.

Application à l’aéraulique : représentation du facteur d’amplification – champ initial (gauche) et champ projeté (droite) (tiré de (Klinyam 2019)).

Résultats sur le temps d’exécution des différents algorithmes de minimisation

On compare ici quelques unes des fonctions proposées la librairie optimize de scipy pour un problème de projection de champ analytique en 2D. De manière intéressante, la méthode L-BFGS-B est la plus performante. À partir de plus de 2000 points sur lesquels on calcule la fonction d’interpolation $Q(x,y)$ , on constate une augmentation linéaire du temps de calcul.