From 087ab69ff351ad0276ab9d6f5012821ea6ed8344 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Fri, 18 Jul 2014 15:54:12 +0200 Subject: [PATCH] NGP fix --- python/_project.pyx | 34 +++++++++++++++---------- python_sample/icgen/borgicgen.py | 32 +++++++++++++---------- python_sample/icgen/gen_ic_from_borg.py | 7 ++--- 3 files changed, 43 insertions(+), 30 deletions(-) diff --git a/python/_project.pyx b/python/_project.pyx index 1261c65..4fea24f 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -1,7 +1,7 @@ from cpython cimport bool from cython cimport view from cython.parallel import prange, parallel -from libc.math cimport sin, cos, abs, floor, sqrt +from libc.math cimport sin, cos, abs, floor, round, sqrt import numpy as np cimport numpy as npx cimport cython @@ -99,9 +99,9 @@ cdef int ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, ry = (inv_delta*y) rz = (inv_delta*z) - ix = int(floor(rx)) - iy = int(floor(ry)) - iz = int(floor(rz)) + ix = int(round(rx)) + iy = int(round(ry)) + iz = int(round(rz)) ix = ix%Ngrid @@ -130,9 +130,9 @@ cdef int ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y, ry = (inv_delta*y) rz = (inv_delta*z) - ix = int(floor(rx)) - iy = int(floor(ry)) - iz = int(floor(rz)) + ix = int(round(rx)) + iy = int(round(ry)) + iz = int(round(rz)) if (ix < 0 or ix >= Ngrid): return -1 @@ -250,7 +250,7 @@ def interp3d(x not None, y not None, in_slice = d Nelt = ax.size with nogil: - if myngp: + if not myngp: if myperiodic: for i in prange(Nelt): if interp3d_INTERNAL_periodic(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i]) < 0: @@ -274,12 +274,20 @@ def interp3d(x not None, y not None, raise ierror return out else: - if periodic: - if interp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: - raise ierror + if not myngp: + if periodic: + if interp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: + raise ierror + else: + if interp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: + raise ierror else: - if interp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: - raise ierror + if periodic: + if ngp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: + raise ierror + else: + if ngp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0: + raise ierror return retval @cython.boundscheck(False) diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index 423d7e8..5e5c029 100644 --- a/python_sample/icgen/borgicgen.py +++ b/python_sample/icgen/borgicgen.py @@ -103,7 +103,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True, @ct.timeit_quiet -def whitify(density, L, cosmo, supergenerate=1, func='HU_WIGGLES'): +def whitify(density, L, cosmo, supergenerate=1, zero_fill=False, func='HU_WIGGLES'): N = density.shape[0] p = ct.CosmologyPower(**cosmo) @@ -141,18 +141,22 @@ def whitify(density, L, cosmo, supergenerate=1, func='HU_WIGGLES'): if supergenerate > 1: cond=np.isnan(density_hat_super) - 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 + if zero_fill: + density_hat_super[cond] = 0 + else: + + 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 - nyquist = density_hat_super[:, :, hNs] - Nplane = nyquist.size - nyquist.flat[:Nplane/2] = np.sqrt(2.0)*nyquist.flat[Nplane:Nplane/2:-1].conj() + hNs = Ns/2 + nyquist = density_hat_super[:, :, hNs] + 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 @@ -166,8 +170,8 @@ def write_icfiles(*generated_ic, **kwargs): """Write the initial conditions from the tuple returned by run_generation""" supergenerate=1 - if 'supergenerate' in kwargs: - supergenerate=kwargs['supergenerate'] + supergenerate=kwargs.get('supergenerate', 1) + zero_fill=kwargs.get('zero_fill', False) posx,vel,density,N,L,a_ic,cosmo = generated_ic ct.simpleWriteGadget("Data/borg.gad", posx, velocities=vel, boxsize=L, Hubble=cosmo['h'], Omega_M=cosmo['omega_M_0'], time=a_ic) @@ -176,7 +180,7 @@ def write_icfiles(*generated_ic, **kwargs): ct.writeGrafic("Data/ic_deltab", density, L, a_ic, **cosmo) - ct.writeWhitePhase("Data/white.dat", whitify(density, L, cosmo, supergenerate=supergenerate)) + ct.writeWhitePhase("Data/white.dat", whitify(density, L, cosmo, supergenerate=supergenerate,zero_fill=zero_fill)) with file("Data/white_params", mode="w") as f: f.write("4\n%lg, %lg, %lg\n" % (cosmo['omega_M_0'], cosmo['omega_lambda_0'], 100*cosmo['h'])) diff --git a/python_sample/icgen/gen_ic_from_borg.py b/python_sample/icgen/gen_ic_from_borg.py index b2a6fe5..9249725 100644 --- a/python_sample/icgen/gen_ic_from_borg.py +++ b/python_sample/icgen/gen_ic_from_borg.py @@ -9,10 +9,11 @@ cosmo['omega_B_0']=0.049 cosmo['SIGMA8']=0.8344 cosmo['ns']=0.9624 -supergen=1 +supergen=2 zstart=50 -astart=1/(1.+zstart) +astart=1/69.#1/(1.+zstart) halfPixelShift=False +zero_fill=True if __name__=="__main__": - bic.write_icfiles(*bic.run_generation("initial_density_988.dat", 0.001, astart, cosmo, supersample=1, shiftPixel=halfPixelShift, do_lpt2=False), supergenerate=supergen) + bic.write_icfiles(*bic.run_generation("initial_density_1380.dat", 0.001, astart, cosmo, supersample=1, shiftPixel=halfPixelShift, do_lpt2=False), supergenerate=supergen, zero_fill=zero_fill)