guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed integrator to stop at exactly the right distance

Browse Source
This commit is contained in:
Guilhem Lavaux 2015-07-14 12:21:20 +02:00
parent b668f9947d
commit 21507e3013
2 changed files with 23 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
python/_project.pyx
View File

@ -600,14 +600,11 @@ def tophat_fourier(x not None):
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) nogil: cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t alpha_max) nogil:
cdef DTYPE_t alpha_max
cdef DTYPE_t tmp_a cdef DTYPE_t tmp_a
cdef DTYPE_t v[3] cdef DTYPE_t v[3]
cdef int i, j cdef int i, j
alpha_max = 10.0 # A big number
for i in xrange(3): for i in xrange(3):
if u[i] == 0.: if u[i] == 0.:
continue continue
@ -630,14 +627,11 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) nogil:
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]) nogil: cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], DTYPE_t alpha_max) nogil:
cdef DTYPE_t alpha_max
cdef DTYPE_t I, tmp_a cdef DTYPE_t I, tmp_a
cdef DTYPE_t v[3], term[4] cdef DTYPE_t v[3], term[4]
cdef int i, j, q cdef int i, j, q
alpha_max = 10.0 # A big number
j = 0 j = 0
for i in range(3): for i in range(3):
if u[i] == 0.: if u[i] == 0.:
@ -651,7 +645,7 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
if tmp_a < alpha_max: if tmp_a < alpha_max:
alpha_max = tmp_a alpha_max = tmp_a
j = i j = i
I = compute_projection(vertex_value, u, u0, alpha_max) I = compute_projection(vertex_value, u, u0, alpha_max)
for i in xrange(3): for i in xrange(3):
@ -665,16 +659,16 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
@cython.boundscheck(False) @cython.boundscheck(False)
cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil: DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil:
cdef DTYPE_t d cdef DTYPE_t d
d = density[iu0[0], iu0[1], iu0[2]] d = density[iu0[0], iu0[1], iu0[2]]
return cube_integral(u, u0, jumper)*d return cube_integral(u, u0, jumper, alpha_max)*d
@cython.boundscheck(False) @cython.boundscheck(False)
cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil: DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil:
cdef DTYPE_t vertex_value[8] cdef DTYPE_t vertex_value[8]
cdef DTYPE_t d cdef DTYPE_t d
cdef int a[3][2], i cdef int a[3][2], i
@ -682,7 +676,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
for i in xrange(3): for i in xrange(3):
a[i][0] = iu0[i] a[i][0] = iu0[i]
a[i][1] = iu0[i]+1 a[i][1] = iu0[i]+1
vertex_value[0 + 2*0 + 4*0] = density[a[0][0], a[1][0], a[2][0]] vertex_value[0 + 2*0 + 4*0] = density[a[0][0], a[1][0], a[2][0]]
vertex_value[1 + 2*0 + 4*0] = density[a[0][1], a[1][0], a[2][0]] vertex_value[1 + 2*0 + 4*0] = density[a[0][1], a[1][0], a[2][0]]
vertex_value[0 + 2*1 + 4*0] = density[a[0][0], a[1][1], a[2][0]] vertex_value[0 + 2*1 + 4*0] = density[a[0][0], a[1][1], a[2][0]]
@ -693,7 +687,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
vertex_value[0 + 2*1 + 4*1] = density[a[0][0], a[1][1], a[2][1]] vertex_value[0 + 2*1 + 4*1] = density[a[0][0], a[1][1], a[2][1]]
vertex_value[1 + 2*1 + 4*1] = density[a[0][1], a[1][1], a[2][1]] vertex_value[1 + 2*1 + 4*1] = density[a[0][1], a[1][1], a[2][1]]
return cube_integral_trilin(u, u0, jumper, vertex_value) return cube_integral_trilin(u, u0, jumper, vertex_value, alpha_max)
@ -714,7 +708,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
cdef int jumper[1] cdef int jumper[1]
cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:], cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil
if integrator_id == 0: if integrator_id == 0:
integrator = integrator0 integrator = integrator0
@ -741,16 +735,17 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
raise RuntimeError("u0[%d] = %g !" % (i,u0[i])) raise RuntimeError("u0[%d] = %g !" % (i,u0[i]))
completed = 0 completed = 0
if ((iu0[0] >= N) or (iu0[0] <= 0) or if ((iu0[0] >= N-1) or (iu0[0] <= 0) or
(iu0[1] >= N) or (iu0[1] <= 0) or (iu0[1] >= N-1) or (iu0[1] <= 0) or
(iu0[2] >= N) or (iu0[2] <= 0)): (iu0[2] >= N-1) or (iu0[2] <= 0)):
completed = 1 completed = 1
I0 = 0 I0 = 0
jumper[0] = 0 jumper[0] = 0
dist2 = 0
while completed == 0: while completed == 0:
I0 += integrator(density, u, u0, u_delta, iu0, jumper) I0 += integrator(density, u, u0, u_delta, iu0, jumper, max_distance-sqrt(dist2))
if u[jumper[0]] < 0: if u[jumper[0]] < 0:
iu0[jumper[0]] -= 1 iu0[jumper[0]] -= 1
@ -760,9 +755,9 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
u0[jumper[0]] = 0 u0[jumper[0]] = 0
if ((iu0[0] >= N) or (iu0[0] <= 0) or if ((iu0[0] >= N-1) or (iu0[0] <= 0) or
(iu0[1] >= N) or (iu0[1] <= 0) or (iu0[1] >= N-1) or (iu0[1] <= 0) or
(iu0[2] >= N) or (iu0[2] <= 0)): (iu0[2] >= N-1) or (iu0[2] <= 0)):
completed = 1 completed = 1
else: else:
dist2 = 0 dist2 = 0
@ -788,7 +783,7 @@ def line_of_sight_projection(DTYPE_t[:,:,:] density not None,
u[1] = a_u[1] u[1] = a_u[1]
u[2] = a_u[2] u[2] = a_u[2]
C_line_of_sight_projection(density, return C_line_of_sight_projection(density,
u, u,
min_distance, min_distance,
max_distance, shifter, integrator_id) max_distance, shifter, integrator_id)
@ -846,10 +841,13 @@ def spherical_projection(int Nside,
shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t)) shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
shifter[:] = 0 shifter[:] = 0
print("allocating map")
outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE) outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE)
print("initializing views")
outm = outm_array outm = outm_array
density_view = density density_view = density
print("progress?")
if progress != 0: if progress != 0:
p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start() p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start()
@ -864,7 +862,7 @@ def spherical_projection(int Nside,
theta,phi = hp.pix2ang(Nside, np.arange(N0)) theta,phi = hp.pix2ang(Nside, np.arange(N0))
with nogil, parallel(): with nogil, parallel():
tid = smp_get_thread_id() tid = smp_get_thread_id()
for i in prange(N0): for i in prange(N0,schedule='dynamic',chunksize=256):
if progress != 0 and (i%booster) == 0: if progress != 0 and (i%booster) == 0:
with gil: with gil:
p.update(_mysum(job_done)) p.update(_mysum(job_done))

2
src/mykdtree.tcc
View File

@ -280,7 +280,7 @@ namespace CosmoTool {
coords tmpBound; coords tmpBound;
NodeIntType nodeId; NodeIntType nodeId;
#pragma omp atomic #pragma omp atomic capture
nodeId = (this->lastNode)++; nodeId = (this->lastNode)++;
node = &nodes[nodeId]; node = &nodes[nodeId];