Added pixel shifting in ICgen

python_sample/icgen/borgadaptor.py

@@ -6,10 +6,11 @@ def fourier_analysis(borg_vol):
     
   return np.fft.rfftn(borg_vol.density)*(L/N)**3, L, N
 
-def half_pixel_shift(borg):
+def half_pixel_shift(borg, doshift=False):
 
     dhat,L,N = fourier_analysis(borg)
-    return dhat, L
+    if not doshift:
+      return dhat, L
 
     ik = np.fft.fftfreq(N,d=L/N)*2*np.pi
     phi = 0.5*L/N*(ik[:,None,None]+ik[None,:,None]+ik[None,None,:(N/2+1)])

python_sample/icgen/borgicgen.py

@@ -50,7 +50,7 @@ def compute_ref_power(L, N, cosmo, bins=10, range=(0,1), func='HU_WIGGLES'):
 
   return bin_power(p.compute(k)*cosmo['h']**3, L, bins=bins, range=range)
 
-def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True):
+def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True, shiftPixel=False):
     """ Generate particles and velocities from a BORG snapshot. Returns a tuple of 
     (positions,velocities,N,BoxSize,scale_factor)."""
 
@@ -59,7 +59,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True)
 
     cgrowth = CosmoGrowth(**cosmo)
 
-    density_hat, L = ba.half_pixel_shift(borg_vol)
+    density_hat, L = ba.half_pixel_shift(borg_vol, doshift=shiftPixel)
 
     lpt = LagrangianPerturbation(density_hat, L, fourier=True, supersample=supersample)
 
@@ -89,16 +89,16 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True)
 
     print("velmul=%lg" % (cosmo['h']*velmul))
 
-    density = np.fft.irfftn(density_hat*D1_0)*(N/L)**3
+    density = np.fft.irfftn(lpt.dhat*D1_0)*(supersample*N/L)**3
 
-    return posx,vel,density,N,L,a_ic
+    return posx,vel,density,N*supersample,L,a_ic
     
 def write_icfiles(*generated_ic, **cosmo):
     """Write the initial conditions from the tuple returned by run_generation"""
     posx,vel,density,N,L,a_ic = generated_ic
 
     ct.simpleWriteGadget("borg.gad", posx, velocities=vel, boxsize=L, Hubble=cosmo['h'], Omega_M=cosmo['omega_M_0'], time=a_ic) 
-    for i,c in enumerate(['x','y','z']):
+    for i,c in enumerate(["x","y","z"]):
         ct.writeGrafic("ic_velc%s" % c, vel[i].reshape((N,N,N)).transpose(), L, a_ic, **cosmo)
 
     ct.writeGrafic("ic_deltab", density, L, a_ic, **cosmo)

python_sample/icgen/gen_ic_from_borg.py

@@ -8,9 +8,9 @@ cosmo['omega_k_0'] = 0
 cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
 
-zstart=200
+zstart=50
 astart=1/(1.+zstart)
-snapid=1
+halfPixelShift=True
 
 if __name__=="__main__":
-    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, do_lpt2=False), **cosmo)
+    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=2, shiftPixel=halfPixelShift, do_lpt2=False), **cosmo)

python_sample/icgen/test_ic_from_borg.py

@@ -1,25 +1,27 @@
 import numpy as np
 import cosmotool as ct
 import borgicgen as bic
+import sys
 
 cosmo={'omega_M_0':0.3175, 'h':0.6711}
 cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
 cosmo['omega_k_0'] = 0
 cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
+N0=128
 
-snap_id=11
+snap_id=int(sys.argv[1])
 
 s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim/", snap_id, doublePrecision=True)
 astart=s.getTime()
 
-pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=1, do_lpt2=True)
+pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=2, do_lpt2=True)
 
-dcic = ct.cicParticles(pos, L, N)
+dcic = ct.cicParticles(pos, L, N0)
 dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
 dcic -= 1
 
-dsim = ct.cicParticles(s.getPositions(), L, N)
+dsim = ct.cicParticles(s.getPositions(), L, N0)
 dsim /= np.average(np.average(np.average(dsim, axis=0), axis=0), axis=0)
 dsim -= 1
 
@@ -27,8 +29,8 @@ dsim -= 1
 dcic_hat = np.fft.rfftn(dcic)*(L/N)**3
 dsim_hat = np.fft.rfftn(dsim)*(L/N)**3
 
-Pcic, bcic = bic.bin_power(np.abs(dcic_hat)**2/L**3, L, bins=50)
-Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, bins=50)
+Pcic, bcic = bic.bin_power(np.abs(dcic_hat)**2/L**3, L, range=(0,1.), bins=150)
+Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, range=(0,1.), bins=150)
 
 borg_evolved = ct.read_borg_vol("final_density_1380.dat")