2D Hybrid Simulation of the Jeans Instability of a Dusty Plasma

Original ver.; 1998 Nov. 03, T. Hada
Rewritten for Python; 2018 May 26, M. Nakanotani

basic equations:

  1. Equation of motion for dust particles $$\frac{d{\bf r}}{dt}={\bf v}, \frac{d{\bf v}}{dt}=-\nabla(\phi({\bf r})+\psi({\bf r}))$$
  2. Electric potential is determined by $$(\nabla^2-\frac{1}{\lambda^2})\phi=-\Sigma_j({\bf r}-{\bf r}_j)$$
  3. Gravity Potential is determined by $$\nabla^2\psi=\Gamma\Sigma_j({\bf r}-{\bf r}_j).$$

Here, $\phi$ is the electric potential, $\psi$ is the gravitational potential, $\Gamma=(\omega_{Jd}/\omega_{pd})^2$ is the squared ratio between the Jeans and the dust plasma frequencies, and $\lambda$ is the Debye length defined from the ambient plasma.

Input parameters:

physical:

  • dlmd = plasma Debye length
  • ggg = $(\omega_{Jd}/\omega_{pd})^2$
  • vth = initial dust thermal velocity

control:

  • ppc = # of particles in a cell
  • nx, ny = # of meshes in the $x$ and $y$ directions
  • dx, dy = grid size in the $x$ and $y$ directions
  • ntmax = maximum time to be computed
  • tmax = maximum time steps to be computed
In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from numba.decorators import jit
from numba import f8,u8
import scipy.fftpack as sf
from IPython import display
In [2]:
def jeans2(ppc,nx,ny,dx,dy,ntmax,tmax,dtmax,dlmd,ggg,vth):
    np.seterr(divide='ignore', invalid='ignore')
    nptcl=ppc*nx*ny
    lx=nx*dx
    ly=ny*dy
    #xx=dx*np.arange(nx)
    #yy=dy*np.arange(ny); X,Y=np.meshgrid(xx,yy)
    t =0.0
    xp=np.random.rand(nptcl)*lx
    yp=np.random.rand(nptcl)*ly
    vx=np.random.normal(0,vth,nptcl)
    vy=np.random.normal(0,vth,nptcl)

    ds =np.zeros((ny,nx))
    epx=np.zeros((ny,nx))
    epy=np.zeros((ny,nx))
    gpx=np.zeros((ny,nx))
    gpy=np.zeros((ny,nx))
    kx=2*mt.pi/lx*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    ky=2*mt.pi/ly*np.r_[np.arange(ny/2),np.arange(-ny/2,0)]
    kkx,kky=np.meshgrid(kx,ky)    
    
    for it in range(ntmax):
        dt=check_dt(dx,dy,vx,vy,dtmax)
        t =t+dt
        if(t>tmax): break
        print(it,dt,t)

        xp=xp+dt*vx
        yp=yp+dt*vy
        xp[xp>lx]=xp[xp>lx]-lx
        xp[xp<0 ]=xp[xp<0 ]+lx
        yp[yp>ly]=yp[yp>ly]-ly
        yp[yp<0 ]=yp[yp<0 ]+ly

        ds=np.zeros((ny,nx))
        ds=dens(dx,dy,nx,ny,nptcl,xp,yp,ds)#ds =np.exp(-0.01*((X-lx/2)**2+(Y-ly/2)**2))
        phi=sf.fft2(ds)/(kkx**2+kky**2+1.0/dlmd**2)
        psi=-ggg*sf.fft2(ds)/(kkx**2+kky**2)
        psi[0,0]=0
        epx=np.real(sf.ifft2(1j*kkx*phi))
        epy=np.real(sf.ifft2(1j*kky*phi))
        gpx=np.real(sf.ifft2(1j*kkx*psi))
        gpy=np.real(sf.ifft2(1j*kky*psi))

        vx,vy=acc(dx,dy,nx,ny,nptcl,xp,yp,vx,vy,epx,epy,gpx,gpy,dt)

        display.clear_output(True)
        plt.figure(figsize=(18,6))
        plt.subplot(131);plt.imshow(np.real(sf.ifft2(phi)),cmap='jet',origin='lower',aspect='auto');plt.title('Electric Potential')
        plt.subplot(132);plt.imshow(np.real(sf.ifft2(psi)),cmap='jet',origin='lower',aspect='auto');plt.title('Gravity Potential')
        plt.subplot(133);plt.plot(xp,yp,'k,');plt.title('Motion of Dust Particles');plt.xlim(0,nx);plt.ylim(0,ny)
        plt.show()

        #plt.plot(xp,vx,',')
        #plt.show()
    return locals()

@jit('f8(f8,f8,f8[:],f8[:],f8)')
def check_dt(dx,dy,vx,vy,dtmax):
    vmaxx=max(vx)
    vmaxy=max(vy)
    ddt  =min(dtmax,dx/vmaxx,dy/vmaxy)
    return ddt

@jit('f8[:](f8,f8,u8,u8,u8,f8[:],f8[:],f8[:],f8[:],f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8)')
def acc(dx,dy,nx,ny,nptcl,xp,yp,vx,vy,epx,epy,gpx,gpy,dt):
    for ip in range(nptcl):
        ixm=mt.floor(xp[ip]/dx); ixp=ixm+1
        iym=mt.floor(yp[ip]/dy); iyp=iym+1
        wxp=xp[ip]/dx-ixm; wxm=1-wxp
        wyp=yp[ip]/dy-iym; wym=1-wyp

        if ixp>nx-1: ixp=ixp-nx
        if iyp>ny-1: iyp=iyp-ny

        vx[ip]=vx[ip]-dt*(wxm*wym*epx[iym,ixm]+wxm*wyp*epx[iyp,ixm]+wxp*wym*epx[iym,ixp]+wxp*wyp*epx[iyp,ixp]\
                         +wxm*wym*gpx[iym,ixm]+wxm*wyp*gpx[iyp,ixm]+wxp*wym*gpx[iym,ixp]+wxp*wyp*gpx[iyp,ixp])
        vy[ip]=vy[ip]-dt*(wxm*wym*epy[iym,ixm]+wxm*wyp*epy[iyp,ixm]+wxp*wym*epy[iym,ixp]+wxp*wyp*epy[iyp,ixp]\
                         +wxm*wym*gpy[iym,ixm]+wxm*wyp*gpy[iyp,ixm]+wxp*wym*gpy[iym,ixp]+wxp*wyp*gpy[iyp,ixp])

    return vx,vy

@jit('f8[:,:](f8,f8,u8,u8,u8,f8[:],f8[:],f8[:,:])')
def dens(dx,dy,nx,ny,nptcl,xp,yp,ds):
    for ip in range(nptcl):
        ixm=mt.floor(xp[ip]/dx); ixp=ixm+1
        iym=mt.floor(yp[ip]/dy); iyp=iym+1
        wxp=xp[ip]/dx-ixm; wxm=1-wxp
        wyp=yp[ip]/dy-iym; wym=1-wyp

        if ixp>nx-1: ixp=ixp-nx
        if iyp>ny-1: iyp=iyp-ny

        ds[iym,ixm]=ds[iym,ixm]+wxm*wym
        ds[iym,ixp]=ds[iym,ixp]+wxp*wym
        ds[iyp,ixm]=ds[iyp,ixm]+wxm*wyp
        ds[iyp,ixp]=ds[iyp,ixp]+wxp*wyp
    ds=ds*nx*ny/nptcl
    return ds
In [3]:
%%time 
data=jeans2(20,128,128,1,1,100,6,0.2,1,1,0.05) #ppc,nx,ny,dx,dy,ntmax,maxt,dtmax,dlmd,ggg,vth
locals().update(data)

CPU times: user 1min 40s, sys: 5.08 s, total: 1min 45s
Wall time: 2min