NGP fix

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))
+    ix = int(round(rx))
-    iy = int(floor(ry))
+    iy = int(round(ry))
-    iz = int(floor(rz))
+    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))
+    ix = int(round(rx))
-    iy = int(floor(ry))
+    iy = int(round(ry))
-    iz = int(floor(rz))
+    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 not myngp:
-        if interp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0:
+        if periodic:
-          raise ierror
+          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:
+        if periodic:
-          raise ierror
+          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)

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
+    if zero_fill:
-    Nz = np.count_nonzero(cond)
+      density_hat_super[cond] = 0
-    density_hat_super.real[cond] = np.random.randn(Nz)
+    else:
-    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
-    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
+      hNs = Ns/2
-    nyquist = density_hat_super[:, :, hNs]
+      nyquist = density_hat_super[:, :, hNs]
-    Nplane = nyquist.size
+      Nplane = nyquist.size
-    nyquist.flat[:Nplane/2] = np.sqrt(2.0)*nyquist.flat[Nplane:Nplane/2:-1].conj()
+      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.get('supergenerate', 1)
-    supergenerate=kwargs['supergenerate']
+  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']))

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)