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 mpl_toolkits.mplot3d import Axes3D


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


    

In [2]:

    

#import pyqtgraph as pg


    

In [3]:

    

%matplotlib inline


    

In [4]:

    

class DefinePotential_Dirac3D( GPU_Dirac3D ):
    def __init__ (self):
        #....................Defining the geometry.....................................    
        gridDIM_X = 256
        gridDIM_Y = 256
        gridDIM_Z = 256
        
        X_amplitude =15
        Y_amplitude =15
        Z_amplitude =15
        
        dt = 0.01

        timeSteps =  1200
        skipFrames =  100 
        #...................Defining the kinematic-dynamical constants.................
        self.mass = 1.
        self.c = 1.
        self.hBar = 1.
        #...................Defining the potential and initial state parameters........
        hBar = 1
    
        self.px = 5.
        self.py = 0.
        self.pz = 0.
        
        self.energy = self.c*np.sqrt( (self.mass*self.c)**2 +  self.px**2 + self.py**2 + self.pz**2  )
        #...................Defining the potential.....................................
       
        V0 = 7.
        w  = 0.3
        
        self.Potential_0_String = '{0}*0.5*( 1. + tanh((x-2.)/{1})   )'.format(V0,w)
        
        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/Dirac3D/Klein.hdf5'

        self.Compute_Ehrenfest_P       = True
        self.antiParticleStepFiltering = False
        #..............................................................................	

        amplitude = (X_amplitude,Y_amplitude,Z_amplitude)
        gridDIM   = (gridDIM_X,gridDIM_Y,gridDIM_Z)

        self.Compute_Ehrenfest_P = True

        GPU_Dirac3D.__init__( self,gridDIM, amplitude ,dt,timeSteps, skipFrames = skipFrames,frameSaveMode='Spinor')

        
    def Set_Initial_Condition_SpinUp (self):

        def gaussian(x,y,z):
            x0 = -3    
            return np.exp( - ( (x-x0)**2 + y**2 + z**2)/2. )
        
        self.Psi_init = self.Spinor_Particle_SpinUp(  (self.px,self.py,self.pz) , gaussian )
        
        norm = self.Norm( self.Psi_init )
        
        #self.FilterElectrons(1) 
        
        self.Psi_init /= norm  
        
        
    """def Set_Initial_Condition_SpinDown (self):

        def gaussian(x,y,z):
            return np.exp( - ( x**2 + y**2 + z**2)/2. )
        
        self.Psi_init = self.Spinor_Particle_SpinDown(  (self.px,self.py,self.pz) , gaussian )
        
        norm = self.Norm_X_CPU( self.Psi_init )
        
        self.Psi_init /= norm  """


    

In [5]:

    

instance = DefinePotential_Dirac3D()


    

    




/home/rcabrera/Documents/source/python/PyWignerCUDA/GPU_Dirac3D.py:753: UserWarning: The CUDA compiler succeeded, but said the following:
kernel.cu(13): warning: variable "dPx" was declared but never referenced

kernel.cu(14): warning: variable "dPy" was declared but never referenced

kernel.cu(14): warning: variable "dPz" was declared but never referenced

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

kernel.cu(14): warning: variable "dPy" was declared but never referenced

kernel.cu(14): warning: variable "dPz" was declared but never referenced


  A3=self.Potential_3_String) ).get_function( "DiracPropagatorA_Kernel" )

In [6]:

    

instance = DefinePotential_Dirac3D()

instance.Set_Initial_Condition_SpinUp()

#instance_Dirac3D.FilterElectrons(1)

print '                                      '
print ' dx = ', instance.dX

Psi_end = instance.Run ()


    

    




                                      
 dx =  0.1171875
--------------------------------------------
              Dirac Propagator 3D           
--------------------------------------------
  save Mode  =   Spinor
         GPU memory Total        5.24945068359 GB
         GPU memory Free         2.55571746826 GB
                                                               
number of steps  =   1200  dt =  0.01
dX =  0.1171875 dY =  0.1171875 dZ =  0.1171875
dPx =  0.209439510239 dPy =  0.209439510239 dPz =  0.209439510239
                                                               
  
 progress  0 %
 norm step 1 =  9786.46696204
 progress  8 %
 progress  16 %
 progress  24 %
 progress  33 %
 progress  41 %
 progress  49 %
 progress  58 %
 progress  66 %
 progress  74 %
 progress  83 %
 progress  91 %
 progress  99 %
 computational time  =  890.230374813

In [7]:

    

instance.Psi_end.shape


    

    
Out[7]:





(4, 256, 256, 256)

In [8]:

    

axis_font = {'size':'24'}

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

ax.plot( instance.timeRange ,  instance.X_average   , 'g',label= '$x^1 $')

ax.plot( instance.timeRange ,  instance.Y_average   , 'r',label= '$x^2 $')

ax.plot( instance.timeRange ,  instance.Z_average   , 'b--',label= '$x^3 $')


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

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


    

    
Out[8]:





<matplotlib.legend.Legend at 0x7f24ac26bf50>






    

In [9]:

    

if instance.Compute_Ehrenfest_P == True:
    axis_font = {'size':'24'}
    
    fig, ax = plt.subplots(figsize=(20, 7))
    
    ax.plot( instance.timeRange ,  instance.Px_average   , 'g',label= '$p^1 $')
    
    ax.plot( instance.timeRange ,  instance.Py_average   , 'r',label= '$p^2 $')
    
    ax.plot( instance.timeRange ,  instance.Pz_average   , 'b--',label= '$p^3 $')
    
    ax.set_ylim(-6,6)
    
    ax.set_xlabel(r'$t$',**axis_font)
    
    ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})


    

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]:





<matplotlib.legend.Legend at 0x7f24b489efd0>






    

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{dy^1}{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]:





<matplotlib.legend.Legend at 0x7f24b4cab590>






    

In [19]:

    

axis_font = {'size':'24'}


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

ax.plot( instance.timeRange ,  np.gradient( instance.Z_average , instance.dt)  , 'g',
        label= '$\\frac{dx^3}{dt} $')

ax.plot( instance.timeRange ,  instance.Alpha3_average ,'r--' ,label='$c \\alpha^3$')

ax.set_ylim(-0.1,0.1)

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

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


    

    
Out[19]:





<matplotlib.legend.Legend at 0x7f24654b3450>






    

In [13]:

    

if instance.Compute_Ehrenfest_P == True:
    
    axis_font = {'size':'24'}
    
    fig, ax = plt.subplots(figsize=(20, 7))
    
    ax.plot( instance.timeRange ,  np.gradient( instance.Px_average , instance.dt)  , 'g',
            label= '$\\frac{dp^1}{dt} $')
    
    
    ax.set_xlabel(r'$t$',**axis_font)
    
    ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})


    

In [14]:

    

if instance.Compute_Ehrenfest_P == True:
    
    axis_font = {'size':'24'}
    
    fig, ax = plt.subplots(figsize=(20, 7))
    
    ax.plot( instance.timeRange ,  np.gradient( instance.Py_average , instance.dt)  , 'g',
            label= '$\\frac{dp^2}{dt} $')
    
    
    ax.set_xlabel(r'$t$',**axis_font)
    
    ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})


    

In [15]:

    

if instance.Compute_Ehrenfest_P == True:
    
    axis_font = {'size':'24'}
    
    fig, ax = plt.subplots(figsize=(20, 7))
    
    ax.plot( instance.timeRange ,  np.gradient( instance.Pz_average , instance.dt)  , 'g',
            label= '$\\frac{dp^3}{dt} $')
    
    
    ax.set_xlabel(r'$t$',**axis_font)
    
    ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})


    

In [16]:

    

if instance.Compute_Ehrenfest_P == True:
    
    axis_font = {'size':'24'}
    
    fig, ax = plt.subplots(figsize=(20, 7))
    
    ax.plot( instance.timeRange ,  instance.K_Energy_average + instance.Potential_0_average , 'g',
            label= '$ Energy$')
    
    ax.set_ylim(0,7)
    
    ax.set_xlabel(r'$t$',**axis_font)
    
    ax.grid('on')
    ax.legend(bbox_to_anchor=(1.05, 0.5), loc=1, prop={'size':22})


    

In [17]:

    

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


ax = fig.gca(projection='3d')

ax.scatter(
            instance.X_average,
            instance.Y_average,
            instance.Z_average
            )

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

ax.set_ylabel(r'$y$',**axis_font)

ax.set_zlabel(r'$z$',**axis_font)

ax.set_xlim(-2 ,6)
ax.set_ylim(-1,1)
ax.set_zlim(-0.1,0.1)


    

    




/usr/lib/pymodules/python2.7/matplotlib/collections.py:548: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == 'face':






    
Out[17]:





(-0.1, 0.1)






    

In [ ]:

    

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


#ax = fig.gca(projection='3d')

ax.scatter(
            instance.X_average,
            instance.Y_average            )

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

ax.set_ylabel(r'$y$',**axis_font)

ax.set_xlim(-1.5 ,1.5)
ax.set_ylim(-1   ,1)
#ax.set_zlim(-5,0)


    

In [ ]:

    

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


#ax = fig.gca(projection='3d')

ax.scatter(
            instance.X_average,
            instance.Y_average            )

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

ax.set_ylabel(r'$y$',**axis_font)

ax.set_xlim(-1.5 ,1.5)
ax.set_ylim(-1   ,1)
#ax.set_zlim(-5,0)