From 7258dd20cb0bb5cbbb03454685711bc4371787d6 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 4 Jun 2014 09:13:34 +0200 Subject: [PATCH] Added pixel shifting in ICgen --- python_sample/icgen/borgadaptor.py | 5 +++-- python_sample/icgen/borgicgen.py | 10 +++++----- python_sample/icgen/gen_ic_from_borg.py | 6 +++--- python_sample/icgen/test_ic_from_borg.py | 14 ++++++++------ 4 files changed, 19 insertions(+), 16 deletions(-) diff --git a/python_sample/icgen/borgadaptor.py b/python_sample/icgen/borgadaptor.py index 6d906a5..2ee6018 100644 --- a/python_sample/icgen/borgadaptor.py +++ b/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)]) diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index df6d151..79ccb0c 100644 --- a/python_sample/icgen/borgicgen.py +++ b/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) diff --git a/python_sample/icgen/gen_ic_from_borg.py b/python_sample/icgen/gen_ic_from_borg.py index de73e4f..9030a99 100644 --- a/python_sample/icgen/gen_ic_from_borg.py +++ b/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) diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py index dcef2a6..cfedb09 100644 --- a/python_sample/icgen/test_ic_from_borg.py +++ b/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")