Projection de champs

1. Présentation du problème

Un écueil courant lors de la résolution de problèmes avec des physiques multiples est que l’utilisation de plusieurs outils entraîne souvent la création de maillages différents (voir la figure ci-dessous). Il s’agit alors de projeter les champs sur une grille commune afin de pouvoir réaliser des calculs en un même point de l’espace pour les deux grandeurs calculées.

Différence de maillage entre l’outil pour le calcul des flux solaires et des vitesses d’air sur une géométrie.

Dans la section qui suit, on présente un outil mathématique permettant de projeter un champ sur un autre, régulier au non.

 

2. Méthodes de projection

Il existe plusieurs algorithmes de complexité variable pour réaliser une projection de champ. On en présente quelques variantes dans cette section.

2.1. Les plus proches voisins

Une technique pour projeter un champ sur l’autre consiste à résoudre un problème de minimisation. On explicite ici la méthode succincte donnée dans code Aster en dimension deux pour l’algorithme des plus proches voisins.

Il s’agit de projeter un nuage de points N portant un champ \varphi sur un ensemble de points M différent de N. Appelons \psi la projection de \varphi sur M.

Deux nuages de points N et M, respectivement initial et de projection

Pour chaque point k appartenant à M de coordonnées (x_k,y_k), on cherche ses n plus proches voisins appartenant à N, ainsi qu’illustré sur la Figure qui suit.

Illustrations des n=5 plus proches voisins du point k

La valeur du champ sur le point projeté \psi_k est calculée avec les coefficients a_k,b_k,c_k de sorte que :

\psi_k = a_k + b_k x_k + c_k y_k

Il s’agit donc de trouver, pour l’ensemble des points k du nuage M, les vecteurs [a_k], [b_k], [c_k] des coefficients qui minimisent la distance entre le champ projeté \psi et le champ \varphi. Ainsi, pour tout point k du nuage M on cherche à minimiser les moindres carrés \varepsilon_k entre la fonction \psi_k et le champ \varphi_i sur les n plus proches voisins du point :

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

Où la fonction de pondération w(d_i) dépend de la distance d_i entre le point (x_i,y_i) du nuage initial et le point k de la grille de projection (x_k,y_k) :

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

Avec :
d_i est la distance (i-k) soit d_i=\sqrt{(x_i-x_k)^2+(y_i-y_k)^2},
d_r est la distance de référence du maillage. Celle-ci est définie comme le rayon du plus petit cercle centré en un point de M contenant trois points du nuage N, comme illustré sur la figure ci-après (on pourra déterminer ce rayon avec un 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}}D est la distance entre les deux points les plus éloignés du nuage.

Illustration de la distance de référence d_r.

 

Le choix du jeu de paramètres (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 suffisants, tout en gardant à l’esprit que la durée de calcul est proportionnelle à M. De plus, s’il l’on choisit \beta>1, l’influence des points lointains est rapidement amortie lorsque leur distance au point de référence est supérieure à 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.

Exemple de champ projeté sur une grille régulière selon la méthode des plus proches voisins.

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}

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}

d_k est la distance entre le point (x,y) et le point (x_k,y_k) soit { 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 tel que :
  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

\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^*], où 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^* et 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

Avec d_i la distance entre les points i et 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) du nuage 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.

Temps d’exécution pour trois méthodes de minimisation (tiré de (Klinyam 2019)).