diff --git a/python/_cosmotool.pyx b/python/_cosmotool.pyx index 706285c..0181b47 100644 --- a/python/_cosmotool.pyx +++ b/python/_cosmotool.pyx @@ -6,6 +6,47 @@ from cpython cimport PyObject, Py_INCREF np.import_array() +cdef extern from "cosmopower.hpp" namespace "CosmoTool": + + cdef enum CosmoFunction "CosmoTool::CosmoPower::CosmoFunction": + POWER_EFSTATHIOU "CosmoTool::CosmoPower::POWER_EFSTATHIOU", + HU_WIGGLES "CosmoTool::CosmoPower::HU_WIGGLES", + HU_BARYON "CosmoTool::CosmoPower::HU_BARYON", + OLD_POWERSPECTRUM, + POWER_BARDEEN "CosmoTool::CosmoPower::POWER_BARDEEN", + POWER_SUGIYAMA "CosmoTool::CosmoPower::POWER_SUGIYAMA", + POWER_BDM, + POWER_TEST + + cdef cppclass CosmoPower: + double n + double K0 + double V_LG_CMB + + double CMB_VECTOR[3] + double h + double SIGMA8 + double OMEGA_B + double OMEGA_C + double omega_B + double omega_C + double Theta_27 + double OMEGA_0 + double Omega + double beta + double OmegaEff + double Gamma0 + double normPower + + + CosmoPower() + void setFunction(CosmoFunction) + void updateCosmology() + void updatePhysicalCosmology() + void normalize() + void setNormalization(double) + double power(double) + cdef extern from "loadSimu.hpp" namespace "CosmoTool": @@ -302,4 +343,60 @@ def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision = return _PySimulationAdaptor(wrap_simudata(data, flags)) -def loadAllRamses(str basepath, int + +cdef class CosmologyPower: + + cdef CosmoPower power + + def __init__(self,**cosmo): + + self.power = CosmoPower() + self.power.OMEGA_B = cosmo['omega_B_0'] + self.power.OMEGA_C = cosmo['omega_M_0']-cosmo['omega_B_0'] + self.power.h = cosmo['h'] + if 'ns' in cosmo: + self.power.n = cosmo['ns'] + + assert self.power.OMEGA_C > 0 + + self.power.updateCosmology() + + def normalize(self,s8): + self.power.SIGMA8 = s8 + self.power.normalize() + + + def setFunction(self,funcname): + cdef CosmoFunction f + + f = POWER_EFSTATHIOU + + if funcname=='EFSTATHIOU': + f = POWER_EFSTATHIOU + elif funcname=='HU_WIGGLES': + f = HU_WIGGLES + elif funcname=='HU_BARYON': + f = HU_BARYON + elif funcname=='BARDEEN': + f = POWER_BARDEEN + elif funcname=='SUGIYAMA': + f = POWER_SUGIYAMA + + self.power.setFunction(f) + + cdef double _compute(self, double k): + k *= self.power.h + return self.power.power(k) + + def compute(self, k): + cdef np.ndarray out + cdef double kval + cdef tuple i + + if isinstance(k, np.ndarray): + out = np.empty(k.shape, dtype=np.float64) + for i,kval in np.ndenumerate(k): + out[i] = self._compute(kval) + return out + else: + return self._compute(k) diff --git a/python/cosmotool/__init__.py b/python/cosmotool/__init__.py index de5d9fa..eebdb9c 100644 --- a/python/cosmotool/__init__.py +++ b/python/cosmotool/__init__.py @@ -1,5 +1,5 @@ from _cosmotool import * -from grafic import writeGrafic, writeWhitePhase +from grafic import writeGrafic, writeWhitePhase, readGrafic from borg import read_borg_vol @@ -13,9 +13,9 @@ class SimulationBare(PySimulationBase): raise TypeError("Simulation object to mirror must be a PySimulationBase") s = args[0] - self.positions = s.getPositions() - self.velocities = s.getVelocities() - self.identifiers = s.getIdentifiers() + self.positions = [q.copy() for q in s.getPositions()] if s.getPositions() is not None else None + self.velocities = [q.copy() for q in s.getVelocities()] if s.getVelocities() is not None else None + self.identifiers = s.getIdentifiers().copy() if s.getIdentifiers() is not None else None self.boxsize = s.getBoxsize() self.time = s.getTime() self.Hubble = s.getHubble() @@ -28,7 +28,7 @@ class SimulationBare(PySimulationBase): def _safe_merge(a, b): if b: if a: - a = [np.append(q, r) for q,r in zip(a,other.positions)] + a = [np.append(q, r) for q,r in zip(a,b)] else: a = b return a @@ -42,15 +42,15 @@ class SimulationBare(PySimulationBase): return a - assert self.time == other.time - assert self.Hubble == other.Hubble - assert self.boxsize == other.boxsize - assert self.Omega_M == other.Omega_M - assert self.Omega_Lambda == other.Omega_Lambda + assert self.time == other.getTime() + assert self.Hubble == other.getHubble() + assert self.boxsize == other.getBoxsize() + assert self.Omega_M == other.getOmega_M() + assert self.Omega_Lambda == other.getOmega_Lambda() - self.positions = _safe_merge(self.positions, other.positions) - self.velocities = _safe_merge(self.velocities, other.velocities) - self.identifiers = _safe_merge0(self.idenfiers, other.identifiers) + self.positions = _safe_merge(self.positions, other.getPositions()) + self.velocities = _safe_merge(self.velocities, other.getVelocities()) + self.identifiers = _safe_merge0(self.identifiers, other.getIdentifiers()) def getPositions(self): return self.positions @@ -106,7 +106,7 @@ def loadRamsesAll(basepath, snapshot_id, **kwargs): cpu_id = 0 output = None while True: - s = loadRamses("%s/output_%05d" % (basepath,snapshot_id), cpu_id, **kwargs) + s = loadRamses("%s/output_%05d" % (basepath,snapshot_id), snapshot_id, cpu_id, **kwargs) if s == None: break if output == None: @@ -115,3 +115,5 @@ def loadRamsesAll(basepath, snapshot_id, **kwargs): output.merge(s) cpu_id += 1 + + return output diff --git a/python/cosmotool/grafic.py b/python/cosmotool/grafic.py index f863167..13d93ac 100644 --- a/python/cosmotool/grafic.py +++ b/python/cosmotool/grafic.py @@ -1,6 +1,32 @@ import struct import numpy as np +def readGrafic(filename): + + with file(filename, mode="rb") as f: + p = struct.unpack("IIIIffffffffI", f.read(4*11 + 2*4)) + checkPoint0, Nx, Ny, Nz, delta, _, _, _, scalefac, omega0, omegalambda0, h, checkPoint1 = p + if checkPoint0 != checkPoint1 or checkPoint0 != 4*11: + raise ValueError("Invalid unformatted access") + + a = np.empty((Nx,Ny,Nz), dtype=np.float32) + + BoxSize = delta * Nx * h + + checkPoint = 4*Ny*Nz + for i in xrange(Nx): + checkPoint = struct.unpack("I", f.read(4))[0] + if checkPoint != 4*Ny*Nz: + raise ValueError("Invalid unformatted access") + + a[i, :, :] = np.fromfile(f, dtype=np.float32, count=Ny*Nz).reshape((Ny, Nz)) + + checkPoint = struct.unpack("I", f.read(4))[0] + if checkPoint != 4*Ny*Nz: + raise ValueError("Invalid unformatted access") + + return a, BoxSize, scalefac, omega0, omegalambda0, h + def writeGrafic(filename, field, BoxSize, scalefac, **cosmo): with file(filename, mode="wb") as f: diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index 865eeb6..4c279db 100644 --- a/python_sample/icgen/borgicgen.py +++ b/python_sample/icgen/borgicgen.py @@ -15,6 +15,40 @@ def gen_posgrid(N, L): return x.flatten(), y.flatten(), z.flatten() +def bin_power(P, L, bins=20, range=(0,1.)): + + N = P.shape[0] + ik = np.fft.fftfreq(N, d=L/N)*2*np.pi + + k = np.sqrt(ik[:,None,None]**2 + ik[None,:,None]**2 + ik[None,None,:(N/2+1)]**2) + + H,b = np.histogram(k, bins=bins, range=range) + Hw,b = np.histogram(k, bins=bins, weights=P, range=range) + + return Hw/H, 0.5*(b[1:]+b[0:bins]) + +def compute_power_from_borg(input_borg, a_borg, cosmo, bins=10, range=(0,1)): + borg_vol = ct.read_borg_vol(input_borg) + N = borg_vol.density.shape[0] + + cgrowth = CosmoGrowth(**cosmo) + D1 = cgrowth.D(1) + D1_0 = D1/cgrowth.D(a_borg) + + density_hat, L = ba.half_pixel_shift(borg_vol) + + return bin_power(D1_0**2*np.abs(density_hat)**2/L**3, L, bins=bins, range=range) + +def compute_ref_power(L, N, cosmo, bins=10, range=(0,1), func='HU_WIGGLES'): + ik = np.fft.fftfreq(N, d=L/N)*2*np.pi + + k = np.sqrt(ik[:,None,None]**2 + ik[None,:,None]**2 + ik[None,None,:(N/2+1)]**2) + p = ct.CosmologyPower(**cosmo) + p.setFunction(func) + p.normalize(cosmo['SIGMA8']) + + return bin_power(p.compute(k)*cosmo['h']**3, L, bins=bins, range=range) + def run_generation(input_borg, a_borg, a_ic, **cosmo): """ Generate particles and velocities from a BORG snapshot. Returns a tuple of (positions,velocities,N,BoxSize,scale_factor).""" @@ -26,7 +60,6 @@ def run_generation(input_borg, a_borg, a_ic, **cosmo): density_hat, L = ba.half_pixel_shift(borg_vol) - lpt = LagrangianPerturbation(density_hat, L, fourier=True) # Generate grid diff --git a/python_sample/icgen/gen_ic_from_borg.py b/python_sample/icgen/gen_ic_from_borg.py index 0c08630..ab45742 100644 --- a/python_sample/icgen/gen_ic_from_borg.py +++ b/python_sample/icgen/gen_ic_from_borg.py @@ -1,11 +1,14 @@ +import cosmotool as ct import borgicgen as bic 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 -zstart=50 +zstart=0 astart=1/(1.+zstart) -if __name__=="__main__": - bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo), **cosmo) +#if __name__=="__main__": +# bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo), **cosmo) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 966d0b0..dc7ee5a 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -88,6 +88,8 @@ if (Boost_FOUND) target_link_libraries(simple3DFilter ${tolink}) endif (Boost_FOUND) -add_executable(gadgetToDensity gadgetToDensity.cpp) -target_link_libraries(gadgetToDensity ${tolink}) +IF (ENABLE_OPENMP) + add_executable(gadgetToDensity gadgetToDensity.cpp) + target_link_libraries(gadgetToDensity ${tolink}) +ENDIF (ENABLE_OPENMP) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4feb000..9737765 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -11,9 +11,12 @@ SET(CosmoTool_SRCS miniargs.cpp growthFactor.cpp cosmopower.cpp - cic.cpp ) +IF (ENABLE_OPENMP) + set(CosmoTool_SRCS ${CosmoTool_SRCS} cic.cpp) +ENDIF (ENABLE_OPENMP) + IF(FOUND_NETCDF3) SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc3.cpp) ELSE(FOUND_NETCDF3) diff --git a/src/cic.cpp b/src/cic.cpp index ba6bac5..bfec994 100644 --- a/src/cic.cpp +++ b/src/cic.cpp @@ -32,7 +32,7 @@ same conditions as regards security. The fact that you are presently reading this means that you have had knowledge of the CeCILL license and that you accept its terms. +*/ - +#include #include #include #include @@ -60,89 +60,29 @@ void CICFilter::resetMesh() void CICFilter::putParticles(CICParticles *particles, uint32_t N) { -#if 0 - uint32_t numCorners = 1 << NUMDIMS; + int threadUsed = omp_get_max_threads(); + double *threadedDensity[threadUsed]; + bool threadActivated[threadUsed]; + uint32_t tUsedMin[threadUsed], tUsedMax[threadUsed]; - for (uint32_t i = 0; i < N; i++) - { - Coordinates xyz; - int32_t ixyz[NUMDIMS]; - int32_t rxyz[NUMDIMS]; - CICType alpha[NUMDIMS]; - CICType beta[NUMDIMS]; - for (int j = 0; j < NUMDIMS; j++) - { - xyz[j] = (particles[i].coords[j] / spatialLen * szGrid); - ixyz[j] = (int32_t)floor(xyz[j] - 0.5); - beta[j] = xyz[j] - ixyz[j] - 0.5; - alpha[j] = 1 - beta[j]; - if (ixyz[j] < 0) - ixyz[j] = szGrid-1; - } + for (int t = 0; t < threadUsed; t++) + { + threadedDensity[t] = new double[totalSize]; + std::fill(threadedDensity[t], threadedDensity[t]+totalSize, 0); + } - CICType tot_mass = 0; - for (int j = 0; j < numCorners; j++) - { - CICType rel_mass = 1; - uint32_t idx = 0; - uint32_t mul = 1; - uint32_t mul2 = 1; + std::fill(threadActivated, threadActivated+threadUsed, false); + std::fill(tUsedMin, tUsedMin+threadUsed, totalSize); + std::fill(tUsedMax, tUsedMax+threadUsed, 0); - for (int k = 0; k < NUMDIMS; k++) - { - uint32_t ipos = ((j & mul2) != 0); - - if (ipos == 1) - { - rel_mass *= beta[k]; - } - else - { - rel_mass *= alpha[k]; - } - - rxyz[k] = ixyz[k] + ipos; - - if (rxyz[k] >= szGrid) - idx += (rxyz[k] - szGrid) * mul; - else - idx += rxyz[k] * mul; - - mul2 *= 2; - mul *= szGrid; - } - - assert(rel_mass > 0); - assert(rel_mass < 1); - assert(idx < totalSize); - densityGrid[idx] += rel_mass * particles[i].mass; - tot_mass += rel_mass; - } - assert(tot_mass < 1.1); - assert(tot_mass > 0.9); - } -#endif -#if 0 - for (uint32_t i = 0; i < N; i++) - { - Coordinates xyz; - int32_t ixyz[NUMDIMS]; - for (int j = 0; j < NUMDIMS; j++) - { - xyz[j] = (particles[i].coords[j] / spatialLen * szGrid); - ixyz[j] = (int32_t)round(xyz[j] - 0.5); - if (ixyz[j] < 0) - ixyz[j] = szGrid-1; - else if (ixyz[j] >= szGrid) - ixyz[j] = 0; - } - - uint32_t idx = ixyz[0] + ixyz[1] * szGrid + ixyz[2] * szGrid * szGrid; - densityGrid[idx] += particles[i].mass; - } - -#endif +#pragma omp parallel + { + int thisThread = omp_get_thread_num(); + double *dg = threadedDensity[thisThread]; + threadActivated[thisThread] = true; + +#pragma omp for schedule(static) for (uint32_t i = 0; i < N; i++) { CICType x, y, z; @@ -193,44 +133,57 @@ void CICFilter::putParticles(CICParticles *particles, uint32_t N) // 000 idx = ix + (iy + iz * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * beta_x * beta_y * beta_z; // 100 idx = ix2 + (iy + iz * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * alpha_x * beta_y * beta_z; // 010 idx = ix + (iy2 + iz * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * beta_x * alpha_y * beta_z; // 110 idx = ix2 + (iy2 + iz * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * alpha_x * alpha_y * beta_z; // 001 idx = ix + (iy + iz2 * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * beta_x * beta_y * alpha_z; // 101 idx = ix2 + (iy + iz2 * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * alpha_x * beta_y * alpha_z; // 011 idx = ix + (iy2 + iz2 * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * beta_x * alpha_y * alpha_z; // 111 idx = ix2 + (iy2 + iz2 * szGrid) * szGrid; - densityGrid[idx] += + dg[idx] += mass * alpha_x * alpha_y * alpha_z; + + tUsedMin[thisThread] = std::min(tUsedMin[thisThread], idx); + tUsedMax[thisThread] = std::max(tUsedMax[thisThread], idx); } + } + + for (int t = 0; t < threadUsed; t++) + { + if (!threadActivated[t]) + continue; +#pragma omp parallel for schedule(static) + for (long p = tUsedMin[t]; p < tUsedMax[t]; p++) + densityGrid[p] += threadedDensity[t][p]; + } } void CICFilter::getDensityField(CICType*& field, uint32_t& res) @@ -238,3 +191,4 @@ void CICFilter::getDensityField(CICType*& field, uint32_t& res) field = densityGrid; res = totalSize; } +