Problème FPU : implémentation CUDA

1. Problème Fermi Pasta Ulam (FPU)

Le problème Fermi-Pasta-Ulam ([1]), consistait à étudier par une simulation informatique le système constitué d'une chaîne d'oscillateurs non linéaires couplés. Cette simuation est connue comme étant l'une des premières simulations (non militaires) de systèmes physiques sur ordinateur.

On reprend les notations définies dans la page chaîne d'oscillateurs linéaires à bords fixes. La force entre deux masses voisines de la chaîne comporte un terme non linéaire. L'équation du mouvement d'une masse s'écrit :

$$\frac{d^2{q}_k}{d{t}^2}=-{\omega }_0^2(2q_k-q_{k-1}-q_{k+1})+b[{(q_{k+1}-q_k)}^n-{(q_k-q_{k-1})}^n]\text{(1)}$$

Le terme non linéaire est soit quadratique (n=2), soit cubique (n=3). Ce dernier cas correspond à une force fonction impaire du déplacement, physiquement plus réaliste que le premier cas. Pour le calcul numérique, on pose ω₀=1.

La chaîne comporte N masses. Les bords de la chaîne sont supposés fixes.

Dans le cas d'une chaîne linéaire (b=0), les N modes propres sont indépendants les uns des autres et conservent leur énergie au cours du mouvement. Par exemple, si la chaîne est mise initialement dans le mode fondamental, elle le reste indéfiniment.

La physique statistique est une théorie qui repose sur le postulat suivant ([2]): pour un système isolé à l'équilibre, tous les états microscopiques accessibles au système ont la même probabilité. Un des résultats de la physique statistique est le théorème d'équipartition de l'énergie. Appliqué à la chaîne d'oscillateur, ce théorème implique une égale répartition de l'énergie entre les N modes propres. Ceci montre que la chaîne d'oscillateurs linéaires, initialement dans le mode fondamental, ne peut évoluer vers un état d'équilibre thermodynamique.

Lorsque des oscillateurs couplés comportent des non-linéarités, celles-ci permettent des échanges d'énergie entre modes propres. Avant le travail de Fermi-Pasta-Ulam, on pensait donc (sans preuve) que pour N grand, une petite non linéarité suffirait à faire évoluer la chaîne vers un état d'équilibre vérifiant le théorème d'équipartition de l'énergie. L'expérience numérique FPU visait à vérifier cette hypothèse.

2. Méthode de calcul

Les équations du mouvement (1) sont intégrés numériquement avec une méthode symplectique, la méthode de méthode de Stormer-verlet. Les méthodes symplectiques s'appliquent aux systèmes conservatifs et permettent d'obtenir une énergie à peu près constante sur les temps longs.

Initialement, la chaîne est placée dans un mode propre, par exemple le mode fondamental.

L'analyse du mouvement effectuée par Fermi-Pasta-Ulam consiste à calculer les coordonnées normales du problème linéaire (b=0), définies par :

$${\xi }_m=\sqrt{2/(N+1)}\sum _{k=1}^N{q}_{k}sin\left(\frac{\pi mk}{N+1}\right)\text{(2)}$$

Les dérivées par rapport au temps sont calculées par la même transformation. L'énergie du mode m est alors calculée par :

$$E_m=\frac{1}{2}({\dot{\xi }}_m^2+{\omega }_m^2{\xi }_m^2)\text{(3)}$$

Lorsque les non-linéarités sont faibles, on s'attend à une variation lente des énergies des modes propres. Il faut donc effectuer la simulation sur des temps très longs pour observer des variations.

Pour effectuer efficacement ce calcul, il est préférable utiliser l'algorithme FFT (Fast Fourier Transform), qui problement n'était pas connu en 1955 (algorithme établi en 1965 par Cooley et Tukey). Dans l'implémentation CUDA proposée ci-dessous, nous effectuons directement le produit matriciel (2). Les tests montrent que le temps de calcul des coordonnées normales avec cette méthode est infime par rapport au temps total de calcul.

3. Programme Python-CUDA

Le programme ci-dessous fpuCuda.py tourne avec une carte graphique NVidia compatible CUDA. Il nécessite le module Python PyCUDA et Numpy.

Chaque masse de la chaîne est traitée par un thread. Tous les threads appartiennent au même bloc. Le nombre N maximal de masses dans la chaîne est donc égal au nombre maximal de threads par bloc (512 ou 1024 sur les cartes les plus récentes).

Le module comporte une classe FpuCuda. Le nombre de variables nvar (c'est-à-dire le nombre N d'oscillateurs de la chaîne), est fourni au constructeur. Celui-ci définit les différentes fonctions (langage C) qui seront exécutées sur le GPU.

La première fonction définit la dérivée de l'énergie potentielle par rapport aux positions. Chaque Thread traite une variable, repérée par l'indice k.

import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
import numpy
import math
import time

class FpuCuda:
    def __init__(self,nvar):
        self.nvar = nvar
        f1 = """__device__ void dPotential(int k, int nvar, float b, float n, float *q, float *dp) {
            float x1,x2;
            if (k==0) {x1=q[0]; x2=q[1]-q[0];}
            else if (k==nvar-1) {x1 = q[k]-q[k-1]; x2=-q[k];}
            else {x1=q[k]-q[k-1]; x2=-q[k]+q[k+1];}
            dp[k] = x1-x2+b*(powf(x1,n)-powf(x2,n));
        }"""
            

La seconde fonction effectue un pas élémentaire de la méthode de Stormer-Verlet. Le pas de temps est h.

        f2 = """__device__ void verletStep(int k, float h, int nvar, float b, float n, float *q, float *p, float *dp) {
            dPotential(k,nvar,b,n,q,dp);
            p[k] -= 0.5*h*dp[k];
            q[k] += h*p[k];
            __syncthreads();
            dPotential(k,nvar,b,n,q,dp);
            p[k] -= 0.5*h*dp[k];
            __syncthreads();
        }"""    
            

La fonction suivante calcule l'énergie de la chaîne. Chaque thread calcule l'énergie d'une masse puis le calcul de la somme est partagé entre les threads.

        f3 = """__device__ void energy(int nvar,float b, float n, float *q, float *p, float *e) {
            float x1,x2;
            int k = threadIdx.x;
            if (k==0) {x1=q[k]; x2=0.0;}
            else if (k==nvar-1) {x1=q[k]-q[k-1]; x2=-q[k];}
            else {x1=q[k]-q[k-1]; x2=0.0;}
            e[k] = 0.5*p[k]*p[k]+0.5*(x1*x1+x2*x2)+1.0/(n+1)*b*(powf(x1,n+1)+powf(x2,n+1));
            __syncthreads();
            int i = blockDim.x/2;
            while (i!=0) {
                if (k < i) e[k] += e[k+i];
                __syncthreads();
                i /= 2;
            }
        }"""
            

La fonction suivante est appelée depuis le CPU (mot clé __global__). Elle effectue plusieurs pas de la méthode de Stormer-verlet (nsteps). L'état initial est fourni dans les tableaux q0 et p0, lesquels servent aussi à récupérer l'état final. Le tableau en (un seul élément) contient à la fin l'énergie du système. Les tableaux utilisés pour les calculs sont dans la mémoire locale du processeur (mot clé __shared__).

        f4="""__global__ void verletMultiSteps(int nsteps, float h, int nvar, float b, float n, float *q0, float *p0, float *en) {
            __shared__ float q[%d];
            __shared__ float p[%d];
            __shared__ float dp[%d];
            __shared__ float e[%d];
            int i;
            int k = threadIdx.x;
            q[k] = q0[k];
            p[k] = p0[k];
            __syncthreads();
            for (i=0; i<nsteps; i++) verletStep(k,h,nvar,b,n,q,p,dp);
            q0[k] = q[k];
            p0[k] = p[k];
            energy(nvar,b,n,q,p,e);
            if (k==0) en[0] = e[0];
        }"""%(nvar,nvar,nvar,nvar)
              

La fonction suivante calcule l'énergie des modes propres. Chaque bloc de threads calcule un mode. Le tableau A contient la matrice des vecteurs propres (matrice de passage des déplacements aux coordonnées normales). Le tableau omeg2 contient les valeurs propres (carré des pulsations propres). Les énergies des modes sont retournées dans le tableau emode.

        f5="""__global__ void modeEnergy(float *q, float *p, float *A, float *omeg2, float *emode) {
            __shared__ float xi[%d];
            __shared__ float dxi[%d];
            int i = blockIdx.x; // mode
            int j = threadIdx.x; // variable
            int index = j+blockDim.x*i;
            xi[j] = A[index]*q[j];
            dxi[j] = A[index]*p[j];
            __syncthreads();
            int j1 = blockDim.x/2;
            while (j1!=0) {
                if (j < j1) {
                    xi[j] += xi[j+j1];
                    dxi[j] += dxi[j+j1];
                }
                __syncthreads();
                j1 /= 2;
            }
            if (j==0) emode[i] = 0.5*(dxi[0]*dxi[0]+omeg2[i]*xi[0]*xi[0]);
        }"""%(nvar,nvar)     
              

Compilation du code CUDA :

        mod=SourceModule(f1+"\n"+f2+"\n"+f3+"\n"+f4+"\n"+f5)
        self.verletMultiSteps = mod.get_function("verletMultiSteps")
        self.modeEnergy = mod.get_function("modeEnergy")
              

Calcul de la matrice de changement de coordonnées (matrice des vecteurs propres) et du tableau des valeurs propres :

        self.A = numpy.empty(shape=(nvar,nvar),dtype=numpy.float32)
        s=math.sqrt(2.0/(nvar+1))
        r=math.pi/(nvar+1)
        for i in range(nvar):
            for j in range(nvar):
                self.A[i,j]=s*math.sin(r*(i+1)*(j+1))
        self.omeg2=numpy.empty(shape=nvar,dtype=numpy.float32)
        r=math.pi/(2*(nvar+1))
        for i in range(nvar):
            self.omeg2[i] = 4*math.pow(math.sin((i+1)*r),2)
              

Création des tableaux des déplacements et des vitesses :

        self.q0 = numpy.empty(shape=nvar,dtype=numpy.float32)
        self.p0 = numpy.empty(shape=nvar,dtype=numpy.float32)
              

La fonction init initialise la chaîne dans un mode, dont l'indice est précisé en argument.

    def init(self,m):
        for k in range(self.nvar):
            self.q0[k] = math.sin((k+1)*m*math.pi/(self.nvar+1))
            self.p0[k] = 0.0
              

La fonction run effectue un calcul complet. h est le pas de temps, n1 est la taille d'un bloc de calcul, c'est-à-dire le nombre de pas élémentaires entre deux calculs des coordonnées normales, n2 le nombre de blocs de calcul. b et n sont les paramètres du terme non linéaire, définis dans l'équation (1). L'argument nm est le nombre de modes à calculer.

    def run(self,h,n1,n2,b,n,nm):
        if nm > self.nvar:
            nm = self.nvar
        self.t = numpy.arange(n2)*h*n1
        self.qmat=numpy.empty(shape=(n2,self.nvar),dtype=numpy.float32)
        self.pmat=numpy.empty(shape=(n2,self.nvar),dtype=numpy.float32)
        en=numpy.empty(shape=(1),dtype=numpy.float32)
        self.emat=numpy.empty(shape=(n2,1),dtype=numpy.float32)
        self.emodemat=numpy.empty(shape=(n2,self.nvar),dtype=numpy.float32)
        #memoire GPU
        dev_q0 = cuda.mem_alloc(self.q0.nbytes)
        dev_p0 = cuda.mem_alloc(self.p0.nbytes)
        dev_en = cuda.mem_alloc(en.nbytes)
        dev_A = cuda.mem_alloc(self.A.nbytes)
        dev_omeg2 = cuda.mem_alloc(self.omeg2.nbytes)
        dev_emode = cuda.mem_alloc(self.q0.nbytes)
        cuda.memcpy_htod(dev_q0,self.q0)
        cuda.memcpy_htod(dev_p0,self.p0)
        cuda.memcpy_htod(dev_A,self.A)
        cuda.memcpy_htod(dev_omeg2,self.omeg2)
        ta = time.time()
        # boucle de calcul
        for i2 in range(n2): 
            self.verletMultiSteps(numpy.uint32(n1),numpy.float32(h),numpy.uint32(self.nvar),numpy.float32(b),numpy.float32(n),dev_q0,dev_p0,dev_en,grid=(1,1),block=(self.nvar,1,1))
            self.modeEnergy(dev_q0,dev_p0,dev_A,dev_omeg2,dev_emode,grid=(nm,1),block=(self.nvar,1,1))
            cuda.memcpy_dtoh(self.qmat[i2],dev_q0)
            cuda.memcpy_dtoh(self.pmat[i2],dev_p0)
            cuda.memcpy_dtoh(self.emat[i2],dev_en)
            cuda.memcpy_dtoh(self.emodemat[i2],dev_emode)
            print "%d sur %d"%(i2,n2)
        cuda.memcpy_dtoh(self.q0,dev_q0)
        cuda.memcpy_dtoh(self.p0,dev_p0)
        print "Duree du calcul (s) : %f"%(time.time()-ta)
        self.deltaE = numpy.std(self.emat)/numpy.average(self.emat)
        print "Ecart type relatif de E : %e"%self.deltaE
              

Les valeurs mémorisées à la fin de chaque bloc de calcul sont disponibles dans les tableaux suivants :

• t : instants
• qmat : déplacements des points de la chaîne.
• pmat : vitesses des points de la chaîne.
• emat : énergie totale de la chaîne.
• emodemat : énergies des modes propres.
• deltaE : écart type relatif de l'énergie totale.

4. Exemple

Dans l'exemple ci-dessous, on reprend un des calculs effectués dans [1]. La chaîne a 64 masses. Une faible non-linéarité quadratique est définie. La chaîne est placée initialement dans le mode 1.

import fpuCuda
from pylab import *
fpu=fpuCuda.FpuCuda(64)
fpu.init(1)
h=0.1
n1=2000
n2=1000
b=0.1
n=2
nm=10
fpu.run(h,n1,n2,b,n,nm)
            

On trace l'énergie des 10 premiers modes en fonction du temps (divisé par n1h) :

figure(0)
for m in range(nm):
    plot(fpu.emodemat[:,m])
axis([0,n2*1.5,0,0.05])
legend(("mode 1","mode 2","mode 3","mode4","mode 5","mode 6","mode 7",),loc='upper right')
xlabel("t/(n1h)")
ylabel("Em")

            

plotA.pdf

Au début du l'évolution, le mode 1 s'atténue rapidement au profit des modes 2,3 et 4. On constate néanmoins que le mode 1 ressurgit un peu plus tard (avec une moindre amplitude), avant de s'atténuer à nouveau. Il se produit une récurrence périodique du mode fondamental, et non pas une évolution vers un état d'équipartition des modes. Comme l'amplitude du mode 1 semble se réduire à chaque récurrence, il faut poursuivre la simulation sur des temps plus longs, ce que Fermi, Pasta et Ulam n'avaient pu faire.

On peut aussi obtenir une représentation de l'état de la chaîne à différents instants. Par exemple, au voisinage du moment où les modes 1,2 et 3 ont la même énergie :

figure(1)
for t in range(10):
    plot(fpu.qmat[100+t,:])
xlabel("k")
ylabel("qk")
            

plotB.pdf