In [1]:
import numpy as np
import scipy.fftpack as fftpack
import pylab as plt
import matplotlib as matplotlib

import pycuda.gpuarray as gpuarray

#-------------------------------------------------------------------------------------
from pywignercuda_path import SetPyWignerCUDA_Path
SetPyWignerCUDA_Path()
from GPU_WignerDirac2D_MassPotential_4x4 import *

In [2]:
%matplotlib inline

In [3]:
class Klein(GPU_WignerDirac2D_4x4):
    def __init__ (self):
    #....................Defining the geometry..................................... 
        X_gridDIM = 512
        P_gridDIM = 512
        
        X_amplitude = 15  
        P_amplitude = 15   

        
        timeSteps  =  600
        dt = 0.01 #dX/c
        
        skipFrames =    100

        #...................Defining the kinematic-dynamical constants.................
        
        mass = 1.
        c = 1.
        
        #self.dt = dX/self.c
        #...................Defining the potential and initial state parameters........
        V0 = 0.
        w  = 0.25
        
        #.........................ODM damping ........................................
        self.gammaDamping = 0.0
        
        #............................................................................
        
        self.D_Theta      = 0.01
        self.D_Lambda     = 0.0
        
        self.grossPitaevskiiCoefficient = 0.0
        
        #self.pX = 9.5
        self.Potential_0_String = '{0}*0.5*( 1. + tanh((x-5.)/{1})   )'.format(V0,w)
        self.Potential_1_String = ' 0.'
        self.Potential_2_String = ' 0.'
        self.Potential_3_String = ' 0.'
        self.Mass_String = '1. + 0.1*x*x'

        #.............................................................................
        GPU_WignerDirac2D_4x4.__init__(self,
            X_gridDIM, P_gridDIM, X_amplitude, P_amplitude, mass, c, dt,
            timeSteps,skipFrames,frameSaveMode='Density',antiParticleNorm = True, computeEnergy=True)
        #.............................................................................
        
          
   
        
    def  Set_Initial_Majorana  (self) :

        #..................Defining the output directory/file ........................

        self.fileName ='/home/rcabrera/DATA/WignerDirac2D_4x4/MassMajorana_D'+str(self.D_Theta)+'.hdf5'
        
        self.W_init = np.empty([4,4,instance.P_gridDIM,instance.X_gridDIM],dtype = np.complex128)
        
        init_x  = 0.
        self.pX = 5.
        
        psiL = self.GaussianSpinor_ParticleUp(  -init_x , self.pX, 1, self.X - 0.5*self.Theta ) + 0j
        psiL = self.MajoranaSpinorPlus(psiL)
        #psiL = self.MajoranaSpinorMinus(psiL)
            
            
        psiR = self.GaussianSpinor_ParticleUp(  -init_x , self.pX, 1, self.X + 0.5*self.Theta )  +0j 
        psiR =  self.MajoranaSpinorPlus(psiR)
        #psiR =  self.MajoranaSpinorMinus(psiR)
        
        for i in range(4):
            for j in range(4):
                self.W_init[i,j][:,:] = psiL[i]*psiR[j].conj()
        
        # To XP       
        self.Fourier_4X4_Theta_To_P(self.W_init)
        
        #instance.FilterElectrons( self.W_init , 1)
        
        norm = self.Wigner_4x4_Norm(self.W_init)
        self.W_init *= 1./ norm        
        
    def  Set_Initial_Cat_State  (self) :
 
        #..................Defining the output directory/file ........................

        self.fileName = '/home/rcabrera/DATA/WignerDiracMass2D_4x4/MassCatState_D'+str(self.D_Theta)+'.hdf5'
        
        self.W_init = np.empty([4,4,instance.P_gridDIM,instance.X_gridDIM],dtype = np.complex128)
        
        init_x  = 0
        self.pX = 5.
        
        #
        
        psiL1  = self.GaussianSpinor_ParticleUp(   init_x , self.pX , 1, self.X - 0.5*self.Theta )       
        psiL1 += self.GaussianSpinor_ParticleUp(   init_x , -self.pX , 1, self.X - 0.5*self.Theta )    
            
            
        psiR1  = self.GaussianSpinor_ParticleUp(   init_x , self.pX, 1, self.X + 0.5*self.Theta )         
        psiR1 += self.GaussianSpinor_ParticleUp(   init_x , -self.pX, 1, self.X + 0.5*self.Theta )         
        #
        
        for i in range(4):
            for j in range(4):
                self.W_init[i,j][:,:] = psiL1[i]*psiR1[j].conj()
        
        # To XP       
        self.Fourier_4X4_Theta_To_P(self.W_init)
        
        #instance.FilterElectrons( self.W_init , 1)
        
        norm = self.Wigner_4x4_Norm(self.W_init)
        self.W_init *= 1./ norm

In [4]:
instance = Klein()


/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1797: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(15): warning: variable "P_gridDIM" was declared but never referenced

kernel.cu(15): warning: variable "P_gridDIM" was declared but never referenced


  transmission_source%(self.CUDA_constants)).get_function("transmission_Kernel")
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1809: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(30): warning: variable "p" was declared but never referenced

kernel.cu(30): warning: variable "p" was declared but never referenced


  self.CUDA_constants,self.Potential_0_String),arch="sm_20").get_function("Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1813: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(30): warning: variable "p" was declared but never referenced

kernel.cu(30): warning: variable "p" was declared but never referenced


  self.CUDA_constants,self.Potential_1_String),arch="sm_20").get_function("Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1817: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(30): warning: variable "p" was declared but never referenced

kernel.cu(30): warning: variable "p" was declared but never referenced


  self.CUDA_constants,self.Potential_2_String),arch="sm_20").get_function("Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1821: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(30): warning: variable "p" was declared but never referenced

kernel.cu(30): warning: variable "p" was declared but never referenced


  self.CUDA_constants,self.Potential_3_String),arch="sm_20").get_function("Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1828: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(17): warning: variable "i" was declared but never referenced

kernel.cu(17): warning: variable "i" was declared but never referenced


  SourceModule( X_Average_source%(self.CUDA_constants),arch="sm_20").get_function( "Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1831: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(18): warning: variable "j" was declared but never referenced

kernel.cu(18): warning: variable "j" was declared but never referenced


  SourceModule( P_Average_source%(self.CUDA_constants),arch="sm_20").get_function( "Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1837: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced

kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced


  SourceModule( XX_Average_source%(self.CUDA_constants),arch="sm_20").get_function( "Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1840: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced

kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced


  SourceModule( Alpha_1_Average_source%(self.CUDA_constants),arch="sm_20").get_function( "Kernel" )
/home/rcabrera/Documents/source/Python/PyWignerCUDA/GPU_WignerDirac2D_MassPotential_4x4.py:1846: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced

kernel.cu(13): warning: variable "P_gridDIM" was declared but never referenced


  SourceModule( X1_Alpha_1_Average_source%(self.CUDA_constants),arch="sm_20").get_function( "Kernel" )
  D_1_Potential_0 =  0 + 0.*x

In [5]:
#instance.Set_Initial_Majorana()
instance.Set_Initial_Cat_State()

instance.Run ()


----------------------------------------------
 Relativistic Wigner-Dirac Propagator:  x-Px  
----------------------------------------------
 dt      =  0.01
 dx      =  0.05859375
 dp      =  0.05859375
 dLambda =  0.209439510239
            
         GPU memory Total        5.17700195312 GB
         GPU memory Free         4.82745361328 GB
 progress  0 %
 cuda grid =   ( (512, 1, 1)  ,  (512, 1, 1) )
 progress  16 %
 progress  33 %
 progress  49 %
 progress  66 %
 progress  83 %
 progress  99 %
 computation time =  47.9556221962  seconds

In [6]:
def PlotWigner(W):
    
    W0 = fftpack.fftshift(  instance.Wigner_4X4__SpinTrace( W ).real  )
    
    x_min = -instance.X_amplitude
    x_max = instance.X_amplitude - instance.dX
    
    p_min = -instance.P_amplitude
    p_max = instance.P_amplitude - instance.dP
    
    global_max = 0.31          #  Maximum value used to select the color range
    global_min = -0.27        # 

    print 'min = ', np.min( W0 ), ' max = ', np.max( W0 )
    print 'normalization = ', np.sum( W0 )*instance.dX*instance.dP

    zero_position =  abs( global_min) / (abs( global_max) + abs(global_min)) 
    wigner_cdict = {'red' 	: 	((0., 0., 0.),
							(zero_position, 1., 1.), 
							(1., 1., 1.)),
					'green' :	((0., 0., 0.),
							(zero_position, 1., 1.),
							(1., 0., 0.)),
					'blue'	:	((0., 1., 1.),
							(zero_position, 1., 1.),
							(1., 0., 0.)) }
    wigner_cmap = matplotlib.colors.LinearSegmentedColormap('wigner_colormap', wigner_cdict, 256)
    #wigner_cmap = plt.colors.LinearSegmentedColormap('wigner_colormap', wigner_cdict, 256)
    

    fig, ax = plt.subplots(figsize=(20, 7))

    
        
    cax = ax.imshow( W0 ,origin='lower',interpolation='nearest',\
    extent=[x_min, x_max, p_min, p_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)

    ax.set_xlabel('x')
    ax.set_ylabel('p')
    #ax.set_xlim((x_min,x_max))
    #ax.set_ylim((-5 , p_max/3.5))
    #ax.set_ylim((-16,16))    
    ax.set_aspect(1)
    ax.grid('on')

In [7]:
def PlotWignerDensity(W):
    
    W0 = fftpack.fftshift( W )
    
    x_min = -instance.X_amplitude
    x_max = instance.X_amplitude - instance.dX
    
    p_min = -instance.P_amplitude
    p_max = instance.P_amplitude - instance.dP
    
    global_max = 0.31          #  Maximum value used to select the color range
    global_min = -0.27        # 

    print 'min = ', np.min( W0 ), ' max = ', np.max( W0 )
    print 'normalization = ', np.sum( W0 )*instance.dX*instance.dP

    zero_position =  abs( global_min) / (abs( global_max) + abs(global_min)) 
    wigner_cdict = {'red' 	: 	((0., 0., 0.),
							(zero_position, 1., 1.), 
							(1., 1., 1.)),
					'green' :	((0., 0., 0.),
							(zero_position, 1., 1.),
							(1., 0., 0.)),
					'blue'	:	((0., 1., 1.),
							(zero_position, 1., 1.),
							(1., 0., 0.)) }
    wigner_cmap = matplotlib.colors.LinearSegmentedColormap('wigner_colormap', wigner_cdict, 256)
    #wigner_cmap = plt.colors.LinearSegmentedColormap('wigner_colormap', wigner_cdict, 256)
    

    fig, ax = plt.subplots(figsize=(20, 7))

    
        
    cax = ax.imshow( W0 ,origin='lower',interpolation='nearest',\
    extent=[x_min, x_max, p_min, p_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)

    ax.set_xlabel('x')
    ax.set_ylabel('p')
    
    ax.set_aspect(1)
    ax.grid('on')

In [8]:
def PlotMarginal_P(instance):
    
    W_0 = fftpack.fftshift(instance.Wigner_4X4__SpinTrace(instance.W_init).real)
        
    print ' norm =  ', np.sum(W_0).real*instance.dX*instance.dP
    
    fig, ax = plt.subplots(figsize=(10, 5))

    prob_P = np.sum(W_0,axis=1)*instance.dX
    ax.plot(instance.P_range, prob_P , label = 'init')
    
    W_0 = fftpack.fftshift( instance.Wigner_4X4__SpinTrace(instance.W_end).real )
    
    print ' norm =  ', np.sum(W_0).real*instance.dX*instance.dP
    
    prob_P = np.sum(W_0,axis=1)*instance.dX
    ax.plot(instance.P_range, prob_P , label = 'final')
    
    ax.set_xlim(-18,18)
    ax.set_xlabel('p')
    ax.set_ylabel('Prob')
    ax.grid('on')
    
    ax.legend(bbox_to_anchor=(0.75, 0.5), loc=2, prop={'size':22})
    
def PlotMarginal_X(instance):
    
    W_0 = fftpack.fftshift( instance.Wigner_4X4__SpinTrace(instance.W_init).real)
        
    
    fig, ax = plt.subplots(figsize=(10, 5))

    prob_X = np.sum(W_0,axis=0)*instance.dP
    ax.plot(instance.X_range, prob_X , label = 'init')
    
    W_0 = fftpack.fftshift( instance.Wigner_4X4__SpinTrace(instance.W_end).real )
    
    
    prob_X = np.sum(W_0,axis=0)*instance.dP
    ax.plot(instance.X_range, prob_X , label = 'final')
    
    ax.set_xlabel('x')
    ax.set_ylabel('Prob')
    ax.grid('on')
    
    ax.legend(bbox_to_anchor=(0.75, 0.5), loc=2, prop={'size':22})

In [9]:
PlotWigner(10*instance.W_init)


min =  -0.560096764418  max =  1.59094241878
normalization =  10.0

In [10]:
PlotWigner( 10*instance.W_end )
print ' time = ', instance.timeRange[-1]


min =  -0.215347446587  max =  1.42291461328
normalization =  10.0
 time =  6.0

In [11]:
PlotMarginal_P( instance )
PlotMarginal_X( instance )


 norm =   1.0
 norm =   1.0

In [12]:
WendFilter = instance.W_end.copy()
instance.FilterElectrons( WendFilter , 1)

In [13]:
PlotWigner( 10*WendFilter )


min =  -0.0842993825887  max =  1.32830795434
normalization =  9.24953068696

Ehrenfest Theorems


In [14]:
axis_font = {'size':'24'}


fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,  np.gradient( instance.X_Average.real , instance.dt)  , 'g',
        label= '$\\frac{dx^1}{dt} $')

ax.plot( instance.timeRange[1:] ,  instance.Alpha_1_Average.real ,'r--' ,label='$c \\alpha^1$')

ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=1, prop={'size':22})


Out[14]:
<matplotlib.legend.Legend at 0x7fcbc0dbde50>

In [15]:
axis_font = {'size':'24'}


fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,  np.gradient( instance.P_Average.real , instance.dt)  , 'g',
        label= '$\\frac{dp^1}{dt} $')

ax.plot( instance.timeRange[1:] , 
        -instance.D_1_Potential_0_Average.real - 2.*instance.mass*instance.gammaDamping*instance.Alpha_1_Average.real ,'r--' ,label='$-c e\, \\partial_1 A_0  $')


ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[15]:
<matplotlib.legend.Legend at 0x7fcbc048ac50>

In [16]:
axis_font = {'size':'24'}

fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,
        2*np.gradient( instance.XP_Average.real , instance.dt)  , 'g',
        label= '$\\frac{d}{dt}( x^1 p_1 ) $')

ax.plot( instance.timeRange[1:] ,
        -2*instance.X1_D_1_Potential_0_Average.real + 2*instance.c*instance.P1_Alpha_1_Average.real -4.*instance.mass*instance.gammaDamping*instance.X1_Alpha_1_Average,
        'r--' ,label='$-2 x^1 \\partial_1 e A_0  + 2 c p_1 \\alpha^1$')


ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[16]:
<matplotlib.legend.Legend at 0x7fcbc0e491d0>

In [17]:
axis_font = {'size':'24'}

fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,
        np.gradient( instance.XX_Average.real , instance.dt)  , 'g',
        label= '$\\frac{d}{dt}( x^1 x_1 ) $')

ax.plot( instance.timeRange[1:] ,
      2*instance.c*instance.X1_Alpha_1_Average.real,
        'r--' ,label='$2 c x_1 \\alpha^1$')


ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[17]:
<matplotlib.legend.Legend at 0x7fcbc0f5a5d0>

In [18]:
axis_font = {'size':'24'}

fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,
        np.gradient( instance.PP_Average.real , instance.dt)  , 'g',
        label= '$\\frac{d}{dt}( x^1 x_1 ) $')

ax.plot( instance.timeRange[1:] ,
       -2*instance.P1_D_1_Potential_0_Average.real - \
       4.*instance.mass*instance.gammaDamping*instance.P1_Alpha_1_Average.real,  
        'r--' ,label='$- c p^1 \\partial_1 A_0 $')


ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[18]:
<matplotlib.legend.Legend at 0x7fcbc14f6d90>

In [19]:
axis_font = {'size':'24'}

fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,
        instance.antiParticle_population.real  , 'g',
        label= 'Antiparticle population')


ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[19]:
<matplotlib.legend.Legend at 0x7fcbc0e96f90>

In [20]:
axis_font = {'size':'24'}

fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange[1:] ,
        instance.Dirac_energy.real  , 'g',
        label= 'Energy')

ax.set_ylim(-1, 1)
ax.set_xlabel(r'$t$',**axis_font)

ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})


Out[20]:
<matplotlib.legend.Legend at 0x7fcbc1076c10>

Loading files


In [21]:
W0 = instance.Load_Density(instance.fileName,0)
PlotWignerDensity( 10*W0 )


min =  -0.560096764418  max =  1.59094241878
normalization =  10.0

In [22]:
W0 = instance.Load_Density(instance.fileName,300)
PlotWignerDensity( 10*W0 )


min =  -0.729437963057  max =  1.49809385506
normalization =  10.0

In [23]:
W0 = instance.Load_Density(instance.fileName,600)
PlotWignerDensity( 10*W0 )


min =  -0.215347446587  max =  1.42291461328
normalization =  10.0

In [23]: