Journal de bord de thèse

Josselin

2019-10-07

Relation de dispersion pour Vlasov-Poisson

Vlasov-Poisson : {tf+vxf+Evf=0xE=Rfdv1 \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)=feq(v)+εg(t,x,v)f(t,x,v) = f_\text{eq}(v) + \varepsilon g(t,x,v) et E(t,x)=εe(t,x)E(t,x) = \varepsilon e(t,x): {tg+vxg+edvfeq=0xe=Rgdv \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) : tiωxik \partial_t\cdot \rightarrow -i\omega\cdot \qquad \partial_x\cdot \rightarrow ik\cdot d’où {ik(vωk)g=edvfeqike=gdv \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=iedvfeqk(vωk) g = \frac{ie\mathrm{d}_vf_\text{eq}}{k\left(v-\frac{\omega}{k}\right)} La seconde donne alors : e=ikiedvfeqk(vωk)dv e = -\frac{i}{k}\int \frac{ie\mathrm{d}_vf_\text{eq}}{k\left(v-\frac{\omega}{k}\right)}\,\mathrm{d}v soit : e(11k2dvfeqvωkdv)=0 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 : D(ω,k)=11k2Rdvfeqvωkdv \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 feq(v)=Mρ,u,T(v)=ρ2πTevu22Tf_\text{eq}(v)=\mathcal{M}_{\rho,u,T}(v) = \frac{\rho}{\sqrt{2\pi T}}e^{-\frac{|v-u|^2}{2T}} : RdvMρ,u,T(v)vβdv=(vu)TMρ,u,T(v)vβdv=1T(vβ+βu)Mρ,u,T(v)vβdv=ρTβuTMρ,u,T(v)vβdv \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 : Mρ,u,T(v)vβdv\int\frac{\mathcal{M}_{\rho,u,T}(v)}{v-\beta}\,\mathrm{d}v : Mρ,u,T(v)vβdv=ρ2πTevu22Tvβdv=ρ2πTe(vu2T)2vu+uβdv \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=vu2Tz=\frac{v-u}{\sqrt{2T}}, dz=dv2T\mathrm{d}z = \frac{\mathrm{d}v}{\sqrt{2T}} : Mρ,u,T(v)vβdv=ρ2πTez22Tz+uβ2Tdz=ρ2πTez2z+uβ2Tdz \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 Z:ξ1πez2zξdz\mathcal{Z}:\xi\mapsto\frac{1}{\sqrt{\pi}}\int\frac{e^{-z^2}}{z-\xi}\,\mathrm{d}z : Mρ,u,T(v)vβdv=ρ2TZ(βu2T) \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ù : dvMρ,u,T(v)vβdv=ρTρ(βu)T2TZ(βu2T) \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 : RdvMρ,u,T(v)vβdv=ρT(1+βu2TZ(βu2T)) \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 feq(v)=δu(v)f_\text{eq}(v)= \delta_u(v) :

Rdvδu(v)vβdv=δu(v)(1vβ)dv=δu(v)1(vβ)2dv \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 : Rdvδu(v)vβdv=1(uβ)2 \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 : feq=M1α,0,Tc+Mα/2,u,1+Mα/2,u,1 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 : D(ω,k)=11k2[1αTc(1+ωk2TcZ(ωk2Tc))α2(1+ωku2Z(ωku2))α2(1+ωk+u2Z(ωk+u2))] \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 :

feq=(1α)δ0+Mα/2,u,1+Mα/2,u,1 f_\text{eq} = (1-\alpha)\delta_0 + \mathcal{M}_{^\alpha/_2,u,1} + \mathcal{M}_{^\alpha/_2,-u,1} Ce qui correspond à Tc0T_c \to 0. La relation de dispersion s’écrit alors : D(ω,k)=11k2[(1α)(kω)2α2(1+ωku2Z(ωku2))α2(1+ωk+u2Z(ωk+u2))] \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+zZ(z)z+2z2 1+z\mathcal{Z}(z) \overset{z\to+\infty}{\approx} \frac{2}{z^2} avec zCz\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) : {tuc=EtE=ρcucvfhdvtf^h=ikvf^hEvfh^ \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 fh=f^hf_h = \hat{f}_h. Ce problème se réécrit : t(ucEfh)U=(010ρc0000ikv)A(ucEfh)+(0vfhdvEvfh^)N(U) \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)} tU=AU+N(U) \partial_t U = AU + N(U) Pour introduire une méthode de Lawson, on utilise le changement de variable suivant V=etAUV=e^{-tA}U, l’exponentielle de la matrice AA peut se calculer et vaut : etA=(cos(ρct)sin(ρct)ρc0ρcsin(ρct)cos(ρct)000eikvt) 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 tV\partial_t V : tV=AetAU+etAtU \partial_t V = -Ae^{-tA}U + e^{-tA}\partial_t U soit : tV=AetAU+etAAU+etAN(U) \partial_t V = -Ae^{-tA}U + e^{-tA}AU + e^{-tA}N(U)

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

d’où : tV=etAN(U)=etAN(etAV)=N~(t,V) \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(ss,nn), ce qui permettra d’exprimer le schéma de Lawson associé. Pour l’exemple, nous utiliserons la méthode RK(33,33) 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 u˙=L(t,u)\dot{u} = L(t,u) : u(1)=un+ΔtL(tn,un)u(2)=34un+14u(1)+14ΔtL(tn+Δt,u(1))un+1=13un+23u(2)+23ΔtL(tn+Δt2,u(2)) \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 tV=etAN(etAV)\partial_t V = e^{-tA}N(e^{tA}V) : V(1)=Vn+ΔtetnAN(etnAVn)V(2)=34Vn+14V(1)+14Δte(tn+Δt)AN(e(tn+Δt)AV(1))Vn+1=13Vn+23V(2)+23Δte(tn+Δt2)AN(e(tn+Δt2)AV(2)) \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)N(U) par : N(U)(00000000λ)=Λ N(U) \approx \begin{pmatrix}0&0&0\\0&0&0\\0&0&\lambda\\\end{pmatrix} = \Lambda Λ\Lambda commute avec etAe^{tA}, donc la méthode RK(33,33) dit de Shu-Osher se réécrit : V(1)=Vn+ΔtΛVn=(I+ΔtΛ)VnV(2)=34Vn+14V(1)+14ΔtΛV(1)=(34I+14(I+ΔtΛ)2)VnVn+1=13Vn+23V(2)+23ΔtΛV(2)=pRK(3,3)(ΔtΛ)Vn \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 pRK(3,3)p_{RK(3,3)} le polynôme défini par : pRK(3,3)(z)=1+z+z22+z36 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 UU : e(tn+Δt)AUn+1=pRK(3,3)(ΔtΛ)etnAUn e^{(t^n+\Delta t)A}U^{n+1} = p_{RK(3,3)}(\Delta t\Lambda)e^{t^nA}U^n soit : Un+1=e(tn+Δt)ApRK(3,3)(ΔtΛ)etnAUn 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 etAe^{tA} et pRK(3,3)(ΔtΛ)p_{RK(3,3)}(\Delta t\Lambda) commutent, on sait déjà que II et etAe^{tA} commutent, Λ\Lambda et etAe^{tA} aussi, calculons Λ2\Lambda^2 : Λ2=(00000000λ2) \Lambda^2 = \begin{pmatrix}0&0&0\\0&0&0\\0&0&\lambda^2\\\end{pmatrix} donc Λ2\Lambda^2 et etAe^{tA} commutent, d’où : Un+1=pRK(3,3)(ΔtΛ)eΔtAUn U^{n+1} = p_{RK(3,3)}(\Delta t\Lambda)e^{\Delta tA}U^n On peut calculer la norme de etAe^{tA} avec notre matrice AA et on obtient etA2=1\left|\left|e^{tA}\right|\right|_2 = 1, t\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(ss,nn) sous-jacente. Le calcul a été effectué pour une méthode RK(33,33), mais se généralise facilement puisque Λ\Lambda commute avec etAe^{tA}, tR\forall t\in\mathbb{R}. De façon plus générale nous obtenons donc le résultat : pLRK(s,n)(z)=pRK(s,n)(z)eΔtA p_{LRK(s,n)}(z) = p_{RK(s,n)}(z)e^{\Delta tA} pLRK(s,n)p_{LRK(s,n)} est la fonction de stabilité de la méthode de Lawson associée à la méthode Runge-Kutta RK(ss,nn) de fonction de stabilité pRK(s,n)p_{RK(s,n)}.

Intéressons nous maintenant à la linéarisation de N(U)N(U) suivante : N(U)(00000ν00λ)=Λ˙ 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 etAe^{tA} ne commutent pas, on a pas l’égalité de fonction de stabilité : pLRK(s,n)(ΔtΛ˙)pRK(s,n)(ΔtΛ˙)eΔtA 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 : JcJh 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 vv des particules chaudes fhf_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)

{tuc=EtE=ρcucvfhdvtf^h+ikvf^h+Evfh^=0 \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 : ikE^(t,k)=f^hdvik\hat{E}(t,k) = \int\hat{f}_h\,\mathrm{d}v, est-elle vérifiée ?

t(ucEf^h)+(010ρc0000ikv)(ucEf^h)+(0fhdvEvfh^)=0 \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 : tU+AU+N(U)=0 \partial_t U + AU + N(U) = 0 une forme propice à l’utilisation des méthodes de Lawson.

Splitting en 3 étapes

  1. φΔt[a]\varphi^{[a]}_{\Delta t} : {tfh+vxfh=0tuc=0tE=vfhdv\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}
  2. φΔt[b]\varphi^{[b]}_{\Delta t} : {tfh+Evfh=0tuc=EtE=0\begin{cases}\partial_t f_h + E\partial_v f_h = 0 \\ \partial_t u_c = E \\ \partial_t E = 0 \end{cases}
  3. φΔt[c]\varphi^{[c]}_{\Delta t} : {tfh=0tuc=0tE=ρcuc\begin{cases} \partial_t f_h = 0 \\ \partial_t u_c = 0 \\ \partial_t E = -\rho_cu_c \end{cases}

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

avec dans l’ordre :

  1. φΔt[c](Un)=(fh(1)uc(1)E(1))=(fhnucnEnΔtρcucn)=U(1) \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)}
  2. φΔt[b](U(1))=(fh(2)uc(2)E(2))=(fh(1)(x,vΔtE(1))uc(1)+ΔtE(1)E(1))=U(2) \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)}
  3. φΔt[a](U(2))=(f^hn+1ucn+1E^n+1)=(eikvΔtf^h(2)uc(2)E^(2)ikf^hn+1dv+ikf^h(2)dv)=Un+1 \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 E^n+1\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 eciΔtAe^{c_i\Delta t A} : eciΔtA=(cos(ciΔt)sin(ciΔt)0sin(ciΔt)cos(ciΔt)000eiciΔtkv) 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 cic_i est un réel quelconque, mais pas de problème lorsque ci>0c_i>0 ou ci<0c_i<0, de même que le cas ci=0c_i=0 qui ne se rencontre jamais. Si cic_i est quelconque il travaille sur ci|c_i| et fait intervenir des cici\frac{|c_i|}{c_i} pour jouer sur la parité des fonctions sinus et cosinus.

Une fonction calculant l’exponentielle eciΔtAe^{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 :

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 : ut+ux=0 u_t + u_x = 0 avec Tf=0.25T_f=0.25, Δx=1100\Delta x=\frac{1}{100}, x[0,1]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 Δt\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 uxu_x par Δt\Delta t. On obtient alors des résultats plus cohérent (correction effectuer jeudi à 14h).

Comparaison des différentes méthodes WENO

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 : ut+23ux=0 u_t + \frac{2}{3}u_x = 0

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

Comparaison des différentes méthodes WENO

Dans ces conditions, WENO-SL fait aussi bien que WENO-Z. Un test similaire a été effectué avec une vitesse négative 23-\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(Δt)\log(\Delta t). Mon erreur dans le code provient (vraisemblablement) du fait que je choisis des Δt\Delta t qui dépendent de Δx\Delta x, ce qui implique un nombre d’itération TfΔt\frac{T_f}{\Delta t} non entier. Je dois donc faire en sorte, pour les 2 mesures d’ordre, d’avoir TfΔtN\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 uu. 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:

Avancement sur les petits trucs à vérifier

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

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 : tU+AU+N(U) \partial_t U + AU + N(U) en linéarisant N(U)N(U) avec une matrice Λ˙\dot{\Lambda} définie par : Λ˙=(00000μ00λ) \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 : Un+1=[10Δtμ(Δt2λ2eiΔtkvsin(Δt2+tn)+ΔtλeiΔtkv2sin(Δt+tn)+2ΔtλeiΔtkvsin(Δt2+tn)+e3iΔtkv2sin(tn)+eiΔtkv2sin(Δt+tn)+4eiΔtkvsin(Δt2+tn))eikv(3Δt2+tn)601Δtμ(Δt2λ2eiΔtkvcos(Δt2+tn)+ΔtλeiΔtkv2cos(Δt+tn)+2ΔtλeiΔtkvcos(Δt2+tn)+e3iΔtkv2cos(tn)+eiΔtkv2cos(Δt+tn)+4eiΔtkvcos(Δt2+tn))eikv(3Δt2+tn)600Δt3λ36+Δt2λ22+Δtλ+1]Un 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 ΔtΛ˙\Delta t\dot{\Lambda} : Δt3eΔtA2tnAΛ˙eΔtA2Λ˙eΔtAΛ˙etnA6+Δt2eΔtAtnAΛ˙eΔtAΛ˙etnA6+Δt2eΔtA2tnAΛ˙eΔtA2Λ˙eΔtA+tnA6+Δt2eΔtA2tnAΛ˙eΔtA2Λ˙etnA6+ΔtetnAΛ˙etnA6+ΔteΔtAtnAΛ˙eΔtA+tnA6+2ΔteΔtA2tnAΛ˙eΔtA2+tnA3+1 \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 Δtλ\Delta t\lambda : Δt3λ2μ\Delta t^3 \lambda^2\mu, Δt2λμ\Delta t^2 \lambda\mu, Δtμ\Delta t\mu. Et pour μ=0\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 etAe^{tA}, et je connais les instants tt où cette fonction sera évaluée dans le schéma (en prennant les notations d’un tableau de Butcher, c’est en e±(tn+ciΔt)Ae^{\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 : u(1)=un+12ΔtL(tn,un)u(2)=3un2u(1)+2ΔtL(tn+12Δt,u(1))un+1=13un+u(1)+13u(2)+16ΔtL(tn+Δt,u(2)) \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 etAe^{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 : tU=AU+N(U) \partial_t U = AU + N(U) avec U=(uvEf^h)TU = \begin{pmatrix}u_v & E & \hat{f}_h\end{pmatrix}^\textsf{T}, et : A=(010ρc0000ikv)N(U)=(0vfhdvEvfh^) 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 : U(1)=eΔtAUn+ΔteΔtAN(Un)U(2)=34eΔt2AUn+14eΔt2AU(1)+14ΔteΔt2AN(U(1))Un+1=13eΔtLUn+23eΔt2AU(2)+23ΔteΔt2AN(U(2)) \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 : (uc(1)E(1)f^h(1))=(cos(ρcΔt)sin(ρcΔt)ρc0ρcsin(ρcΔt)cos(ρcΔt)000eikvΔt)(ucnEnf^hn)+Δt(cos(ρcΔt)sin(ρcΔt)ρc0ρcsin(ρcΔt)cos(ρcΔt)000eikvΔt)(0vfhndvEnvfhn^) \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 (ucni)i({u_c^n}_i)_i, (Ein)i(E^n_i)_i, et (f^hnι,k)ι,k({\hat{f}_h^n}_{\iota,k})_{\iota,k} avec ii indice d’espace xi=iΔxx_i = i\Delta x, kk indice de vitesse vk=kΔvv_k = k\Delta v et ι\iota indice de phase κι=2πιL\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 xi=iΔxx_i = i\Delta x et la phase κι=2πιL\kappa_\iota = \frac{2\pi\iota}{L}.

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

  1. Calcul de (Envfhn)j,k\left(E^n\partial_vf_h^n\right)_{j,k} avec WENO ou autre, et calcul de (vfhndv)j(\int vf_h^n\,\mathrm{d}v)_j (ceci nécessite d’effectuer une FFT inverse de f^hn\hat{f}_h^n, et il s’agit finalement le la seule étape où fhf_h est nécessaire).
  2. Calcul des (uc(1))j(u_c^{(1)})_j et (E(1))j(E^{(1)})_j : ucj(1)=ucjnρccos(ρcΔt)Ejnsin(ρcΔt)ρc+(vfhndv)jsin(ρcΔt)Ej(1)=ucjnρcsin(ρcΔt)+Ejncos(ρcΔt)+(vfhndv)jcos(ρcΔt),j,xj=jΔx \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
  3. Boucle pour tout vkv_k : effectuer une FFT pour obtenir (Envfhn^)ι,k\left(\widehat{E^n\partial_vf_h^n}\right)_{\iota,k} à kk fixé, et calcul des (f^h(1))ι,k(\hat{f}_h^{(1)})_{\iota,k} : f^hι,k(1)=eiκιvkΔtf^hι,kneiκιvkΔt(Envfhn^)ι,k,ι,κι=2πιL {\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(Un)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 vkv_k il n’est pas possible d’effectuer le calcul des (uc(1)i,Ei(1))({u_c^{(1)}}_i,E^{(1)}_i) en même temps que f^h(1)ι,k{\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 (X˙\dot{X}) indiquent des variables temporaires ne servant qu’à un étage, ou dans une boucle.

  1. étage 1 :
    • (f˙)j,kIFFT(f^hnι,k)(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^n}_{\iota,k})
    • (J˙)jkvkf˙j,kΔv(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v
    • (Eδvf˙)j,kWENO((Envf˙)j,k)(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^n\partial_v\dot{f})_{j,k})
    • pour tout xjx_j :
      • uc(1)jucnjcos(ρcΔt)+Ejnsin(ρcΔt)ρcΔtJ˙jsin(Δt){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)
      • Ej(1)ucnjρcsin(ρcΔt)+Ejncos(ρcΔt)ΔtJ˙jcos(Δ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 vjv_j :
      • (δ˙)ιFFT(Eδvf˙j,k)(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})
      • pour tout κι\kappa_\iota :
        • f^h(1)ι,kf^hnι,keiκιvkΔtΔtδ˙ιeiκιvkΔt{\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}
  2. étage 2:
    • (f˙)j,kIFFT(f^h(1)ι,k)(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^{(1)}}_{\iota,k})
    • (J˙)jkvkf˙j,kΔv(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v
    • (Eδvf˙)j,kWENO((E(1)vf˙)j,k)(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^{(1)}\partial_v\dot{f})_{j,k})
    • pour tout xjx_j :
      • uc(1)j34ucnjcos(ρcΔt2)+34Ejnsin(Δt2)ρc+14uc(1)jcos(ρcΔt2)14Ej(1)sin(ρcΔt2)ρc+14ΔtJ˙jsin(ρcΔt2){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})
      • Ej(1)34ucnjρcsin(ρcΔt2)+34Ejncos(ρcΔt2)+14uc(1)jρcsin(ρcΔt2)+14Ej(1)cos(ρcΔt2)14ΔtJ˙jcos(ρcΔt2)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 vjv_j :
      • (δ˙)ιFFT(Eδvf˙j,k)(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})
      • pour tout κι\kappa_\iota :
        • f^h(1)ι,k34f^hnι,keiκιvkΔt2+14f^h(1)ι,keiκιvkΔt214Δtδ˙ιeiκιvkΔt2{\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}}
  3. étage 3:
    • (f˙)j,kIFFT(f^h(2)ι,k)(\dot{f})_{j,k} \gets \text{IFFT}({\hat{f}_h^{(2)}}_{\iota,k})
    • (J˙)jkvkf˙j,kΔv(\dot{J})_j \gets \sum_k v_k{\dot{f}}_{j,k}\Delta v
    • (Eδvf˙)j,kWENO((E(2)vf˙)j,k)(\dot{E\delta_vf})_{j,k} \gets \text{WENO}((E^{(2)}\partial_v\dot{f})_{j,k})
    • pour tout xjx_j :
      • ucn+1j13ucnjcos(ρcΔt)+13Ejnsin(ρcΔt)ρc+23uc(2)jcos(ρcΔt2)+23Ej(2)sin(ρcΔt2)ρc23ΔtJ˙jsin(ρcΔt2){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})
      • Ejn+113ucnjρcsin(Δt)+13Ejncos(ρcΔt)23uc(2)jρcsin(ρcΔt2)+23Ej(2)cos(ρcΔt2)23ΔtJ˙jcos(ρcΔt2)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 vjv_j :
      • (δ˙)ιFFT(Eδvf˙j,k)(\dot{\delta})_\iota \gets \text{FFT}(\dot{E\delta_vf}_{j,k})
      • pour tout κι\kappa_\iota :
        • f^hn+1ι,k13f^hnι,keiκιvkΔt+23f^h(2)ι,keiκιvkΔt223Δtδ˙ιeiκιvkΔt2{\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 :

tU=AU+N(U) \partial_t U = AU + N(U) avec : U=(ucEf^h)A=(01010000eikv)N(U)=(0vfhdvEvfh^) 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α)fc+αfh×(1+ϵcos(kx))f = (1-\alpha)f_c + \alpha f_h\times(1 + \epsilon\cos(kx)). Pour fhf_h nous prenons : fh=M[α/2,u,1]+M[α/2,u,1] f_h = \mathcal{M}_{[\alpha/2,u,1]} + \mathcal{M}_{[\alpha/2,-u,1]} il est également possible de prendre : fh=M[1,0,1]v10S f_h = \mathcal{M}_{[1,0,1]}\frac{v^{10}}{\mathcal{S}} S\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 ρc\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=(010ρc0000eikv)N(U)=(0vfhdvEvfh^) 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 ρc=1α\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 :

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 :

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 : ut+(au)x=0 u_t + (au)_x = 0 On commence par une intégration sur [tn,tn+1][t^n,t^{n+1}] : u(tn+1,x)=u(tn,x)([tn,tn+1]au(τ,x)dτ)x 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 : uin+1=uinFxx=x1 u_i^{n+1} = u_i^n - \mathcal{F}_{x|x=x_1} F\mathcal{F} est défini par : F(x)=[tn,tn+1]au(τ,x)dτ=1Δx[xΔx2,x+Δx2]H(ξ)dξ \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 à xx de F\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) : Fx=1Δx(H(x+Δx2)H(xΔx2)) \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 : uin+1=uin1Δx(H(xi+12)H(xi12)) 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 Ωi\Omega_i défini par les trois points (xi,tn+1)(x_i,t^{n+1}), (xi,tn)(x_i,t^n) et (xi,tn)(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

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

Ωiut+(au)x=0 \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 (xi,tn)(x_i^\star,t^n),(xi,tn+1)(x_i,t^{n+1})) : Ωiut+(au)x=xixiu(tn,ξ)dξ+tntn+1au(τ,xi)dτ \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 : xixiu(tn,ξ)dξ=tntn+1au(τ,xi)dτ=F(xi) \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) F(xi)=Hˉi\mathcal{F}(x_i) = \bar{\mathcal{H}}_i (valeur moyenne de H\mathcal{H} sur la cellule [xi12,xi+12][x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]). Donc on utilise le pied exacte de la caractéristique xix_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 {Hˉi}i=1N\{\bar{\mathcal{H}}_i\}_{i=1}^N à partir des {uin}i=1N\{u_i^n\}_{i=1}^N, plus précisemment les points nécessaires sont : Hˉi=R[xi,xi](uip1n,,ui+q1n) \bar{\mathcal{H}}_i = \mathcal{R}[x_i^\star,x_i](u_{i-p_1}^n,\dots,u_{i+q_1}^n) R[a,b]\mathcal{R}[a,b] est la reconstruction de [a,b]u(t,ξ)dξ\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 xix_i^\star et xix_i, aussi loin soit xix_i^\star. Dans la pratique, lors de l’écriture du schéma, seul le stencil autour de xix_i^\star est nécessaire, mais pour conserver la formulation volumes-finis, il est nécessaire de considérer la reconstruction jusqu’à xix_i.

La suite de l’articule est tout un jeu de construction de l’opérateur de reconstruction R[a,b]\mathcal{R}[a,b] en utilisant une méthode WENO, dans un premier temps avec xix_i^\star dans la cellule [xi1,xi][x_{i-1},x_i] (volumes finis classiques), puis en dehors de cette cellule (avec la prise en compte de ξ=xixiΔx\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ù xix_i^\star est dans la cellule [xi1,xi][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 : ut+ux=0 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 :

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 δ0(v)\delta_0(v) par une gaussienne. Pour cela on simule à tt fixé (5\sim 5 ou 1010). La gaussienne doit être au moins représentée par une vingtaine de points pour être correctement représentée, donc : ΔvTc20 \Delta v \sim \frac{\sqrt{T_c}}{20} avec Tc=104,103,102,101T_c = 10^{-4},\,10^{-3},\,10^{-2},\,10^{-1}. Donc pour s’assurer d’avoir toujours le même Δv\Delta v (et donc le même Δt\Delta t), on prend Δv10420=5.0104\Delta v \sim \frac{\sqrt{10^{-4}}}{20} = 5.0\cdot 10^{-4}. On prend le Δt\Delta t qui CFLise ce Δv\Delta v : Δt=σΔvEn \Delta t = \sigma\frac{\Delta v}{||E^n||_\infty} Pour le moment le code VHL avec Lawson utilise un RK(33,33), donc σ=1.433\sigma = 1.433 et les premières simus donnent En0.12||E^n||_\infty \leq 0.12, on obtient alors pour Δv\Delta v, NvN_v et Δt\Delta t : Δv5.0104Nv=16Δv32000Δt5.9103 \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 NvN_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.
TcT_c ΔvTc20\Delta v \sim \frac{\sqrt{T_c}}{20} Nv16ΔvN_v \sim \frac{16}{\Delta v} Δt1.433Δv0.12\Delta t \sim 1.433\frac{\Delta v}{0.12}
10410^{-4} 5.001045.00\cdot 10^{-4} 3200032\,000 5.971035.97\cdot 10^{-3}
10310^{-3} 1.581031.58\cdot 10^{-3} 1011910\,119 1.891021.89\cdot 10^{-2}
10210^{-2} 5.001035.00\cdot 10^{-3} 32003\,200 5.971025.97\cdot 10^{-2}
10110^{-1} 1.581021.58\cdot 10^{-2} 10121\,012 1.891011.89\cdot 10^{-1}

Il est intéressant de noter que les modélisations hybrides supposent Tc0T_c\to 0 mais n’imposent pas de raffinement de maillage. Il peut être intéressant de comparer les résultats avec Tc=104T_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 :

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

2019-12-19

À propos des relations de dispersion

Considérons le modèle de Vlasov-Poisson : {tf+vxf+Evf=0xE=Rfdv1 \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)=feq(v)+εg(t,x,v)f(t,x,v) = f_\text{eq}(v) + \varepsilon g(t,x,v) : {εtg+vεxg+Edvfeq+εEvg=0xE=Rfeqdv+εRgdv1 \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 11, par conséquent, l’équation de Poisson implique que EE est de taille ε\varepsilon, on notera alors E(t,x)=εe(t,x)E(t,x) = \varepsilon e(t,x) : {εtg+vεxg+eεdvfeq+ε2evg=0εxe=εRgdv \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 : {tg+vxg+edvfeq=0xe=Rgdv \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 d\mathrm{d} droit la dérivée en vv de feqf_\text{eq} car feqf_\text{eq} ne dépend que de vv, 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 : xiktt=0iω \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(ω,k)D(\omega,k).

{ik(vωk)g=edvfeqgt=0ike=gdv \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 gg : g=iedvfeqk(vωk)+igt=0k(vωk) 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=ikiedvfeqk(vωk)dvikigt=0k(vωk)dv 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=1k2gt=0vωkdv11k2dvfeqvωkdv 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=N(ω,k)D(ω,k)e = \frac{N(\omega,k)}{D(\omega,k)} avec : N(ω,k)=1k2gt=0vωkdvD(ω,k)=11k2dvfeqvωkdv 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(ω,k)D(\omega,k).

Calcul d’intégrales

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

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

Cas d’une distribution maxwellienne

RdvMρ,u,T(v)vβdv=R(vu)TMρ,u,T(v)vβdv=1TR(vβ+βu)Mρ,u,T(v)vβdv=ρTβuTRMρ,u,T(v)vβdv \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 : Mvβdv\int\frac{\mathcal{M}}{v-\beta}\,\mathrm{d}v : RMρ,u,T(v)vβdv=ρ2πTRevu22Tvβdv=ρ2πTRe(vu2T)2vu+uβdv \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=vu2Tz=\frac{v-u}{\sqrt{2T}}, dz=dv2T\mathrm{d}z=\frac{\mathrm{d}v}{\sqrt{2T}} : RMρ,u,T(v)vβdv=ρ2πTRez22Tz+uβ2Tdz=ρ2πTRez2z+uβ2Tdz \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 : Z:ζ1π+ez2zζdz \mathcal{Z}: \zeta \mapsto \frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty} \frac{e^{-z^2}}{z-\zeta}\,\mathrm{d}z D’où : RMρ,u,T(v)vβdv=ρ2TZ(βu2T) \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 : RdvMρ,u,T(v)vβdv=ρT(1+βu2TZ(βu2T)) \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

Rdvδu(v)vβdv=Rδu(v)dv(1vβ)dv=Rδu(v)1(vβ)2dv \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 : Rdvδu(v)vβdv=1(uβ)2 \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)=M(1α),0,Tc(v)+(Mα ⁣/ ⁣2,u,1(v)+Mα ⁣/ ⁣2,u,1(v))(1+ϵcos(kx)) 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, uu la vitesse caractéristique des particules chaudes, et TcT_c un paramètre pris arbitrairement bas correpondant à la température des particules froides (Tc1T_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 Tc0T_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α)δ0(v)+(Mα ⁣/ ⁣2,u,1(v)+Mα ⁣/ ⁣2,u,1(v))(1+ϵcos(kx)) 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 feqf_\text{eq} depuis ces conditions initiales dépendant de xx et vv il suffit de supprimer la perturbation en xx (c’est-à-dire prendre ϵ=0\epsilon=0). On remarque alors que feqf_\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 : DK(ω,k)=11k2[1αTc(1+ωk2TcZ(ωk2Tc))α2(1+ωku2Z(ωku2))α2(1+ωk+u2Z(ωk+u2))] \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é : DH(ω,k)=11k2[(1α)(kω)2α2(1+ωku2Z(ωku2))α2(1+ωk+u2Z(ωk+u2))] \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 : E(t)=v2fdvdx+E2dx=E(0) \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)=M1α,0,Tc(v)+(Mα/2,u,1(v)+Mα/2,u,1(v))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 : EK(t)=[(1α)Tc+αu2+α]L \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 xx cela donnera 00.

Dans le cas hybride linéarisé, f(t=0,x,v)=(1α)δ0(v)+(Mα/2,u,1(v)+Mα/2,u,1(v))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 : EHL=[αu2+α]L \mathcal{E}_{HL} = \left[ \alpha u^2 + \alpha \right]L

D’où : EK(t)EHL(t)=(1α)LTc \mathcal{E}_K(t) - \mathcal{E}_{HL}(t) = (1-\alpha)LT_c ce qui donne en version discrète : EKn+O(Δtp)EHLn+O(Δtq)=(1α)LTc \mathcal{E}^n_K + \mathcal{O}(\Delta t^p) - \mathcal{E}^n_{HL} + \mathcal{O}(\Delta t^q) = (1-\alpha)LT_c avec pp et qq respectivement les ordres en temps de codes de résolution du modèle cinétique et du modèle hybride linéarisé. Donc pour pp et qq suffisamment grand on ne trace que (1α)LTc(1-\alpha)LT_c. Il est intéressant de constater que l’erreur sur l’énergie est linéaire par rapport au produit ρcTc\rho_c T_c, avec ρc\rho_c la densité des particules froides dans la condition initiale (1α)(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 fcf_c et les partiules chaudes fhf_h : {tfc+vxfc+Evfc=0tfh+vxfh+Evfh=0tE=v(fh+fc)dv \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 fcf_c pour l’écrire comme un modèle fluide à partir de grandeurs intégrales, on calcule donc les moments de la première équation : Rt(1v)fcdv+Rvx(1v)fcdv+REv(1v)fcdv=0 \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 fc(t,x,v)=ρc(t,x)δuc(t,x)(v)f_c(t,x,v) = \rho_c(t,x)\delta_{u_c(t,x)}(v), ce qui nous permet d’obtenir : {tρc+x(ρcuc)=0t(ρcuc)+x(ρcuc2)ρcE=0 \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 : {tρc+x(ρcuc)=0t(ρcuc)+x(ρcuc2)ρcE=0tE=vfhdvρcucxE=fhdv+ρc1tfh+vxfh+Evfh=0 \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 EE à 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 : ρc(t,x)=ρc(0)(x)+ερc(1)(t,x)uc(t,x)=εuc(1)(t,x)E(t,x)=εE(1)(t,x) \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) : {εtρc(1)+εx(ρc(0)u(1))=0ερc(0)tuc(1)ερc(0)E(1)=0εtE(1)=vfhdvερc(0)uc(1)tfh+vxfh+E(1)vfh=0 \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 ρc(1)\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 (fhf_h étant supposé petit par rapport à fcf_c), on peut donc simplifier par ε\varepsilon pour obtenir le modèle de Vlasov hybride linéarisé : {tuc(1)=E(1)tE(1)=vfhdvρc(0)uc(1)tfh+vxfh+E(1)vfh=0 \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

Résultats

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

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

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

É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

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

Attention : pour bien observer une droite la température TcT_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

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 : tU=AU+F(U) \partial_t U = AU+F(U)

avec : U=[jc,xjc,yBxByExEyf^] 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=[0100Ωpe20010000Ωpe200000000000000010000000100000000000ikvz] ,F(U)=[00iEykiExkiByk+R vx(f^)iBxk+R vy(f^)(f^)] 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 :

U(1)=eΔtAUn+ΔteΔtAF(Un)U(2)=34eΔt2AUn+14eΔt2AU(1)+Δt4eΔt2AF(U1)Un+1=13eΔtAUn+23eΔt2AU(2)+23ΔteΔt2AF(U(2)) \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 FF 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ΔtAe^{\Delta tA} peut se faire de manière exacte (pour la valeur de Ωpe=2\Omega_{pe}=2) :

eΔtA=[0.378732187481834cos(2Δt9172)+0.621267812518167cos(2Δt17+92)0.378732187481834sin(2Δt9172)0.621267812518167sin(2Δt17+92)000.970142500145332sin(2Δt9172)+0.970142500145332sin(2Δt17+92)0.970142500145332cos(2Δt9172)+0.970142500145332cos(2Δt17+92)00.378732187481834sin(2Δt9172)+0.621267812518167sin(2Δt17+92)0.378732187481834cos(2Δt9172)+0.621267812518167cos(2Δt17+92)000.970142500145332cos(2Δt9172)0.970142500145332cos(2Δt17+92)0.970142500145332sin(2Δt9172)+0.970142500145332sin(2Δt17+92)0001.000000001.00000.242535625036333sin(2Δt9172)0.242535625036333sin(2Δt17+92)0.242535625036333cos(2Δt9172)0.242535625036333cos(2Δt17+92)000.621267812518167cos(2Δt9172)+0.378732187481834cos(2Δt17+92)0.621267812518166sin(2Δt9172)0.378732187481833sin(2Δt17+92)00.242535625036333cos(2Δt9172)+0.242535625036333cos(2Δt17+92)0.242535625036333sin(2Δt9172)0.242535625036333sin(2Δt17+92)000.621267812518167sin(2Δt9172)+0.378732187481834sin(2Δt17+92)0.621267812518166cos(2Δt9172)+0.378732187481833cos(2Δt17+92)0000000eiΔtkvz] 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 Ωpe\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 Rvxyfdv\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 ff, 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 AA laisse sous-entendre que jc,xj_{c,x} et jc,yj_{c,y} seront résolus dans l’espace réel, mais le brassage dû à eAe^{A} implique que ces variables sont des variables de Fourier. Puisque la transformée de Fourier est définie par :

f^(ξ,v)=Rf(z,v)ei2πξzdz \hat{f}(\xi,v) = \int_\mathbb{R} f(z,v)e^{-i2\pi\xi z}\,\mathrm{d}z on peut écrire que : Rvxyf(z,v)dv^=RRvxyf(z,v)dvei2πξzdz=RvxyRf(z,v)ei2πξzdzdv=Rvxyf^dv \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 vv (caché dans (f^)\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 jc,xyj_{c,x|y}, ExyE_{x|y} ou BxyB_{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.

  1. Calcul des variables j^c,x(1)\hat{j}_{c,x}^{(1)}, j^c,y(1)\hat{j}_{c,y}^{(1)}, E^x(1)\hat{E}_x^{(1)}, E^y(1)\hat{E}_y^{(1)}, B^x(1)\hat{B}_x^{(1)} et B^y(1)\hat{B}_y^{(1)} :
    • calcul des intégrales vxf^dv\int v_x \hat{f}\,\mathrm{d}v et vyf^dv\int v_y \hat{f}\,\mathrm{d}v :
      • (j^h,x)ikx,ky,kzvkxf^i,kx,ky,kzΔ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
      • (j^h,y)ikx,ky,kzvkyf^i,kx,ky,kzΔ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 xix_i (code donné par sympy) :
      • j^c,x(1)\hat{j}_{c,x}^{(1)}\gets \cdots
      • j^c,y(1)\hat{j}_{c,y}^{(1)}\gets \cdots
      • E^c,x(1)\hat{E}_{c,x}^{(1)}\gets \cdots
      • E^c,y(1)\hat{E}_{c,y}^{(1)}\gets \cdots
      • B^c,x(1)\hat{B}_{c,x}^{(1)}\gets \cdots
      • B^c,y(1)\hat{B}_{c,y}^{(1)}\gets \cdots
  2. Calcul de la variable f^(1)\hat{f}^{(1)} :
    • transformée inverse de f^n\hat{f}^{n} :
      • (f)i,kx,kz,kziFFT(f^,kx,kz,kzn)i\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 vxv_x, vyv_y et vzv_z)
    • calcul de l’approximation de (Ex+vyB0+vzBy)vxf+(EyvxB0+vzBx)vyf+(vxByvyBx)vzf(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 vkxv_{k_x} :
        • pour tout vkyv_{k_y} :
          • pour tout vkyv_{k_y} :
            • pour tout vkzv_{k_z} :
              • pour tout kik_i :
                • velocity_vxEx,i+vkyB0+vkzBy,i\texttt{velocity\_vx} \gets E_{x,i}+v_{k_y}B_0 + v_{k_z}B_{y,i}
                • velocity_vyEy,ivkxB0+vkzBx,i\texttt{velocity\_vy} \gets E_{y,i}-v_{k_x}B_0 + v_{k_z}B_{x,i}
                • velocity_vzvkxBy,ivkyBx,i\texttt{velocity\_vz} \gets v_{k_x}B_{y,i} - v_{k_y}B_{x,i}
                • vfi,kx,ky,kzWENO(velocity_vx,fi,kx3:kx+3,ky,kz)+WENO(velocity_vy,fi,kx,ky3:ky+3,kz)+WENO(velocity_vz,fi,kx,ky,kz3:kz+3)\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 f^(1)\hat{f}^{(1)} :
      • pour tout vkxv_{k_x} :
        • pour tout vkyv_{k_y} :
          • pour tout vkyv_{k_y} :
            • pour tout vkzv_{k_z} :
              • (vf^)iFFT(vf,kx,ky,kz)i\left(\widehat{\partial_vf}\right)_i\gets\text{FFT}(\partial_vf_{\cdot,k_x,k_y,k_z})_i
              • pour tout kik_i :
                • f^i,kx,ky,kz(1)f^i,kx,ky,kzn+Δtvf^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

WENO\text{WENO} est une fonction qui calcule une approximation de axu(xi)a\partial_xu(x_i) en prenant en argument aa et le stencil ui3:i+3u_{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 f^i+12±\hat{f}^\pm_{i+\frac{1}{2}}. C’est-à-dire WENO(a,ui3,ui2,ui1,ui,ui+1,ui+2,ui+3)\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+ui+12+ui12+Δx+aui+12ui12Δxa^+\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×v)vf(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 B0B_0 (composante en zz du vecteur B=(Bx,By,B0)TB = (B_x,B_y,B_0)^\textsf{T}).

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

avec :

v=(vxvy)E=(ExEy)J=(0110) 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,etJw)g(t,w) = f(t,e^{tJ}w), soit g(t,etJv)=f(t,v)g(t,e^{-tJ}v) = f(t,v), avec w=(w1,w2)w=(w_1,w_2). On obtient alors :

tg(t,w)=t(f(t,eetJw))=Evf(t,etJw) \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 : tg(t,w)=Evf(t,etJw) \partial_tg(t,w) = -E\cdot\nabla_vf(t,e^{tJ}w)

On veut maintenir construire notre équation que vérifie gg :

etJwg(t,w)=(vf)(t,etJw)EetJwg(t,w)=E(vf)(t,etJw)=tg(t,w) \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ù :

tg(t,w)+EetJwg(t,w)=0 \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)=f0(v)g(t=0,w) = f_0(v). Après avoir faire un schéma sur ggavec un schéma Runge-Kutta, om retrouve la solution sur ff par : f(t,v)=g(t,etJv)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)(x,y) (plutôt qu’un (vx,vy)(v_x,v_y) correspondant plus au cas 1dz-3dv).

Donc avec ces notations le toy-model devient :

{tu+(E+JX)Xf=0u(t=0,x,y)=u0(x,y) \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=(xy)E=(e1e2)J=(0110) 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 : tg(t,z)=EetJzg(t,z)\partial_t g(t,z) = E\cdot e^{tJ}\nabla_z g(t,z) avec z=(z1,z2)z=(z_1,z_2).

L’équation sur uu est résolu avec une méthode Runge-Kutta, un RK(3,3) suffira pour le test, et l’équation sur gg aussi. La CFL sur uu est du type : ΔtσΔxE+JX\Delta t \leq \sigma\frac{\Delta x}{\|E+JX\|_\infty}, ce qui est grosso-modo une CFL en ΔtσΔxE+xmax\Delta t \leq \sigma\frac{\Delta x}{\|E\|_\infty + x_\text{max}} (on considère que domaine symétrie en xx et yy avec (xy)min=(xy)max(x|y)_\text{min} = -(x|y)_\text{max}, la notation (xy)(x|y) indique indifféremment la variable xx ou yy). La CFL de l’équation filtrée est, quant à elle, en ΔtσΔxE\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 EE 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 :

ut+aux+buy=0 u_t + a u_x + b u_y = 0

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

tuyxu+xyu=0 \partial_t u - y\partial_x u + x\partial_y u = 0

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

{tuyxu=0{tu+xyu=0 \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 Δtmin(σΔxymax,σΔyxmax)\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 :

tu+J(xy)u=0 \partial_t u + J\begin{pmatrix}x\\y\end{pmatrix}\cdot\nabla u = 0

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

ΔtσΔxJ(xy)? \Delta t \leq \sigma \frac{\Delta x}{\left\|J\begin{pmatrix}x\\y\end{pmatrix}\right\|_\text{?}}

On a remarqué que prendre xmaxx_\text{max}, ce qui revient à prendre la norme J(xy)\left\|J\begin{pmatrix}x\\y\end{pmatrix}\right\|_\infty ne fonctionnait pas. Lorsque j’essaie de prendre la norme 2, donc 2xmax\sqrt{2}x_\text{max} ça ne fonctionne pas non plus, je semble obligé de prendre 2xmax2x_\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 : ut+ux=0u_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 aa et bb avec :

aij=e1yj,bij=e2+xi a_{ij} = e_1 - y_j ,\quad b_{ij} = e_2 + x_i

puis une autre version avec le filtrage :

a(t)=e1cos(t)e2sin(t),b(t)=e1sin(t)+e2cos(t) 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,etJv)f(t,v) = g(t,e^{-tJ}v) Nicolas m’a conseillé le splitting en 3 étapes de Joackim. Effectuer une rotation de tt revient à calculer la solution hh de l’équation th+Jvvh=0\partial_t h + Jv\cdot\partial_v h=0, h(0)=g(t)h(0)=g(t) (dans notre cas).

ϕ1={th+v2v1h=0ϕ2={thv1v2h=0ϕ3={th+v2v1h=0 \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)=ϕ3[tan(t/2)]ϕ2[sin(t)]ϕ1[tan(t/2)](g(t))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 :

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

ΔtσΔxE+Jv10.00637255\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 gg :

ΔtσΔxE10.325 \Delta t \leq \sigma\frac{\Delta x}{\left\|E\right\|_1} \approx 0.325

Dans le cadre de WENO-RK(3,3) σ1.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.325

La solution calculée sur ff dans le cas où Δt=0.325\Delta t = 0.325 ne donne que des nan, car on est complètement hors CFL, la solution sur gg 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 B0=0B_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 Δt=5103\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 B0=0B_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,vz)=f(t,z,etB0Jv,vz)g(t,z,w,v_z) = f(t,z,e^{-tB_0J}v,v_z) pour obtenir l’équation sur gg qui est :

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

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

tg=(vzExcos(B0t)+Eysin(B0t)+vzBxsin(B0t)vzBycos(B0t)Exsin(B0t)+Eycos(B0t)+vzBxcos(B0t)+vzBysin(B0t)Bx(w1sin(B0t)+w2cos(B0t))+By(w1cos(B0t)w2sin(B0t)))(zw1w2vz) \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 : tu=vzzu+velocity_v1v1u+velocity_v2v2u+velocity_vzvzu\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 vxfdv\int v_xf\,\mathrm{d}v et vyfdv\int v_yf\,\mathrm{d}v qui deviennent respectivement (w1cos(B0t)w2sin(B0t))gdwdvz\int (w_1\cos(B_0t) - w_2\sin(B_0t))g\mathrm{d}w\mathrm{d}v_z et (w1sin(B0t)+w2cos(B0t))gdwdvz\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 :

tf+vzzf(Ex+vyB0vzBy)vxf(EyvxB0+vzBx)vyf(vxByvyBx)vzf=0 \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 :

tf+vzzf(E+B0Jv)vf(v×(BxBy0))vf=0 \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,vz),v=(vx,vy)v = (v_\perp,v_z), \qquad v_\perp = (v_x,v_y)

On pose le changement de variable : v=etB0Jwv_\perp = e^{-tB_0J}w pour avoir une fonction gg :

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

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

tg=t(f(t,etB0Jw))=[tf](t,etB0Jw)[B0JetB0Jwvf](t,etB0Jw)=[vzzf+(E+B0Jv)vf+((vperpvz)×(BxBy0))vf](v=etB0Jw)[B0JetB0Jwvf](v=etB0Jw)=vzzf(etB0Jw)+E[vf](etB0Jw)+B0JetB0Jw[vf](etB0Jw)+(etB0Jwvz)×(BxBy0)[vf](etB0Jw)B0JetB0Jwvf(etB0Jw) \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 :

wg(t,z,w,vz)=wf(etB0Jw)=etB0J[vf](etB0Jw)etB0Jwg(w)=[vf](etB0Jw) \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ù :

tg=vzzg(t,z,w,vz)+etB0JEwg+((etB0Jwvz)×(BxBy0))(etB0Jwgvzg) \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 :

etB0J=(cssc)avec c=cos(B0t),s=sin(B0t) 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 :

etB0JEwg=(cssc)(ExEy)wg=(Exc+EysExs+Eyc)wg \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 :

((etB0Jwvz)×(BxBy0))(etB0Jwgvzg)=((cw1sw2sw1+cw2vz)×(BxBy0))(cw1gsw2gsw1g+cw2gvzg)=(vzByvzBx(cw1sw2)By(sw1+cw2)Bx)(cw1gsw2gsw1g+cw2gvzg) \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 :

tg=vzzg+(Exc+Eys+vzBxsvzBycExs+Eyc+vzBxc+vzBysBx(sw1+cw2)+By(cw1sw2))(w1gw2gvzg) \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 velocity_vx,y,z\texttt{velocity\_v}_{x,y,z} : velocity_vx=Exc+Eys+vzBxsvzBycvelocity_vy=Exs+Eyc+vzBxc+vzBysvelocity_vz=Bx(sw1+cw2)+By(cw1sw2) \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 : u˙=Lu+N(u) \dot{u} = Lu + N(u)

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

u(1)=eΔtLun+ΔteΔtLN(un)u(2)=34eΔt2Lun+14eΔt2Lu(1)+14ΔteΔt2LN(u(1))un+1=13eΔtLun+23eΔt2Lu(2)+23ΔteΔt2LN(u(2)) \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:uNuN:u\mapsto N\cdot u. Si LL et NN commutent, alors on trouve un+1=eΔtL[I+ΔtN+Δt22N2+Δt36N3]unu^{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)=[eΔtL+ΔteΔtLN]un 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+ΔtN)(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.

u(2)=[34eΔt2L+14eΔt2L(eΔtL+ΔteΔtLN)+14ΔteΔt2LN(eΔtL+ΔteΔtLN)]un=[eΔt2L+14Δt(eΔt2LN+eΔt2LNeΔtL)+14Δt2eΔt2LNeΔtLN]un \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 un+1u^{n+1} (après calcul) :

un+1=[eΔtL+Δt(23eΔt2LNeΔt2L+16eΔtLN+16NeΔtL)+Δt22(13NeΔtLN+13eΔt2LNeΔt2LN+13eΔt2LNeΔt2LNeΔtL)+Δt36eΔt2LNeΔt2LNeΔtLN]un \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 LL et NN 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=NLLN=NL :

un+1=eΔtL[I+ΔtN+Δt22N2+Δt36N3]un=eΔtL[eΔtN+O(Δt4N4)]un \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 LL et NN commutent donc eLeN=eL+Ne^{L}e^{N} = e^{L+N}.

un+1=eΔt(L+N)un+eΔtLO(Δt4N4)un 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 LL et NN 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 un+1=eΔtL[I+ΔtN+Δt22N2+Δt36N3]unu^{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^nLL et NN commutent. Cela revient donc à étudier un+1=[eΔtL+O((ΔtL))][eΔtN+O(Δt4N4)]unu^{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 LL et NN commutant cela revient à un+1=eΔt(L+N)+O()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 : u˙=Lu+N(u) \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(ss,nn), 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) :

u(1)=eΔt2Lun+12ΔteΔt2LN(un)u(2)=eΔt2Lun+12ΔteΔt2LN(u(1))u(3)=eΔtLun+ΔteΔt2LN(u(2))un+1=13eΔtLun+13eΔt2Lu(1)+23eΔt2Lu(2)+13u(3)+16ΔtN(u(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αΔtLe^{\alpha \Delta t L} est compliqué formellement (i.e. avec sympy), nous nous proposons de réécrire ce schéma en substituant la fonction MeMM\mapsto e^{M} par MPp,q(M)M\mapsto \mathrm{P}_{p,q}(M) pour l’approximant de Padé d’ordre (p,q)(p,q) de eMe^{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) : u(1)=Pp,q(Δt2L)un+12ΔtPp,q(Δt2L)N(un)u(2)=Pp,q(Δt2L)un+12ΔtPp,q(Δt2L)N(u(1))u(3)=Pp,q(ΔtL)un+ΔtPp,q(Δt2L)N(u(2))un+1=13Pp,q(ΔtL)un+13Pp,q(Δt2L)u(1)+23Pp,q(Δt2L)u(2)+13u(3)+16ΔtN(u(3)) \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 Pp,q\mathrm{P}_{p,q} est définie comme :

hp.q(z)=i=0pp!(pi)!(p+q)!(p+qi)!zii!kp.q(z)=j=0q(1)jq!(qj)!(p+q)!(p+qj)!zjj!ezPp,q(z)=hp,q(z)kp,q(z) \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 MM, on a : eMPp,q(M)=hp,q(M)(kp,q(M))1 e^{M} \approx \mathrm{P}_{p,q}(M) = h_{p,q}(M)\cdot \left(k_{p,q}(M)\right)^{-1}

Le calcul des Pp,q(αΔtL)\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 tP(tL)t\mapsto\mathrm{P}(t\cdot L) et de l’évaluer en certaines valeurs pour t=αΔtt=\alpha \Delta t. L’inversion n’est alors effectuée qu’une seule fois (lors de la construction de la fonction Pp,q[L]:tPp,q(tL)\mathrm{P}^{[L]}_{p,q}:t\mapsto\mathrm{P}_{p,q}(t\cdot L). L’évaluation de cette fonction en une valeur quelconque de tt (même un symbol comme Δt\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 :

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) :

{tjc,x=Ωpe2Exjc,yB0tjc,y=Ωpe2Ey+jc,xB0tBx=zEytBy=zExtEx=zByjc,x+vxfhdvtEy=zBxjc,y+vyfhdvtfh=vzzfh+(Ex+vyB0vzBy)vxfh+(EyvxB0+vzBx)vyfh+(vxByvyBx)vzfh \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 fhf_h et le calcul des courants jh,=vfhdvj_{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 : U˙=(A00ikvz)LU+N(U) \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 LL à différents temps. La matrice LL est une matrice diagonale par blocs, et l’exponentielle d’une matrice diagonale est l’exponentielle des termes diagonaux, c’est-à-dire : etL=(etA00eikvz) 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 etAe^{tA}, c’est-à-dire : etL(Pp,q(tA)00eikvz) 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 Pp.q(M)\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 fhf_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 LL est : L=(0B000Ωpe200B00000Ωpe200000000000000010000000100000000000vzz) 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 etLe^{tL} pour tout tt et tout kk (kk intervenant dans la transformée de Fourier en zz de la dérivée spatiale dans l’équation de Vlasov). Nous allons surtout étudier, pour B0=1B_0=1 et Ωpe2=4\Omega_{pe}^2=4, la matrice réduite AA, sans l’équation de Vlasov : A=(010040100004000000000000100000010000) 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 etAe^{tA} pour différentes valeurs de tt et comparer Pp,p(tA)\mathrm{P}_{p,p}(tA) pour ces mêmes valeurs de tt et différentes valeurs de pNp\in\mathbb{N}.

Tout d’abord on regarde etAe^{t A} : (0.378732187481834cos(2t9172)+0.621267812518167cos(2t17+92)0.378732187481834sin(2t9172)0.621267812518167sin(2t17+92)000.970142500145332sin(2t9172)+0.970142500145332sin(2t17+92)0.970142500145332cos(2t9172)+0.970142500145332cos(2t17+92)0.378732187481834sin(2t9172)+0.621267812518166sin(2t17+92)0.378732187481834cos(2t9172)+0.621267812518167cos(2t17+92)000.970142500145332cos(2t9172)0.970142500145332cos(2t17+92)0.970142500145332sin(2t9172)+0.970142500145332sin(2t17+92)001.00000001.0000.242535625036333sin(2t9172)0.242535625036333sin(2t17+92)0.242535625036333cos(2t9172)0.242535625036333cos(2t17+92)000.621267812518166cos(2t9172)+0.378732187481834cos(2t17+92)0.621267812518167sin(2t9172)0.378732187481834sin(2t17+92)0.242535625036333cos(2t9172)+0.242535625036333cos(2t17+92)0.242535625036333sin(2t9172)0.242535625036333sin(2t17+92)000.621267812518166sin(2t9172)+0.378732187481834sin(2t17+92)0.621267812518167cos(2t9172)+0.378732187481834cos(2t17+92)) \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 P3,3(tA)\mathrm{P}_{3,3}(t A) : (2t2(3240t85985t648600t4486000t2)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(600t910170t732400t5189000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340t2)(480t916200t764800t5243000t3+1620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+32400002t2(600t910170t732400t5189000t31620000t)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(480t9+16200t7+64800t5+243000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(3240t8+5985t6+48600t4+486000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340t2)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000002t2(1080t926370t797200t5432000t3)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(10440t82475t6102600t4486000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340t2)(12960t8+23940t6+194400t4+1944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+32400002t2(10440t82475t6102600t4486000t2+3240000)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(12960t823940t6194400t41944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(1080t9+26370t7+97200t5+432000t3)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340t2)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+32400002t2(600t910170t732400t5189000t31620000t)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(480t916200t764800t5243000t3+1620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(3240t85985t648600t4486000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340+t2)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+32400002t2(3240t8+5985t6+48600t4+486000t2)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(600t910170t732400t5189000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340+t2)(480t9+16200t7+64800t5+243000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000002t2(10440t82475t6102600t4486000t2+3240000)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(12960t8+23940t6+194400t4+1944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(1080t926370t797200t5432000t3)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340+t2)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+32400002t2(1080t9+26370t7+97200t5+432000t3)5(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(1t22)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t36+2t)(10440t82475t6102600t4486000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(3t340+t2)(12960t823940t6194400t41944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000001000000100t3(3240t85985t648600t4486000t2)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+t2(480t916200t764800t5243000t3+1620000t)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(600t910170t732400t5189000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000t3(600t910170t732400t5189000t31620000t)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+t2(7200t88460t6151200t4972000t2+3240000)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(3240t8+5985t6+48600t4+486000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(480t9+16200t7+64800t5+243000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+324000000t3(1080t926370t797200t5432000t3)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+t2(12960t8+23940t6+194400t4+1944000t2)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(10440t82475t6102600t4486000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000t3(10440t82475t6102600t4486000t2+3240000)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+t2(2400t9+40680t7+129600t5+756000t3+6480000t)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(1080t9+26370t7+97200t5+432000t3)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(12960t823940t6194400t41944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000t3(600t910170t732400t5189000t31620000t)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)t2(7200t88460t6151200t4972000t2+3240000)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(3240t85985t648600t4486000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(480t916200t764800t5243000t3+1620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000t3(3240t8+5985t6+48600t4+486000t2)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)t2(480t9+16200t7+64800t5+243000t31620000t)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(600t910170t732400t5189000t31620000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(7200t88460t6151200t4972000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+324000000t3(10440t82475t6102600t4486000t2+3240000)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)t2(2400t9+40680t7+129600t5+756000t3+6480000t)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(1080t926370t797200t5432000t3)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(12960t8+23940t6+194400t4+1944000t2)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000t3(1080t9+26370t7+97200t5+432000t3)30(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)t2(12960t823940t6194400t41944000t2)10(64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000)+(12t25)(10440t82475t6102600t4486000t2+3240000)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000+(t324t2)(2400t9+40680t7+129600t5+756000t3+6480000t)64t12+864t10+11124t8+105705t6+394200t4+1458000t2+3240000) \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
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 Pp,p(A)\mathrm{P}_{p,p}(A) convergent bien vers celles de eAe^A, et qu’elles sont très proches numériquement à partir de n3n\geq 3. La norme (ici de Frobenius) de Pp,p(tA)\mathrm{P}_{p,p}(tA) converge bien vers celle de etAe^{tA}, et à partir de n10n\geq 10 l’erreur entre les deux est inférieure à l’erreur machine. Le paramètre t[1,1]t\in[-1,1] représente le pas de temps Δt\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ΔtAe^{\Delta t A}, eΔt2Ae^{\frac{\Delta t}{2}A} et eΔt2Ae^{-\frac{\Delta t}{2}A} (en tout cas dans les méthodes que je considère), de plus le pas de temps Δt\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[1,1]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 (M1=1detMtcomAM^{-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 Nz×Nvx×Nvy×Nvz=15×20×20×41N_z\times N_{v_x}\times N_{v_y}\times N_{v_z} = 15\times 20\times 20\times 41, et Δt=0.1\Delta t=0.1 avec une méthode Lawson-RK(4,4), et un approximant de Padé P2,2\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 10710^{-7} jusqu’à 10510^{-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 : xexact LawsonnxPadeˊ Lawsonn,tn x^n_{\textrm{exact Lawson}} - x^n_{\textrm{Padé Lawson}}\,,\quad \forall t^n

avec xnx^n_{\star} correspondant à une énergie (énergie électrique, énergie magnétique ou énergie cinétique des particules froides) au temps tnt^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 Nz×Nvx×Nvy×Nvz=15×20×20×41N_z\times N_{v_x}\times N_{v_y}\times N_{v_z} = 15\times 20\times 20\times 41, et Δt=0.1\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=(0B000Ωpe200B00000Ωpe2000000z00000z00100z00001z0000000000vzz) 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 AA : A=(01004010000400000ik0000ik0100ik0001ik000) 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 ikik proviennent de la transformée de Fourier dans la direction zz des dérivées spatiales.

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

Pour différentes valeurs de kk on trace dans un premier temps les valeurs propre de eAe^{A} ainsi que celles de P3,3(A)\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 p3p\geq 3, les valeurs propres de eAe^A et Pn.n(A)\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 \mathrm{P}_{3,3}(A) pour différentes valeurs de k

Comparaison des valeurs propres de eAe^{A} et de P3,3(A)\mathrm{P}_{3,3}(A) pour différentes valeurs de k[ ⁣[Nz,Nz] ⁣]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 P3,3\mathrm{P}_{3,3}.

On peut aussi tracer, comme pour le cas sans les équations de Maxwell, l’erreur entre etAe^{tA} et Pn,n(tA)\mathrm{P}_{n,n}(tA) pour t[1,1]t\in[-1,1] et p[ ⁣[1,14] ⁣]p\in[\![1,14]\!], j’ai fixé regardé ici 2 modes, k=2,15k=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=15

Erreur entre etAe^{tA} et Pp,p(tA)\mathrm{P}_{p,p}(tA) pour deux modes kk

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 etAe^{tA} et P2,2(tA)\mathrm{P}_{2,2}(tA) pour différents modes de Fourier kk.

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))

Numériquement dans ce cas, toutes les valeurs propres sont de module 1 (l’échelle est en 1012+110^{-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))

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é P3,2\mathrm{P}_{3,2}, et on peut, sans problème généralisé pour tout nn et mm tels que n>mn>m.

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

Dans le cas de l’approximant de Padé P2,3(A(k))\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 nn et mm tels que m>nm>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}_{2,3}(A(k))

Valeurs prorpres de P3,2(A(k))\mathrm{P}_{3,2}(A(k)) (à gauche) et P2,3(A(k))\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 T5(A)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.

Les énergies au cours du temps

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 Δt\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 Δt=0.5\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 Δt=0.05\Delta t=0.05 que dans le papier, et un pas de temps plus grand Δt=0.1\Delta t=0.1 qui devient intéressant.

Les énergies au cours du temps

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…)

Les énergies au cours du temps

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)(1,1) à (2,2)(2,2), et il n’y a que le Padé (2,1)(2,1) qui explose, ce qui s’explique parce que le numérateur a un degré supérieur au dénominateur.

Les énergies au cours du temps

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).

Les énergies au cours du temps

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).

Les énergies au cours du temps

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 EE et de fhf_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.

Les énergies au cours du temps

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

Les énergies au cours du temps

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 Nz×Nvx×Nvy×Nvz=27×32×32×41N_z\times N_{v_x} \times N_{v_y} \times N_{v_z}=27\times 32\times 32 \times 41, jusqu’au temps final Tf=200T_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) :

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 Δt=0.1\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 : Nz×Nvx×Nvy×Nvz=27×32×32×41N_z \times N_{v_x} \times N_{v_y} \times N_{v_z}=27\times32\times32\times41 et Δt0=0.05\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 :

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.

Erreur relative des différentes simulations

2021-07-29

Valeurs propres de Taylor

On souhaite calculer une approximation de l’exponentielle de la partie linéaire AA de notre problème, où AA est définie par : A=(0100Ωpe2010000Ωpe200000iκ0000iκ0100iκ0001iκ000) 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 κ[ ⁣[Nz,Nz] ⁣]\kappa\in[\![-N_z,N_z]\!] et Nz=15N_z=15, et Ωpe=2\Omega_{pe}=2. On sait déjà que AA étant semblable à une matrice anti-hermitienne, ses valeurs propres sont imaginaires pures.

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

Les valeurs propres de AA étant imaginaire pures, on sait que les valeurs propres de etAe^{tA} sont sur le cercle unité (tR\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 Tm(A)T_m(A), avec m=1,,6m=1,\dots,6.

Valeur propre de Taylor 1

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

Valeur propre de Taylor 2

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 T2(A)T_2(A) sur une parabole. On a aucune valeur propre dans le cercle unité pour Taylor 1 et 2.

Valeur propre de Taylor 3

Pour T3(A)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é à T3T_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é à T1T_1, T2T_2 et T3T_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é
Valeur propre de Taylor 4

Pour T4T_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.

Valeur propre de Taylor 5

Un peu la même chose pour T5T_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).

Valeur propre de Taylor 6

Et là, je pense que cela pourrait être intéressant d’ajouter une simu avec T6T_6 car on remarque que les valeurs propres ne sont pas dans le cercle unité (elles en sont pas loin, comme pour T2T_2 et T3T_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 T4T_4 et T5T_5 sont stables (on allume suffisamment peu de modes, et pour ces quelques modes on est stable ?), mais illustrer que c’est instable avec T6T_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 zz 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].

Le nième point de Leja d’un intervalle KK est défini par : zn=argmaxzK(V(z)=i=0n1zzi) z_n = \arg\max_{z\in K}\left( V(z) = \prod_{i=0}^{n-1} |z-z_i| \right) et z0z_0 est pris tel que z0=maxzKz|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 civjκΔtc_iv_j\kappa\Delta t, avec les coefficients de Butcher ci[0,1]c_i\in[0,1] (généralement), la vitesse vj[2,2]Rv_j\in[-2,2]\approx\mathbb{R}, les modes de Fourier κ[ ⁣[13,13] ⁣]\kappa\in [\![-13,13]\!] et le pas de temps Δt[0,4]\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]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][-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 imaginaire 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]

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)(2,2).)

Partie réelle 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]
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é