Journal de bord de thèse

Josselin

2019-10-07

Relation de dispersion pour Vlasov-Poisson

Vlasov-Poisson : $\begin{cases} \partial_t f + v\partial_x f + E\partial_v f = 0 \\ \partial_x E = \int_{\mathbb{R}} f\,\mathrm{d}v - 1 \end{cases}$

Linéarisation autour d’un état d’équilibre $f(t,x,v) = f_\text{eq}(v) + \varepsilon g(t,x,v)$ et $E(t,x) = \varepsilon e(t,x)$ : $\begin{cases} \partial_t g + v\partial_x g + e\mathrm{d}_vf_\text{eq} = 0 \\ \partial_x e = \int_{\mathbb{R}}g\,\mathrm{d}v \end{cases}$ Fourier-Laplace (condition initiale nulle) : $\partial_t\cdot \rightarrow -i\omega\cdot \qquad \partial_x\cdot \rightarrow ik\cdot$ d’où $\begin{cases} ik\left(v-\frac{\omega}{k}\right)g = - e\mathrm{d}_vf_\text{eq} \\ ike = \int g\,\mathrm{d}v \end{cases}$ La première équation nous donne : $g = \frac{ie\mathrm{d}_vf_\text{eq}}{k\left(v-\frac{\omega}{k}\right)}$ La seconde donne alors : $e = -\frac{i}{k}\int \frac{ie\mathrm{d}_vf_\text{eq}}{k\left(v-\frac{\omega}{k}\right)}\,\mathrm{d}v$ soit : $e\left(1-\frac{1}{k^2}\int\frac{\mathrm{d}_vf_\text{eq}}{v-\frac{\omega}{k}}\,\mathrm{d}v\right) = 0$ On note alors : $\boxed{ D(\omega,k) = 1-\frac{1}{k^2}\int_{\mathbb{R}}\frac{\mathrm{d}_vf_\text{eq}}{v-\frac{\omega}{k}}\,\mathrm{d}v }$

Calcul de l’intégrale pour une maxwellienne :

Étude de l’intégrale dans le cas $f_\text{eq}(v)=\mathcal{M}_{\rho,u,T}(v) = \frac{\rho}{\sqrt{2\pi T}}e^{-\frac{|v-u|^2}{2T}}$ : $\begin{aligned} \int_{\mathbb{R}}\frac{\mathrm{d}_v\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \int\frac{-\frac{(v-u)}{T}\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v \\ &= -\frac{1}{T}\int\frac{(v-\beta+\beta-u)\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v \\ &= -\frac{\rho}{T}-\frac{\beta-u}{T}\int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v \end{aligned}$ Étudions plus particulièrement l’intégrale : $\int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v$ : $\begin{aligned} \int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \frac{\rho}{\sqrt{2\pi T}}\int\frac{e^{-\frac{|v-u|^2}{2T}}}{v-\beta}\,\mathrm{d}v \\ &= \frac{\rho}{\sqrt{2\pi T}}\int\frac{e^{-\left(\frac{v-u}{\sqrt{2T}}\right)^2}}{v-u+u-\beta}\,\mathrm{d}v \end{aligned}$ on effectue le changement de variable $z=\frac{v-u}{\sqrt{2T}}$ , $\mathrm{d}z = \frac{\mathrm{d}v}{\sqrt{2T}}$ : $\begin{aligned} \int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \frac{\rho}{\sqrt{2\pi T}}\int\frac{e^{-z^2}}{\sqrt{2T}z+u-\beta}\sqrt{2T}\,\mathrm{d}z\\ &= \frac{\rho}{\sqrt{2\pi T}}\int\frac{e^{-z^2}}{z+\frac{u-\beta}{\sqrt{2T}}}\mathrm{d}z \end{aligned}$ On pose $\mathcal{Z}:\xi\mapsto\frac{1}{\sqrt{\pi}}\int\frac{e^{-z^2}}{z-\xi}\,\mathrm{d}z$ : $\int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v = \frac{\rho}{\sqrt{2T}}\mathcal{Z}\left(\frac{\beta-u}{\sqrt{2T}}\right)$ D’où : $\int\frac{\mathrm{d}_v\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v = -\frac{\rho}{T} - \frac{\rho(\beta-u)}{T\sqrt{2T}}\mathcal{Z}\left(\frac{\beta-u}{\sqrt{2T}}\right)$ soit : $\boxed{\int_{\mathbb{R}}\frac{\mathrm{d}_v\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v = -\frac{\rho}{T}\left(1+\frac{\beta-u}{\sqrt{2T}}\mathcal{Z}\left(\frac{\beta-u}{\sqrt{2T}}\right)\right) }$

Calcul de l’intégrale pour une distribution de Dirac

Étude de l’intégrale dans le cas $f_\text{eq}(v)= \delta_u(v)$ :

$\begin{aligned} \int_{\mathbb{R}}\frac{\mathrm{d}_v\delta_u(v)}{v-\beta}\,\mathrm{d}v &= -\int\delta_u(v)\left(\frac{1}{v-\beta}\right)'\,\mathrm{d}v \\ &= -\int\delta_u(v)\frac{1}{(v-\beta)^2} \,\mathrm{d}v \end{aligned}$ soit : $\boxed{ \int_{\mathbb{R}}\frac{\mathrm{d}_v\delta_u(v)}{v-\beta}\,\mathrm{d}v = -\frac{1}{(u-\beta)^2} }$

Test numérique

On étudie le cas d’un triple bump : $f_\text{eq} = \mathcal{M}_{1-\alpha,0,T_c} + \mathcal{M}_{^\alpha/_2,u,1} + \mathcal{M}_{^\alpha/_2,-u,1}$ la relation de dispersion s’écrit alors : $\boxed{ \begin{aligned} D(\omega,k) = 1-\frac{1}{k^2}\left[\vphantom{\frac{\frac{}{}}{\sqrt{}}}\right. &-\frac{1-\alpha}{T_c}\left(1+\frac{\omega}{k\sqrt{2T_c}}\mathcal{Z}\left(\frac{\omega}{k\sqrt{2T_c}}\right) \right) \\ & -\frac{\alpha}{2}\left(1+\frac{\frac{\omega}{k}-u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}-u}{\sqrt{2}}\right)\right) \\ & -\frac{\alpha}{2}\left(1+\frac{\frac{\omega}{k}+u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}+u}{\sqrt{2}}\right)\right) \left.\vphantom{\frac{\frac{}{}}{\sqrt{}}}\right] \end{aligned} }$

On souhaite comparer ce test à celui de la modélisation hybride avec une distribution de Dirac pour représenter les particules froides :

$f_\text{eq} = (1-\alpha)\delta_0 + \mathcal{M}_{^\alpha/_2,u,1} + \mathcal{M}_{^\alpha/_2,-u,1}$ Ce qui correspond à $T_c \to 0$ . La relation de dispersion s’écrit alors : $\boxed{ \begin{aligned} D(\omega,k) = 1 - \frac{1}{k^2}\left[\vphantom{\frac{\frac{}{}}{\sqrt{}}}\right. & -(1-\alpha)\left(\frac{k}{\omega}\right)^2 \\ & -\frac{\alpha}{2}\left(1+\frac{\frac{\omega}{k}-u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}-u}{\sqrt{2}}\right)\right) \\ & -\frac{\alpha}{2}\left(1+\frac{\frac{\omega}{k}+u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}+u}{\sqrt{2}}\right)\right) \left.\vphantom{\frac{\frac{}{}}{\sqrt{}}}\right] \end{aligned} }$ Ce qui revient à dire que : $1+z\mathcal{Z}(z) \overset{z\to+\infty}{\approx} \frac{2}{z^2}$ avec $z\in\mathbb{C}$ n’évoluant que selon un axe parallèle à la droite des réelles, c’est-à-dire que seulement sa partie réelle tend vers l’infini.

2019-10-25

Lawson pour VPHL

On s’intéresse au problème de Vlasov-Poisson hybride linéarisé (VPHL) : $\begin{cases} \partial_t u_c = E \\ \partial_t E = -\rho_cu_c - \int vf_h\,\mathrm{d}v \\ \partial_t \hat{f}_h = -ikv\hat{f}_h - \widehat{E\partial_v f_h} \end{cases}$ par commodité, nous noterons $f_h = \hat{f}_h$ . Ce problème se réécrit : $\partial_t\underbrace{\begin{pmatrix}u_c\\E\\f_h\end{pmatrix}}_{U} = \underbrace{\begin{pmatrix}0 & 1 & 0 \\ -\rho_c & 0 & 0 \\ 0 & 0 & -ikv \end{pmatrix}}_{A}\begin{pmatrix}u_c\\E\\f_h\end{pmatrix} + \underbrace{\begin{pmatrix}0 \\ -\int vf_h\,\mathrm{d}v \\ -\widehat{E\partial_v f_h} \end{pmatrix}}_{N(U)}$ $\partial_t U = AU + N(U)$ Pour introduire une méthode de Lawson, on utilise le changement de variable suivant $V=e^{-tA}U$ , l’exponentielle de la matrice $A$ peut se calculer et vaut : $e^{-tA} = \begin{pmatrix} \cos(\sqrt{\rho_c}t) & -\frac{\sin(\sqrt{\rho_c}t)}{\sqrt{\rho_c}} & 0\\ \sqrt{\rho_c}\sin(\sqrt{\rho_c}t) & \cos(\sqrt{\rho_c}t) & 0 \\ 0 & 0 & e^{ikvt} \end{pmatrix}$ On calcule maintenant la valeur de $\partial_t V$ : $\partial_t V = -Ae^{-tA}U + e^{-tA}\partial_t U$ soit : $\partial_t V = -Ae^{-tA}U + e^{-tA}AU + e^{-tA}N(U)$

En écrivant $e^{tA}$ sous forme de série : $e^{tA }= \sum_{k\in\mathbb{N}}\frac{t^kA^k}{k!}$ il devient évident que $Ae^{tA} = e^{tA}A$ , i.e. que $A$ et $e^{tA}$ commutent.

d’où : $\partial_t V = e^{-tA}N(U) = e^{-tA}N(e^{tA}V) = \tilde{N}(t,V)$ On obtient une forme propice à l’utilisation d’une méthode RK( $s$ , $n$ ), ce qui permettra d’exprimer le schéma de Lawson associé. Pour l’exemple, nous utiliserons la méthode RK( $3$ , $3$ ) dit de Shu-Osher. Les méthodes dites de Shu-Osher sont des méthodes de type Runge-Kutta réécrite pour réutiliser les résultats précédents de n’effectuer qu’un seul appel à la fonction non linéaire par étage, ces méthodes sont intéressantes lorsque le coût de calcul est plus long que le temps d’accès à la mémoire (et que cette dernière n’est pas limitée). Dans ce cas, ce qui nous intéresse est l’obtention de la fonction de stabilité, et celle-ci est plus simple à obtenir à la main dans le cas d’une méthode dite de Shu-Osher. Pour les notations, nous écrivons cette méthode pour résoudre un problème du type $\dot{u} = L(t,u)$ : $\begin{aligned} u^{(1)} &= u^n + \Delta t L(t^n,u^n) \\ u^{(2)} &= \frac{3}{4}u^n + \frac{1}{4}u^{(1)} + \frac{1}{4}\Delta t L(t^n+\Delta t,u^{(1)}) \\ u^{n+1} &= \frac{1}{3}u^n + \frac{2}{3}u^{(2)} + \frac{2}{3}\Delta t L(t^n+\frac{\Delta t}{2},u^{(2)}) \\ \end{aligned}$ Nous écrivons maintenant cette méthode pour le problème $\partial_t V = e^{-tA}N(e^{tA}V)$ : $\begin{aligned} V^{(1)} &= V^n + \Delta t e^{-t^nA}N(e^{t^nA}V^n) \\ V^{(2)} &= \frac{3}{4}V^n + \frac{1}{4}V^{(1)} + \frac{1}{4}\Delta t e^{-(t^n+\Delta t)A}N(e^{(t^n+\Delta t)A}V^{(1)}) \\ V^{n+1} &= \frac{1}{3}V^n + \frac{2}{3}V^{(2)} + \frac{2}{3}\Delta t e^{-(t^n+\frac{\Delta t}{2})A}N(e^{(t^n+\frac{\Delta t}{2})A}V^{(2)}) \\ \end{aligned}$ Dans un premier temps nous approximons la non linéarité $N(U)$ par : $N(U) \approx \begin{pmatrix}0&0&0\\0&0&0\\0&0&\lambda\\\end{pmatrix} = \Lambda$ $\Lambda$ commute avec $e^{tA}$ , donc la méthode RK( $3$ , $3$ ) dit de Shu-Osher se réécrit : $\begin{aligned} V^{(1)} &= V^n + \Delta t \Lambda V^n = (I+\Delta t \Lambda)V^n\\ V^{(2)} &= \frac{3}{4}V^n + \frac{1}{4}V^{(1)} + \frac{1}{4}\Delta t \Lambda V^{(1)} = \left( \frac{3}{4}I + \frac{1}{4}(I+\Delta t \Lambda)^2 \right)V^n\\ V^{n+1} &= \frac{1}{3}V^n + \frac{2}{3}V^{(2)} + \frac{2}{3}\Delta t \Lambda V^{(2)} = p_{RK(3,3)}(\Delta t\Lambda)V^n \\ \end{aligned}$ avec $p_{RK(3,3)}$ le polynôme défini par : $p_{RK(3,3)}(z) = 1+z+\frac{z^2}{2}+\frac{z^3}{6}$ On cherche maintenant à exprimer ceci en fonction de notre inconnue $U$ : $e^{(t^n+\Delta t)A}U^{n+1} = p_{RK(3,3)}(\Delta t\Lambda)e^{t^nA}U^n$ soit : $U^{n+1} = e^{-(t^n+\Delta t)A}p_{RK(3,3)}(\Delta t\Lambda)e^{t^nA}U^n$ Maintenant il reste à savoir si $e^{tA}$ et $p_{RK(3,3)}(\Delta t\Lambda)$ commutent, on sait déjà que $I$ et $e^{tA}$ commutent, $\Lambda$ et $e^{tA}$ aussi, calculons $\Lambda^2$ : $\Lambda^2 = \begin{pmatrix}0&0&0\\0&0&0\\0&0&\lambda^2\\\end{pmatrix}$ donc $\Lambda^2$ et $e^{tA}$ commutent, d’où : $U^{n+1} = p_{RK(3,3)}(\Delta t\Lambda)e^{\Delta tA}U^n$ On peut calculer la norme de $e^{tA}$ avec notre matrice $A$ et on obtient $\left|\left|e^{tA}\right|\right|_2 = 1$ , $\forall t$ . Comme pour le cas cinétique K, l’étude de stabilité de la méthode de Lawson revient à l’étude de la stabilité de la méthode RK( $s$ , $n$ ) sous-jacente. Le calcul a été effectué pour une méthode RK( $3$ , $3$ ), mais se généralise facilement puisque $\Lambda$ commute avec $e^{tA}$ , $\forall t\in\mathbb{R}$ . De façon plus générale nous obtenons donc le résultat : $p_{LRK(s,n)}(z) = p_{RK(s,n)}(z)e^{\Delta tA}$ où $p_{LRK(s,n)}$ est la fonction de stabilité de la méthode de Lawson associée à la méthode Runge-Kutta RK( $s$ , $n$ ) de fonction de stabilité $p_{RK(s,n)}$ .

Intéressons nous maintenant à la linéarisation de $N(U)$ suivante : $N(U) \approx \begin{pmatrix}0&0&0\\0&0&\nu\\0&0&\lambda\\\end{pmatrix} = \dot{\Lambda}$ On remarque alors que $\dot{\Lambda}$ et $e^{tA}$ ne commutent pas, on a pas l’égalité de fonction de stabilité : $p_{LRK(s,n)}(\Delta t\dot{\Lambda}) \neq p_{RK(s,n)}(\Delta t\dot{\Lambda})e^{\Delta tA}$

La linéarisation permettant d’obtenir le même résultat que dans le cas scalaire peut s’interpréter d’un point de vue physique comme : $J_c \gg J_h$ i.e. le courant des particules chaudes est négligeable devant celui des particules froides ; ou, comme le cas scalaire, cela revient à ne considérer que le transport 2D pour étudier la stabilité des méthodes numériques sur l’équation de Vlasov.

On considère donc que la seule source de CFL est le transport en $v$ des particules chaudes $f_h$ .

Dépôt sur HAL

J’ai effectué le dépôt sur HAL de l’article avec Lukas et Nicolas Exponential methods for solving hyperbolic problems with application to kinetic equations. J’ai enfin compris ce qui bloquait dans le dépôt simultané sur arXiv, toutes les images doivent être au format jpeg, png ou pdf ; une image était au format eps (autorisé uniquement pour TeX ou LaTeX et non pdfLaTeX…) J’ai effectué les modifications requises et j’attends la validation du dépôt par HAL. Nicolas évoquait de le soumettre à JCP.

Actualisation : le dépôt est validé : hal-02321916, version 1.

2019-10-30

Résolution de VPHL

Vlasov-Poisson hybride linéarisé (VPHL)

$\begin{cases} \partial_tu_c = E \\ \partial_tE = -\rho_cu_c - \int vf_h\,\mathrm{d}v \\ \partial_t\hat{f}_h + ikv\hat{f}_h + \widehat{E\partial_vf_h} = 0 \\ \end{cases}$

L’équation de Poisson : $ik\hat{E}(t,k) = \int\hat{f}_h\,\mathrm{d}v$ , est-elle vérifiée ?

$\partial_t\begin{pmatrix}u_c \\ E \\ \hat{f}_h \end{pmatrix} + \begin{pmatrix} 0 & -1 & 0 \\ \rho_c & 0 & 0 \\ 0 & 0 & ikv \end{pmatrix}\begin{pmatrix}u_c \\ E \\ \hat{f}_h\end{pmatrix} + \begin{pmatrix} 0 \\ \int f_h\,\mathrm{d}v \\ \widehat{E\partial_vf_h} \end{pmatrix} = 0$ de la forme : $\partial_t U + AU + N(U) = 0$ une forme propice à l’utilisation des méthodes de Lawson.

Splitting en 3 étapes

$\varphi^{[a]}_{\Delta t}$ : $\begin{cases}\partial_t f_h + v\partial_x f_h =0 \\ \partial_t u_c = 0 \\ \partial_t E = -\int vf_h\,\mathrm{d}v\end{cases}$
$\varphi^{[b]}_{\Delta t}$ : $\begin{cases}\partial_t f_h + E\partial_v f_h = 0 \\ \partial_t u_c = E \\ \partial_t E = 0 \end{cases}$
$\varphi^{[c]}_{\Delta t}$ : $\begin{cases} \partial_t f_h = 0 \\ \partial_t u_c = 0 \\ \partial_t E = -\rho_cu_c \end{cases}$

$U^{n+1} = \varphi^{[a]}_{\Delta t} \circ \varphi^{[b]}_{\Delta t} \circ \varphi^{[c]}_{\Delta t} (U^n)$

avec dans l’ordre :

$\varphi^{[c]}_{\Delta t}(U^n) = \begin{pmatrix} f_h^{(1)} \\ u_c^{(1)} \\ E^{(1)} \end{pmatrix} = \begin{pmatrix} f_h^n \\ u_c^n \\ E^n - \Delta t \rho_cu_c^n \end{pmatrix} = U^{(1)}$
$\varphi^{[b]}_{\Delta t}(U^{(1)}) = \begin{pmatrix} f_h^{(2)} \\ u_c^{(2)} \\ E^{(2)} \end{pmatrix} = \begin{pmatrix}f_h^{(1)}\left(x,v-\Delta t E^{(1)}\right) \\ u_c^{(1)} + \Delta t E^{(1)} \\ E^{(1)} \end{pmatrix} = U^{(2)}$
$\varphi^{[a]}_{\Delta t}(U^{(2)}) = \begin{pmatrix} \hat{f}_h^{n+1} \\ u_c^{n+1} \\ \hat{E}^{n+1} \end{pmatrix} = \begin{pmatrix}e^{-ikv\Delta t}\hat{f}_h^{(2)} \\ u_c^{(2)} \\ \hat{E}^{(2)} - \frac{i}{k}\int \hat{f}_h^{n+1}\,\mathrm{d}v + \frac{i}{k}\int \hat{f}_h^{(2)}\,\mathrm{d}v \end{pmatrix} = U^{n+1}$ La formulation de la dernière équation d’Ampère, le calcul de $\hat{E}^{n+1}$ , permet de vérifier simultanément l’équation de Poisson.

2019-11-04

Le schéma de Lawson pour la résolution de VPHL n’a pas lieu de changer. Le schéma numérique fait intervenir des exponentielles qui ne dépendent que du pas de temps. Dans le cas scalaire, il était simple de recalculer ces valeurs, dans le cas matriciel il est intéressant d’assembler les matrices en amont. Le calcul du schéma fait intervenir la matrice $e^{c_i\Delta t A}$ : $e^{c_i\Delta t A} = \begin{pmatrix} \cos(c_i\Delta t) & -\sin(c_i\Delta t) & 0 \\ \sin(c_i\Delta t) & \cos(c_i\Delta t) & 0 \\ 0 & 0 & e^{i c_{i} \Delta t k v} \end{pmatrix}$

sympy ne simplifie pas comme il faut les termes lorsque $c_i$ est un réel quelconque, mais pas de problème lorsque $c_i>0$ ou $c_i<0$ , de même que le cas $c_i=0$ qui ne se rencontre jamais. Si $c_i$ est quelconque il travaille sur $|c_i|$ et fait intervenir des $\frac{|c_i|}{c_i}$ pour jouer sur la parité des fonctions sinus et cosinus.

Une fonction calculant l’exponentielle $e^{c_i\Delta t A}$ est intéressant (plutôt qu’espérer que ublas de boost me calcule à chaque itération l’exponentielle), il n’est pas envisageable d’avoir un assemblage définitif de la matrice en début de simulation car :

$e^{i c_i\Delta t kv}$ peut prendre beaucoup de valeurs, mais ce n’est qu’un terme diagonal à modifier.
On souhaite faire évoluer nos simulations vers des méthodes à pas de temps adaptatif, donc $\Delta t$ n’est pas constant, donc le bloc diagonal ne peut pas s’assembler définitivement en début de simulation.

La matrice ne peut être assemblée définitivement, mais on peut conserver la même zone mémoire tout au long de la simulation. La librairie ublas de boost:: en particulier d’effectue pas d’affectation par défaut des valeurs lors de l’allocation de la matrice ; conserver la même zone mémoire permet de ne pas affecter 4 zéros à chaque fois que cette matrice est nécessaire.

Dans la pratique il est même inutile d’effectuer l’assemblage de la matrice, il sera bien plus performant d’effectue directement l’assemblage du produit matrice-vecteur dans une fonction.

2019-11-06

J’ai implémenté le WENO-SL (semi-Lagrangien) dans un code de test 1Dx, de transport à vitesse à vitesse constante, etc. Le WENO-SL est décrit dans Qiu:2011 ou Qiu:2010.

Simulation de : $u_t + u_x = 0$ avec $T_f=0.25$ , $\Delta x=\frac{1}{100}$ , $x\in[0,1]$ (périodique), RK(4,4) comme intégrateur en temps et 5 condtions initiales différentes : une période de c

Comparaison des différentes méthodes WENO

On remarque que le schéma WENO3 diffuse un peu plus que les autres WENO d’ordre 5 : WENO5 (qui est le WENO standard WENO-JS), WENO-M et WENO-Z, on remarque que le WENO-B de Banks oscille (constatation déjà faite). Le WENO-SL diffuse encore plus que les autres. Je trouve ce résultat relativement étrange, je n’épargne donc pas une erreur dans mon implémentation.

La présentation du WENO-SL est faite avec comme intégrateur en temps du Euler explicite. Puisqu’il est nécessaire de faire entrer du $\Delta t$ dans la méthode WENO-SL, je ne sais pas si cela se mêle correctement avec un autre intégrateur en temps (comme RK(4,4) utilisé dans cette simulation).

Correction/compréhension de WENO-SL

En réalité, dans la construction de la méthode WENO-SL, la discrétisation s’effectue simultanément en temps et en espace, il n’y a pas d’intégrateur en temps à proprement parlé, dans mon code de simulation de test 1D cela revient à utiliser une méthode d’Euler explicite et de diviser l’approximation $u_x$ par $\Delta t$ . On obtient alors des résultats plus cohérent (correction effectuer jeudi à 14h).

On remarque alors que le WENO-SL donne un très bon résultat sur une discontinuité, bien meilleur que WENO-Z (qui est une optimisation du WENO5 classique WENO-JS pour minimiser la diminution d’ordre lors d’une discontinuité). Ce résultat est trompeur car le transport s’effectue à vitesse 1, donc une méthode semi-lagrangienne est considérée comme exacte. Je n’ai pas effectué de mesure d’ordre sur ce schéma (ordre en espace et en temps du coup ?).

Résolution de : $u_t + \frac{2}{3}u_x = 0$

Avec les mêmes paramètres numériques.

Dans ces conditions, WENO-SL fait aussi bien que WENO-Z. Un test similaire a été effectué avec une vitesse négative $-\frac{2}{3}$ , et on retrouve bien les mêmes résultats. Cela permet de vérifier l’implémentation des flux, puisque dépendant du signe de la vitesse.

2019-11-08

Compte-rendu avec Nicolas

J’ai commencé à étudier l’ordre en temps (normalement exact) et en espace (normalement d’ordre 5) du WENO-SL, mais mon résultat préliminaire me donnait un schéma d’ordre 1 en temps et d’ordre 5 en espace (avec une erreur qui sature très vite à $\log(\Delta t)$ . Mon erreur dans le code provient (vraisemblablement) du fait que je choisis des $\Delta t$ qui dépendent de $\Delta x$ , ce qui implique un nombre d’itération $\frac{T_f}{\Delta t}$ non entier. Je dois donc faire en sorte, pour les 2 mesures d’ordre, d’avoir $\frac{T_f}{\Delta t}\in\mathbb{N}$ .

On a également discuté du schéma de Lawson pour VPHL. J’avais laissé tombé les calculs pour la non-linéarité $\dot{\Lambda}$ avec un terme extra-diagonal, car je ne trouvais pas de forme simple car les matrices ne commutent pas. En détaillant chaque terme des matrices avec Nicolas, il semblerait possible d’avoir une forme exploitable pour l’étude généralisée des ordres supérieur à 1.

On a longuement discuté de WENO-SL et de son interprétation comme une interpolation semi-lagrangienne. La version de WENO-SL par Qiu continue de faire apparaître des flux et non une interpolation. La version semi-lagrangienne par Francis Filbet semble plus sympathique, mais utilise des polynômes de Hermite, et donc à un moment une approximation de la dérivée de $u$ . Je dois me replonger un peu plus dans les papiers de Shu sur la construction de WENO, pour être certain de la construction des flux et l’interprétation que Nicolas veut en faire comme une interpolation (pour retrouver une méthode semi-lagrangienne).

TODO:

corriger le script pour obtenir l’ordre en $t$ et $x$

calculer les valeurs propres de $e^{tA}$ pour être certain du résultat donné par sympy : $\left|\left|e^{tA} \right|\right|_2 = 1$

comparer les polynômes d’interpolation de WENO (et comprendre comment on les calcule), et comparer la version linéarisé avec le polynôme d’interpolation du semi-lagrangien. Voir papier de Francis sur du HWENO (pour WENO semi-lagrangien avec polynôme de Hermite).

Avancement sur les petits trucs à vérifier

L’erreur est bien constante avec $\Delta t$ , et on atteint l’erreur numérique. En $x$ on retrouve bien la pente de 5, et n’est pas saturé par $\Delta t$ .
Les valeurs propres de $e^{-tA}$ sont $e^{-itkv}$ et $\cos(t) \pm i\sin(t)$ . Ceci nous donne bien des valeurs propres de module 1.

Ordre de WENO-SL en \Delta t, le schéma est exacte en temps — Ordre de WENO-SL en $\Delta t$ , le schéma est *exacte* en temps

Ordre de WENO-SL en \Delta x, le schéma est d’ordre 5 en espace — Ordre de WENO-SL en $\Delta x$ , le schéma est d’ordre 5 en espace

2019-11-15

Cette semaine j’ai surtout effectué de l’enseignement et de la correction de copies. J’ai aussi été confronté à sympy et la difficulté de lui faire effectuer des simplifications qui me semblaient évidentes. J’ai réussi à obtenir le polynôme de stabilité d’une méthode de Lawson basée sur RK(3,3) pour le problème suivant : $\partial_t U + AU + N(U)$ en linéarisant $N(U)$ avec une matrice $\dot{\Lambda}$ définie par : $\dot{\Lambda} = \begin{pmatrix}0&0&0 \\ 0&0&\mu \\ 0&0&\lambda \end{pmatrix}$

On obtient alors de manière un peu brute avec sympy : $U^{n+1} = \begin{bmatrix}1 & 0 & - \frac{\Delta t \mu \left(\Delta t^{2} \lambda^{2} e^{i \Delta t k v} \sin{\left(\frac{\Delta t}{2} + t^{n} \right)} + \Delta t \lambda e^{\frac{i \Delta t k v}{2}} \sin{\left(\Delta t + t^{n} \right)} + 2 \Delta t \lambda e^{i \Delta t k v} \sin{\left(\frac{\Delta t}{2} + t^{n} \right)} + e^{\frac{3 i \Delta t k v}{2}} \sin{\left(t^{n} \right)} + e^{\frac{i \Delta t k v}{2}} \sin{\left(\Delta t + t^{n} \right)} + 4 e^{i \Delta t k v} \sin{\left(\frac{\Delta t}{2} + t^{n} \right)}\right) e^{- i k v \left(\frac{3 \Delta t}{2} + t^{n}\right)}}{6}\\0 & 1 & \frac{\Delta t \mu \left(\Delta t^{2} \lambda^{2} e^{i \Delta t k v} \cos{\left(\frac{\Delta t}{2} + t^{n} \right)} + \Delta t \lambda e^{\frac{i \Delta t k v}{2}} \cos{\left(\Delta t + t^{n} \right)} + 2 \Delta t \lambda e^{i \Delta t k v} \cos{\left(\frac{\Delta t}{2} + t^{n} \right)} + e^{\frac{3 i \Delta t k v}{2}} \cos{\left(t^{n} \right)} + e^{\frac{i \Delta t k v}{2}} \cos{\left(\Delta t + t^{n} \right)} + 4 e^{i \Delta t k v} \cos{\left(\frac{\Delta t}{2} + t^{n} \right)}\right) e^{- i k v \left(\frac{3 \Delta t}{2} + t^{n}\right)}}{6}\\0 & 0 & \frac{\Delta t^{3} \lambda^{3}}{6} + \frac{\Delta t^{2} \lambda^{2}}{2} + \Delta t \lambda + 1\end{bmatrix}U^n$ Cette matrice est le fruit du polynôme en $\Delta t\dot{\Lambda}$ : $\frac{\Delta t^{3} e^{- \frac{\Delta t A}{2} - t^{n} A} \dot{\Lambda} e^{- \frac{\Delta t A}{2}} \dot{\Lambda} e^{\Delta t A} \dot{\Lambda} e^{t^{n} A}}{6} + \frac{\Delta t^{2} e^{- \Delta t A - t^{n} A} \dot{\Lambda} e^{\Delta t A} \dot{\Lambda} e^{t^{n} A}}{6} + \frac{\Delta t^{2} e^{- \frac{\Delta t A}{2} - t^{n} A} \dot{\Lambda} e^{- \frac{\Delta t A}{2}} \dot{\Lambda} e^{\Delta t A + t^{n} A}}{6} + \frac{\Delta t^{2} e^{- \frac{\Delta t A}{2} - t^{n} A} \dot{\Lambda} e^{\frac{\Delta t A}{2}} \dot{\Lambda} e^{t^{n} A}}{6} + \frac{\Delta t e^{- t^{n} A} \dot{\Lambda} e^{t^{n} A}}{6} + \frac{\Delta t e^{- \Delta t A - t^{n} A} \dot{\Lambda} e^{\Delta t A + t^{n} A}}{6} + \frac{2 \Delta t e^{- \frac{\Delta t A}{2} - t^{n} A} \dot{\Lambda} e^{\frac{\Delta t A}{2} + t^{n} A}}{3} + 1$ Ce qui est important c’est que uniquement la dernière colonne est remplie. On a à chaque fois quelque chose de homogène à une puissance de $\Delta t\lambda$ : $\Delta t^3 \lambda^2\mu$ , $\Delta t^2 \lambda\mu$ , $\Delta t\mu$ . Et pour $\mu=0$ on retrouve bien la même matrice que dans le cas déjà traité précédemment.

J’ai passé du temps avec sympy pour conserver des cosinus et sinus dans la forme finale (ceux-ci étaient écrits en complexe sous leur forme polaire), et ainsi qu’écrire ceci de manière à pouvoir tester le plus simplement possible avec d’autres schémas. Puisque je ne sais pas trop quoi conclure avec cette matrice, je n’ai pas cherché encore à effectuer ce travail sur d’autres schémas.

Pour information, pour obtenir ce résultat, j’ai défini une fonction permettant de calcul $e^{tA}$ , et je connais les instants $t$ où cette fonction sera évaluée dans le schéma (en prennant les notations d’un tableau de Butcher, c’est en $e^{\pm(t^n+c_i\Delta t)A}$ . Puis je demande à sympy de me substituer toutes ces exponentielles de matrices et $\dot{\Lambda}$ simultanément pour obtenir ce résultat.

def m_exp(x):
  t = (x/A).simplify()
  M = sp.Matrix([[sp.cos(t),sp.sin(t),0],[-sp.sin(t),sp.cos(t),0],[0,0,sp.exp(-sp.I*k*v*t)]])
  return M
list_t = sum([[t,-t] for t in [tn,tn+dt,tn+dt/2,dt,dt/2]],[])
list_subs = [(sp.exp(t*A).expand(),m_exp(t*A)) for t in list_t]+[(L,mL)]

p.subs(list_subs,simultaneous=True) # p est notre polynôme en $\Delta t\dot{\Lambda}$

J’ai pu testé avec un autre schéma RK(3,3) que le schéma RK(3,3) de Shu-Osher, défini par : $\begin{aligned} u^{(1)} & = u^n + \frac{1}{2}\Delta t L(t^n,u^n) \\ u^{(2)} & = 3u^n - 2u^{(1)} + 2\Delta t L(t^n+\frac{1}{2}\Delta t, u^{(1)}) \\ u^{n+1} & = -\frac{1}{3}u^n + u^{(1)} + \frac{1}{3}u^{(2)} + \frac{1}{6}\Delta t L(t^n+\Delta t, u^{(2)}) \end{aligned}$ Et la fonction de stabilité du schéma de Lawson induit, avec une non-linéarité $\dot{\Lambda}$ (qui ne commute pas avec $e^{tA}$ ) est différente, et donc la matrice de stabilité varie aussi.

2019-11-20

Prenons le schéma de Lawson induit par le schéma RK(3,3) de Shu-Osher, pour la résolution de l’équation : $\partial_t U = AU + N(U)$ avec $U = \begin{pmatrix}u_v & E & \hat{f}_h\end{pmatrix}^\textsf{T}$ , et : $A = \begin{pmatrix}0 & 1 & 0 \\ -\rho_c & 0 & 0 \\ 0 & 0 & -ikv \end{pmatrix} \qquad N(U) = \begin{pmatrix}0 \\ -\int vf_h\,\mathrm{d}v \\ - \widehat{E\partial_v f_h} \end{pmatrix}$ Avec ces notations, le schéma s’écrit : $\begin{aligned} U^{(1)} & = e^{\Delta t A}U^n + \Delta t e^{\Delta t A} N(U^n) \\ U^{(2)} & = \frac{3}{4}e^{\frac{\Delta t}{2}A}U^n + \frac{1}{4}e^{-\frac{\Delta t}{2}A}U^{(1)} + \frac{1}{4}\Delta te^{-\frac{\Delta t}{2}A}N(U^{(1)}) \\ U^{n+1} & = \frac{1}{3}e^{\Delta t L}U^n + \frac{2}{3}e^{\frac{\Delta t}{2}A}U^{(2)} + \frac{2}{3}\Delta te^{\frac{\Delta t}{2}A}N(U^{(2)}) \end{aligned}$

Détaillons maintenant le premier étage : $\begin{pmatrix} u_c^{(1)} \\ E^{(1)} \\ \hat{f}_h^{(1)} \end{pmatrix} = \begin{pmatrix} \cos(\sqrt{\rho_c}\Delta t) & \frac{\sin(\sqrt{\rho_c}\Delta t)}{\sqrt{\rho_c}} & 0 \\ -\sqrt{\rho_c}\sin(\sqrt{\rho_c}\Delta t) & \cos(\sqrt{\rho_c}\Delta t) & 0 \\ 0 & 0 & e^{-ikv\Delta t} \end{pmatrix} \begin{pmatrix} u_c^n \\ E^n \\ \hat{f}_h^n \end{pmatrix} + \Delta t \begin{pmatrix} \cos(\sqrt{\rho_c}\Delta t) & \frac{\sin(\sqrt{\rho_c}\Delta t)}{\sqrt{\rho_c}} & 0 \\ -\sqrt{\rho_c}\sin(\sqrt{\rho_c}\Delta t) & \cos(\sqrt{\rho_c}\Delta t) & 0 \\ 0 & 0 & e^{-ikv\Delta t} \end{pmatrix} \begin{pmatrix} 0 \\ -\int vf_h^n\,\mathrm{d}v \\ -\widehat{E^n\partial_vf_h^n} \end{pmatrix}$ avec $({u_c^n}_i)_i$ , $(E^n_i)_i$ , et $({\hat{f}_h^n}_{\iota,k})_{\iota,k}$ avec $i$ indice d’espace $x_i = i\Delta x$ , $k$ indice de vitesse $v_k = k\Delta v$ et $\iota$ indice de phase $\kappa_\iota = \frac{2\pi \iota}{L}$ . Pour la transformée de Fourier nous utilisons une FFT (Fast Fourier Transform), et l’algorithme impose un même nombre de cellules en espace que de nombres d’onde, il est donc possible d’utiliser le même indice pour parcourir l’espace $x_i = i\Delta x$ et la phase $\kappa_\iota = \frac{2\pi\iota}{L}$ .

Nous obtenons alors par exemple pour le premier étage l’algorithme suivant :

Calcul de $\left(E^n\partial_vf_h^n\right)_{j,k}$ avec WENO ou autre, et calcul de $(\int vf_h^n\,\mathrm{d}v)_j$ (ceci nécessite d’effectuer une FFT inverse de $\hat{f}_h^n$ , et il s’agit finalement le la seule étape où $f_h$ est nécessaire).
Calcul des $(u_c^{(1)})_j$ et $(E^{(1)})_j$ : $\begin{aligned} {u_c}^{(1)}_j & = {u_c}^n_j\sqrt{\rho_c}\cos(\sqrt{\rho_c}\Delta t) - E^n_j\frac{\sin(\sqrt{\rho_c}\Delta t)}{\sqrt{\rho_c}} + \left(\int vf_h^n\,\mathrm{d}v\right)_j\sin(\sqrt{\rho_c}\Delta t) \\ E^{(1)}_j & = -{u_c}^n_j\sqrt{\rho_c}\sin(\sqrt{\rho_c}\Delta t) + E^n_j\cos(\sqrt{\rho_c}\Delta t) + \left(\int vf_h^n\,\mathrm{d}v\right)_j\cos(\sqrt{\rho_c}\Delta t) \\ \end{aligned}\qquad,\forall j, x_j = j\Delta x$
Boucle pour tout $v_k$ : effectuer une FFT pour obtenir $\left(\widehat{E^n\partial_vf_h^n}\right)_{\iota,k}$ à $k$ fixé, et calcul des $(\hat{f}_h^{(1)})_{\iota,k}$ : ${\hat{f}_h}^{(1)}_{\iota,k} = e^{-i\kappa_\iota v_k\Delta t}{\hat{f}_h}^{n}_{\iota,k} - e^{-i\kappa_\iota v_k\Delta t}\left(\widehat{E^n\partial_vf_h^n}\right)_{\iota,k}\quad,\forall \iota, \kappa_\iota = \frac{2\pi\iota}{L}$

L’étape 1 est le calcul de $N(U^n)$ , puis les étapes 2 et 3 sont indépendantes (donc parallélisables). Les autres étages du schéma de Lawson (quelque soit le schéma RK sous-jacent) suivent la même logique. La boucle sur les $\kappa_\iota$ étant à l’intérieur de la boucle en $v_k$ il n’est pas possible d’effectuer le calcul des $({u_c^{(1)}}_i,E^{(1)}_i)$ en même temps que ${\hat{f}_h^{(1)}}_{\iota,k}$ .

Dans l’algorithme suivant, qui détaille les 3 étages de la méthode, toutes les variables avec des points ( $\dot{X}$ ) indiquent des variables temporaires ne servant qu’à un étage, ou dans une boucle.

étage 1 :
- $(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^n}_{\iota,k})$
- $(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v$
- $(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^n\partial_v\dot{f})_{j,k})$
- pour tout $x_j$ :
  - ${u_c^{(1)}}_j \gets {u_c^n}_j\cos(\sqrt{\rho_c}\Delta t) + E^n_j\frac{\sin(\sqrt{\rho_c}\Delta t)}{\sqrt{\rho_c}} - \Delta t\dot{J}_j\sin(\Delta t)$
  - $E^{(1)}_j \gets -{u_c^n}_j\sqrt{\rho_c}\sin(\sqrt{\rho_c}\Delta t) + E^n_j\cos(\sqrt{\rho_c}\Delta t) - \Delta t\dot{J}_j\cos(\Delta t)$
- pour tout $v_j$ :
  - $(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})$
  - pour tout $\kappa_\iota$ :
    - ${\hat{f}_h^{(1)}}_{\iota,k} \gets {\hat{f}_h^{n}}_{\iota,k}e^{-i\kappa_\iota v_k\Delta t} - \Delta t\dot{\delta}_\iota e^{-i\kappa_\iota v_k \Delta t}$
étage 2:
- $(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^{(1)}}_{\iota,k})$
- $(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v$
- $(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^{(1)}\partial_v\dot{f})_{j,k})$
- pour tout $x_j$ :
  - ${u_c^{(1)}}_j \gets \frac{3}{4}{u_c^n}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{3}{4}E^n_j\frac{\sin(\frac{\Delta t}{2})}{\sqrt{\rho_c}} + \frac{1}{4}{u_c^{(1)}}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) - \frac{1}{4}E^{(1)}_j\frac{\sin(\sqrt{\rho_c}\frac{\Delta t}{2})}{\sqrt{\rho_c}} + \frac{1}{4}\Delta t\dot{J}_j\sin(\frac{\sqrt{\rho_c}\Delta t}{2})$
  - $E^{(1)}_j \gets -\frac{3}{4}{u_c^n}_j\sqrt{\rho_c}\sin(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{3}{4}E^n_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{1}{4}{u_c^{(1)}}_j\sqrt{\rho_c}\sin(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{1}{4}E^{(1)}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) - \frac{1}{4}\Delta t\dot{J}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2})$
- pour tout $v_j$ :
  - $(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})$
  - pour tout $\kappa_\iota$ :
    - ${\hat{f}_h^{(1)}}_{\iota,k} \gets \frac{3}{4}{\hat{f}_h^{n}}_{\iota,k}e^{-i\kappa_\iota v_k\frac{\Delta t}{2}} + \frac{1}{4}{\hat{f}_h^{(1)}}_{\iota,k}e^{i\kappa_\iota v_k \frac{\Delta t}{2}} - \frac{1}{4}\Delta t\dot{\delta}_\iota e^{i\kappa_\iota v_k \frac{\Delta t}{2}}$
étage 3:
- $(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^{(2)}}_{\iota,k})$
- $(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v$
- $(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^{(2)}\partial_v\dot{f})_{j,k})$
- pour tout $x_j$ :
  - ${u_c^{n+1}}_j \gets \frac{1}{3}{u_c^n}_j\cos(\sqrt{\rho_c}\Delta t) + \frac{1}{3}E^n_j\frac{\sin(\sqrt{\rho_c}\Delta t)}{\sqrt{\rho_c}} + \frac{2}{3}{u_c^{(2)}}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{2}{3}E^{(2)}_j\frac{\sin(\sqrt{\rho_c}\frac{\Delta t}{2})}{\sqrt{\rho_c}} - \frac{2}{3}\Delta t\dot{J}_j\sin(\sqrt{\rho_c}\frac{\Delta t}{2})$
  - $E^{n+1}_j \gets -\frac{1}{3}{u_c^n}_j\sqrt{\rho_c}\sin(\Delta t) + \frac{1}{3}E^n_j\cos(\sqrt{\rho_c}\Delta t) - \frac{2}{3}{u_c^{(2)}}_j\sqrt{\rho_c}\sin(\sqrt{\rho_c}\frac{\Delta t}{2}) + \frac{2}{3}E^{(2)}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2}) - \frac{2}{3}\Delta t\dot{J}_j\cos(\sqrt{\rho_c}\frac{\Delta t}{2})$
- pour tout $v_j$ :
  - $(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})$
  - pour tout $\kappa_\iota$ :
    - ${\hat{f}_h^{n+1}}_{\iota,k} \gets \frac{1}{3}{\hat{f}_h^{n}}_{\iota,k}e^{-i\kappa_\iota v_k \Delta t} + \frac{2}{3}{\hat{f}_h^{(2)}}_{\iota,k}e^{-i\kappa_\iota v_k \frac{\Delta t}{2}} - \frac{2}{3}\Delta t\dot{\delta}_\iota e^{-i\kappa_\iota v_k \frac{\Delta t}{2}}$

2019-12-05

Semaine à Nantes

Correctif

Cette semaine j’ai corrigé le modèle de Vlasov-Poisson Hybride Linéarisé que je considérais comme :

$\partial_t U = AU + N(U)$ avec : $U = \begin{pmatrix}u_c\\E\\\hat{f}_h\end{pmatrix} \qquad A = \begin{pmatrix}0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & e^{-ikv}\end{pmatrix} \qquad N(U) = \begin{pmatrix} 0 \\ \int vf_h\,\mathrm{d}v \\ -\widehat{E\partial_vf_h} \end{pmatrix}$ avec comme condition initiale : $f = (1-\alpha)f_c + \alpha f_h\times(1 + \epsilon\cos(kx))$ . Pour $f_h$ nous prenons : $f_h = \mathcal{M}_{[\alpha/2,u,1]} + \mathcal{M}_{[\alpha/2,-u,1]}$ il est également possible de prendre : $f_h = \mathcal{M}_{[1,0,1]}\frac{v^{10}}{\mathcal{S}}$ où $\mathcal{S}$ est une constante de normalisation. Dans la pratique il est nécessaire dans le modèle de considérer la densité des particules froides $\rho_c$ même si lors de la linéarisation celle-ci est considérer comme constante, et il y a une erreur de signe. Le modèle devient alors : $A = \begin{pmatrix}0 & 1 & 0 \\ -\rho_c & 0 & 0 \\ 0 & 0 & e^{-ikv}\end{pmatrix} \qquad N(U) = \begin{pmatrix} 0 \\ -\int vf_h\,\mathrm{d}v \\ -\widehat{E\partial_vf_h} \end{pmatrix}$ avec $\rho_c = 1-\alpha$ en utilisant les notations des conditions initiales proposées.

Tout le document a été corrigé conformément à cette erreur.

Simulation

Nous avons opté pour la comparaison entre 2 méthodes de simulation pour étudier ce modèle :

Méthode de Lawson, avec résolution du terme $E\partial_vf_h$ par WENO5
Méthode de splitting à 3 étages, avec résolution de l’équation $\partial_tf_h+E\partial_vf_h=0$ par une méthode semi-lagrangienne d’ordre 5 dans un premier temps, avant d’utiliser une méthode WENO-SL.

Il existe un code opérationnelle pour la première méthode ainsi que la deuxième (avec un polynôme de Lagrange de degré 5).

J’ai enfin effectué un nouveau commit sur GitHub pour avoir une sauvegarde de mon code. Ces codes s’y trouvent.

J’ai donc maintenant 3 codes de simulation de Vlasov-Poisson :

un code full cinétique avec Lawson+WENO
un code hybride avec Lawson+WENO
un code hybride avec splitting à 3 étages

2019-12-12

Étudions le WENO-SL présenté par Qiu et Shu dans [Qiu:2011]. Il s’agit d’une méthode volumes finis semi-Lagrangienne. Pour cela on s’intéresse à l’équation : $u_t + (au)_x = 0$ On commence par une intégration sur $[t^n,t^{n+1}]$ : $u(t^{n+1},x) = u(t^n,x) - \left(\int_{[t^n,t^{n+1}]}au(\tau,x)\,\mathrm{d}\tau\right)_x$ soit : $u_i^{n+1} = u_i^n - \mathcal{F}_{x|x=x_1}$ où $\mathcal{F}$ est défini par : $\mathcal{F}(x) = \int_{[t^n,t^{n+1}]} au(\tau,x)\,\mathrm{d}\tau = \frac{1}{\Delta x}\int_{[x-\frac{\Delta x}{2},x+\frac{\Delta x}{2}]}\mathcal{H}(\xi)\,\mathrm{d}\xi$ Cela nous permet dans un premier temps de réécrire la dérivée par rapport à $x$ de $\mathcal{F}$ (ce que l’on doit calculer pour écrire notre schéma en temps) comme une différence de flux (méthode des volumes finis) : $\mathcal{F}_x = \frac{1}{\Delta x}\left(\mathcal{H}\left(x+\frac{\Delta x}{2}\right) - \mathcal{H}\left(x-\frac{\Delta x}{2}\right) \right)$ D’où l’écriture de : $u_i^{n+1} = u_i^n - \frac{1}{\Delta x}\left(\mathcal{H}(x_{i+\frac{1}{2}}) - \mathcal{H}(x_{i-\frac{1}{2}})\right)$

D’un autre côté, on intègre sur le domaine triangulaire $\Omega_i$ défini par les trois points $(x_i,t^{n+1})$ , $(x_i,t^n)$ et $(x_i^\star,t^n)$ (voir figure pour plus d’informations), l’équation de départ.

Domaine d’intégration \Omega pour la construction de la méthode semi-lagrangienne — Domaine d’intégration $\Omega$ pour la construction de la méthode semi-lagrangienne

Le point $x_i^\star$ est défini comme le pied de la caractérisitique issue de $x_i$ . Par intégration de l’équation sur $\Omega_i$ :

$\int_{\Omega_i} u_t + (au)_x = 0$

par le théorème de Green, on obtient directement (dans l’article ils utilisent le théorème de la divergence, mais je ne comprends pas alors comment est gérer le bord $(x_i^\star,t^n)$ , $(x_i,t^{n+1})$ ) : $\int_{\Omega_i} u_t + (au)_x = -\int_{x_i^\star}^{x_i} u(t^n,\xi)\,\mathrm{d}\xi + \int_{t^n}^{t^{n+1}} au(\tau,x_i)\,\mathrm{d}\tau$ ce qui nous donne : $\int_{x_i^\star}^{x_i} u(t^n,\xi)\,\mathrm{d}\xi = \int_{t^n}^{t^{n+1}} au(\tau,x_i)\,\mathrm{d}\tau = \mathcal{F}(x_i)$ avec (méthode des volumes finis oblige) $\mathcal{F}(x_i) = \bar{\mathcal{H}}_i$ (valeur moyenne de $\mathcal{H}$ sur la cellule $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ ). Donc on utilise le pied exacte de la caractéristique $x_i^\star$ pour faire une estimation du flux, on construit donc une méthode volumes finis semi-lagrangienne.

Il est nécessaire de reconstruire les $\{\bar{\mathcal{H}}_i\}_{i=1}^N$ à partir des $\{u_i^n\}_{i=1}^N$ , plus précisemment les points nécessaires sont : $\bar{\mathcal{H}}_i = \mathcal{R}[x_i^\star,x_i](u_{i-p_1}^n,\dots,u_{i+q_1}^n)$ où $\mathcal{R}[a,b]$ est la reconstruction de $\int_{[a,b]}u(t,\xi)\,\mathrm{d}\xi$ , en utilisant les points indiqués dans le stencil entre parenthèses. Le problème que l’on peut remarqué avec cette méthode volumes-finis/semi-lagrangienne est qu’il est nécessaire d’avoir tous les points entre $x_i^\star$ et $x_i$ , aussi loin soit $x_i^\star$ . Dans la pratique, lors de l’écriture du schéma, seul le stencil autour de $x_i^\star$ est nécessaire, mais pour conserver la formulation volumes-finis, il est nécessaire de considérer la reconstruction jusqu’à $x_i$ .

La suite de l’articule est tout un jeu de construction de l’opérateur de reconstruction $\mathcal{R}[a,b]$ en utilisant une méthode WENO, dans un premier temps avec $x_i^\star$ dans la cellule $[x_{i-1},x_i]$ (volumes finis classiques), puis en dehors de cette cellule (avec la prise en compte de $\xi=\frac{x_{i^\star}-x_i^\star}{\Delta x}$ , donc un vrai WENO-SL).

Ce que j’ai pour le moment implémenté le 2019-11-06 n’est que la version où $x_i^\star$ est dans la cellule $[x_{i-1},x_i]$ .

2019-12-17

Ça y est, le projet ponio (Python Objects for Numerical IntegratOr) est partagé sur pypi, il suffit donc maintenant d’effectuer :

  pip install ponio

pour bénéficier des meilleurs outils d’analyse de méthodes Runge-Kutta !

J’ajouterai prochainement tous mes outils d’analyse de WENO (5 et 3) linéarisé, et potentiellement d’autres intégrateur en espace. Il s’agit d’une bibliothèque d’outils pour l’étude de la stabilité de l’équation d’advection : $u_t + u_x = 0$ Je pense faire une mise-à-jour majeur (ajout des outils d’analyse de von Neumann) pendant les prochaines vacances.

Préparation de simus de comparaisons

Maintenant que j’ai 3 codes de simulations de Vlasov :

Full-kinetic avec Lawson
Vlasov hybrid linearized avec Lawson
Vlasov hybrid linearized avec splitting

il faut comparer les résultats. Nicolas propose de faire des simulations avec un maillage très fin (surtout en vitesse) pour comparer le modèle full-kinetic avec les modélisations hybrides, et ne voir que l’erreur d’approximation du $\delta_0(v)$ par une gaussienne. Pour cela on simule à $t$ fixé ( $\sim 5$ ou $10$ ). La gaussienne doit être au moins représentée par une vingtaine de points pour être correctement représentée, donc : $\Delta v \sim \frac{\sqrt{T_c}}{20}$ avec $T_c = 10^{-4},\,10^{-3},\,10^{-2},\,10^{-1}$ . Donc pour s’assurer d’avoir toujours le même $\Delta v$ (et donc le même $\Delta t$ ), on prend $\Delta v \sim \frac{\sqrt{10^{-4}}}{20} = 5.0\cdot 10^{-4}$ . On prend le $\Delta t$ qui CFLise ce $\Delta v$ : $\Delta t = \sigma\frac{\Delta v}{||E^n||_\infty}$ Pour le moment le code VHL avec Lawson utilise un RK( $3$ , $3$ ), donc $\sigma = 1.433$ et les premières simus donnent $||E^n||_\infty \leq 0.12$ , on obtient alors pour $\Delta v$ , $N_v$ et $\Delta t$ : $\Delta v \approx 5.0\cdot 10^{-4} \qquad N_v = \frac{16}{\Delta v}\approx 32\,000 \qquad \Delta t \approx 5.9\cdot 10^{-3}$ Ce qui nous donne environ 1700 itérations pour atteindre le temps final 10, mais avec un $N_v$ très grand (le maximum que j’avais pris jusque là était 2048).

Comparaison des pas de temps maximaux possibles avec différentes températures et taille de la grille en vitesse.
$T_c$	$\Delta v \sim \frac{\sqrt{T_c}}{20}$	$N_v \sim \frac{16}{\Delta v}$	$\Delta t \sim 1.433\frac{\Delta v}{0.12}$
$10^{-4}$	$5.00\cdot 10^{-4}$	$32\,000$	$5.97\cdot 10^{-3}$
$10^{-3}$	$1.58\cdot 10^{-3}$	$10\,119$	$1.89\cdot 10^{-2}$
$10^{-2}$	$5.00\cdot 10^{-3}$	$3\,200$	$5.97\cdot 10^{-2}$
$10^{-1}$	$1.58\cdot 10^{-2}$	$1\,012$	$1.89\cdot 10^{-1}$

Il est intéressant de noter que les modélisations hybrides supposent $T_c\to 0$ mais n’imposent pas de raffinement de maillage. Il peut être intéressant de comparer les résultats avec $T_c = 10^{-4}$ (et la grille en vitesse que cela impose) et une grille bien plus grossière avec une modélisation hybride).

On souhaite, en plus de ces 4 simulations en full-kinetic, avoir les résultats sur au moins 1 des 2 codes de modélisation hybride, pour comparer les quantités :

le champ électrique $\left(\sum_i |E^{\varepsilon}_i(t) - E_i(t)|^2\right)^{^1/_2}$ , il est donc nécessaire d’avoir des simulations avec les mêmes pas de temps.
la densité $\rho^{\varepsilon}(t)$ contre $\rho_c(t)+\rho_h(t)$ , sachant que $\rho_h$ est directement donnée par $f_h$ , mais $\rho_c$ est une donnée perdue du problème qu’il faut récupérer par l’équation de Poisson (avec Fourier).
le courant $J=\int vf\,\mathrm{d}v$ contre $\rho_cu_c + \int vf_h\,\mathrm{d}v$
l’énergie totale $H=\iint v^2f\,\mathrm{d}v\mathrm{d}x + \int E^2\,\mathrm{d}x$ , contre $H=\int \rho_cu_c^2\,\mathrm{d}x + \iint v^2f_h\,\mathrm{d}x + \int E^2\,\mathrm{d}x$

Bref il faut stocker avec les mêmes pas de temps :

l’énergie électrique (une seule valeure par itéation)
l’énergie totale (une seule valeure par itération)
le champ électrique en fin de simu
la densité en fin de simu
le courant en fin de simu

2019-12-19

À propos des relations de dispersion

Considérons le modèle de Vlasov-Poisson : $\begin{cases} \partial_t f + v\partial_x f + E\partial_v f = 0 \\ \partial_x E = \int_{\mathbb{R}} f\,\mathrm{d}v - 1 \end{cases}$ On souhaite linéariser le systène autour d’un état d’équilibre $f(t,x,v) = f_\text{eq}(v) + \varepsilon g(t,x,v)$ : $\begin{cases} \varepsilon\partial_t g + v\varepsilon\partial_x g + E\mathrm{d}_v f_\text{eq} + \varepsilon E\partial_v g = 0 \\ \partial_x E = \int_{\mathbb{R}} f_\text{eq}\,\mathrm{d}v + \varepsilon\int_\mathbb{R} g\,\mathrm{d}v - 1 \end{cases}$ On peut normalisé en considérant que la densité de particules à l’équilibre est de masse $1$ , par conséquent, l’équation de Poisson implique que $E$ est de taille $\varepsilon$ , on notera alors $E(t,x) = \varepsilon e(t,x)$ : $\begin{cases} \varepsilon\partial_t g + v\varepsilon\partial_x g + e\varepsilon\mathrm{d}_v f_\text{eq} + \varepsilon^2 e\partial_v g = 0 \\ \varepsilon\partial_x e = \varepsilon\int_\mathbb{R} g\,\mathrm{d}v \end{cases}$ en ne conservant que les termes d’ordre $\varepsilon$ on obtient finalement : $\begin{cases} \partial_t g + v\partial_x g + e\mathrm{d}_v f_\text{eq} = 0 \\ \partial_x e = \int_\mathbb{R} g\,\mathrm{d}v \end{cases}$

On note avec un $\mathrm{d}$ droit la dérivée en $v$ de $f_\text{eq}$ car $f_\text{eq}$ ne dépend que de $v$ , c’est une notation usuelle en physique.

On peut maintenant appliquer les transformées de Fourier et de Laplace, qui peuvent se résumer comme suit : $\partial_x\cdot \to ik\cdot \qquad \partial_t\cdot \to \cdot_{|t=0} - i\omega\cdot$

On continuera à utiliser le même nom de variable après les transformées successives de Fourier et de Laplace.

On peut simplifier le calcul en prennant une condition initiale nulle dans la transformée de Laplace. Dans tous les cas cela mènera à la même fonction $D(\omega,k)$ .

$\begin{cases} ik\left(v-\frac{\omega}{k}\right) g = e\mathrm{d}_v f_\text{eq} - g_{|t=0} \\ ik e = \int g\,\mathrm{d}v \end{cases}$ La première équation nous permet d’expliciter $g$ : $g = \frac{ie\mathrm{d}_v f_\text{eq}}{k\left(v-\frac{\omega}{k}\right)} + \frac{i g_{|t=0}}{k\left(v-\frac{\omega}{k}\right)}$ La seconde équation donne alors : $e = -\frac{i}{k}\int \frac{ie\mathrm{d}_v f_\text{eq}}{k\left(v-\frac{\omega}{k}\right)}\,\mathrm{d}v - \frac{i}{k}\int \frac{i g_{|t=0}}{k\left(v-\frac{\omega}{k}\right)}\,\mathrm{d}v$ soit : $e = \frac{\frac{1}{k^2}\int \frac{g_{|t=0}}{v-\frac{\omega}{k}}\,\mathrm{d}v }{ 1 - \frac{1}{k^2}\int\frac{\mathrm{d}_v f_\text{eq}}{v-\frac{\omega}{k}}\,\mathrm{d}v }$ On notera alors $e = \frac{N(\omega,k)}{D(\omega,k)}$ avec : $N(\omega,k) = \frac{1}{k^2}\int \frac{g_{|t=0}}{v-\frac{\omega}{k}}\,\mathrm{d}v \qquad D(\omega,k) = 1 - \frac{1}{k^2}\int\frac{\mathrm{d}_v f_\text{eq}}{v-\frac{\omega}{k}}\,\mathrm{d}v$ Bien entendu on s’intéresse aux pôles de cette fraction, c’est-à-dire aux zéros de $D(\omega,k)$ .

Calcul d’intégrales

Il est intéressant de calculer pour des cas classiques la valeur de $D(\omega,k)$ , considérons alors 2 cas classiques :

$f_\text{eq}(v) = \mathcal{M}_{\rho,u,T}(v) = \frac{\rho}{\sqrt{2\pi T}}e^{-\frac{|v-u|^2}{2T}}$
$f_\text{eq}(v) = \delta_u(v)$

Dans ces 2 cas, nous calculerons $\int_\mathbb{R} \frac{\mathrm{d}_v f_\text{eq}}{v-\beta}\,\mathrm{d}v$ .

Cas d’une distribution maxwellienne

$\begin{aligned} \int_\mathbb{R}\frac{\mathrm{d}_v\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \int_\mathbb{R} \frac{ -\frac{(v-u)}{T}\mathcal{M}_{\rho,u,T}(v) }{v-\beta}\,\mathrm{d}v \\ &= -\frac{1}{T}\int_\mathbb{R} \frac{(v-\beta+\beta-u)\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v \\ &= -\frac{\rho}{T} - \frac{\beta-u}{T}\int_\mathbb{R}\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v \end{aligned}$ On s’intéresse maintenant plus particulièrement à l’intégrale : $\int\frac{\mathcal{M}}{v-\beta}\,\mathrm{d}v$ : $\begin{aligned} \int_\mathbb{R} \frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \frac{\rho}{\sqrt{2\pi T}}\int_\mathbb{R}\frac{ e^{-\frac{|v-u|^2}{2T}} }{v-\beta}\,\mathrm{d}v \\ &= \frac{\rho}{\sqrt{2\pi T}}\int_\mathbb{R}\frac{ e^{-\left(\frac{v-u}{\sqrt{2T}}\right)^2} }{v-u + u-\beta}\,\mathrm{d}v \end{aligned}$ on effectue alors le changement de variable $z=\frac{v-u}{\sqrt{2T}}$ , $\mathrm{d}z=\frac{\mathrm{d}v}{\sqrt{2T}}$ : $\begin{aligned} \int_\mathbb{R} \frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v &= \frac{\rho}{\sqrt{2\pi T}}\int_\mathbb{R} \frac{e^{-z^2}}{\sqrt{2T}z + u-\beta}\sqrt{2T}\mathrm{d}z \\ &= \frac{\rho}{\sqrt{2\pi T}}\int_\mathbb{R} \frac{e^{-z^2}}{z+\frac{u-\beta}{\sqrt{2T}}}\,\mathrm{d}z \end{aligned}$ On introduit alors la fonction de Fried-Conte : $\mathcal{Z}: \zeta \mapsto \frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty} \frac{e^{-z^2}}{z-\zeta}\,\mathrm{d}z$ D’où : $\int_\mathbb{R}\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v = \frac{\rho}{\sqrt{2T}}\mathcal{Z}\left(\frac{\beta-u}{\sqrt{2T}}\right)$ ce qui nous permet de donner une expression de notre première intégrale : $\boxed{ \int_\mathbb{R} \frac{\mathrm{d}_v\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v = -\frac{\rho}{T}\left( 1 + \frac{\beta-u}{\sqrt{2T}}\mathcal{Z}\left(\frac{\beta-u}{\sqrt{2T}}\right) \right) }$

Cas d’une distribution de Dirac

$\begin{aligned} \int_\mathbb{R}\frac{\mathrm{d}_v\delta_u(v)}{v-\beta}\,\mathrm{d}v &= -\int_\mathbb{R} \delta_u(v)\mathrm{d}_v\left(\frac{1}{v-\beta}\right)\,\mathrm{d}v \\ &= \int_\mathbb{R} \delta_u(v)\frac{1}{(v-\beta)^2}\,\mathrm{d}v \end{aligned}$ soit : $\boxed{ \int_\mathbb{R} \frac{\mathrm{d}_v\delta_u(v)}{v-\beta}\,\mathrm{d}v = \frac{1}{(u-\beta)^2} }$

Tests numériques

On souhaite calculer quelques relations de dispersion pour différentes conditions initiales. En particulier nous souhaitons étudier le cas triple beams qui correspond à la condition initiale : $f(t=0,x,v) = \mathcal{M}_{(1-\alpha),0,T_c}(v) + \left(\mathcal{M}_{^\alpha\!/\!_2,u,1}(v) + \mathcal{M}_{^\alpha\!/\!_2,-u,1}(v)\right)(1+\epsilon\cos(k x))$ avec $\alpha$ la densité de particules chaudes, $u$ la vitesse caractéristique des particules chaudes, et $T_c$ un paramètre pris arbitrairement bas correpondant à la température des particules froides ( $T_c \ll 1$ ). Il est possible de dérivée de l’équation de Vlasov-Poisson un modèle hybride linéarisé qui permet de réaliser des simulations avec le cas critique $T_c\to 0$ qui nécessite, dans le cadre cinétique, une grille en vitesse infiniment fine. Dans ce cas la condition initiale se réécrit : $f(t=0,x,v) = (1-\alpha)\delta_0(v) + \left(\mathcal{M}_{^\alpha\!/\!_2,u,1}(v) + \mathcal{M}_{^\alpha\!/\!_2,-u,1}(v)\right)(1+\epsilon\cos(k x))$

Pour obtenir $f_\text{eq}$ depuis ces conditions initiales dépendant de $x$ et $v$ il suffit de supprimer la perturbation en $x$ (c’est-à-dire prendre $\epsilon=0$ ). On remarque alors que $f_\text{eq}$ n’est qu’une somme de distribution maxwellienne et de Dirac, il suffit alors, par linéarité de l’intégrale, de reprendre les calculs d’intégrales précédents.

On obtient alors pour le cas cinétique la relation de dispersion : $\boxed{ \begin{aligned} D^K(\omega,k) = 1 - \frac{1}{k^2} \left[\vphantom{\frac{}{\sqrt{}}}\right.&-\frac{1-\alpha}{T_c}\left( 1 + \frac{\omega}{k\sqrt{2T_c}}\mathcal{Z}\left(\frac{\omega}{k\sqrt{2T_c}}\right) \right) \\ &-\frac{\alpha}{2}\left( 1 + \frac{\frac{\omega}{k}-u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}-u}{\sqrt{2}}\right) \right) \\ &-\frac{\alpha}{2}\left( 1 + \frac{\frac{\omega}{k}+u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}+u}{\sqrt{2}}\right) \right)\left. \vphantom{\frac{}{\sqrt{}}}\right] \end{aligned} }$ Et dans le cas hybride linéarisé : $\boxed{ \begin{aligned} D^{H}(\omega,k) = 1 - \frac{1}{k^2}\left[\vphantom{\frac{}{\sqrt{}}}\right.&(1-\alpha)\left(\frac{k}{\omega}\right)^2 \\ &-\frac{\alpha}{2}\left( 1 + \frac{\frac{\omega}{k}-u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}-u}{\sqrt{2}}\right) \right) \\ &-\frac{\alpha}{2}\left( 1 + \frac{\frac{\omega}{k}+u}{\sqrt{2}}\mathcal{Z}\left(\frac{\frac{\omega}{k}+u}{\sqrt{2}}\right) \right)\left. \vphantom{\frac{}{\sqrt{}}}\right] \end{aligned} }$

2020-01-10

Calcul de l’énergie au cours du temps : $\begin{aligned} \mathcal{E}(t) & = \iint v^2f\,\mathrm{d}v\mathrm{d}x + \int E^2\,\mathrm{d}x \\ & = \mathcal{E}(0) \end{aligned}$

Dans le cas cinétique, $f(t=0,x,v) = \mathcal{M}_{1-\alpha,0,T_c}(v) + \left( \mathcal{M}_{\alpha/2,u,1}(v) + \mathcal{M}_{\alpha/2,-u,1}(v) \right)$ , ce qui donne : $\mathcal{E}_K(t) = \left[ (1-\alpha)T_c + \alpha u^2 + \alpha \right] L$

La perturbation avec le cosinus passe à la trape car par intégration sur tout le domaine en $x$ cela donnera $0$ .

Dans le cas hybride linéarisé, $f(t=0,x,v) = (1-\alpha)\delta_0 (v) + \left( \mathcal{M}_{\alpha/2,u,1}(v) + \mathcal{M}_{\alpha/2,-u,1}(v) \right)$ , ce qui donne : $\mathcal{E}_{HL} = \left[ \alpha u^2 + \alpha \right]L$

D’où : $\mathcal{E}_K(t) - \mathcal{E}_{HL}(t) = (1-\alpha)LT_c$ ce qui donne en version discrète : $\mathcal{E}^n_K + \mathcal{O}(\Delta t^p) - \mathcal{E}^n_{HL} + \mathcal{O}(\Delta t^q) = (1-\alpha)LT_c$ avec $p$ et $q$ respectivement les ordres en temps de codes de résolution du modèle cinétique et du modèle hybride linéarisé. Donc pour $p$ et $q$ suffisamment grand on ne trace que $(1-\alpha)LT_c$ . Il est intéressant de constater que l’erreur sur l’énergie est linéaire par rapport au produit $\rho_c T_c$ , avec $\rho_c$ la densité des particules froides dans la condition initiale $(1-\alpha)$ .

Obtention du modèle de Vlasov hybride linéarisé

Parce qu’on a refait le calcul avec Nicolas, car j’ai eu un problème pendant les vacances pour réobtenir ce modèle, et pour être certain de l’avoir quelque part.

On part du modèle de Vlasov-Ampère, avec 2 espères, les particules froides $f_c$ et les partiules chaudes $f_h$ : $\begin{cases} \partial_t f_c + v\partial_x f_c + E \partial_v f_c = 0 \\ \partial_t f_h + v\partial_x f_h + E \partial_v f_h = 0 \\ \partial_t E = - \int v\left(f_h+f_c\right)\,\mathrm{d}v \end{cases}$ Le travail s’effectue essentiellement sur $f_c$ pour l’écrire comme un modèle fluide à partir de grandeurs intégrales, on calcule donc les moments de la première équation : $\int_\mathbb{R} \partial_t \begin{pmatrix}1\\ v\end{pmatrix} f_c\,\mathrm{d}v + \int_\mathbb{R} v\partial_x \begin{pmatrix}1\\ v\end{pmatrix} f_c\,\mathrm{d}v + \int_\mathbb{R} E\partial_v \begin{pmatrix}1\\ v\end{pmatrix} f_c\,\mathrm{d}v = 0$ de plus, on suppose que $f_c(t,x,v) = \rho_c(t,x)\delta_{u_c(t,x)}(v)$ , ce qui nous permet d’obtenir : $\begin{cases} \partial_t \rho_c + \partial_x (\rho_cu_c) = 0 \\ \partial_t (\rho_cu_c) + \partial_x (\rho_cu_c^2)- \rho_cE = 0 \end{cases}$ On obtient ainsi le modèle de Vlasov hybride non-linéaire : $\begin{cases} \partial_t \rho_c + \partial_x(\rho_cu_c) = 0 \\ \partial_t (\rho_cu_c) + \partial_x (\rho_cu_c^2) - \rho_c E = 0 \\ \partial_t E = - \int vf_h\,\mathrm{d}v - \rho_cu_c \\ \partial_x E = \int f_h\,\mathrm{d}v + \rho_c - 1 \\ \partial_t f_h + v\partial_x f_h + E \partial_v f_h = 0 \\ \end{cases}$ Il est a noter que l’équation de Poisson (4e équation) est toujours vérifiée si on calcule le champ électrique $E$ à l’aide de l’équation d’Ampère (3e équation) et que l’équation de Poisson est vérifiée par la condition initiale, on ne l’écrira plus par la suite et on supposera seulement que notre condition initiale vérifie Poisson.

On souhaite maintenant linéarisé ce problème autour de : $\begin{aligned} \rho_c(t,x) & = \rho_c^{(0)}(x) + & \varepsilon\rho_c^{(1)}(t,x) \\ u_c(t,x) & = & \varepsilon u_c^{(1)}(t,x) \\ E(t,x) & = & \varepsilon E^{(1)}(t,x) \end{aligned}$ On obtient alors (en négligeant tous les termes quadratiques ou plus) : $\begin{cases} \varepsilon\partial_t \rho_c^{(1)} + \varepsilon\partial_x(\rho_c^{(0)}u^{(1)}) = 0 \\ \varepsilon\rho_c^{(0)}\partial_t u_c^{(1)} - \varepsilon \rho_c^{(0)}E^{(1)} = 0 \\ \varepsilon\partial_t E^{(1)} = - \int vf_h\,\mathrm{d}v - \varepsilon\rho_c^{(0)}u_c^{(1)} \\ \partial_t f_h + v\partial_x f_h + E^{(1)}\partial_v f_h = 0 \end{cases}$ La première équation est la seule faisant intervenir $\rho_c^{(1)}$ , il ne sert pas pour le calcul des autres variables et il peut être réobtenu par l’équation de Poisson (qui est toujours vérifiée si la condition initiale la vérifie), on retire donc de notre modèle cette première équation ainsi que cette variable. On remarque que les autres équations ne font intervenir que des termes d’ordre $\varepsilon$ ( $f_h$ étant supposé petit par rapport à $f_c$ ), on peut donc simplifier par $\varepsilon$ pour obtenir le modèle de Vlasov hybride linéarisé : $\begin{cases} \partial_t u_c^{(1)} = E^{(1)} \\ \partial_t E^{(1)} = -\int vf_h\,\mathrm{d}v - \rho_c^{(0)}u_c^{(1)} \\ \partial_t f_h + v\partial_x f_h + E^{(1)}\partial_v f_h = 0 \end{cases}$

To Do list

S’inscrire à Modern Scientific C++ du 20 au 23 avril à la maison de la simulation.
Tracer $\mathcal{E}_K(T_c)$ (pour voir une jolie droite). Pour cela choisir un temps arbitraire et faire une interpolation de Lagrange d’ordre 5 pour avoir pour chaque $T_c$ le même temps.
Faire la même chose avec les autres grandeurs que l’on suit au cours du temps $E_\text{max}$ et l’énergie électrique.
Faire des tests à produit $\rho_cT_c$ constant (donc différents $\alpha$ ).
Vérifier avec les solutions des relations de dispersion et tracer $\omega_K(T_c) - \omega_{HL}$ en fonction de $T_c$
Implémenter DP5 pour le VHLL (Vlasov Hybride Linéarisé avec Lawson), et commencer à regarder le pas de temps adaptatif.
Comparaison splitting - Lawson VHL sur l’énergie totale et tester avec des grands $\Delta t$ (test de stabilité). Il est à noter que dans ce cas, $\Delta v$ n’est pas imposé par $T_c$ et peut donc être pris grand.

Résultats

Les courbes suivantes sont issus de l’interpolation de données au temps $t_1 = 9.5$ .

Évolution de l’énergie totale H en fonction de la température T_c — Évolution de l’énergie totale $H$ en fonction de la température $T_c$

La pente de 10 concorde avec ce qui est attendu, $\alpha = 0.2$ et $L = 4\pi$ , donc $(1-\alpha)L \approx 10.05$ .

Évolution de E_\text{max} en fonction de la température T_c — Évolution de $E_\text{max}$ en fonction de la température $T_c$

Évolution de l’énergie électrique en fonction de la température T_c — Évolution de l’énergie électrique en fonction de la température $T_c$

On peut aussi comparer les zéros des relations de dispersion entre le modèle cinétique (pour différentes valeurs de $T_c$ ) et celle du modèle hybride, on trace ainsi $\Delta \omega_c = \omega_K(T_c) - \omega_{HL}$ .

Attention : pour bien observer une droite la température $T_c$ a été prise sur un domaine plus grand.

Différence des zéros de la relation de dispersion en fonction de la température T_c — Différence des zéros de la relation de dispersion en fonction de la température $T_c$

On représente aussi l’erreur relative faite sur l’énergie totale. Le code cinétique donne, quelque soit la température, une erreur relative sur l’énergie totale nulle, à l’inverse du code hybride (variation très faible).

Évolution temporelle de l’erreur relative

2020-09-03

Je travaille sur la résolution numérique du modèle : $\partial_t U = AU+F(U)$

avec : $U = \left[\begin{matrix}j_{c,x}\\j_{c,y}\\B_{x}\\B_{y}\\E_{x}\\E_{y}\\\hat{f}\end{matrix}\right]$ et : $A = \left[\begin{matrix}0 & -1 & 0 & 0 & \Omega_{pe}^{2} & 0 & 0\\1 & 0 & 0 & 0 & 0 & \Omega_{pe}^{2} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0\\-1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & -1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - i k v_{z}\end{matrix}\right] \ , \qquad F(U) = \left[\begin{matrix}0\\0\\i E_{y} k\\- i E_{x} k\\- i B_{y} k + \int_\mathbb{R}\ v_x{\left(\hat{f} \right)}\\i B_{x} k + \int_\mathbb{R}\ v_y{\left(\hat{f} \right)}\\\partial_{\dots}{\left(\hat{f} \right)}\end{matrix}\right]$

On souhaite écrire un schéma exponentiel (Lawson(RK(3,3)) pour résoudre ce problème, celui-ci s’écrit :

$\begin{aligned} U^{(1)} &= e^{\Delta t A}U^n + \Delta te^{\Delta t A}F\left(U^n\right) \\ U^{(2)} &= \frac{3}{4}e^{\frac{\Delta t}{2}A}U^n + \frac{1}{4}e^{-\frac{\Delta t}{2}A}U^{(1)} + \frac{\Delta t}{4}e^{-\frac{\Delta t}{2}A}F\left(U^{1}\right) \\ U^{n+1} &= \frac{1}{3}e^{\Delta t A}U^n + \frac{2}{3}e^{\frac{\Delta t}{2}A}U^{(2)} + \frac{2}{3}\Delta te^{\frac{\Delta t}{2}A}F\left(U^{(2)}\right) \end{aligned}$

Il s’agit de la méthode Lawson induite par la méthode RK(3,3) de Shu-Osher (qui fait donc qu’un seul appel à la fonction non linéaire $F$ par étage, ce qui permet de réduire la complexité temporelle du schéma, au détriment de la complexité spatiale).

Le calcul de $e^{\Delta tA}$ peut se faire de manière exacte (pour la valeur de $\Omega_{pe}=2$ ) :

$e^{\Delta tA} = \left[\begin{matrix}0.378732187481834 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518167 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.378732187481834 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.621267812518167 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.970142500145332 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & - 0.970142500145332 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0\\- 0.378732187481834 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518167 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.378732187481834 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518167 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.970142500145332 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.970142500145332 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.970142500145332 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0\\0 & 0 & 1.0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1.0 & 0 & 0 & 0\\- 0.242535625036333 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.242535625036333 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.621267812518167 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481834 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.621267812518166 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.378732187481833 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0\\- 0.242535625036333 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.242535625036333 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & - 0.242535625036333 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & - 0.621267812518167 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481834 \sin{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.621267812518166 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481833 \cos{\left(\frac{\sqrt{2} \Delta t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0\\0 & 0 & 0 & 0 & 0 & 0 & e^{- i \Delta t k v_{z}}\end{matrix}\right]$

Et plus généralement je n’ai pas de problème pour calculer l’exponentielle de cette matrice pour des valeurs particulière de $\Omega_{pe}$ .

J’ai réussi à générer automatiquement quelque chose de proche du code C++ que je vais écrire à partir de ces expressions de schéma et matricielle, par exemple pour le premier étage :

jx1[i] = 0.970142500145332*Exn[i]*(std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 0.970142500145332*Eyn[i]*(-1.*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)) + std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.))) - 0.970142500145332*dt*(-1.*jcx[i] + I*Byn[i]*k)*(std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 0.970142500145332*dt*(-1.*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)) + std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)))*(I*Bxn[i]*k + jcy[i]) + 1.0*jxn[i]*(0.378732187481834*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.621267812518167*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*jyn[i]*(0.378732187481834*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.621267812518167*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)));
//---
jy1[i] = 1.0*Exn[i]*(0.970142500145332*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.970142500145332*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*Eyn[i]*(0.970142500145332*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.970142500145332*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*dt*(-1.*jcx[i] + I*Byn[i]*k)*(0.970142500145332*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.970142500145332*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*dt*(I*Bxn[i]*k + jcy[i])*(0.970142500145332*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.970142500145332*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*jxn[i]*(0.378732187481834*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.621267812518167*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*jyn[i]*(0.378732187481834*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.621267812518167*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)));
//---
Bx1[i] = 1.0*Bxn[i] + 1.0*I*Eyn[i]*dt*k;
//---
By1[i] = 1.0*Byn[i] - 1.0*I*Exn[i]*dt*k;
//---
Ex1[i] = 1.0*Exn[i]*(0.621267812518167*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.378732187481834*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*Eyn[i]*(0.621267812518166*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.378732187481833*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*dt*(-1.*jcx[i] + I*Byn[i]*k)*(0.621267812518167*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.378732187481834*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*dt*(I*Bxn[i]*k + jcy[i])*(0.621267812518166*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.378732187481833*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*jxn[i]*(0.242535625036333*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.242535625036333*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*jyn[i]*(0.242535625036333*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.242535625036333*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)));
//---
Ey1[i] = -1.0*Exn[i]*(0.621267812518167*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.378732187481834*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*Eyn[i]*(0.621267812518166*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.378732187481833*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*dt*(-1.*jcx[i] + I*Byn[i]*k)*(0.621267812518167*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.378732187481834*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) + 1.0*dt*(I*Bxn[i]*k + jcy[i])*(0.621267812518166*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.378732187481833*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*jxn[i]*(0.242535625036333*std::cos(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) - 0.242535625036333*std::cos(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.))) - 1.0*jyn[i]*(0.242535625036333*std::sin(1.5*dt*std::sqrt(-0.222222222222222*std::sqrt(17.) + 2.)) + 0.242535625036333*std::sin(1.5*dt*std::sqrt(0.222222222222222*std::sqrt(17.) + 2.)));
//---
hf1[k_x][k_y][k_z][i] = 1.0*(dt*\partial_\dots(1.0*hfn[k_x][k_y][k_z][i]) + hfn[k_x][k_y][k_z][i])*std::exp(-1.0*I*dt*k*v_z);

Je n’ai pas réussi à retirer de l’expression sympy les 1.0*..., normalement le compilateur devrait les supprimer lors de la première étage de compilation, au pire je pourrais les retirer à la main. Les autres étages sont générés avec sympy de la même manière.

Vous aurez remarqué la présence des variables jcx[i] et jcy[i] non définies ailleurs, il s’agit de l’approximation de $\int_\mathbb{R}v_{x|y}\,f\,\mathrm{d}v$ , ainsi que de \partial_\dots(1.0*hfn[k_x][k_y][k_z][i] qui est toute la partie non linéaire sur $f$ , faisant intervenir WENO, que je ne sais pas encore comment écrire. Mais l’autre chose que j’ai remarqué c’est que l’écriture de notre matrice $A$ laisse sous-entendre que $j_{c,x}$ et $j_{c,y}$ seront résolus dans l’espace réel, mais le brassage dû à $e^{A}$ implique que ces variables sont des variables de Fourier. Puisque la transformée de Fourier est définie par :

$\hat{f}(\xi,v) = \int_\mathbb{R} f(z,v)e^{-i2\pi\xi z}\,\mathrm{d}z$ on peut écrire que : $\begin{aligned} \widehat{\int_\mathbb{R}v_{x|y}\,f(z,v)\,\mathrm{d}v} &= \int_\mathbb{R} \int_\mathbb{R}v_{x|y}\,f(z,v)\,\mathrm{d}ve^{-i2\pi\xi z}\,\mathrm{d}z = \int_\mathbb{R}v_{x|y} \int_\mathbb{R}f(z,v)e^{-i2\pi\xi z}\,\mathrm{d}z\,\mathrm{d}v \\ &= \int_\mathbb{R} v_{x|y}\,\hat{f}\,\mathrm{d}v \end{aligned}$

Bref tout pour se faire dans l’espace de Fourier sauf le transport en $v$ (caché dans $\partial_{\dots}\left(\hat{f} \right)$ ). N’y a-t-il pas moyen d’effectuer ce transport dans l’espace de Fourier, ce qui éviterait des transformées de Fourier et transformées inverses, ainsi que le dédoublement de la mémoire nécessaire pour effectuer ces opérations ? En tout cas, vu que la seule opération nécessitant des variables dans l’espace réel est une opération ne modifiant pas $j_{c,x|y}$ , $E_{x|y}$ ou $B_{x|y}$ , je pense que lors de la simulation je vais directement travailler sur ces variables dans l’espace de Fourier (et faire les transformées inverses lorsque c’est nécessaire).

L’étape restante la plus problématique est l’appel au WENO en 3D, donc je pense que dans le courant de la semaine prochaine j’aurais mes premières simulations qui tourneront (ou à débuger…).

2020-09-07

Chaque étage se décompose de manière similaire, voici le détail pour le premier étage en pseudo-code.

Calcul des variables $\hat{j}_{c,x}^{(1)}$ , $\hat{j}_{c,y}^{(1)}$ , $\hat{E}_x^{(1)}$ , $\hat{E}_y^{(1)}$ , $\hat{B}_x^{(1)}$ et $\hat{B}_y^{(1)}$ :
- calcul des intégrales $\int v_x \hat{f}\,\mathrm{d}v$ et $\int v_y \hat{f}\,\mathrm{d}v$ :
  - $\left(\hat{j}_{h,x}\right)_i \gets \sum_{k_x,k_y,k_z} v_{k_x}\,\hat{f}_{i,k_x,k_y,k_z}\,\Delta v$
  - $\left(\hat{j}_{h,y}\right)_i \gets \sum_{k_x,k_y,k_z} v_{k_y}\,\hat{f}_{i,k_x,k_y,k_z}\,\Delta v$
- pour tout $x_i$ (code donné par sympy) :
  - $\hat{j}_{c,x}^{(1)}\gets \cdots$
  - $\hat{j}_{c,y}^{(1)}\gets \cdots$
  - $\hat{E}_{c,x}^{(1)}\gets \cdots$
  - $\hat{E}_{c,y}^{(1)}\gets \cdots$
  - $\hat{B}_{c,x}^{(1)}\gets \cdots$
  - $\hat{B}_{c,y}^{(1)}\gets \cdots$
Calcul de la variable $\hat{f}^{(1)}$ :
- transformée inverse de $\hat{f}^{n}$ :
  - $\left(f\right)_{i,k_x,k_z,k_z} \gets \text{iFFT}(\hat{f}^n_{\cdot,k_x,k_z,k_z})_i$ (boucle en $v_x$ , $v_y$ et $v_z$ )
- calcul de l’approximation de $(E_x+v_yB_0 + v_zB_y)\partial_{v_x}f + (E_y-v_xB_0+v_zB_x)\partial_{v_y}f +(v_xB_y - v_yB_x)\partial_{v_z}f$ :
  - pour tout $v_{k_x}$ :
    - pour tout $v_{k_y}$ :
      - pour tout $v_{k_y}$ :
        
        pour tout $v_{k_z}$ :
        
        pour tout $k_i$ :
        
        $\texttt{velocity\_vx} \gets E_{x,i}+v_{k_y}B_0 + v_{k_z}B_{y,i}$
        
        $\texttt{velocity\_vy} \gets E_{y,i}-v_{k_x}B_0 + v_{k_z}B_{x,i}$
        
        $\texttt{velocity\_vz} \gets v_{k_x}B_{y,i} - v_{k_y}B_{x,i}$
        
        $\partial_vf_{i,k_x,k_y,k_z}\gets \text{WENO}(\texttt{velocity\_vx},f_{i,k_x-3:k_x+3,k_y,k_z}) + \text{WENO}(\texttt{velocity\_vy},f_{i,k_x,k_y-3:k_y+3,k_z}) + \text{WENO}(\texttt{velocity\_vz},f_{i,k_x,k_y,k_z-3:k_z+3})$
- incrémentation de $\hat{f}^{(1)}$ :
  - pour tout $v_{k_x}$ :
    - pour tout $v_{k_y}$ :
      - pour tout $v_{k_y}$ :
        
        pour tout $v_{k_z}$ :
        
        $\left(\widehat{\partial_vf}\right)_i\gets\text{FFT}(\partial_vf_{\cdot,k_x,k_y,k_z})_i$
        
        pour tout $k_i$ :
        
        $\hat{f}^{(1)}_{i,k_x,k_y,k_z} \gets \hat{f}^n_{i,k_x,k_y,k_z} + \Delta t\widehat{\partial_vf}_i$

où $\text{WENO}$ est une fonction qui calcule une approximation de $a\partial_xu(x_i)$ en prenant en argument $a$ et le stencil $u_{i-3:i+3}$ . Cette fonction WENO diffère de ce que j’avais déjà pu implémenter avant, car il s’agissait alors d’une fonction renvoyant le flux numérique $\hat{f}^\pm_{i+\frac{1}{2}}$ . C’est-à-dire $\text{WENO}(a,u_{i-3},u_{i-2},u_{i-1},u_{i},u_{i+1},u_{i+2},u_{i+3})$ est une fonction renvoyant l’approximation : $a^+\frac{u^+_{i+\frac{1}{2}}-u^+_{i-\frac{1}{2}}}{\Delta x} + a^-\frac{u^-_{i+\frac{1}{2}}-u^-_{i-\frac{1}{2}}}{\Delta x}$ donc les flux numériques sont calculés 2 fois, ce qui engendre des pertes en temps de calcul, mais un gain important en utilisation de la mémoire (puisqu’il n’est pas nécessaire de sauvergarder tous les flux numériques avant calcul des différences).

2020-09-26

Toy model

On travaille sur le toy-model de Nicolas. Dans la résolution du modèle hybride-linéarisé de Vlasov-Maxwell en 3dv-1dz on se retrouve dans l’équation de Vlasov avec un terme en $(E+B\times v)\cdot\nabla_v f$ , et l’idée de ce toy-model est d’inclure le plus possible de terme dans la partie linéaire du schéma de Lawson, ici le terme en $B_0$ (composante en $z$ du vecteur $B = (B_x,B_y,B_0)^\textsf{T}$ ).

$\begin{cases} \partial_t f + (E+Jv)\cdot\nabla_v f = 0 \\ f(t=0,v) = f_0(v) \end{cases}$

avec :

$v = \begin{pmatrix}v_x \\ v_y\end{pmatrix} \quad E = \begin{pmatrix}E_x \\ E_y\end{pmatrix} \quad J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$

On pose $g(t,w) = f(t,e^{tJ}w)$ , soit $g(t,e^{-tJ}v) = f(t,v)$ , avec $w=(w_1,w_2)$ . On obtient alors :

$\partial_t g(t,w) = \partial_t\left(f(t,e^{e^{tJ}w}) \right) = -E\cdot\nabla_vf(t,e^{tJ}w)$ (pour le détail des calculs voir le document de Nicolas) et : $\partial_tg(t,w) = -E\cdot\nabla_vf(t,e^{tJ}w)$

On veut maintenir construire notre équation que vérifie $g$ :

$\begin{aligned} e^{-tJ}\nabla_w g(t,w) &= (\nabla_vf)(t,e^{tJ}w) \\ -E\cdot e^{-tJ}\nabla_wg(t,w) &= -E\cdot(\nabla_vf)(t,e^{tJ}w) \\ & = \partial_tg(t,w) \end{aligned}$

d’où :

$\partial_t g(t,w) + E\cdot e^{tJ}\nabla_wg(t,w) = 0$

on retrouve bien l’équation du document de Nicolas, avec comme condition initiale $g(t=0,w) = f_0(v)$ . Après avoir faire un schéma sur $g$ avec un schéma Runge-Kutta, om retrouve la solution sur $f$ par : $f(t,v) = g(t,e^{-tJ}v)$ .

Implémentation du toy-model

Je souhaite effectuer différents tests sur une équation de transport (simulation avec arithmétique stochastique, ou comparaison entre une méthode de splitting avec de la lazy evaluation et sans, etc.) donc le code pour simuler ce toy-model je change le repère de coordonnées en $(x,y)$ (plutôt qu’un $(v_x,v_y)$ correspondant plus au cas 1dz-3dv).

Donc avec ces notations le toy-model devient :

$\begin{cases} \partial_t u + (E+JX)\cdot\nabla_X f = 0 \\ u(t=0,x,y) = u_0(x,y) \end{cases}$

avec :

$X = \begin{pmatrix}x\\ y \end{pmatrix} \quad E = \begin{pmatrix}e_1\\ e_2 \end{pmatrix} \quad J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$

Ce qui nous donne la version filtrée suivante : $\partial_t g(t,z) = E\cdot e^{tJ}\nabla_z g(t,z)$ avec $z=(z_1,z_2)$ .

L’équation sur $u$ est résolu avec une méthode Runge-Kutta, un RK(3,3) suffira pour le test, et l’équation sur $g$ aussi. La CFL sur $u$ est du type : $\Delta t \leq \sigma\frac{\Delta x}{\|E+JX\|_\infty}$ , ce qui est grosso-modo une CFL en $\Delta t \leq \sigma\frac{\Delta x}{\|E\|_\infty + x_\text{max}}$ (on considère que domaine symétrie en $x$ et $y$ avec $(x|y)_\text{min} = -(x|y)_\text{max}$ , la notation $(x|y)$ indique indifféremment la variable $x$ ou $y$ ). La CFL de l’équation filtrée est, quant à elle, en $\Delta t\leq \sigma\frac{\Delta x}{\|E\|_\infty}$ , la taille du domaine n’est plus impactée dans la CFL.

Dans le cadre de la simulation du modèle 1dz-3dv, on sait que $E$ peut-être considéré comme petit.

2020-09-28

Questionnement de CFL

Pour mettre en place le toy-model j’ai commencé par coder un petit programme pour simuler l’équation :

$u_t + a u_x + b u_y = 0$

de manière suffisamment souple pour que $a$ et $b$ puissent être des vecteurs, ou des scalaires. Donc en particulier j’ai regardé le cas $a=-y$ et $b=x$ , pour obtenir une rotation (dans le sens trigonométrique si je ne dis pas trop de bêtises) :

$\partial_t u - y\partial_x u + x\partial_y u = 0$

J’effectue la simulation sur un domaine carré, centré en $(0,0)$ , avec le même maillage dans la direction $x$ et $y$ . Sur cette équation, on s’était toujours dit que la CFL devait être en $\Delta t \leq \sigma \frac{\Delta x}{x_\text{max}}$ , par analogie à une équation d’advection du type $u_t + au_x = 0$ où la CFL est en $\Delta t \leq \sigma \frac{\Delta x}{|a|}$ . Plus exactement la CFL est équivalente à considérer le splitting :

$\begin{aligned} \begin{cases}\partial_t u - y\partial_x u = 0 \end{cases} \\ \begin{cases}\partial_t u + x\partial_y u = 0 \end{cases} \end{aligned}$

où la CFL serait (si on résout chaque système avec une méthode de Runge-Kutta) en $\Delta t \leq \min\left(\sigma\frac{\Delta x}{y_\text{max}},\sigma\frac{\Delta y}{x_\text{max}}\right)$ .

Or… lorsque j’effectue une simulation avec cette CFL, cela explose. Puisque dans ma méthode Runge-Kutta je ne considère pas ce splitting, je résous l’équation suivante :

$\partial_t u + J\begin{pmatrix}x\\y\end{pmatrix}\cdot\nabla u = 0$

où la CFL sera (par analogie au cas 1d) en :

$\Delta t \leq \sigma \frac{\Delta x}{\left\|J\begin{pmatrix}x\\y\end{pmatrix}\right\|_\text{?}}$

On a remarqué que prendre $x_\text{max}$ , ce qui revient à prendre la norme $\left\|J\begin{pmatrix}x\\y\end{pmatrix}\right\|_\infty$ ne fonctionnait pas. Lorsque j’essaie de prendre la norme 2, donc $\sqrt{2}x_\text{max}$ ça ne fonctionne pas non plus, je semble obligé de prendre $2x_\text{max}$ ce qui revient à prendre la norme 1.

Je n’ai pas réfléchi à pourquoi la norme 1, autrement qu’en remarquant que c’est la seule qui semble fonctionner dans mon problème. Je n’ai pas remis en question le calcul de mon nombre de CFL $\sigma$ , calculé pour une équation 1d : $u_t + u_x = 0$ .

La question derrière tout cela étant aussi, comment cela se fait-il que je n’ai pas remarqué ce problème plus tôt ?

J’ai essayé d’implémenter le splitting en utilisant pour résoudre chaque sous-système une méthode Runge-Kutta RK(3,3), je m’attendais à ce que ça ne donne pas de bons résultats (j’ai un étirement de la solution plus qu’une rotation), mais je suis bien sous CFL en utilisant la condition CFL avec la norme infinie.

Sur cela, j’implémente la version toy-model, cela demande de préciser les variables $a$ et $b$ avec :

$a_{ij} = e_1 - y_j ,\quad b_{ij} = e_2 + x_i$

puis une autre version avec le filtrage :

$a(t) = e_1\cos(t) - e_2\sin(t),\quad b(t) = e_1\sin(t) + e_2\cos(t)$

2020-09-30

Nouvelles du filtrage

Le toy-model a été implémenté, et pour effectuer la rotation de fin : $f(t,v) = g(t,e^{-tJ}v)$ Nicolas m’a conseillé le splitting en 3 étapes de Joackim. Effectuer une rotation de $t$ revient à calculer la solution $h$ de l’équation $\partial_t h + Jv\cdot\partial_v h=0$ , $h(0)=g(t)$ (dans notre cas).

$\begin{aligned} \phi_1 &= \begin{cases} \partial_t h + v_2\partial_{v_1}h = 0 \end{cases} \\ \phi_2 &= \begin{cases} \partial_t h - v_1\partial_{v_2}h = 0 \end{cases}\\ \phi_3 &= \begin{cases} \partial_t h + v_2\partial_{v_1}h = 0 \end{cases} \end{aligned}$ et on a $f(t) = \phi_3^{[\tan(t/2)]}\circ\phi_2^{[\sin(t)]}\circ\phi_1^{[\tan(t/2)]}(g(t))$ . Chaque système est résolu par une méthode de Lagrange d’ordre 5.

On effectue une simulation du toy-model avec les paramètres numériques suivants :

$N_x=N_y=200$
$T_f = 20$
$\Omega = [-5,5]\times [-5,5]$
$f(t=0,x,y) = \exp\left( -\frac{x^2}{0.75} -\frac{(y-1)^2}{0.25} \right)$

Une première simulation sera effectuée sous la condition CFL de l’équation sur $f$ , ce qui nous donne :

$\Delta t \leq \sigma\frac{\Delta x}{\left\|E+Jv\right\|_1} \approx 0.00637255$

et une seconde sous la condition CFL de l’équation sur $g$ :

$\Delta t \leq \sigma\frac{\Delta x}{\left\|E\right\|_1} \approx 0.325$

Dans le cadre de WENO-RK(3,3) $\sigma \approx 1.3$ .

Résultats de simulation pour f, g et f après filtrage, pour \Delta t = 0.00637255 — Résultats de simulation pour $f$ , $g$ et $f$ après filtrage, pour $\Delta t = 0.00637255$

Résultats de simulation pour f, g et f après filtrage, pour \Delta t = 0.325 — Résultats de simulation pour $f$ , $g$ et $f$ après filtrage, pour $\Delta t = 0.325$

La solution calculée sur $f$ dans le cas où $\Delta t = 0.325$ ne donne que des nan, car on est complètement hors CFL, la solution sur $g$ reste bonne en norme de l’oeil.

Par contre… la simulation hybrid1dx3dv_lawson.cc explose toujours alors que je semble être largement sous CFL. Actuellement tourne un run avec $B_0=0$ . Les simulations de jeudi dernier qui concordaient bien entre la méthode de Lawson et la méthode de splitting étaient sur un maillage grossier, sur un maillage plus fin je semble toujours hors CFL alors que $\Delta t = 5\cdot 10^{-3}$ .

Est-ce que implémenter directement le filtrage pour le modèle hybride avec la méthode de Lawson ne permettrait pas de résoudre ce problème de CFL ?

Filtrage du hybride 1dz-3dv

J’ai lancé une simulation avec $B_0 = 0$ pour vérifier qu’il n’y avait pas d’autres erreurs dans le code de Lawson, et comparer cela avec la méthode de splitting, cela tourne très bien sans explosé à cause d’une contrainte de CFL. Donc d’après Nicolas, puisque de toute façon on voulait effectuer ce filtrage, je passe à l’implémentation de filtrage sans que la méthode de Lawson ne tourne (ou plus exactement sans comprendre d’où vient ce problème de CFL).

Je reprends les calculs de Nicolas, on pose $g(t,z,w,v_z) = f(t,z,e^{-tB_0J}v,v_z)$ pour obtenir l’équation sur $g$ qui est :

$\partial_t g + v_z\partial_z g - e^{-tB_0J}E\cdot\nabla_w g - \mathcal{B}g = 0$

avec $\mathcal{B}g = (\textrm{v}\times B)\cdot\nabla_\textrm{v}f$ ( $\textrm{v} = (v_x,v_y,v_z)^T$ ). J’obtiens finalement, après calcul :

$\partial_t g = \begin{pmatrix} -v_z \\ E_x\cos(B_0t) + E_y\sin(B_0t) + v_zB_x\sin(B_0t) - v_zB_y\cos(B_0t) \\ -E_x\sin(B_0t) + E_y\cos(B_0t) + v_zB_x\cos(B_0t) + v_zB_y\sin(B_0t) \\ -B_x\left( w_1\sin(B_0t) + w_2\cos(B_0t) \right) + B_y\left( w_1\cos(B_0t) - w_2\sin(B_0t) \right) \end{pmatrix} \cdot \begin{pmatrix} \partial_{z} \\ \partial_{w_1} \\ \partial_{w_2} \\ \partial_{v_z} \\ \end{pmatrix}$

Modif à faire dans le code

On a toujours une équation du type : $\partial_t u = -v_z\partial_z u + \texttt{velocity\_v1}\partial_{v_1} u + \texttt{velocity\_v2}\partial_{v_2} u + \texttt{velocity\_vz}\partial_{v_z} u$ , donc en reprenant le pseudo-code du 7 septembre, cela revient à simplement modifier les valeurs des variables velocity_vi.

Attention : les velocity_vi dépendent maintenant du temps, il est donc important de prendre en compte les pas de temps de chaque étage dans la méthode Runge-Kutta.

L’autre modification à effectuer est dans le calcul des flux $\int v_xf\,\mathrm{d}v$ et $\int v_yf\,\mathrm{d}v$ qui deviennent respectivement $\int (w_1\cos(B_0t) - w_2\sin(B_0t))g\mathrm{d}w\mathrm{d}v_z$ et $\int (w_1\sin(B_0t) + w_2\cos(B_0t))g\mathrm{d}w\mathrm{d}v_z$ . Ceux-ci dépendent également du temps.

2020-10-02

Reprise du calcul de filtrage pour le modèle 1dz3dv

On rappelle l’équation des particules chaudes dans le modèle hybride 1dz-3dv :

$\partial_tf + v_z\partial_zf - \left( E_x + v_yB_0 - v_zB_y \right)\partial_{v_x}f - \left( E_y - v_xB_0 + v_zB_x \right)\partial_{v_y}f - \left( v_xB_y - v_yB_x \right)\partial_{v_z}f = 0$

Équation que l’on réécrit comme :

$\partial_tf + v_z\partial_zf - \left( E + B_0Jv_\perp \right)\cdot\nabla_{v_\perp}f - \left( v\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix} \right)\cdot\nabla_vf = 0$

Les notations sont : $v = (v_\perp,v_z), \qquad v_\perp = (v_x,v_y)$

On pose le changement de variable : $v_\perp = e^{-tB_0J}w$ pour avoir une fonction $g$ :

$g(t,z,w,v_z) = f(t,z,e^{-tB_0J}w,v_z)$

On souhaite calculer $\partial_t g$ et connaître l’équation que vérifie $g$ :

$\begin{aligned} \partial_t g &= \partial_t\left(f\left(t,e^{-tB_0J}w\right)\right) = \left[\partial_tf\right](t,e^{-tB_0J}w) - \left[ B_0Je^{-tB_0J}w\cdot\nabla_{v_\perp}f \right](t,e^{-tB_0J}w) \\ % &= -v_z\partial_zf(e^{-tB_0J}w) + \left[ (E+B_0Jv_\perp)\cdot\nabla_{v_\perp}f \right](v_\perp=e^{-tB_0J}w) + \left[ \left(\begin{pmatrix}v_\perp\\ v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix}\right)\cdot\nabla_vf \right](v_\perp=e^{-tB_0J}w) - B_0Je^{-tB_0J}w\cdot\nabla_{v_\perp}f(e^{-tB_0J}w) \\ &= \left[ -v_z\partial_z f + (E+B_0Jv_{\perp})\cdot\nabla_{v_\perp}f + \left( \begin{pmatrix}v_perp \\ v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix} \right)\cdot\nabla_v f \right]( v_\perp = e^{-tB_0J}w ) - \left[ B_0Je^{-tB_0J}w\cdot\nabla_{v_\perp}f \right](v_\perp = e^{-tB_0J}w) \\ &= -v_z\partial_zf(e^{-tB_0J}w) + E\cdot\left[\nabla_{v_\perp}f\right](e^{-tB_0J}w) + B_0Je^{-tB_0J}w\cdot\left[\nabla_{v_\perp}f\right](e^{-tB_0J}w) + \begin{pmatrix}e^{-tB_0J}w\\v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix}\cdot\left[\nabla_vf\right](e^{-tB_0J}w) - B_0Je^{-tB_0J}w\cdot\nabla_{v_\perp}f(e^{-tB_0J}w) \end{aligned}$

Or :

$\begin{aligned} \nabla_w g(t,z,w,v_z) &= \nabla_wf(e^{-tB_0J}w) \\ &= e^{tB_0J}\left[\nabla_{v_\perp}f\right](e^{-tB_0J}w) \\ e^{-tB_0J}\nabla_wg(w) &= \left[ \nabla_{v_\perp}f \right](e^{-tB_0J}w) \end{aligned}$

d’où :

$\partial_t g = -v_z\partial_zg(t,z,w,v_z) + e^{tB_0J}E\cdot\nabla_wg + \left(\begin{pmatrix}e^{-tB_0J}w\\v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix}\right)\cdot\begin{pmatrix}e^{-tB_0J}\nabla_wg\\ \partial_{v_z}g\end{pmatrix}$

On souhaite avoir une expression plus précise du dernier terme, pour cela on note :

$e^{-tB_0J} = \begin{pmatrix}c&-s\\ s & c \end{pmatrix} \qquad \text{avec }\quad c = \cos(B_0t),\quad s=\sin(B_0t)$

On a donc :

$\begin{aligned} e^{tB_0J}E\cdot\nabla_wg &= \begin{pmatrix}c&s\\-s&c\end{pmatrix}\begin{pmatrix}E_x\\E_y\end{pmatrix}\cdot\nabla_wg \\ &= \begin{pmatrix}E_x c + E_y s \\ -E_x s + E_y c\end{pmatrix}\cdot\nabla_wg \end{aligned}$

et d’autre part :

$\begin{aligned} &\left(\begin{pmatrix}e^{-tB_0J}w\\v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix}\right)\cdot\begin{pmatrix}e^{-tB_0J}\nabla_wg\\ \partial_{v_z}g\end{pmatrix}\\ =&\left(\begin{pmatrix}cw_1 - sw_2\\sw_1 + cw_2\\v_z\end{pmatrix}\times\begin{pmatrix}B_x\\B_y\\0\end{pmatrix}\right)\cdot\begin{pmatrix}c\partial_{w_1}g - s\partial_{w_2}g\\ s\partial_{w_1}g + c\partial_{w_2}g \\ \partial_{v_z}g\end{pmatrix}\\ =& \begin{pmatrix} -v_zB_y \\ v_zB_x \\ (c w_1 - sw_2)B_y - (sw_1 + cw_2)B_x \end{pmatrix} \cdot \begin{pmatrix} c\partial_{w_1}g - s\partial_{w_2}g \\ s\partial_{w_1}g + c\partial_{w_2}g \\ \partial_{v_z}g \end{pmatrix} \end{aligned}$

soit :

$\partial_t g = -v_z\partial_z g + \begin{pmatrix} E_x c + E_y s + v_zB_xs - v_zB_yc \\ -E_x s + E_y c + v_zB_xc + v_zB_ys \\ -B_x(sw_1 + cw_2) + B_y(cw_1 - sw_2) \end{pmatrix}\cdot\begin{pmatrix} \partial_{w_1}g \\ \partial_{w_2}g \\ \partial_{v_z}g \end{pmatrix}$

Ce qui nous donne les valeurs de $\texttt{velocity\_v}_{x,y,z}$ : $\begin{aligned} \texttt{velocity\_vx} &= E_x c + E_y s + v_zB_xs - v_zB_yc \\ \texttt{velocity\_vy} &= -E_x s + E_y c + v_zB_xc + v_zB_ys \\ \texttt{velocity\_vz} &= -B_x(sw_1 + cw_2) + B_y(cw_1 - sw_2) \end{aligned}$

2021-03-18

des calculs dans la fonction de stabilité d’une méthode de Lawson matricielle

On étudie la résolution numérique de l’équation : $\dot{u} = Lu + N(u)$

Soit le schéma Lawson induit par la méthode RK(3,3) de Shu-Osher :

$\begin{aligned} u^{(1)} &= e^{\Delta t L}u^n + \Delta t e^{\Delta t L}N(u^n) \\ u^{(2)} &= \frac{3}{4}e^{\frac{\Delta t}{2}L}u^n + \frac{1}{4}e^{-\frac{\Delta t}{2}L}u^{(1)} + \frac{1}{4}\Delta t e^{-\frac{\Delta t}{2}L}N(u^{(1)}) \\ u^{n+1} &= \frac{1}{3}e^{\Delta t L}u^n + \frac{2}{3}e^{\frac{\Delta t}{2}L}u^{(2)} + \frac{2}{3}\Delta t e^{\frac{\Delta t}{2}L}N(u^{(2)}) \end{aligned}$

linéarisons (pour étudier l’ordre du schéma), c’est-à-dire considérons : $N:u\mapsto N\cdot u$ . Si $L$ et $N$ commutent, alors on trouve $u^{n+1} = e^{\Delta t L}\left[I+\Delta tN + \frac{\Delta t^2}{2}N^2 + \frac{\Delta t^3}{6}N^3\right]u^n$ , mais ne supposons pas cela ici.

$u^{(1)} = \left[ e^{\Delta t L} + \Delta t e^{\Delta t L}N \right]u^n$

Plusieurs fois dans ce calcul il est possible d’effectuer des factorisations par $(I+\Delta tN)$ , cela ne sera pas effectuer ici, car nous souhaitons développer le résultat, donc il ne sert finalement à rien d’effectuer ces jolies factorisation.

$\begin{aligned} u^{(2)} &= \left[ \frac{3}{4}e^{\frac{\Delta t}{2}L} + \frac{1}{4}e^{-\frac{\Delta t}{2}L}\left( e^{\Delta tL}+\Delta te^{\Delta tL}N \right) + \frac{1}{4}\Delta t e^{-\frac{\Delta t}{2}L}N\left( e^{\Delta tL}+\Delta te^{\Delta tL}N \right) \right]u^n \\ &= \left[ e^{\frac{\Delta t}{2}L} + \frac{1}{4}\Delta t\left(e^{\frac{\Delta t}{2}L}N + e^{-\frac{\Delta t}{2}L}Ne^{\Delta tL}\right) + \frac{1}{4}\Delta t^2 e^{-\frac{\Delta t}{2}L}Ne^{\Delta tL}N \right]u^n \end{aligned}$

Et enfin $u^{n+1}$ (après calcul) :

$\begin{aligned} u^{n+1} = \Big[ e^{\Delta tL} &+ \Delta t\left(\frac{2}{3}e^{\frac{\Delta t}{2}L}Ne^{\frac{\Delta t}{2}L}+\frac{1}{6}e^{\Delta tL}N + \frac{1}{6}Ne^{\Delta tL}\right) \\ & + \frac{\Delta t^2}{2}\left(\frac{1}{3}Ne^{\Delta tL}N + \frac{1}{3}e^{\frac{\Delta t}{2}L}Ne^{\frac{\Delta t}{2}L}N + \frac{1}{3} e^{\frac{\Delta t}{2}L}Ne^{-\frac{\Delta t}{2}L}Ne^{\Delta tL} \right) \\ & + \frac{\Delta t^3}{6}e^{\frac{\Delta t}{2}L}Ne^{-\frac{\Delta t}{2}L}Ne^{\Delta tL}N \Big]u^n \end{aligned}$

Si $L$ et $N$ commutent, on retrouve bien la forme souhaitée, ce qui laisse sous-entendre qu’il n’y a pas d’erreur dans le calcul.

Je ne sais comment étudier ce truc (mais au moins il est calculé) donc regardons le cas où $LN=NL$ :

$\begin{aligned} u^{n+1} &= e^{\Delta tL}\left[I + \Delta tN + \frac{\Delta t^2}{2}N^2 + \frac{\Delta t^3}{6}N^3\right]u^n \\ &= e^{\Delta tL}\left[e^{\Delta tN} + \mathcal{O}\left(\Delta t^4N^4\right)\right]u^n \end{aligned}$

On a supposé que $L$ et $N$ commutent donc $e^{L}e^{N} = e^{L+N}$ .

$u^{n+1} = e^{\Delta t(L+N)}u^n + e^{\Delta tL}\mathcal{O}\left(\Delta t^4N^4\right)u^n$

L’ordre est une notion asymptotique, que l’on obtient en linéarisant notre problème, il me semble donc normal de supposer que $L$ et $N$ commutent pour trouver formellement l’ordre 3 de la méthode.

Maintenant si on souhaite étudier l’erreur avec un approximant de Padé (ou un développement de Taylor peut-être dans un premier temps), je pense que l’on peut partir de la forme $u^{n+1}=e^{\Delta tL}\left[I + \Delta tN + \frac{\Delta t^2}{2}N^2 + \frac{\Delta t^3}{6}N^3\right]u^n$ où $L$ et $N$ commutent. Cela revient donc à étudier $u^{n+1} = \left[e^{\Delta tL} + \mathcal{O}\left((\Delta tL)^{\cdots}\right)\right]\left[e^{\Delta tN} + \mathcal{O}\left(\Delta t^4N^4\right)\right]u^n$ . On remarque alors que l’on a fait disparaître tout le problème, car $L$ et $N$ commutant cela revient à $u^{n+1} = e^{\Delta t(L+N)}+\mathcal{O}(\dots)$ .

Le problème est que je ne sais pas comment analyser le cas non-commutatif dans le cadre général, qui semble évidemment plus intéressant.

2021-04-17

Grands principes de Padé-Lawson-Runge-Kutta

On s’intéresse à la résolution d’un problème du type : $\dot{u} = Lu + N(u)$ Pour cela on se propose d’utiliser une méthode de Lawson induite par une méthode de type Runge-Kutta RK( $s$ , $n$ ), ici on se concentrera sur une méthode RK(4,4) qui forme les premiers étages de la méthode à pas de temps adaptatif DP4(3) :

$\begin{aligned} u^{(1)} & = e^{\frac{\Delta t}{2}L}u^n + \frac{1}{2}\Delta t e^{\frac{\Delta t}{2}L}N(u^n) \\ u^{(2)} & = e^{\frac{\Delta t}{2}L}u^n + \frac{1}{2}\Delta t e^{\frac{\Delta t}{2}L}N(u^{(1)}) \\ u^{(3)} & = e^{\Delta t L}u^n + \Delta t e^{\frac{\Delta t}{2}L} N(u^{(2)}) \\ u^{n+1} & = -\frac{1}{3}e^{\Delta t L}u^n + \frac{1}{3}e^{\frac{\Delta t}{2}L}u^{(1)} + \frac{2}{3}e^{\frac{\Delta t}{2}L}u^{(2)} + \frac{1}{3}u^{(3)} + \frac{1}{6}\Delta t N(u^{(3)}) \\ \end{aligned}$

Nous souhaitons utiliser cette méthode dans un contexte où le calcul de $e^{\alpha \Delta t L}$ est compliqué formellement (i.e. avec sympy), nous nous proposons de réécrire ce schéma en substituant la fonction $M\mapsto e^{M}$ par $M\mapsto \mathrm{P}_{p,q}(M)$ pour l’approximant de Padé d’ordre $(p,q)$ de $e^{M}$ .

Il est a noter que cette étape, d’un point de vue informatique, signifie simplement que l’on propose une autre implémentation de la fonction exp.

Le schéma devient (simple réécriture) : $\begin{aligned} u^{(1)} & = \mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)u^n + \frac{1}{2}\Delta t \mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)N(u^n) \\ u^{(2)} & = \mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)u^n + \frac{1}{2}\Delta t \mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)N(u^{(1)}) \\ u^{(3)} & = \mathrm{P}_{p,q}\left(\Delta tL\right)u^n + \Delta t \mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right) N(u^{(2)}) \\ u^{n+1} & = -\frac{1}{3}\mathrm{P}_{p,q}\left(\Delta t L\right)u^n + \frac{1}{3}\mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)u^{(1)} + \frac{2}{3}\mathrm{P}_{p,q}\left(\frac{\Delta t}{2}L\right)u^{(2)} + \frac{1}{3}u^{(3)} + \frac{1}{6}\Delta t N(u^{(3)}) \\ \end{aligned}$

où la fonction $\mathrm{P}_{p,q}$ est définie comme :

$\begin{aligned} h_{p.q}(z) = \sum_{i=0}^p \frac{\frac{p!}{(p-i)!}}{\frac{(p+q)!}{(p+q-i)!}}\frac{z^i}{i!} \\ k_{p.q}(z) = \sum_{j=0}^q (-1)^j \frac{\frac{q!}{(q-j)!}}{\frac{(p+q)!}{(p+q-j)!}} \frac{z^j}{j!} \\ e^{z} \approx \mathrm{P}_{p,q}(z) = \frac{h_{p,q}(z)}{k_{p,q}(z)} \end{aligned}$ pour une matrice quelconque $M$ , on a : $e^{M} \approx \mathrm{P}_{p,q}(M) = h_{p,q}(M)\cdot \left(k_{p,q}(M)\right)^{-1}$

Le calcul des $\mathrm{P}_{p,q}(\alpha\Delta t L)$ nécessite à chaque fois une inversion de matrice, ce qui a un coût numérique élevé. Une astuce consiste donc à considérer la fonction $t\mapsto\mathrm{P}(t\cdot L)$ et de l’évaluer en certaines valeurs pour $t=\alpha \Delta t$ . L’inversion n’est alors effectuée qu’une seule fois (lors de la construction de la fonction $\mathrm{P}^{[L]}_{p,q}:t\mapsto\mathrm{P}_{p,q}(t\cdot L)$ . L’évaluation de cette fonction en une valeur quelconque de $t$ (même un symbol comme $\Delta t$ pour construire le schéma) est quasi-immédiate.

Maintenant pour générer un code C++ il suffit de convertir les expressions sympy ainsi obtenues en des chaînes de caractères représentant du code C++.

Cette étape se fait en plusieurs petites sous-étapes qui sont simplement techniques :

Conversion des symbols en des noms de variables C++, et conversion des fonctions mathématiques en des fonctions ayant pour nom leur équivalent C++ (sp.sin $\mapsto$ std::sin par exemple)

Évaluation des fractions (pour être sûr de ne pas avoir de division entre entiers qui donneraient 0 en C++ (1/2 == 0 et 1./2==0.5 par exemple).

Conversion en chaînes de caractères de ces expressions sympy (normalement le plus gros des manipulations sur les expressions est fait).

Mise en forme de ces expressions dans un squelette de code (utilisation de jinja2 qui est un moteur de template).

Et contrairement à ce que je m’attendais, ce sont finalement ces étapes qui prennent le plus de temps dans la génération de code, le calcul de l’approximant de Padé (ne s’effectuant qu’une seule fois) est relativement rapide.

Étude de l’approximant de Padé pour Vlasov-Maxwell hybride linéarisé

Dans cette section on se propose d’étudier le comportement de l’approximant de Padé pour l’équation de Vlasov-Maxwell hybride linéarisé (VMHL) :

$\begin{cases} \partial_t j_{c,x} & = \Omega_{pe}^2E_x - j_{c,y}B_0 \\ \partial_t j_{c,y} & = \Omega_{pe}^2E_y + j_{c,x}B_0 \\ \partial_t B_{x} & = \partial_zE_y \\ \partial_t B_{y} & = -\partial_zE_x \\ \partial_t E_{x} & = -\partial_zB_y - j_{c,x} + \int v_xf_h\,\mathrm{d}v \\ \partial_t E_{y} & = \partial_zB_x - j_{c,y} + \int v_yf_h\,\mathrm{d}v \\ \partial_t f_h & = -v_z\partial_zf_h + (E_x+v_yB_0-v_zB_y)\partial_{v_x}f_h + (E_y-v_xB_0+v_zB_x)\partial_{v_y}f_h + (v_xB_y - v_yB_x)\partial_{v_z}f_h \end{cases}$

L’opération de filtrage n’implique pas de grandes modifications dans l’implémentation, de plus celle-ci touche seulement l’équation de Vlasov sur la variable $f_h$ et le calcul des courants $j_{h,\star} = \int v_{\star}f_h\,\mathrm{d}v$ , et ce que nous allons présenter ici touche très peu ces termes.

Ce système peut s’écrire sous la forme : $\dot{U} = \underbrace{\begin{pmatrix}A & 0 \\ 0 & -ikv_z \end{pmatrix}}_{L} U + N(U)$ Il est nécessaire dans un schéma de Lawson de calculer l’exponentielle de la partie linéaire $L$ à différents temps. La matrice $L$ est une matrice diagonale par blocs, et l’exponentielle d’une matrice diagonale est l’exponentielle des termes diagonaux, c’est-à-dire : $e^{tL} = \begin{pmatrix}e^{tA} & 0 \\ 0 & e^{-ikv_z} \end{pmatrix}$ Si on souhaite utiliser un approximant de Padé, il est possible de le faire seulement pour approximer $e^{tA}$ , c’est-à-dire : $e^{tL} \approx \begin{pmatrix}\mathrm{P}_{p,q}(tA) & 0 \\ 0 & e^{-ikv_z} \end{pmatrix}$ Ceci permet de minimiser la taille de la matrice à inverser dans le calcul de $\mathrm{P}_{p.q}(M)$ , minimiser la taille des expressions à manipuler avec sympy (espace mémoire et temps de calcul sur les arbres d’expression), et conserver exacte l’équation sur $f_h$ .

Maintenant nous nous proposons d’étudier le comportement de l’approximant de Padé, tout d’abord avec le choix de partie linéaire sans les équations de Maxwell, et ensuite avec les équations de Maxwell.

Partie linéaire sans les équations de Maxwell

Lorsque la partie linéaire ne contient pas les équations de Maxwell, la matrice de la partie linéaire $L$ est : $L = \begin{pmatrix} 0 & -B_0 & 0 & 0 & \Omega_{pe}^2 & 0 & 0 \\ B_0 & 0 & 0 & 0 & 0 & \Omega_{pe}^2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -v_z\partial_z \\ \end{pmatrix}$

Pour cette matrice sympy permet de calculer $e^{tL}$ pour tout $t$ et tout $k$ ( $k$ intervenant dans la transformée de Fourier en $z$ de la dérivée spatiale dans l’équation de Vlasov). Nous allons surtout étudier, pour $B_0=1$ et $\Omega_{pe}^2=4$ , la matrice réduite $A$ , sans l’équation de Vlasov : $A = \begin{pmatrix} 0 & -1 & 0 & 0 & 4 & 0 \\ 1 & 0 & 0 & 0 & 0 & 4 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 0 \\ \end{pmatrix}$ et regarder $e^{tA}$ pour différentes valeurs de $t$ et comparer $\mathrm{P}_{p,p}(tA)$ pour ces mêmes valeurs de $t$ et différentes valeurs de $p\in\mathbb{N}$ .

Tout d’abord on regarde $e^{t A}$ : $\begin{pmatrix} 0.378732187481834 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518167 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.378732187481834 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.621267812518167 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.970142500145332 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & - 0.970142500145332 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)}\\- 0.378732187481834 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518166 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.378732187481834 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.621267812518167 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.970142500145332 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.970142500145332 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.970142500145332 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.970142500145332 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)}\\0 & 0 & 1.0 & 0 & 0 & 0\\0 & 0 & 0 & 1.0 & 0 & 0\\- 0.242535625036333 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.242535625036333 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & 0.621267812518166 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481834 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.621267812518167 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.378732187481834 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)}\\- 0.242535625036333 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.242535625036333 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & - 0.242535625036333 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} - 0.242535625036333 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0 & 0 & - 0.621267812518166 \sin{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481834 \sin{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} & 0.621267812518167 \cos{\left(\frac{\sqrt{2} t \sqrt{9 - \sqrt{17}}}{2} \right)} + 0.378732187481834 \cos{\left(\frac{\sqrt{2} t \sqrt{\sqrt{17} + 9}}{2} \right)} \end{pmatrix}$ à comparer avec $\mathrm{P}_{3,3}(t A)$ : $\begin{pmatrix} - \frac{2 t^{2} \left(3240 t^{8} - 5985 t^{6} - 48600 t^{4} - 486000 t^{2}\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{3 t^{3}}{40} - \frac{t}{2}\right) \left(480 t^{9} - 16200 t^{7} - 64800 t^{5} - 243000 t^{3} + 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & - \frac{2 t^{2} \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 480 t^{9} + 16200 t^{7} + 64800 t^{5} + 243000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(- 3240 t^{8} + 5985 t^{6} + 48600 t^{4} + 486000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{3 t^{3}}{40} - \frac{t}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & 0 & 0 & - \frac{2 t^{2} \left(1080 t^{9} - 26370 t^{7} - 97200 t^{5} - 432000 t^{3}\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{3 t^{3}}{40} - \frac{t}{2}\right) \left(- 12960 t^{8} + 23940 t^{6} + 194400 t^{4} + 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & - \frac{2 t^{2} \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(12960 t^{8} - 23940 t^{6} - 194400 t^{4} - 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(- 1080 t^{9} + 26370 t^{7} + 97200 t^{5} + 432000 t^{3}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{3 t^{3}}{40} - \frac{t}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000}\\ \frac{2 t^{2} \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(480 t^{9} - 16200 t^{7} - 64800 t^{5} - 243000 t^{3} + 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(3240 t^{8} - 5985 t^{6} - 48600 t^{4} - 486000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{3 t^{3}}{40} + \frac{t}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & \frac{2 t^{2} \left(- 3240 t^{8} + 5985 t^{6} + 48600 t^{4} + 486000 t^{2}\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{3 t^{3}}{40} + \frac{t}{2}\right) \left(- 480 t^{9} + 16200 t^{7} + 64800 t^{5} + 243000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & 0 & 0 & \frac{2 t^{2} \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 12960 t^{8} + 23940 t^{6} + 194400 t^{4} + 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(1080 t^{9} - 26370 t^{7} - 97200 t^{5} - 432000 t^{3}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{3 t^{3}}{40} + \frac{t}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & \frac{2 t^{2} \left(- 1080 t^{9} + 26370 t^{7} + 97200 t^{5} + 432000 t^{3}\right)}{5 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{t^{2}}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{t^{3}}{6} + 2 t\right) \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(- \frac{3 t^{3}}{40} + \frac{t}{2}\right) \left(12960 t^{8} - 23940 t^{6} - 194400 t^{4} - 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000}\\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0\\ \frac{t^{3} \left(3240 t^{8} - 5985 t^{6} - 48600 t^{4} - 486000 t^{2}\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{t^{2} \left(480 t^{9} - 16200 t^{7} - 64800 t^{5} - 243000 t^{3} + 1620000 t\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & \frac{t^{3} \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{t^{2} \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(- 3240 t^{8} + 5985 t^{6} + 48600 t^{4} + 486000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 480 t^{9} + 16200 t^{7} + 64800 t^{5} + 243000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & 0 & 0 & \frac{t^{3} \left(1080 t^{9} - 26370 t^{7} - 97200 t^{5} - 432000 t^{3}\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{t^{2} \left(- 12960 t^{8} + 23940 t^{6} + 194400 t^{4} + 1944000 t^{2}\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & \frac{t^{3} \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{t^{2} \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(- 1080 t^{9} + 26370 t^{7} + 97200 t^{5} + 432000 t^{3}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(12960 t^{8} - 23940 t^{6} - 194400 t^{4} - 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000}\\ - \frac{t^{3} \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} - \frac{t^{2} \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(3240 t^{8} - 5985 t^{6} - 48600 t^{4} - 486000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(480 t^{9} - 16200 t^{7} - 64800 t^{5} - 243000 t^{3} + 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & - \frac{t^{3} \left(- 3240 t^{8} + 5985 t^{6} + 48600 t^{4} + 486000 t^{2}\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} - \frac{t^{2} \left(- 480 t^{9} + 16200 t^{7} + 64800 t^{5} + 243000 t^{3} - 1620000 t\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(600 t^{9} - 10170 t^{7} - 32400 t^{5} - 189000 t^{3} - 1620000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 7200 t^{8} - 8460 t^{6} - 151200 t^{4} - 972000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & 0 & 0 & - \frac{t^{3} \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} - \frac{t^{2} \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(1080 t^{9} - 26370 t^{7} - 97200 t^{5} - 432000 t^{3}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 12960 t^{8} + 23940 t^{6} + 194400 t^{4} + 1944000 t^{2}\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} & - \frac{t^{3} \left(- 1080 t^{9} + 26370 t^{7} + 97200 t^{5} + 432000 t^{3}\right)}{30 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} - \frac{t^{2} \left(12960 t^{8} - 23940 t^{6} - 194400 t^{4} - 1944000 t^{2}\right)}{10 \left(64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000\right)} + \frac{\left(1 - \frac{2 t^{2}}{5}\right) \left(- 10440 t^{8} - 2475 t^{6} - 102600 t^{4} - 486000 t^{2} + 3240000\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} + \frac{\left(\frac{t^{3}}{24} - \frac{t}{2}\right) \left(- 2400 t^{9} + 40680 t^{7} + 129600 t^{5} + 756000 t^{3} + 6480000 t\right)}{64 t^{12} + 864 t^{10} + 11124 t^{8} + 105705 t^{6} + 394200 t^{4} + 1458000 t^{2} + 3240000} \end{pmatrix}$

Valeurs propres de e^A et comparaison avec les valeurs propre de \mathrm{P}_{p,p}(A) pour p=1,\dots,14 — Valeurs propres de $e^A$ et comparaison avec les valeurs propre de $\mathrm{P}_{p,p}(A)$ pour $p=1,\dots,14$

Comparaison de e^{tA} et \mathrm{P}_{p,p}(tA) pour p=1,\dots,14 pour t\in[-1,1]. Ceci correspond à l’évaluation de e^{\alpha\Delta tA} pour des valeurs de \alpha parfois négatives et un \Delta t<1, bref un cas classique dans un schéma de Lawson. — Comparaison de $e^{tA}$ et $\mathrm{P}_{p,p}(tA)$ pour $p=1,\dots,14$ pour $t\in[-1,1]$ . Ceci correspond à l’évaluation de $e^{\alpha\Delta tA}$ pour des valeurs de $\alpha$ parfois négatives et un $\Delta t<1$ , bref un cas classique dans un schéma de Lawson.

Étude préliminaire de l’approximant de Padé de la matrice réduite de la partie linéaire sans les équations de Maxwell

On remarque que les valeurs propres de la matrice $\mathrm{P}_{p,p}(A)$ convergent bien vers celles de $e^A$ , et qu’elles sont très proches numériquement à partir de $n\geq 3$ . La norme (ici de Frobenius) de $\mathrm{P}_{p,p}(tA)$ converge bien vers celle de $e^{tA}$ , et à partir de $n\geq 10$ l’erreur entre les deux est inférieure à l’erreur machine. Le paramètre $t\in[-1,1]$ représente le pas de temps $\Delta t$ . Les différents étages d’une méthode de Lawson basée sur RK(3,3) ou RK(4,4) demandent le calcul de $e^{\Delta t A}$ , $e^{\frac{\Delta t}{2}A}$ et $e^{-\frac{\Delta t}{2}A}$ (en tout cas dans les méthodes que je considère), de plus le pas de temps $\Delta t$ est généralement pris petit (dans ce cas à cause de la CFL induit par les équations de Maxwell dans la partie non-linéaire), ce qui justifie $t\in[-1,1]$ .

La génération dans ce cadre là fonctionne sans nécessité d’assistance particulière à partir des fonctions de bases de sympy. Mais puisque j’ai rencontré des problèmes dans le contexte (qui est celui qui nous intéresse le plus) avec les équations de Maxwell dans la partie linéaire, j’ai re-coder une fonction de calcul de déterminant et d’inverse de matrice. Ces nouvelles fonction n’effectuent aucune vérification sur l’inversibilité de la matrice, pas de recherche de la meilleure ligne ou colonne pour effectuer le développement du déterminant, etc. Il s’agit d’implémentation naïves ( $M^{-1} = \frac{1}{\det{M}}^t\mathrm{com}A$ ) mais qui tournent bien plus rapidement que leur équivalent sympy (qui font des vérifications sur la validité du calcul). Les expressions que j’obtiens ainsi sont très longues, donc leur traitement avec sympy est relativement long aussi.

Résultats numériques

Et maintenant des résultats de simulations sur le maillage $N_z\times N_{v_x}\times N_{v_y}\times N_{v_z} = 15\times 20\times 20\times 41$ , et $\Delta t=0.1$ avec une méthode Lawson-RK(4,4), et un approximant de Padé $\mathrm{P}_{2,2}$ , bref normalement des conditions pas terribles pour faire tourner le schéma.

L’évolution temporelle de l’énergie électrique, magnétique et cinétique des particules froides

On voit que les résultats sont très semblables entre la méthode Lawson-RK(4,4) et la méthode Lawson-RK(4,4) avec approximant de Padé (2,2). La pente est bien capturée malgré le maillage. Ce résultat est confirmée avec l’erreur relative sur l’énergie totale, qui sont aussi très similaires.

L’évolution temporelle de l’erreur relative sur l’énergie totale

L’évolution temporelle de l’erreur entre les deux méthodes sur les différentes énergies

On peut tout de même mesurer une différence dans la résolution des deux méthodes, qui est de l’ordre de $10^{-7}$ jusqu’à $10^{-5}$ selon l’énergie considérée (ce ne sont pas des erreurs relatives !). Cette différence est simplement calculée en faisant la différence à chaque temps entre les deux énergies : $x^n_{\textrm{exact Lawson}} - x^n_{\textrm{Padé Lawson}}\,,\quad \forall t^n$

avec $x^n_{\star}$ correspondant à une énergie (énergie électrique, énergie magnétique ou énergie cinétique des particules froides) au temps $t^n$ pour la méthode $\star$ (soit la méthode de Lawson exacte ou la méthode de Lawson avec approximant de Padé (2,2)).

J’ai un doute sur l’exactitude du résultat de la différence, car je me demande si je n’ai pas un problème de décalage sur la première itération. Dans tous les cas l’ordre de grandeur reste bon, et on observe bien la phase linéaire (où l’erreur est très faible), et la phase non-linéaire, où les champs sont plus important, ce qui pourrait expliquer la différence entre les méthodes.

Temps de calcul des simulation avec le maillage $N_z\times N_{v_x}\times N_{v_y}\times N_{v_z} = 15\times 20\times 20\times 41$ , et $\Delta t=0.1$ .
Simulation	Temps de calul
LRK(4,4)	14h 30min 33s
P(2,2)-LRK(4,4)	14h 52min 28s

La différence de temps, sur une exécution, n’est pas significative, donc il ne semble pas y avoir de surcoût particulier à utiliser un approximant de Padé (en tout cas pour l’ordre faible).

Partie linéaire avec les équations de Maxwell

La matrice de la partie linéaire avec les équations de Maxwell est : $L = \begin{pmatrix} 0 & -B_0 & 0 & 0 & \Omega_{pe}^2 & 0 & 0 \\ B_0 & 0 & 0 & 0 & 0 & \Omega_{pe}^2 & 0 \\ 0 & 0 & 0 & 0 & 0 & \partial_z & 0 \\ 0 & 0 & 0 & 0 & -\partial_z & 0 & 0 \\ -1 & 0 & 0 & -\partial_z & 0 & 0 & 0 \\ 0 & -1 & \partial_z & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -v_z\partial_z \\ \end{pmatrix}$

De même pour ce cas nous allons considérer seulement la matrice réduite $A$ : $A = \begin{pmatrix} 0 & -1 & 0 & 0 & 4 & 0 \\ 1 & 0 & 0 & 0 & 0 & 4 \\ 0 & 0 & 0 & 0 & 0 & ik \\ 0 & 0 & 0 & 0 & -ik & 0 \\ -1 & 0 & 0 & -ik & 0 & 0 \\ 0 & -1 & ik & 0 & 0 & 0 \\ \end{pmatrix}$ où les termes $ik$ proviennent de la transformée de Fourier dans la direction $z$ des dérivées spatiales.

sympy étant incapable de caluler $e^{tA}$ pour une valeur de $k$ donnée, tout les calculs seront entièrement réalisé avec numpy.

Pour différentes valeurs de $k$ on trace dans un premier temps les valeurs propre de $e^{A}$ ainsi que celles de $\mathrm{P}_{3,3}(A)$ . On choisit l’approximant d’ordre 3, car les résultats sans les équations de Maxwell nous indiquent qu’à partir de $p\geq 3$ , les valeurs propres de $e^A$ et $\mathrm{P}_{n.n}(A)$ sont numériquement très proches. J’aurai pu, par sécurité, prendre une valeur plus grande.

Valeurs propre de e^{A} pour différentes valeurs de k — Valeurs propre de $e^{A}$ pour différentes valeurs de $k$

Valeurs propre de \mathrm{P}_{3,3}(A) pour différentes valeurs de k — Valeurs propre de $\mathrm{P}_{3,3}(A)$ pour différentes valeurs de $k$

Comparaison des valeurs propres de $e^{A}$ et de $\mathrm{P}_{3,3}(A)$ pour différentes valeurs de $k\in[\![-N_z,N_z]\!]$ .

On remarque que les valeurs propres sont bien sur le cercle unité, mais leur répartition est légèrement différente, cela est peut-être dû à l’ordre trop faible de $\mathrm{P}_{3,3}$ .

On peut aussi tracer, comme pour le cas sans les équations de Maxwell, l’erreur entre $e^{tA}$ et $\mathrm{P}_{n,n}(tA)$ pour $t\in[-1,1]$ et $p\in[\![1,14]\!]$ , j’ai fixé regardé ici 2 modes, $k=2,15$ (pas de raison particulière).

Erreur de Frobenius entre e^{tA} et \mathrm{P}_{p,p}(tA) pour k=2 — Erreur de Frobenius entre $e^{tA}$ et $\mathrm{P}_{p,p}(tA)$ pour $k=2$

Erreur de Frobenius entre e^{tA} et \mathrm{P}_{p,p}(tA) pour k=15 — Erreur de Frobenius entre $e^{tA}$ et $\mathrm{P}_{p,p}(tA)$ pour $k=15$

Erreur entre $e^{tA}$ et $\mathrm{P}_{p,p}(tA)$ pour deux modes $k$

Dans la pratique on choisit un approximant de Padé d’ordre particulier, ici ça sera 2 (choix purement arbitraire), et on regarde l’erreur entre $e^{tA}$ et $\mathrm{P}_{2,2}(tA)$ pour différents modes de Fourier $k$ .

Erreur de Frobenius entre e^{tA} et \mathrm{P}_{2,2}(tA) — Erreur de Frobenius entre $e^{tA}$ et $\mathrm{P}_{2,2}(tA)$

2021-04-21

À propos des valeurs propres de Padé

On souhaite s’assurer de la stabilité d’une méthode Padé-Lawson-Runge-Kutta, pour cela il semble nécessaire que les valeurs propres de Padé soient de norme inférieure ou égale à 1, pour se retrouver dans un cas similaire à la stabilité d’une méthode de type Runge-Kutta.

On commence l’étude par un approximant de Padé (2,2).

Module des valeurs propres de \mathrm{P}_{2,2}(A(k)) — Module des valeurs propres de $\mathrm{P}_{2,2}(A(k))$

Numériquement dans ce cas, toutes les valeurs propres sont de module 1 (l’échelle est en $10^{-12}+1$ ), la monté en ordre dans l’approximant de Padé donne les mêmes résultats.

Module des valeurs propres de \mathrm{P}_{3,2}(A(k)) — Module des valeurs propres de $\mathrm{P}_{3,2}(A(k))$

On remarque que si le numérateur est de degré plus élevé, l’approximant de Padé toutes les valeurs propres sont de module supérieur à 1 (la figure de droite n’est qu’un zoom proche de 1 pour observer le comportement de la ligne de valeurs propres visible à gauche). Cela implique l’instabilité d’une méthode induite par l’approximant de Padé $\mathrm{P}_{3,2}$ , et on peut, sans problème généralisé pour tout $n$ et $m$ tels que $n>m$ .

Module des valeurs propres de \mathrm{P}_{2,3}(A(k)) — Module des valeurs propres de $\mathrm{P}_{2,3}(A(k))$

Dans le cas de l’approximant de Padé $\mathrm{P}_{2,3}(A(k))$ , les valeurs propres de l’approximant de Padé sont toutes de module inférieur à 1, ce qui laisse sous-entendre la stabilité de la méthode de Lawson induite par cet approximant. On peut très certainement généralisé pour tout $n$ et $m$ tels que $m>n$ .

Et parce que les figures sont jolies on peut aussi regarder la disposition des valeurs propres.

Valeurs propre de \mathrm{P}_{3,2}(A(k)) — Valeurs propre de $\mathrm{P}_{3,2}(A(k))$

Valeurs propre de \mathrm{P}_{2,3}(A(k)) — Valeurs propre de $\mathrm{P}_{2,3}(A(k))$

Valeurs prorpres de $\mathrm{P}_{3,2}(A(k))$ (à gauche) et $\mathrm{P}_{2,3}(A(k))$ (à droite).

2021-07-09

Résultats Lawson approximé par Taylor ou Padé

Je fais ça rapidement, histoire de discuter des résultats numériques pour savoir ce que je garde vraiment ou non. À savoir que précédemment on a déjà illustré que les valeurs propres de la troncature de la série de Taylor ne sont (pas du tout) sur le cercle unité, alors que c’est le cas pour l’apporixmant de Padé. Cela fait dire à Nicolas qu’il faut tout bonnement le virer des résultats numérique, je trouve qu’il donne de meilleurs résultats qu’attendu (surtout lorsqu’on regarde les valeurs propres de $T_5(A)$ ).

Je n’ai pas cherché à avoir de jolies explosions avec une rotation (même code que pour obtenir l’ordre des Lawson-Padé-Taylor mais en mode rotation). En effet ce code engendre des modes élevés (car multiplié par la vitesse qui dépend des bornes du domaine) et il est possible de faire exploser Taylor sans faire exploser Padé. Piste à creuser peut-être.

Vérification

Premier lot de runs, on utilise la même partie linéaire que dans le papier (donc les équations de Maxwell ne sont pas dans la partie linéaire). Il s’agit juste de vérifier la possibilité de la méthode.

Pour le moment rien de passionant, on cherche juste à vérifier que l’on peut approcher l’exponentielle de la partie linéaire sans encore en profiter (une proof of concept).

Recherche d’un plus grand $\Delta t$

Maintenant on met les équations de Maxwell dans la partie linéaire, on regarde si on peut prendre des pas de temps plus grand (des tests précédent m’ont fait remarquer que $\Delta t=0.5$ était un peu trop grand, on le confirmera après avec le pas de temps adaptatif). On prend le même $\Delta t=0.05$ que dans le papier, et un pas de temps plus grand $\Delta t=0.1$ qui devient intéressant.

Run pas très passionant non plus, ça marche, c’est cool !

Taylor

Exploration de différentes troncatures de la série de Taylor, de l’ordre 1 à 4, avec un RK(3,3). (Les runs précédent étaient avec un RK(4,4), il s’agit ici de faire des runs normalement plus rapide pour tester des trucs…)

On remarque que pour les Taylor d’ordre plus faible ça explose, je ne sais pas exactement pourquoi.

Padé

Exploration de différents approximants de Padé, de l’ordre $(1,1)$ à $(2,2)$ , et il n’y a que le Padé $(2,1)$ qui explose, ce qui s’explique parce que le numérateur a un degré supérieur au dénominateur.

Dormand-Prince

On regarde maintenant avec la méthode DP4(3), entre Taylor et Padé. On regarde tout d’abord les énergies, Padé prend de très grands pas de temps en début de simulation, ce qui explique le shift dans les résultats (la bonne solution est celle donnée par Taylor).

On regarde Maintenant l’historique de la taille du pas de temps, et on remarque que la méthode de Padé est nettement plus intéressante. Je présume que c’est le fait que Taylor soit, en théorie, instable qui justifie qu’il n’essaie pas de prendre de grands pas de temps, et se restreint à peu près à la CFL de la partie non-linéaire. On peut aussi regarder l’erreur locale, comme dans le cas du papier on a beaucoup d’itérations rejettées (je n’ai pas encore regardé en détail combien).

On peut regarder la contribution de chaque erreur, j’ai rassembler les courants froids, les champs magnétiques et électriques pour ne pas voir 7 contributions. On remarque que l’erreur provient très largement de $E$ et de $f_h$ , ce qui correspond à ce qui est résolu dans la partie non-linéaire. Les blancs sont dûs aux itérations rejettées ou d’autres problèmes (purement techniques) de mon fichier de sauvegarde.

Par curiosité j’ai aussi regardé l’évolution du temps courant au cours des itérations.

2021-07-22

À propos des temps de calcul

Après un petit script bash pour récupérer les temps des simus, voici les résultats :

cas test	simu	temps
pade	vmhllf_mp11rk33	0d 05h 38m 50s
pade	vmhllf_mp21rk33	0d 05h 40m 19s
pade	vmhllf_mp22rk33	0d 05h 38m 29s
pade	vmhllf_mp12rk33	0d 05h 37m 47s
dtn	vmhllf_mp22dp43	0d 03h 49m 01s
dtn	vmhllf_mt5dp43	0d 06h 27m 39s
max/dt0p05	vmhllf_mp22rk44	0d 14h 13m 56s
max/dt0p05	vmhllf_mt5rk44	0d 14h 18m 13s
max/dt0p1	vmhllf_mp22rk44	0d 07h 09m 54s
max/dt0p1	vmhllf_mt5rk44	0d 07h 06m 30s
max/dt0p1	vmhllf_p22rk44	0d 07h 10m 59s
max/dt0p1	vmhllf_t5rk44	0d 07h 10m 07s
max/dt0p12	vmhllf_t5rk44	0d 05h 56m 46s
max/dt0p12	vmhllf_p22rk44	0d 06h 06m 08s
max/dt0p12	vmhllf_mt5rk44	0d 05h 58m 18s
max/dt0p12	vmhllf_mp22rk44	0d 05h 59m 38s
verif	vmhllf_rk44	0d 14h 10m 14s
verif	vmhllf_p22rk44	0d 14h 26m 49s
verif	vmhllf_t5rk44	0d 14h 21m 53s
taylor	vmhllf_mt4rk33	0d 05h 38m 29s
taylor	vmhllf_mt2rk33	0d 05h 37m 54s
taylor	vmhllf_mt3rk33	0d 05h 39m 09s
taylor	vmhllf_mt1rk33	0d 05h 38m 45s

Toutes les simus sont faites sur le maillage $N_z\times N_{v_x} \times N_{v_y} \times N_{v_z}=27\times 32\times 32 \times 41$ , jusqu’au temps final $T_f=200$ . Les paramètres spécifiques aux cas tests sont (la valeur dt0 est la valeur du pas de temps initial, qui est le pas de temps pour toute la simu si la méthode n’est pas à pas de temps adaptatif) :

pade : dt0=0.1
dtn : dt0=0.1 (pas de temps adaptatif)
max/dt0p05 : dt0=0.05
max/dt0p1 : dt0=0.1
max/dt0p12 : dt0=0.12
verif : dt0=0.05
taylor : dt0=0.1

Le nom des simulations est formé de la manière suivante : vmhllf pour Vlasov Maxwell Hybride linéarisé Lawson filtré puis un m pour indiquer si Maxwell est dans la partie linéaire, puis la méthode d’approximation avec son ordre et enfin la méthode en temps rk33, rk44 ou dp43.

On compare que les simulations où le pas de temps était $\Delta t = 0.1$ (pade, max/dt0p1/ et taylor), dans ce cas on remarque que le coup moyen par étage est de 1h45 (les simulations pade et taylor sont avec du RK(3,3) contre du RK(4,4) sur max/dt0p1). La seule simulation avec une méthode de Lawson standard est la simu verif avec vmhllf_rk44, et on remarque qu’il n’y a pas de surcoût significatif dans l’usage d’approximants de Padé ou de troncature de la série de Taylor.

J’ai lancé d’autres runs juste pour mesurer les temps de calcul sur le même cas test avec différents schémas, voici les résultats.

cas test	simu	(s)	temps	last time
chap3/time	vmhllf_mt5rk44	50403	0d 14h 00m 03s	200,00
chap3/time	vmhllf_mp22rk33	39326	0d 10h 55m 26s	200,00
chap3/time	vmhllf_mp22dp43	14984	0d 04h 09m 44s	200,00
chap3/time	vmhllf_mp11rk33	39251	0d 10h 54m 11s	200,00
chap3/time	vmhllf_mp22rk44	50399	0d 13h 59m 59s	200,00
chap3/time	vmhls_lie	48310	0d 13h 25m 10s	200,00
chap3/time	vmhls_suzuki	259173	2d 23h 59m 33s	191,75
chap3/time	vmhls_suzuki	270324	3d 03h 05m 24s	200,00
chap3/time	vmhllf_rk44	50775	0d 14h 06m 15s	200,00
chap3/time	vmhllf_dp43	42244	0d 11h 44m 04s	200,00
chap3/time	vmhls_strang	61794	0d 17h 09m 54s	200,00
chap3/time	vmhllf	41349	0d 11h 29m 09s	200,00
chap3/time	vmhllf_mt4rk33	39220	0d 10h 53m 40s	200,00

Le cas test est : $N_z \times N_{v_x} \times N_{v_y} \times N_{v_z}=27\times32\times32\times41$ et $\Delta t_0 = 0.05$ (pas de temps constant pour la majorité des simus, sauf la simu DP4(3) avec P(2,2) et la simu DP4(3) sans Maxwell dans la partie linéaire).

Pour la simulation avec Suzuki, le serveur de calcul PlaFRIM ne permet pas par défaut de lancer des jobs de plus de 3 jours, donc le temps indiqué pour la simulation Suzuki qui va jusqu’au temps 200 est une estimation.

Juste pour info j’ajoute aussi l’erreur relative sur ce cas test (les énergies sont toutes similaires et conforment aux relations de dispersion). On avait déjà remarqué que les méthodes de splittings n’ont pas encore convergées et donc n’ont pas le comportement souhaité en temps long. On observe 3 groupes de courbes :

en haut il y a la méthode de Suzuki (toute seule) qui, sans doute à cause des nombres FFT, a un comportement pire que Lie et Strang
~~au milieu il y a presque toutes les simulations~~
en bas, les simulations RK(3,3)-T(4) et RK(4,4)-P(2,2). Je n’ai pas d’explication du pourquoi. Ce sont deux simulations avec Maxwell dans la partie linéaire (ce qui augmente la conservation de quantités par la résolution exacte de nombreux termes).

~~Erreur relative des différentes simulations~~

Erratum : petite erreur dans l’ordre des données et leur label respectif dans la légende. Ce ne sont pas RK(3,3)-T(4) et RK(4,4)-P(2,2) qui sont en dessous, mais Lie et Strang ! (Suzuki reste bien au dessus) Toutes les méthodes de Lawson, avec ou sans Taylor ou Padé, avec ou sans pas de temps adaptatif, donnent la même erreur relative ou presque.

2021-07-29

Valeurs propres de Taylor

On souhaite calculer une approximation de l’exponentielle de la partie linéaire $A$ de notre problème, où $A$ est définie par : $A = \begin{pmatrix} 0 & -1 & 0 & 0 & \Omega_{pe}^2 & 0 \\ 1 & 0 & 0 & 0 & 0 & \Omega_{pe}^2 \\ 0 & 0 & 0 & 0 & 0 & i\kappa \\ 0 & 0 & 0 & 0 & -i\kappa & 0 \\ -1 & 0 & 0 & -i\kappa & 0 & 0 \\ 0 & -1 & i\kappa & 0 & 0 & 0 \\ \end{pmatrix}$ avec les modes de Fourier $\kappa\in[\![-N_z,N_z]\!]$ et $N_z=15$ , et $\Omega_{pe}=2$ . On sait déjà que $A$ étant semblable à une matrice anti-hermitienne, ses valeurs propres sont imaginaires pures.

On approche $e^{M}$ par une troncature de la série de Taylor à l’ordre $m$ : $e^{M}\approx T_m(M) = \sum_{k=0}^{m}\frac{M^k}{k!}.$

Les valeurs propres de $A$ étant imaginaire pures, on sait que les valeurs propres de $e^{tA}$ sont sur le cercle unité ( $\forall t\in\mathbb{R}$ ).

Maintenant que tout est bien défini, regardons les valeurs propres de la troncature de la série de Taylor $T_m(A)$ , avec $m=1,\dots,6$ .

Pour $T_1(A)=I+A$ , rien d’anormal puisque les valeurs propres de $A$ sont imaginaires pures et que $A$ est diagonalisable, on a : $I+A = P(I+D)P^{-1}$ où $P$ est la matrice de diagonalisation de $A$ . On trouve bien que les valeurs propres de $I+A$ sont simplement translatées de 1.

Je n’y connais rien en géométrie et autres courbes quadratiques, mais ça me me semble pas étonnant de retrouver les valeurs propres $T_2(A)$ sur une parabole. On a aucune valeur propre dans le cercle unité pour Taylor 1 et 2.

Pour $T_3(A)$ , on remarque que l’on a quelques valeurs propres (pour les modes les plus bas) dans le cercle unité, ce qui laisse présager d’une certaine stabilité du schéma de Lawson couplé à $T_3$ , ce qui… n’est pas ce que je retrouvais dans mes résultats numériques, où l’on remarque que le schéma de Lawson LRK(3,3) couplé à $T_1$ , $T_2$ et $T_3$ est instable. Après il y a assez peu de modes sur ou dans le cercle (en particulier le mode 2 (orange) qui risque d’être excité relativement tôt dans la simulation).

Énergie électrique, magnétique et cinétique des particules froides dans une simulation de Vlasov-Maxwell hybride linéarisé

Pour $T_4$ on a toutes les valeurs propres des bas modes dans le cercle unité, cela semble justifier de sa stabilité numérique observée dans l’exemple précédent.

Un peu la même chose pour $T_5$ , qui tourne sur de plus grosses simus (LRK(4,4) et avec un pas de temps plus grand que celui pris dans l’exemple précédent).

Et là, je pense que cela pourrait être intéressant d’ajouter une simu avec $T_6$ car on remarque que les valeurs propres ne sont pas dans le cercle unité (elles en sont pas loin, comme pour $T_2$ et $T_3$ ), il semble donc possible d’avoir un schéma instable avec une méthode d’ordre plus élevé. Cela donnerait sans doute un dernier coup fatal contre l’utilisation de Taylor, actuellement on ne comprend pas pourquoi $T_4$ et $T_5$ sont stables (on allume suffisamment peu de modes, et pour ces quelques modes on est stable ?), mais illustrer que c’est instable avec $T_6$ alors que l’on a une erreur de troncature plus faible me semble intéressant.

Une autre solution pour faire exploser toutes les simu avec Taylor est de rafiner un peu en $z$ et utiliser une condition initiale excitée par un bain de modes. Sachant qu’on a pas besoin de capturer la pente, cela peut se faire avec un maillage grossier en vitesse.

2021-11-09

Lorsqu’il est question des points de Leja

Les points de Leja sont des points sélectionnés pour être pas trop mal pour mailler un domaine (le conditionnement du maillage est pas mauvais). Ce maillage n’est pas régulier et propose plus de points sur les bords du maillage qu’au milieu et on représente ici pour 11 points.

les 11 points de Leja sur l’intervalle [-100,100]. — les 11 points de Leja sur l’intervalle $[-100,100]$ .

Le nième point de Leja d’un intervalle $K$ est défini par : $z_n = \arg\max_{z\in K}\left( V(z) = \prod_{i=0}^{n-1} |z-z_i| \right)$ et $z_0$ est pris tel que $|z_0|=\max_{z\in K}|z|$ . Ceci fonctionne donc sur le plan complexe ou un axe. Dans notre cas nous souhaitons principalement approcher la fonction exponentielle sur l’axe imaginaire, plus particulièrement en $c_iv_j\kappa\Delta t$ , avec les coefficients de Butcher $c_i\in[0,1]$ (généralement), la vitesse $v_j\in[-2,2]\approx\mathbb{R}$ , les modes de Fourier $\kappa\in [\![-13,13]\!]$ et le pas de temps $\Delta t\in[0,4]$ (d’après mes simus à pas de temps adaptatif LDP4(3)-P(2,2)), soit l’évaluation de la fonction exponentielle sur l’axe imaginaire entre $i[-100,100]$ environ. On construit donc le polynôme de Lagrange (de degré 10, car on a 11 points), et on va l’évaluer sur l’intervalle d’interpolation $[-100,100]$ et regarder sa partie réelle (qui devrait être un cosinus) sa partie imaginaire (un sinus), son module (qui devrait être égal à 1), ainsi que placer ces points (qui devraient être sur le cercle unité.

Partie réelle de L_{10}(iy) pour y\in[-100,100] — Partie réelle de $L_{10}(iy)$ pour $y\in[-100,100]$

Partie imaginaire de L_{10}(iy) pour y\in[-100,100] — Partie imaginaire de $L_{10}(iy)$ pour $y\in[-100,100]$

Module de L_{10}(iy) pour y\in[-100,100] — Module de $L_{10}(iy)$ pour $y\in[-100,100]$

Points sur le plan complexe L_{10}(iy) pour y\in[-100,100] — Points sur le plan complexe $L_{10}(iy)$ pour $y\in[-100,100]$

Tests avec 11 points de Leja

Et on effectue la même chose pour faire du Lagrange d’ordre 5 (en théorie je devrais faire du Lagrange d’ordre 4 car je dois comparer ça avec du Padé $(2,2)$ .)

Partie réelle de L_{5}(iy) pour y\in[-100,100] — Partie réelle de $L_{5}(iy)$ pour $y\in[-100,100]$

Partie imaginaire de L_{5}(iy) pour y\in[-100,100] — Partie imaginaire de $L_{5}(iy)$ pour $y\in[-100,100]$

Module de L_{5}(iy) pour y\in[-100,100] — Module de $L_{5}(iy)$ pour $y\in[-100,100]$

Points sur le plan complexe L_{5}(iy) pour y\in[-100,100] — Points sur le plan complexe $L_{5}(iy)$ pour $y\in[-100,100]$

Tests avec 6 points de Leja

Forcément il y a quelque chose qui s’appelle le théorème de Shannon et le signal est ici très largement sous échantillonné