In [1]:
import numpy as np
import pycuda.gpuarray as gpuarray
from pycuda.tools import make_default_context
import matplotlib as matplotlib
import pylab as plt
#-------------------------------------------------------------------------------------
from pywignercuda_path import SetPyWignerCUDA_Path
SetPyWignerCUDA_Path()
from GPU_Dirac2D import *
In [2]:
%matplotlib inline
In [3]:
class DefinePotential_Dirac2D( GPU_Dirac2D ):
def __init__ (self):
#....................Defining the geometry.....................................
gridDIM_X = 2*256
gridDIM_Y = 2*256
X_amplitude =12
Y_amplitude =12
dt = 0.01
timeSteps = 1000
skipFrames = 100
#...................Defining the kinematic-dynamical constants.................
self.mass = 1.
self.c = 1.
self.hBar = 1.
#...................Defining the potential and initial state parameters........
self.B = 1
self.py = 1
#...................Defining the potential.....................................
#self.Potential_0_String = ' 2 '
#self.Potential_1_String = '-0.5*{qB}*y'.format( qB = qB )
#self.Potential_2_String = ' 0.5*{qB}*x'.format( qB = qB )
#self.Potential_3_String = ' 0. '
p0 = np.sqrt( (self.mass*self.c)**2 + self.py**2 )
self.Potential_0_String = ' {K}*x '.format(K= self.B*self.py/(2*self.c*self.mass))
self.Potential_1_String = ' {K1}*( {K2}*t - y/{p0} )'.format(
K1=0.5*self.B*self.c*self.mass , K2=self.c*self.py**2/(p0**2), p0=p0 )
self.Potential_2_String = ' {K}*x '.format( K=self.B*p0/(2*self.c*self.mass) )
self.Potential_3_String = ' 0. '
#..................Defining the output directory/file .........................
self.fileName = '/home/rcabrera/DATA/Dirac2D/LinearDisplacement.hdf5'
#..............................................................................
amplitude = (X_amplitude,Y_amplitude)
gridDIM = (gridDIM_X,gridDIM_Y)
GPU_Dirac2D.__init__( self,gridDIM, amplitude ,dt,timeSteps, skipFrames = skipFrames,frameSaveMode='Spinor')
def Set_Initial_Condition (self):
def gaussian(x,y):
return np.exp( - self.B/(4*self.c*self.hBar)*( x**2 + y**2) )
self.Psi_init = self.Spinor_Particle_SpinUp( (0.,self.py) , gaussian )
norm = self.Norm( self.Psi_init )
self.Psi_init /= norm
In [4]:
instance = DefinePotential_Dirac2D()
instance.Set_Initial_Condition()
#instance_Dirac3D.FilterElectrons(1)
print ' '
print ' dx = ', instance.dX
instance.Run ()
Out[4]:
In [5]:
Rho_end = instance.Density_From_Spinor( instance.Psi_end )
In [6]:
Rho_end.shape
Out[6]:
In [7]:
def PlotDensity( Psi ):
global_max = 0.12 # Maximum value used to select the color range
global_min = -0.01
Rho = instance.Density_From_Spinor(Psi)
Rho = fftpack.fftshift( Rho )
x1_min = -instance.X_amplitude
x1_max = instance.X_amplitude - instance.dX
p1_min = -instance.Y_amplitude
p1_max = instance.Y_amplitude - instance.dY
print ' max Rho = ', np.max(Rho)
#print ' min Rho = ', np.min(Rho)
print ' normalization = ', np.sum( Rho )*instance.dX*instance.dY
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, 512)
fig, ax = plt.subplots(figsize=(20, 7))
ax.grid('on')
cax = ax.imshow( Rho \
,origin='lower',interpolation='nearest',\
extent=[x1_min, x1_max, p1_min, p1_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)
In [8]:
PlotDensity( instance.Psi_init )
In [9]:
PlotDensity( instance.load_Spinor(100) )
PlotDensity( instance.load_Spinor(200) )
PlotDensity( instance.load_Spinor(300) )
PlotDensity( instance.Psi_end )
In [10]:
instance.X_average.shape, instance.timeRange.shape
Out[10]:
In [11]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , np.gradient(instance.X_average , instance.dt) , 'g',
label= '$\\frac{dx^1}{dt} $')
ax.plot( instance.timeRange , instance.Alpha1_average ,'r--' ,label='$c \\alpha^1$')
ax.set_xlabel(r'$t$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[11]:
In [12]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , np.gradient(instance.Y_average , instance.dt) , 'g',
label= '$\\frac{dx^2}{dt} $')
ax.plot( instance.timeRange , instance.Alpha2_average ,'r--' ,label='$c \\alpha^2$')
ax.set_xlabel(r'$t$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[12]:
In [13]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , instance.Beta_average , 'g',
label= '$\\langle\\gamma^0\\rangle $')
ax.set_xlabel(r'$t$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[13]:
In [14]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , instance.Potential_0_average + instance.KEnergy_average, 'b-',
label= '$ K + A_0 $')
ax.set_ylim(-0.,2.5)
ax.set_xlabel(r'$t$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[14]:
In [15]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.X_average , instance.Y_average , 'r-')
#ax.set_ylim(-8,5)
#ax.set_xlim(-3,3)
ax.set_xlabel(r'$x$',**axis_font)
ax.set_ylabel(r'$y$',**axis_font)
#ax.legend(bbox_to_anchor=(1.05, 0.9), loc=1, prop={'size':22})
Out[15]:
In [16]:
Psi_end_filtered = instance.FilterElectrons(-1, instance.Psi_end)
In [17]:
PlotDensity( 2*Psi_end_filtered )
In [17]: