1# distutils: include_dirs = LIB_DIR 2# distutils: libraries = STD_LIBS 3""" 4Particle smoothing in cells 5 6 7 8 9""" 10 11 12cimport numpy as np 13 14import numpy as np 15 16cimport cython 17from cpython.exc cimport PyErr_CheckSignals 18from libc.math cimport cos, fabs, sin, sqrt 19from libc.stdlib cimport free, malloc, realloc 20from libc.string cimport memmove 21from oct_container cimport Oct, OctInfo, OctreeContainer 22 23 24cdef void spherical_coord_setup(np.float64_t ipos[3], np.float64_t opos[3]): 25 opos[0] = ipos[0] * sin(ipos[1]) * cos(ipos[2]) 26 opos[1] = ipos[0] * sin(ipos[1]) * sin(ipos[2]) 27 opos[2] = ipos[0] * cos(ipos[1]) 28 29cdef void cart_coord_setup(np.float64_t ipos[3], np.float64_t opos[3]): 30 opos[0] = ipos[0] 31 opos[1] = ipos[1] 32 opos[2] = ipos[2] 33 34cdef class ParticleSmoothOperation: 35 def __init__(self, nvals, nfields, max_neighbors, kernel_name): 36 # This is the set of cells, in grids, blocks or octs, we are handling. 37 cdef int i 38 self.nvals = nvals 39 self.nfields = nfields 40 self.maxn = max_neighbors 41 self.sph_kernel = get_kernel_func(kernel_name) 42 43 def initialize(self, *args): 44 raise NotImplementedError 45 46 def finalize(self, *args): 47 raise NotImplementedError 48 49 @cython.cdivision(True) 50 @cython.boundscheck(False) 51 @cython.wraparound(False) 52 @cython.initializedcheck(False) 53 def process_octree(self, OctreeContainer mesh_octree, 54 np.int64_t [:] mdom_ind, 55 np.float64_t[:,:] positions, 56 np.float64_t[:,:] oct_positions, 57 fields = None, int domain_id = -1, 58 int domain_offset = 0, 59 periodicity = (True, True, True), 60 index_fields = None, 61 OctreeContainer particle_octree = None, 62 np.int64_t [:] pdom_ind = None, 63 geometry = "cartesian"): 64 # This will be a several-step operation. 65 # 66 # We first take all of our particles and assign them to Octs. If they 67 # are not in an Oct, we will assume they are out of bounds. Note that 68 # this means that if we have loaded neighbor particles for which an Oct 69 # does not exist, we are going to be discarding them -- so sparse 70 # octrees will need to ensure that neighbor octs *exist*. Particles 71 # will be assigned in a new NumPy array. Note that this incurs 72 # overhead, but reduces complexity as we will now be able to use 73 # argsort. 74 # 75 # After the particles have been assigned to Octs, we process each Oct 76 # individually. We will do this by calling "get" for the *first* 77 # particle in each set of Octs in the sorted list. After this, we get 78 # neighbors for each Oct. 79 # 80 # Now, with the set of neighbors (and thus their indices) we allocate 81 # an array of particles and their fields, fill these in, and call our 82 # process function. 83 # 84 # This is not terribly efficient -- for starters, the neighbor function 85 # is not the most efficient yet. We will also need to handle some 86 # mechanism of an expandable array for holding pointers to Octs, so 87 # that we can deal with >27 neighbors. 88 if particle_octree is None: 89 particle_octree = mesh_octree 90 pdom_ind = mdom_ind 91 cdef int nf, i, j, n 92 cdef int dims[3] 93 cdef np.float64_t[:] *field_check 94 cdef np.float64_t **field_pointers 95 cdef np.float64_t *field_vals 96 cdef np.float64_t pos[3] 97 cdef np.float64_t dds[3] 98 cdef np.float64_t **octree_field_pointers 99 cdef int nsize = 0 100 cdef np.int64_t *nind = NULL 101 cdef OctInfo moi, poi 102 cdef Oct *oct 103 cdef np.int64_t numpart, offset, local_ind, poff 104 cdef np.int64_t moff_p, moff_m 105 cdef np.int64_t[:] pind, doff, pdoms, pcount 106 cdef np.int64_t[:,:] doff_m 107 cdef np.ndarray[np.float64_t, ndim=1] tarr 108 cdef np.ndarray[np.float64_t, ndim=4] iarr 109 cdef np.float64_t[:,:] cart_positions 110 cdef np.float64_t[:,:] oct_left_edges, oct_dds 111 cdef OctInfo oinfo 112 if geometry == "cartesian": 113 self.pos_setup = cart_coord_setup 114 cart_positions = positions 115 elif geometry == "spherical": 116 self.pos_setup = spherical_coord_setup 117 cart_positions = np.empty((positions.shape[0], 3), dtype="float64") 118 119 cart_positions[:,0] = positions[:,0] * \ 120 np.sin(positions[:,1]) * \ 121 np.cos(positions[:,2]) 122 cart_positions[:,1] = positions[:,0] * \ 123 np.sin(positions[:,1]) * \ 124 np.sin(positions[:,2]) 125 cart_positions[:,2] = positions[:,0] * \ 126 np.cos(positions[:,1]) 127 periodicity = (False, False, False) 128 else: 129 raise NotImplementedError 130 dims[0] = dims[1] = dims[2] = (1 << mesh_octree.oref) 131 cdef int nz = dims[0] * dims[1] * dims[2] 132 numpart = positions.shape[0] 133 # pcount is the number of particles per oct. 134 pcount = np.zeros_like(pdom_ind) 135 oct_left_edges = np.zeros((pdom_ind.shape[0], 3), dtype='float64') 136 oct_dds = np.zeros_like(oct_left_edges) 137 # doff is the offset to a given oct in the sorted particles. 138 doff = np.zeros_like(pdom_ind) - 1 139 doff_m = np.zeros((mdom_ind.shape[0], 2), dtype="int64") 140 moff_p = particle_octree.get_domain_offset(domain_id + domain_offset) 141 moff_m = mesh_octree.get_domain_offset(domain_id + domain_offset) 142 # pdoms points particles at their octs. So the value in this array, for 143 # a given index, is the local oct index. 144 pdoms = np.zeros(positions.shape[0], dtype="int64") - 1 145 nf = len(fields) 146 if fields is None: 147 fields = [] 148 field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf) 149 for i in range(nf): 150 tarr = fields[i] 151 field_pointers[i] = <np.float64_t *> tarr.data 152 if index_fields is None: 153 index_fields = [] 154 nf = len(index_fields) 155 index_field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf) 156 for i in range(nf): 157 iarr = index_fields[i] 158 index_field_pointers[i] = <np.float64_t *> iarr.data 159 for i in range(3): 160 self.DW[i] = (mesh_octree.DRE[i] - mesh_octree.DLE[i]) 161 self.periodicity[i] = periodicity[i] 162 cdef np.float64_t factor = (1 << (particle_octree.oref)) 163 for i in range(positions.shape[0]): 164 for j in range(3): 165 pos[j] = positions[i, j] 166 oct = particle_octree.get(pos, &oinfo) 167 if oct == NULL or (domain_id > 0 and oct.domain != domain_id): 168 continue 169 # Note that this has to be our local index, not our in-file index. 170 # This is the particle count, which we'll use once we have sorted 171 # the particles to calculate the offsets into each oct's particles. 172 offset = oct.domain_ind - moff_p 173 pcount[offset] += 1 174 pdoms[i] = offset # We store the *actual* offset. 175 # store oct positions and dds to avoid searching for neighbors 176 # in octs that we know are too far away 177 for j in range(3): 178 oct_left_edges[offset, j] = oinfo.left_edge[j] 179 oct_dds[offset, j] = oinfo.dds[j] * factor 180 # Now we have oct assignments. Let's sort them. 181 # Note that what we will be providing to our processing functions will 182 # actually be indirectly-sorted fields. This preserves memory at the 183 # expense of additional pointer lookups. 184 pind = np.asarray(np.argsort(pdoms), dtype='int64', order='C') 185 # So what this means is that we now have all the oct-0 particle indices 186 # in order, then the oct-1, etc etc. 187 # This now gives us the indices to the particles for each domain. 188 for i in range(positions.shape[0]): 189 # This value, poff, is the index of the particle in the *unsorted* 190 # arrays. 191 poff = pind[i] 192 offset = pdoms[poff] 193 # If we have yet to assign the starting index to this oct, we do so 194 # now. 195 if doff[offset] < 0: doff[offset] = i 196 #print(domain_id, domain_offset, moff_p, moff_m) 197 #raise RuntimeError 198 # Now doff is full of offsets to the first entry in the pind that 199 # refers to that oct's particles. 200 cdef np.ndarray[np.uint8_t, ndim=1] visited 201 visited = np.zeros(mdom_ind.shape[0], dtype="uint8") 202 cdef int nproc = 0 203 # This should be thread-private if we ever go to OpenMP 204 cdef DistanceQueue dist_queue = DistanceQueue(self.maxn) 205 dist_queue._setup(self.DW, self.periodicity) 206 for i in range(oct_positions.shape[0]): 207 if (i % 10000) == 0: 208 PyErr_CheckSignals() 209 for j in range(3): 210 pos[j] = oct_positions[i, j] 211 oct = mesh_octree.get(pos, &moi) 212 offset = mdom_ind[oct.domain_ind - moff_m] * nz 213 if visited[oct.domain_ind - moff_m] == 1: continue 214 visited[oct.domain_ind - moff_m] = 1 215 if offset < 0: continue 216 nproc += 1 217 self.neighbor_process( 218 dims, moi.left_edge, moi.dds, cart_positions, field_pointers, doff, 219 &nind, pind, pcount, offset, index_field_pointers, 220 particle_octree, domain_id, &nsize, oct_left_edges, 221 oct_dds, dist_queue) 222 #print("VISITED", visited.sum(), visited.size,) 223 #print(100.0*float(visited.sum())/visited.size) 224 if nind != NULL: 225 free(nind) 226 227 @cython.cdivision(True) 228 @cython.boundscheck(False) 229 @cython.wraparound(False) 230 @cython.initializedcheck(False) 231 def process_particles(self, OctreeContainer particle_octree, 232 np.ndarray[np.int64_t, ndim=1] pdom_ind, 233 np.ndarray[np.float64_t, ndim=2] positions, 234 fields = None, int domain_id = -1, 235 int domain_offset = 0, 236 periodicity = (True, True, True), 237 geometry = "cartesian"): 238 # The other functions in this base class process particles in a way 239 # that results in a modification to the *mesh*. This function is 240 # designed to process neighboring particles in such a way that a new 241 # *particle* field is defined -- this means that new particle 242 # attributes (*not* mesh attributes) can be created that rely on the 243 # values of nearby particles. For instance, a smoothing kernel, or a 244 # nearest-neighbor field. 245 cdef int nf, i, j, k, n 246 cdef int dims[3] 247 cdef np.float64_t **field_pointers 248 cdef np.float64_t *field_vals 249 cdef np.float64_t dds[3] 250 cdef np.float64_t pos[3] 251 cdef np.float64_t **octree_field_pointers 252 cdef int nsize = 0 253 cdef np.int64_t *nind = NULL 254 cdef OctInfo moi, poi 255 cdef Oct *oct 256 cdef Oct **neighbors = NULL 257 cdef np.int64_t nneighbors, numpart, offset, local_ind 258 cdef np.int64_t moff_p, moff_m, pind0, poff 259 cdef np.int64_t[:] pind, doff, pdoms, pcount 260 cdef np.ndarray[np.float64_t, ndim=1] tarr 261 cdef np.ndarray[np.float64_t, ndim=2] cart_positions 262 if geometry == "cartesian": 263 self.pos_setup = cart_coord_setup 264 cart_positions = positions 265 elif geometry == "spherical": 266 self.pos_setup = spherical_coord_setup 267 cart_positions = np.empty((positions.shape[0], 3), dtype="float64") 268 269 cart_positions[:,0] = positions[:,0] * \ 270 np.sin(positions[:,1]) * \ 271 np.cos(positions[:,2]) 272 cart_positions[:,1] = positions[:,0] * \ 273 np.sin(positions[:,1]) * \ 274 np.sin(positions[:,2]) 275 cart_positions[:,2] = positions[:,0] * \ 276 np.cos(positions[:,1]) 277 periodicity = (False, False, False) 278 else: 279 raise NotImplementedError 280 numpart = positions.shape[0] 281 pcount = np.zeros_like(pdom_ind) 282 doff = np.zeros_like(pdom_ind) - 1 283 moff_p = particle_octree.get_domain_offset(domain_id + domain_offset) 284 pdoms = np.zeros(positions.shape[0], dtype="int64") - 1 285 nf = len(fields) 286 if fields is None: 287 fields = [] 288 field_pointers = <np.float64_t**> alloca(sizeof(np.float64_t *) * nf) 289 for i in range(nf): 290 tarr = fields[i] 291 field_pointers[i] = <np.float64_t *> tarr.data 292 for i in range(3): 293 self.DW[i] = (particle_octree.DRE[i] - particle_octree.DLE[i]) 294 self.periodicity[i] = periodicity[i] 295 for i in range(positions.shape[0]): 296 for j in range(3): 297 pos[j] = positions[i, j] 298 oct = particle_octree.get(pos) 299 if oct == NULL or (domain_id > 0 and oct.domain != domain_id): 300 continue 301 # Note that this has to be our local index, not our in-file index. 302 # This is the particle count, which we'll use once we have sorted 303 # the particles to calculate the offsets into each oct's particles. 304 offset = oct.domain_ind - moff_p 305 pcount[offset] += 1 306 pdoms[i] = offset # We store the *actual* offset. 307 # Now we have oct assignments. Let's sort them. 308 # Note that what we will be providing to our processing functions will 309 # actually be indirectly-sorted fields. This preserves memory at the 310 # expense of additional pointer lookups. 311 pind = np.asarray(np.argsort(pdoms), dtype='int64', order='C') 312 # So what this means is that we now have all the oct-0 particle indices 313 # in order, then the oct-1, etc etc. 314 # This now gives us the indices to the particles for each domain. 315 for i in range(positions.shape[0]): 316 # This value, poff, is the index of the particle in the *unsorted* 317 # arrays. 318 poff = pind[i] 319 offset = pdoms[poff] 320 # If we have yet to assign the starting index to this oct, we do so 321 # now. 322 if doff[offset] < 0: doff[offset] = i 323 #print(domain_id, domain_offset, moff_p, moff_m) 324 #raise RuntimeError 325 # Now doff is full of offsets to the first entry in the pind that 326 # refers to that oct's particles. 327 cdef int maxnei = 0 328 cdef int nproc = 0 329 # This should be thread-private if we ever go to OpenMP 330 cdef DistanceQueue dist_queue = DistanceQueue(self.maxn) 331 dist_queue._setup(self.DW, self.periodicity) 332 for i in range(doff.shape[0]): 333 if doff[i] < 0: continue 334 offset = pind[doff[i]] 335 for j in range(3): 336 pos[j] = positions[offset, j] 337 for j in range(pcount[i]): 338 pind0 = pind[doff[i] + j] 339 for k in range(3): 340 pos[k] = positions[pind0, k] 341 self.neighbor_process_particle(pos, cart_positions, field_pointers, 342 doff, &nind, pind, pcount, pind0, 343 NULL, particle_octree, domain_id, &nsize, 344 dist_queue) 345 #print("VISITED", visited.sum(), visited.size,) 346 #print(100.0*float(visited.sum())/visited.size) 347 if nind != NULL: 348 free(nind) 349 350 cdef int neighbor_search(self, np.float64_t pos[3], OctreeContainer octree, 351 np.int64_t **nind, int *nsize, 352 np.int64_t nneighbors, np.int64_t domain_id, 353 Oct **oct = NULL, int extra_layer = 0): 354 cdef OctInfo oi 355 cdef Oct *ooct 356 cdef Oct **neighbors 357 cdef Oct **first_layer 358 cdef int j, total_neighbors = 0, initial_layer = 0 359 cdef int layer_ind = 0 360 cdef np.int64_t moff = octree.get_domain_offset(domain_id) 361 ooct = octree.get(pos, &oi) 362 if oct != NULL and ooct == oct[0]: 363 return nneighbors 364 oct[0] = ooct 365 if nind[0] == NULL: 366 nsize[0] = 27 367 nind[0] = <np.int64_t *> malloc(sizeof(np.int64_t)*nsize[0]) 368 # This is our "seed" set of neighbors. If we are asked to, we will 369 # create a master list of neighbors that is much bigger and includes 370 # everything. 371 layer_ind = 0 372 first_layer = NULL 373 while 1: 374 neighbors = octree.neighbors(&oi, &nneighbors, ooct, self.periodicity) 375 # Now we have all our neighbors. And, we should be set for what 376 # else we need to do. 377 if total_neighbors + nneighbors > nsize[0]: 378 nind[0] = <np.int64_t *> realloc( 379 nind[0], sizeof(np.int64_t)*(nneighbors + total_neighbors)) 380 nsize[0] = nneighbors + total_neighbors 381 for j in range(nneighbors): 382 # Particle octree neighbor indices 383 nind[0][j + total_neighbors] = neighbors[j].domain_ind - moff 384 total_neighbors += nneighbors 385 if extra_layer == 0: 386 # Not adding on any additional layers here. 387 free(neighbors) 388 neighbors = NULL 389 break 390 if initial_layer == 0: 391 initial_layer = nneighbors 392 first_layer = neighbors 393 else: 394 # Allocated internally; we free this in the loops if we aren't 395 # tracking it 396 free(neighbors) 397 neighbors = NULL 398 ooct = first_layer[layer_ind] 399 layer_ind += 1 400 if layer_ind == initial_layer: 401 neighbors 402 break 403 404 405 for j in range(total_neighbors): 406 # Particle octree neighbor indices 407 if nind[0][j] == -1: continue 408 for n in range(j): 409 if nind[0][j] == nind[0][n]: 410 nind[0][j] = -1 411 # This is allocated by the neighbors function, so we deallocate it. 412 if first_layer != NULL: 413 free(first_layer) 414 return total_neighbors 415 416 @cython.cdivision(True) 417 @cython.boundscheck(False) 418 @cython.wraparound(False) 419 @cython.initializedcheck(False) 420 def process_grid(self, gobj, 421 np.ndarray[np.float64_t, ndim=2] positions, 422 fields = None): 423 raise NotImplementedError 424 425 cdef void process(self, np.int64_t offset, int i, int j, int k, 426 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 427 np.float64_t **ifields, DistanceQueue dq): 428 raise NotImplementedError 429 430 @cython.cdivision(True) 431 @cython.boundscheck(False) 432 @cython.wraparound(False) 433 @cython.initializedcheck(False) 434 cdef void neighbor_find(self, 435 np.int64_t nneighbors, 436 np.int64_t *nind, 437 np.int64_t[:] doffs, 438 np.int64_t[:] pcounts, 439 np.int64_t[:] pinds, 440 np.float64_t[:,:] ppos, 441 np.float64_t cpos[3], 442 np.float64_t[:,:] oct_left_edges, 443 np.float64_t[:,:] oct_dds, 444 DistanceQueue dq 445 ): 446 # We are now given the number of neighbors, the indices into the 447 # domains for them, and the number of particles for each. 448 cdef int ni, i, j, k 449 cdef np.int64_t offset, pn, pc 450 cdef np.float64_t pos[3] 451 cdef np.float64_t ex[2] 452 cdef np.float64_t DR[2] 453 cdef np.float64_t cp, r2_trunc, r2, dist 454 dq.neighbor_reset() 455 for ni in range(nneighbors): 456 if nind[ni] == -1: continue 457 # terminate early if all 8 corners of oct are farther away than 458 # most distant currently known neighbor 459 if oct_left_edges != None and dq.curn == dq.maxn: 460 r2_trunc = dq.neighbors[dq.curn - 1].r2 461 # iterate over each dimension in the outer loop so we can 462 # consolidate temporary storage 463 # What this next bit does is figure out which component is the 464 # closest, of each possible permutation. 465 # k here is the dimension 466 r2 = 0.0 467 for k in range(3): 468 # We start at left edge, then do halfway, then right edge. 469 ex[0] = oct_left_edges[nind[ni], k] 470 ex[1] = ex[0] + oct_dds[nind[ni], k] 471 # There are three possibilities; we are between, left-of, 472 # or right-of the extrema. Thanks to 473 # http://stackoverflow.com/questions/5254838/calculating-distance-between-a-point-and-a-rectangular-box-nearest-point 474 # for some help. This has been modified to account for 475 # periodicity. 476 dist = 0.0 477 DR[0] = (ex[0] - cpos[k]) 478 DR[1] = (cpos[k] - ex[1]) 479 for j in range(2): 480 if not self.periodicity[k]: 481 pass 482 elif (DR[j] > self.DW[k]/2.0): 483 DR[j] -= self.DW[k] 484 elif (DR[j] < -self.DW[k]/2.0): 485 DR[j] += self.DW[k] 486 dist = fmax(dist, DR[j]) 487 r2 += dist*dist 488 if r2 > r2_trunc: 489 continue 490 offset = doffs[nind[ni]] 491 pc = pcounts[nind[ni]] 492 for i in range(pc): 493 pn = pinds[offset + i] 494 for j in range(3): 495 pos[j] = ppos[pn, j] 496 dq.neighbor_eval(pn, pos, cpos) 497 498 @cython.cdivision(True) 499 @cython.boundscheck(False) 500 @cython.wraparound(False) 501 @cython.initializedcheck(False) 502 cdef void neighbor_process(self, int dim[3], np.float64_t left_edge[3], 503 np.float64_t dds[3], np.float64_t[:,:] ppos, 504 np.float64_t **fields, 505 np.int64_t [:] doffs, np.int64_t **nind, 506 np.int64_t [:] pinds, np.int64_t[:] pcounts, 507 np.int64_t offset, 508 np.float64_t **index_fields, 509 OctreeContainer octree, np.int64_t domain_id, 510 int *nsize, np.float64_t[:,:] oct_left_edges, 511 np.float64_t[:,:] oct_dds, 512 DistanceQueue dq): 513 # Note that we assume that fields[0] == smoothing length in the native 514 # units supplied. We can now iterate over every cell in the block and 515 # every particle to find the nearest. We will use a priority heap. 516 cdef int i, j, k, ntot, nntot, m, nneighbors 517 cdef np.float64_t cpos[3] 518 cdef np.float64_t opos[3] 519 cdef Oct* oct = NULL 520 cpos[0] = left_edge[0] + 0.5*dds[0] 521 for i in range(dim[0]): 522 cpos[1] = left_edge[1] + 0.5*dds[1] 523 for j in range(dim[1]): 524 cpos[2] = left_edge[2] + 0.5*dds[2] 525 for k in range(dim[2]): 526 self.pos_setup(cpos, opos) 527 nneighbors = self.neighbor_search(opos, octree, 528 nind, nsize, nneighbors, domain_id, &oct, 0) 529 self.neighbor_find(nneighbors, nind[0], doffs, pcounts, 530 pinds, ppos, opos, oct_left_edges, 531 oct_dds, dq) 532 # Now we have all our neighbors in our neighbor list. 533 if dq.curn <-1*dq.maxn: 534 ntot = nntot = 0 535 for m in range(nneighbors): 536 if nind[0][m] < 0: continue 537 nntot += 1 538 ntot += pcounts[nind[0][m]] 539 print("SOMETHING WRONG", dq.curn, nneighbors, ntot, nntot) 540 self.process(offset, i, j, k, dim, opos, fields, 541 index_fields, dq) 542 cpos[2] += dds[2] 543 cpos[1] += dds[1] 544 cpos[0] += dds[0] 545 546 @cython.cdivision(True) 547 @cython.boundscheck(False) 548 @cython.wraparound(False) 549 @cython.initializedcheck(False) 550 cdef void neighbor_process_particle(self, np.float64_t cpos[3], 551 np.float64_t[:,:] ppos, 552 np.float64_t **fields, 553 np.int64_t[:] doffs, np.int64_t **nind, 554 np.int64_t[:] pinds, np.int64_t[:] pcounts, 555 np.int64_t offset, 556 np.float64_t **index_fields, 557 OctreeContainer octree, 558 np.int64_t domain_id, int *nsize, 559 DistanceQueue dq): 560 # Note that we assume that fields[0] == smoothing length in the native 561 # units supplied. We can now iterate over every cell in the block and 562 # every particle to find the nearest. We will use a priority heap. 563 cdef int i, j, k, ntot, nntot, m 564 cdef int dim[3] 565 cdef Oct *oct = NULL 566 cdef np.int64_t nneighbors = 0 567 i = j = k = 0 568 dim[0] = dim[1] = dim[2] = 1 569 cdef np.float64_t opos[3] 570 self.pos_setup(cpos, opos) 571 nneighbors = self.neighbor_search(opos, octree, 572 nind, nsize, nneighbors, domain_id, &oct, 0) 573 self.neighbor_find(nneighbors, nind[0], doffs, pcounts, pinds, ppos, 574 opos, None, None, dq) 575 self.process(offset, i, j, k, dim, opos, fields, index_fields, dq) 576 577cdef class VolumeWeightedSmooth(ParticleSmoothOperation): 578 # This smoothing function evaluates the field, *without* normalization, at 579 # every point in the *mesh*. Applying a normalization results in 580 # non-conservation of mass when smoothing density; to avoid this, we do not 581 # apply this normalization factor. The SPLASH paper 582 # (http://arxiv.org/abs/0709.0832v1) discusses this in detail; what we are 583 # applying here is equation 6, with variable smoothing lengths (eq 13). 584 cdef np.float64_t **fp 585 cdef public object vals 586 def initialize(self): 587 cdef int i 588 if self.nfields < 4: 589 # We need four fields -- the mass should be the first, then the 590 # smoothing length for particles, the normalization factor to 591 # ensure mass conservation, then the field we're smoothing. 592 raise RuntimeError 593 cdef np.ndarray tarr 594 self.fp = <np.float64_t **> malloc( 595 sizeof(np.float64_t *) * (self.nfields - 3)) 596 self.vals = [] 597 # We usually only allocate one field; if we are doing multiple field, 598 # single-pass smoothing, then we might have more. 599 for i in range(self.nfields - 3): 600 tarr = np.zeros(self.nvals, dtype="float64", order="F") 601 self.vals.append(tarr) 602 self.fp[i] = <np.float64_t *> tarr.data 603 604 def finalize(self): 605 free(self.fp) 606 return self.vals 607 608 @cython.cdivision(True) 609 @cython.boundscheck(False) 610 @cython.wraparound(False) 611 @cython.initializedcheck(False) 612 cdef void process(self, np.int64_t offset, int i, int j, int k, 613 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 614 np.float64_t **index_fields, DistanceQueue dq): 615 # We have our i, j, k for our cell, as well as the cell position. 616 # We also have a list of neighboring particles with particle numbers. 617 cdef int n, fi 618 cdef np.float64_t weight, r2, val, hsml, dens, mass, coeff, max_r 619 cdef np.float64_t max_hsml, ihsml, ihsml3, kern 620 coeff = 0.0 621 cdef np.int64_t pn 622 # We get back our mass 623 # rho_i = sum(j = 1 .. n) m_j * W_ij 624 max_r = sqrt(dq.neighbors[dq.curn-1].r2) 625 max_hsml = index_fields[0][gind(i,j,k,dim) + offset] 626 for n in range(dq.curn): 627 # No normalization for the moment. 628 # fields[0] is the smoothing length. 629 r2 = dq.neighbors[n].r2 630 pn = dq.neighbors[n].pn 631 # Smoothing kernel weight function 632 mass = fields[0][pn] 633 hsml = fields[1][pn] 634 dens = fields[2][pn] 635 if hsml < 0: 636 hsml = max_r 637 if hsml == 0: continue 638 ihsml = 1.0/hsml 639 hsml = fmax(max_hsml/2.0, hsml) 640 ihsml3 = ihsml*ihsml*ihsml 641 # Usually this density has been computed 642 if dens == 0.0: continue 643 weight = (mass / dens) * ihsml3 644 kern = self.sph_kernel(sqrt(r2) * ihsml) 645 weight *= kern 646 # Mass of the particle times the value 647 for fi in range(self.nfields - 3): 648 val = fields[fi + 3][pn] 649 self.fp[fi][gind(i,j,k,dim) + offset] += val * weight 650 return 651 652volume_weighted_smooth = VolumeWeightedSmooth 653 654cdef class NearestNeighborSmooth(ParticleSmoothOperation): 655 cdef np.float64_t *fp 656 cdef public object vals 657 def initialize(self): 658 cdef np.ndarray tarr 659 assert(self.nfields == 1) 660 tarr = np.zeros(self.nvals, dtype="float64", order="F") 661 self.vals = tarr 662 self.fp = <np.float64_t *> tarr.data 663 664 def finalize(self): 665 return self.vals 666 667 @cython.cdivision(True) 668 @cython.boundscheck(False) 669 @cython.wraparound(False) 670 @cython.initializedcheck(False) 671 cdef void process(self, np.int64_t offset, int i, int j, int k, 672 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 673 np.float64_t **index_fields, DistanceQueue dq): 674 # We have our i, j, k for our cell, as well as the cell position. 675 # We also have a list of neighboring particles with particle numbers. 676 cdef np.int64_t pn 677 # We get back our mass 678 # rho_i = sum(j = 1 .. n) m_j * W_ij 679 pn = dq.neighbors[0].pn 680 self.fp[gind(i,j,k,dim) + offset] = fields[0][pn] 681 #self.fp[gind(i,j,k,dim) + offset] = dq.neighbors[0].r2 682 return 683 684nearest_smooth = NearestNeighborSmooth 685 686cdef class IDWInterpolationSmooth(ParticleSmoothOperation): 687 cdef np.float64_t *fp 688 cdef public int p2 689 cdef public object vals 690 def initialize(self): 691 cdef np.ndarray tarr 692 assert(self.nfields == 1) 693 tarr = np.zeros(self.nvals, dtype="float64", order="F") 694 self.vals = tarr 695 self.fp = <np.float64_t *> tarr.data 696 self.p2 = 2 # Power, for IDW, in units of 2. So we only do even p's. 697 698 def finalize(self): 699 return self.vals 700 701 @cython.cdivision(True) 702 @cython.boundscheck(False) 703 @cython.wraparound(False) 704 @cython.initializedcheck(False) 705 cdef void process(self, np.int64_t offset, int i, int j, int k, 706 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 707 np.float64_t **index_fields, DistanceQueue dq): 708 # We have our i, j, k for our cell, as well as the cell position. 709 # We also have a list of neighboring particles with particle numbers. 710 cdef np.int64_t pn, ni, di 711 cdef np.float64_t total_weight = 0.0, total_value = 0.0, r2, val, w 712 # We're going to do a very simple IDW average 713 if dq.neighbors[0].r2 == 0.0: 714 pn = dq.neighbors[0].pn 715 self.fp[gind(i,j,k,dim) + offset] = fields[0][pn] 716 for ni in range(dq.curn): 717 r2 = dq.neighbors[ni].r2 718 val = fields[0][dq.neighbors[ni].pn] 719 w = r2 720 for di in range(self.p2 - 1): 721 w *= r2 722 total_value += w * val 723 total_weight += w 724 self.fp[gind(i,j,k,dim) + offset] = total_value / total_weight 725 return 726 727idw_smooth = IDWInterpolationSmooth 728 729cdef class NthNeighborDistanceSmooth(ParticleSmoothOperation): 730 731 def initialize(self): 732 return 733 734 def finalize(self): 735 return 736 737 @cython.cdivision(True) 738 @cython.boundscheck(False) 739 @cython.wraparound(False) 740 @cython.initializedcheck(False) 741 cdef void process(self, np.int64_t offset, int i, int j, int k, 742 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 743 np.float64_t **index_fields, DistanceQueue dq): 744 cdef np.float64_t max_r 745 # We assume "offset" here is the particle index. 746 max_r = sqrt(dq.neighbors[dq.curn-1].r2) 747 fields[0][offset] = max_r 748 749nth_neighbor_smooth = NthNeighborDistanceSmooth 750 751cdef class SmoothedDensityEstimate(ParticleSmoothOperation): 752 def initialize(self): 753 return 754 755 def finalize(self): 756 return 757 758 @cython.cdivision(True) 759 @cython.boundscheck(False) 760 @cython.wraparound(False) 761 @cython.initializedcheck(False) 762 cdef void process(self, np.int64_t offset, int i, int j, int k, 763 int dim[3], np.float64_t cpos[3], np.float64_t **fields, 764 np.float64_t **index_fields, DistanceQueue dq): 765 cdef np.float64_t r2, hsml, dens, mass, weight, lw 766 cdef int pn 767 # We assume "offset" here is the particle index. 768 hsml = sqrt(dq.neighbors[dq.curn-1].r2) 769 dens = 0.0 770 weight = 0.0 771 for pn in range(dq.curn): 772 mass = fields[0][dq.neighbors[pn].pn] 773 r2 = dq.neighbors[pn].r2 774 lw = self.sph_kernel(sqrt(r2) / hsml) 775 dens += mass * lw 776 weight = (4.0/3.0) * 3.1415926 * hsml**3 777 fields[1][offset] = dens/weight 778 779density_smooth = SmoothedDensityEstimate 780