diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index f2d173f..2efa9ba 100644 --- a/python_sample/icgen/borgicgen.py +++ b/python_sample/icgen/borgicgen.py @@ -136,8 +136,12 @@ def whitify(density, L, cosmo, supergenerate=1, func='HU_WIGGLES'): if supergenerate > 1: cond=np.isnan(density_hat_super) - x = np.random.randn(np.count_nonzero(cond),2)/np.sqrt(2.0) - density_hat_super[cond] = x[:,0] + 1j * x[:,1] + print np.where(np.isnan(density_hat_super))[0].size + Nz = np.count_nonzero(cond) + density_hat_super.real[cond] = np.random.randn(Nz) + density_hat_super.imag[cond] = np.random.randn(Nz) + density_hat_super[cond] /= np.sqrt(2.0) + print np.where(np.isnan(density_hat_super))[0].size # Now we have to fix the Nyquist plane hNs = Ns/2 @@ -145,6 +149,8 @@ def whitify(density, L, cosmo, supergenerate=1, func='HU_WIGGLES'): Nplane = nyquist.size nyquist.flat[:Nplane/2] = np.sqrt(2.0)*nyquist.flat[Nplane:Nplane/2:-1].conj() + print np.where(np.isnan(density_hat_super))[0].size + return np.fft.irfftn(density_hat_super)*Ns**1.5 diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py index 2287340..0ce037e 100644 --- a/python_sample/icgen/test_ic_from_borg.py +++ b/python_sample/icgen/test_ic_from_borg.py @@ -10,26 +10,34 @@ cosmo['omega_k_0'] = 0 cosmo['omega_B_0']=0.049 cosmo['SIGMA8']=0.8344 cosmo['ns']=0.9624 -N0=128 +N0=256 doSimulation=True +simShift=True snap_id=int(sys.argv[1]) astart=1/100. if doSimulation: - s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim/", snap_id, doublePrecision=True) + s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim2/", snap_id, doublePrecision=True) astart=s.getTime() L = s.getBoxsize() - dsim = ct.cicParticles(s.getPositions(), L, N0) + p = s.getPositions() + Nsim = int( np.round( p[0].size**(1./3)) ) + print("Nsim = %d" % Nsim) + + if simShift: + p = [(q-0.5*L/Nsim)%L for q in p] + + dsim = ct.cicParticles(p[::-1], L, N0) dsim /= np.average(np.average(np.average(dsim, axis=0), axis=0), axis=0) dsim -= 1 dsim_hat = np.fft.rfftn(dsim)*(L/N0)**3 Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, range=(0,1.), bins=150) -pos,_,density,N,L,_,_ = bic.run_generation("initial_density_1380.dat", 0.001, astart, cosmo, supersample=2, do_lpt2=True) +pos,_,density,N,L,_,_ = bic.run_generation("initial_density_2588.dat", 0.001, astart, cosmo, supersample=2, do_lpt2=True) dcic = ct.cicParticles(pos, L, N0) dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0) @@ -50,7 +58,7 @@ Pref, bref = bic.compute_ref_power(L, N0, cosmo, range=(0,1.), bins=150) Pcic /= D1_0**2 Pdens /= D1_0**2 -borg_evolved = ct.read_borg_vol("final_density_1380.dat") -dborg_hat = np.fft.rfftn(borg_evolved.density)*L**3/borg_evolved.density.size +#borg_evolved = ct.read_borg_vol("final_density_1380.dat") +#dborg_hat = np.fft.rfftn(borg_evolved.density)*L**3/borg_evolved.density.size -Pborg, bborg = bic.bin_power(np.abs(dborg_hat)**2/L**3, L, range=(0,1.),bins=150) +#Pborg, bborg = bic.bin_power(np.abs(dborg_hat)**2/L**3, L, range=(0,1.),bins=150)