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_DiracDaviau2D import *
In [2]:
%matplotlib inline
In [3]:
class DefinePotential_DiracDaviau2D( GPU_DiracDaviau2D ):
def __init__ (self):
#....................Defining the geometry.....................................
gridDIM_X = 512
gridDIM_Y = 512
X_amplitude =18
Y_amplitude =18
dt = 0.01
timeSteps = 500
skipFrames = 10
#...................Defining the kinematic-dynamical constants.................
self.mass = 1.
self.c = 1.
#...................Defining the potential and initial state parameters........
hBar = 1
#self.a = 2.
#self.r = 5.
self.px = 1.
self.py = 0.
self.energy = np.sqrt( (self.mass*self.c)**2 + self.px**2 + self.py**2 )
#Q = 1./( self.mass + self.energy )
#...................Defining the potential.....................................
self.Potential_0_String = '0.'
self.Potential_1_String = '0.'
self.Potential_2_String = '0.'
self.Potential_3_String = '0.'
#..................Defining the output directory/file .........................
self.fileName = '/home/rcabrera/DATA/DiracDaviau2D/Free.hdf5'
#..............................................................................
amplitude = (X_amplitude,Y_amplitude)
gridDIM = (gridDIM_X,gridDIM_Y)
#method = 'Hybrid'
#method = 'Bandrauk'
GPU_DiracDaviau2D.__init__( self,gridDIM, amplitude ,dt,timeSteps, skipFrames = skipFrames,frameSaveMode='Spinor')
def Set_Initial_Condition (self):
def gaussian(x,y):
return np.exp( - ( x**2 + y**2)/4. )
self.Psi_init = self.Spinor_Particle_SpinUp( (self.px,self.py) , gaussian )
#self.Psi_init = self.Spinor_Particle_SpinDown( (self.px,self.py) , gaussian )
#self.Psi_init = self.Spinor_AntiParticle_SpinUp( (self.px,self.py) , gaussian )
#self.Psi_init = self.Spinor_AntiParticle_SpinDown( (self.px,self.py) , gaussian )
norm = self.Norm( self.Psi_init )
self.Psi_init /= norm
In [4]:
instance = DefinePotential_DiracDaviau2D()
instance.Set_Initial_Condition()
#instance_Dirac3D.FilterElectrons(1)
print ' '
print ' dx = ', instance.dX
Psi_end = instance.Run ()
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))
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 )
#PlotDensity( instance.load_Spinor(100) )
#PlotDensity( instance.load_Spinor(200) )
#PlotDensity( instance.load_Spinor(500) )
In [9]:
PlotDensity( instance.Psi_end );
In [9]:
In [10]:
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[10]:
In [11]:
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[11]:
In [12]:
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[12]:
In [13]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , instance.KEnergy_average , 'r',
label= '$Kinetic Energy $')
ax.plot( instance.timeRange , instance.Potential_0_average , 'g',
label= '$Potential Energy $')
ax.plot( instance.timeRange , instance.Potential_0_average + instance.KEnergy_average, 'b--',
label= '$ K + A_0 $')
ax.set_ylim(-0.1,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]:
instance.Px;
In [15]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.X_average , instance.Px_average , 'r-', label= '$ x $')
ax.set_ylim(0,2)
ax.set_xlim(-0.1,3)
ax.set_xlabel(r'$x$',**axis_font)
ax.set_ylabel(r'$p_x$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.9), loc=1, prop={'size':22})
Out[15]:
In [16]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.Y_average , instance.Py_average , 'r-', label= '$ y $')
#ax.set_ylim(-10,10)
#ax.set_xlim(-100000,10000)
ax.set_xlabel(r'$y$',**axis_font)
ax.set_ylabel(r'$p_y$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[16]:
In [17]:
axis_font = {'size':'24'}
fig, ax = plt.subplots(figsize=(20, 7))
ax.plot( instance.timeRange , instance.TotalProbabilityX_average , 'g',
label= '$ Probability\, X $')
ax.plot( instance.timeRange , instance.TotalProbabilityP_average , 'r--',
label= '$ Probability\, P $')
ax.set_ylim(0.,1.1)
ax.set_xlabel(r'$t$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})
Out[17]:
In [18]:
Psi_init_filtered = instance.FilterElectrons(-1, instance.Psi_init)
Psi_end_filtered = instance.FilterElectrons(-1, instance.Psi_end)
In [21]:
PlotDensity( 2*Psi_init_filtered )
In [22]:
PlotDensity( 2*Psi_end_filtered )
In [20]: