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
import matplotlib.pyplot as mplt
#-------------------------------------------------------------------------------------
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 = 512
gridDIM_Y = 512
X_amplitude =12
Y_amplitude =12
dt = 0.005
timeSteps = 8000
skipFrames = 500
#...................Defining the kinematic-dynamical constants.................
self.mass = 1.
self.c = 1.
self.hBar = 1.
#...................Defining the potential and initial state parameters........
self.B = 1
# Parameters of the circular trajectory
self.omega = 0.1
self.r = 6.
#...................Defining the potential.....................................
#p0 = np.sqrt( (self.mass*self.c)**2 + self.py**2 )
W = np.sqrt( self.c**2 - (self.r*self.omega)**2 )
self.Potential_0_String = '({omega}*pow(pow({c},2) - pow({omega},2)*pow({r},2),-0.5)*({c}*{hBar} + ({B} + 2*{c}*{m}*{omega})*{r}*x*cos({omega}*t) - {B}*pow({r},2) + ({B} + 2*{c}*{m}*{omega})*{r}*y*sin({omega}*t)))/2.'.format(
m=self.mass,c=self.c,hBar=self.hBar,B=self.B,omega=self.omega,r=self.r,W=W)
self.Potential_1_String = '(pow({c},-1)*pow({W},-1)*(2*{B}*y*pow({omega},2)*pow( {r},2)*pow(cos({omega}*t),2) + 2*{c}*{r}*({B}*{c} - {hBar}*pow({omega},2))*sin({omega}*t) - {B}*(2*y*pow({c},2) + x*pow({omega},2)*pow({r},2)*sin(2*{omega}*t))))/4.'.format(
m=self.mass,c=self.c,hBar=self.hBar,B=self.B,omega=self.omega,r=self.r,W=W)
self.Potential_2_String = '(pow({c},-1)*pow({W},-1)*(-2*{c}*{r}*cos({omega}*t)*({B}*{c} - {hBar}*pow({omega},2)) + {B}*(2*x*pow({c},2) - 2*x*pow({omega},2)*pow({r},2)*pow(sin({omega}*t),2) + y*pow({omega},2)*pow({r},2)*sin(2*{omega}*t))))/4.'.format(
m=self.mass,c=self.c,hBar=self.hBar,B=self.B,omega=self.omega,r=self.r,W=W)
self.Potential_3_String = ' 0. '
#..................Defining the output directory/file .........................
self.fileName = '/home/rcabrera/DATA/Dirac2D/Packet_Circular_Displacement.hdf5'
#..............................................................................
amplitude = (X_amplitude,Y_amplitude)
gridDIM = (gridDIM_X,gridDIM_Y)
GPU_Dirac2D.__init__( self,gridDIM, amplitude ,dt,timeSteps, skipFrames = skipFrames,frameSaveMode='Spinor')
#GPU_Dirac2D.__init__( self,gridDIM, amplitude ,dt,timeSteps, skipFrames = skipFrames,frameSaveMode='Density')
def Set_Initial_Condition (self):
def gaussian(x,y):
return np.exp( - self.B/(4*self.c*self.hBar)*( (x-self.r)**2 + y**2) )
self.Psi_init = self.Spinor_Particle_SpinUp( (0.,self.r*self.omega) , 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='none',\
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 [10]:
#PlotDensity( instance.load_Spinor(200) )
#PlotDensity( instance.load_Spinor(400) )
PlotDensity( instance.load_Spinor(1000) )
PlotDensity( instance.Psi_end )
In [11]:
instance.X_average.shape, instance.timeRange.shape
Out[11]:
In [12]:
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[12]:
In [13]:
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[13]:
In [14]:
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[14]:
In [15]:
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[15]:
In [16]:
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.set_aspect(1)
#ax.legend(bbox_to_anchor=(1.05, 0.9), loc=1, prop={'size':22})
In [17]:
def PosterPlotDensity( Psi , Psi2, Psi3):
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='none',\
extent=[x1_min, x1_max, p1_min, p1_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)
#
Rho2 = instance.Density_From_Spinor(Psi2)
Rho2 = fftpack.fftshift( Rho2 )
masked_data = np.ma.masked_where(Rho2 < 0.002, Rho2)
ax.imshow( masked_data \
,origin='lower',interpolation='none',\
extent=[x1_min, x1_max, p1_min, p1_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)
#
Rho3 = instance.Density_From_Spinor(Psi3)
Rho3 = fftpack.fftshift( Rho3 )
masked_data = np.ma.masked_where(Rho3 < 0.002, Rho3)
ax.imshow( masked_data \
,origin='lower',interpolation='none',\
extent=[x1_min, x1_max, p1_min, p1_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)
plt.show()
In [18]:
PosterPlotDensity( instance.load_Spinor(0), instance.load_Spinor(1000), instance.load_Spinor(2000) )
In [19]:
Psi_end_filtered = instance.FilterElectrons(-1, instance.Psi_end)
In [20]:
PlotDensity( 2*Psi_end_filtered )
In [ ]: