Add ELBDM/UniformGranule test problem

example/test_problem/ELBDM/UniformGranule/Input__TestProb

@@ -0,0 +1,20 @@
+# problem-specific runtime parameters
+ComputeCorrelation      1                        # compute correlation function, must set OPT__RECORD_USER=1 [0]
+ReComputeCorrelation    0                        # use the simulation time of RESTART as initial time for computing time correlation; only available for RESTART [0]
+FilePath_corr           ./Record__Correlation    # output path for correlation function text files
+MinLv                   0                        # do statistics from MinLv to MaxLv [0]
+MaxLv                   0                        # do statistics from MinLv to MaxLv [MAX_LEVEL]
+StepInitial             0                        # inital step for recording correlation function [0]
+StepInterval            1                        # interval for recording correlation function (only for OutputCorrelationMode=0) [1
+StepEnd                 200                      # end step for recording correlation function (only for OutputCorrelationMode=0) [NoMax_int]
+
+# parameters for correlation function profile (only useful when ComputeCorrelation=1)
+RadiusMax_corr          0.04                     # maximum radius of correlation profile [0.5*BoxSize]
+dr_min_corr             0.00234224               # minimum bin size [1e-2*0.5*BoxSize]
+LogBin_corr             1                        # use log bin [0]
+LogBinRatio_corr        1.1940617521534918       # exp((log(RadiusMax_corr)-log(dr_min_corr))/NBin) [2]
+RemoveEmpty_corr        0                        # remove 0 sample bins; false: Data[empty_bin]=Weight[empty_bin]=NCell[empty_bin]=0 [0]
+
+# parameters for initial density profile (only useful when ComputeCorrelation=1)
+dr_min_prof             0.00058556               # minimum bin size [0.25*dr_min_corr]
+

example/test_problem/ELBDM/UniformGranule/README

@@ -0,0 +1,32 @@
+Compilation flags:
+========================================
+Enable  : MODEL=ELBDM, GRAVITY, SUPPORT_HDF5, NCOMP_PASSIVE_USER=1
+Disable : COMOVING
+
+
+Default setup:
+=======================================
+0. Copy "generate_make.sh" to the directory "src/," execute "sh generate_make.sh" to generate "Makefile,"
+   and then execute "make clean" and "make -j 4" to generate the executable "gamer"
+
+1. Code units
+   (1) UNIT_L = kpc
+   (2) UNIT_V = km/s
+   (3) UNIT_D = 1.0e10 Msun/kpc^3
+
+2. ELBDM_MASS              3.0e-19
+   ELBDM_TAYLOR3_AUTO      0
+   ELBDM_BASE_SPECTRAL     1
+
+3. OPT__BC_FLU_*           1  (periodic)
+   OPT__BC_POT             1  (periodic)
+
+
+Default UM_IC info:
+=======================================
+BoxSize                 =  0.08 kpc
+dimension               =  128^3
+average density         =  0.25e10 Msun/kpc^3
+1-D velocity dispersion =  5.1 km/s
+
+

example/test_problem/ELBDM/UniformGranule/clean.sh

@@ -0,0 +1,3 @@
+rm -rf Record__* Data_* PowerSpec_* log
+mkdir  Record__Correlation
+touch  Record__Correlation/.empty

example/test_problem/ELBDM/UniformGranule/download_ic.sh

@@ -0,0 +1,5 @@
+filename=uniform-granule-ic
+
+curl https://girder.hub.yt/api/v1/item/66e0028e32f323dee1b801dc/download -o ${filename}.tgz
+tar -zxvf ${filename}.tgz
+rm ${filename}.tgz

example/test_problem/ELBDM/UniformGranule/generate_make.sh

@@ -0,0 +1,6 @@
+# This script should run in the same directory as configure.py
+
+PYTHON=python3
+
+${PYTHON} configure.py --machine=eureka_intel --mpi=true --hdf5=true --fftw=FFTW3 --gpu=true \
+                       --model=ELBDM --gravity=true --passive=1

example/test_problem/ELBDM/UniformGranule/python_script/get_vel_dispersion.py

@@ -0,0 +1,85 @@
+import yt
+import numpy as np
+import argparse
+import sys
+
+#-------------------------------------------------------------------------------------------------------------------------
+# load the command-line parameters
+parser = argparse.ArgumentParser( description='Get average velocity dispersion of the entire box' )
+
+parser.add_argument( '-s', action='store', required=True,  type=int, dest='idx_start',
+                     help='first data index' )
+parser.add_argument( '-e', action='store', required=True,  type=int, dest='idx_end',
+                     help='last data index' )
+parser.add_argument( '-d', action='store', required=False, type=int, dest='didx',
+                     help='delta data index [%(default)d]', default=1 )
+
+args=parser.parse_args()
+
+idx_start   = args.idx_start
+idx_end     = args.idx_end
+didx        = args.didx
+
+# print command-line parameters
+print( '\nCommand-line arguments:' )
+print( '-------------------------------------------------------------------' )
+for t in range( len(sys.argv) ):
+   print( str(sys.argv[t]))
+print( '' )
+print( '-------------------------------------------------------------------\n' )
+
+for idx in range(idx_start, idx_end+1, didx):
+   ds = yt.load("../Data_%06d"%idx)
+   N = ds.parameters["NX0"]
+   N_tot = N[0]*N[1]*N[2]
+   UNIT_L = ds.parameters["Unit_L"]
+   ma   = ds.parameters["ELBDM_Mass"]
+   hbar = ds.parameters["ELBDM_PlanckConst"]
+   fac = hbar/ma
+
+   grad_Real = ds.add_gradient_fields(("Real"))
+   grad_Imag = ds.add_gradient_fields(("Imag"))
+
+   dd = ds.all_data()
+
+   if ( N_tot != len(dd["Dens"])):
+      print('file size not matched!')
+      sys.exit(1)
+
+   dens = np.array(dd["Dens"])
+   real = np.array(dd["Real"])
+   imag = np.array(dd["Imag"])
+   avedens = np.mean(dens)
+   print("average density = " + str(avedens))
+   print("minimum density = " + str(min(dens)))
+
+   grad_real = np.array([dd["Real_gradient_x"], dd["Real_gradient_y"], dd["Real_gradient_z"]])*UNIT_L
+   grad_imag = np.array([dd["Imag_gradient_x"], dd["Imag_gradient_y"], dd["Imag_gradient_z"]])*UNIT_L
+
+   vx_bk = fac*(imag*grad_real[0] - real*grad_imag[0])/dens
+   vy_bk = fac*(imag*grad_real[1] - real*grad_imag[1])/dens
+   vz_bk = fac*(imag*grad_real[2] - real*grad_imag[2])/dens
+
+   vx_qp = fac*(imag*grad_imag[0] + real*grad_real[0])/dens
+   vy_qp = fac*(imag*grad_imag[1] + real*grad_real[1])/dens
+   vz_qp = fac*(imag*grad_imag[2] + real*grad_real[2])/dens
+
+   sigma_x_square_bk = np.average(dens*vx_bk**2)/avedens - (np.average(dens*vx_bk))**2/avedens**2
+   sigma_y_square_bk = np.average(dens*vy_bk**2)/avedens - (np.average(dens*vy_bk))**2/avedens**2
+   sigma_z_square_bk = np.average(dens*vz_bk**2)/avedens - (np.average(dens*vz_bk))**2/avedens**2
+
+   sigma_x_square_qp = np.average(dens*vx_qp**2)/avedens - (np.average(dens*vx_qp))**2/avedens**2
+   sigma_y_square_qp = np.average(dens*vy_qp**2)/avedens - (np.average(dens*vy_qp))**2/avedens**2
+   sigma_z_square_qp = np.average(dens*vz_qp**2)/avedens - (np.average(dens*vz_qp))**2/avedens**2
+
+
+   sigma_bk = ((sigma_x_square_bk + sigma_y_square_bk +sigma_z_square_bk)/3.)**0.5
+   sigma_qp = ((sigma_x_square_qp + sigma_y_square_qp +sigma_z_square_qp)/3.)**0.5
+
+   sigma_total = (sigma_bk**2+sigma_qp**2)**0.5
+   print('bulk velocity dispersion = %g'%(sigma_bk))
+   print('quantum pressure velocity dispersion = %g'%(sigma_qp))
+   print('total velocity dispersion = %g'%(sigma_total))
+
+   d = 0.35*2*np.pi*hbar/(ma*sigma_total)
+   print('estimated granule diameter = %g\n'%(d))

example/test_problem/ELBDM/UniformGranule/python_script/plot_correlation.py

@@ -0,0 +1,73 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import glob
+import re
+import sys
+import yt
+
+# -------------------------------------------------------------------------------------------------------------------------
+# user-specified parameters
+sigma  = 5.1                                  # velocity dispersion in km/s
+# -------------------------------------------------------------------------------------------------------------------------
+ds     = yt.load( '../Data_000000' )
+UNIT_L = ds.parameters["Unit_L"]              # code length unit (default is kpc)
+UNIT_V = ds.parameters["Unit_V"]              # code velocity unit (default is km/s)
+UNIT_T = ds.parameters["Unit_T"]
+ma     = ds.parameters["ELBDM_Mass"]          # FDM particle mass in code unit
+h_bar  = ds.parameters["ELBDM_PlanckConst"]   # reduced planck constant in code unit
+d      = 0.35*2.0*np.pi*h_bar/ma/sigma        # granule diameter in code uniti (kpc)
+
+# C(t) decay time scale, C(t_corr)=(1+2*0.35**2)**(-1.5)*C(t=0) = 0.72*C(t=0)
+t_corr = d/(2**0.5*np.pi*sigma)*UNIT_T/(3600.*24.*365.*1e6) # expected correlation time scale in Myr
+
+print("velocity dispersion             = %g km/s"%sigma)
+print("estimated granule size          = %g kpc"%d)
+print("expected correlation time scale = %g Myr"%t_corr)
+print("")
+
+file_dir = '../Record__Correlation/'
+files = glob.glob(file_dir + 'correlation_function_t=*.txt')
+
+filename = np.array(files)
+
+r = []
+C = []
+t = []
+for f in filename:
+   Corr = np.loadtxt(f, skiprows=1, dtype=float)
+   if not r:
+      r.append(Corr[:,0])
+   else:
+      if not np.array_equal(r[0], Corr[:,0]):
+         print('radius bin not matched!! filename=\"%s\"'%f)
+         sys.exit(1)
+
+   C.append(Corr[:,1])
+   match = re.search(r'correlation_function_t=(\d+\.\d+e[+-]\d+)', f)
+   if match:
+      time = float(match.group(1))
+      t.append(time)
+   else:
+      print('time pattern not matched!! filename=\"%s\"'%f)
+      sys.exit(1)
+
+if not r:
+   print("no file loaded, check ../Record__Correlation/ !!")
+
+t = np.array(t)*UNIT_T/(3600.*24.*365.*1e6)
+C = np.array(C)
+
+color = plt.cm.turbo(np.linspace(0.1, 0.9, len(r[0])))
+
+plt.figure()
+for i in range(len(r[0])):
+   plt.plot(t, C[:,i], c=color[i], label = r'$r$ = %1.3e kpc'%(r[0][i]))
+plt.axvline(x=t_corr, ls='--', lw = '0.9', color = '0.7', label = 'expected correlation time scale')
+plt.xlabel('$t$ (Myr)')
+plt.ylabel(r'$C(t)$')
+plt.xlim(0, t[-1])
+plt.legend(bbox_to_anchor=(1.03,0.03), loc='lower left')
+plt.savefig('fig_correlation.png', dpi = 150, bbox_inches="tight")
+plt.close()
+
+

example/test_problem/ELBDM/UniformGranule/python_script/plot_powerspectrum.py

@@ -0,0 +1,84 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import yt
+import argparse
+import sys
+
+#-------------------------------------------------------------------------------------------------------------------------
+# load the command-line parameters
+parser = argparse.ArgumentParser( description='Get power spectrum' )
+
+parser.add_argument( '-s', action='store', required=True,  type=int, dest='idx_start',
+                     help='first data index' )
+parser.add_argument( '-e', action='store', required=True,  type=int, dest='idx_end',
+                     help='last data index' )
+parser.add_argument( '-d', action='store', required=False, type=int, dest='didx',
+                     help='delta data index [%(default)d]', default=1 )
+
+args=parser.parse_args()
+
+idx_start   = args.idx_start
+idx_end     = args.idx_end
+didx        = args.didx
+
+# print command-line parameters
+print( '\nCommand-line arguments:' )
+print( '-------------------------------------------------------------------' )
+for t in range( len(sys.argv) ):
+   print( str(sys.argv[t]))
+print( '' )
+print( '-------------------------------------------------------------------\n' )
+
+# plot GAMER generated power spectrum as reference
+GAMER_PS = True
+print('GAMER_PS='+str(GAMER_PS)+'\n')
+
+for idx in range(idx_start, idx_end+1, didx):
+   ds = yt.load("../Data_%06d"%idx)
+   N = ds.parameters["NX0"][0]
+   BoxSize = ds.parameters["BoxSize"][0]
+   dh = ds.parameters["CellSize"][0]
+
+   dd = ds.covering_grid(level=0, left_edge=[0, 0, 0], dims=ds.domain_dimensions)
+   rho = dd["Dens"].d
+   rhok = np.fft.rfftn( rho )
+
+   kx = 2*np.pi*np.fft.fftfreq( N, dh )
+   ky = 2*np.pi*np.fft.fftfreq( N, dh )
+   kz = 2*np.pi*np.fft.rfftfreq( N, dh )
+
+   Pk3d = abs(rhok)**2
+   Pk_total = np.zeros(N//2+1)
+   count = np.zeros(N//2+1)
+
+   for i in range( N ):
+      for j in range( N ):
+         for k in range( N//2+1 ):
+            l = round((kx[i]**2+ky[j]**2+kz[k]**2)**0.5*BoxSize/(2.0*np.pi))
+            if (l < N//2+1):
+               Pk_total[l] = Pk_total[l] + Pk3d[i,j,k]
+               count[l] = count[l] + 1
+   Pk1d = np.divide(Pk_total, count, out=np.zeros_like(Pk_total), where=count!=0)
+   k3Pk = Pk1d*kz**3
+   d = 2*np.pi/kz[np.argmax(k3Pk)]
+   print('estimated granule diameter = %g\n'%(d))
+
+   if(GAMER_PS):
+      gamer_ps = np.loadtxt('../PowerSpec_%06d'%idx, skiprows=1, dtype=float)
+      gamer_k    = gamer_ps[:,0]
+      gamer_Pk1d = gamer_ps[:,1]
+      gamer_k3Pk = gamer_k**3*gamer_Pk1d
+
+   plt.figure()
+   plt.title("Dimensionless Power Spectrum")
+   plt.plot(kz[1:], k3Pk[1:]/max(k3Pk), label = 'numpy')
+   if(GAMER_PS):
+      plt.plot(gamer_k, gamer_k3Pk/max(gamer_k3Pk), label = 'gamer')
+   plt.xlabel('$k$')
+   plt.ylabel(r'$k^3P(k)$')
+   plt.yscale('log')
+   plt.legend(loc='lower right')
+   plt.savefig('fig_powerspectrum_dimensionless_%06d.png'%idx, dpi = 150, bbox_inches="tight")
+   plt.close()
+
+

example/test_problem/ELBDM/UniformGranule/python_script/plot_slice.py

@@ -0,0 +1,67 @@
+import yt
+import numpy as np
+import argparse
+import sys
+
+#-------------------------------------------------------------------------------------------------------------------------
+# load the command-line parameters
+parser = argparse.ArgumentParser( description='Get disk properties' )
+
+parser.add_argument( '-s', action='store', required=True,  type=int, dest='idx_start',
+                     help='first data index' )
+parser.add_argument( '-e', action='store', required=True,  type=int, dest='idx_end',
+                     help='last data index' )
+parser.add_argument( '-d', action='store', required=False, type=int, dest='didx',
+                     help='delta data index [%(default)d]', default=1 )
+
+args=parser.parse_args()
+
+idx_start   = args.idx_start
+idx_end     = args.idx_end
+didx        = args.didx
+
+# print command-line parameters
+print( '\nCommand-line arguments:' )
+print( '-------------------------------------------------------------------' )
+for t in range( len(sys.argv) ):
+   print( str(sys.argv[t]))
+print( '' )
+print( '-------------------------------------------------------------------\n' )
+
+field     = 'Dens'
+colormap  = 'algae'
+#colormap  = 'magma'
+dpi       = 300
+
+
+for idx in range(idx_start, idx_end+1, didx):
+   ds = yt.load("../Data_%06d"%idx)
+   dd = ds.all_data()
+   dens = np.array(dd["Dens"])
+   avedens = np.mean(dens)
+
+   plt = yt.SlicePlot( ds, 0, fields = field, center = 'c')
+   plt.set_zlim( field, avedens*1.0e-4, avedens*1.0e+1, dynamic_range=None)
+   plt.set_axes_unit( 'kpc' )
+   plt.set_unit( field, 'Msun/kpc**3')
+   plt.set_cmap( field, colormap )
+   plt.annotate_timestamp( time_unit='Myr', corner='upper_right', text_args={'color':'k'} )
+   plt.save( mpl_kwargs={"dpi":dpi} )
+
+   plt = yt.SlicePlot( ds, 1, fields = field, center = 'c')
+   plt.set_zlim( field, avedens*1.0e-4, avedens*1.0e+1, dynamic_range=None)
+   plt.set_axes_unit( 'kpc' )
+   plt.set_unit( field, 'Msun/kpc**3')
+   plt.set_cmap( field, colormap )
+   plt.annotate_timestamp( time_unit='Myr', corner='upper_right', text_args={'color':'k'} )
+   plt.save( mpl_kwargs={"dpi":dpi} )
+
+   plt = yt.SlicePlot( ds, 2, fields = field, center = 'c')
+   plt.set_zlim( field, avedens*1.0e-4, avedens*1.0e+1, dynamic_range=None)
+   plt.set_axes_unit( 'kpc' )
+   plt.set_unit( field, 'Msun/kpc**3')
+   plt.set_cmap( field, colormap )
+   plt.annotate_timestamp( time_unit='Myr', corner='upper_right', text_args={'color':'k'} )
+   plt.save( mpl_kwargs={"dpi":dpi} )
+
+

include/Typedef.h

@@ -93,7 +93,9 @@ const TestProbID_t
    TESTPROB_ELBDM_PLANE_WAVE                   = 1010,
    TESTPROB_ELBDM_PERTURBATION                 = 1011,
    TESTPROB_ELBDM_HALO_MERGER                  = 1012,
-   TESTPROB_ELBDM_DISK_HEATING                 = 1013;
+   TESTPROB_ELBDM_DISK_HEATING                 = 1013,
+   TESTPROB_ELBDM_UNIFORM_GRANULE              = 1014;
+
 
 // program initialization options
 typedef int OptInit_t;