| /*========================================================================= |
| |
| Copyright (c) 2007, Los Alamos National Security, LLC |
| |
| All rights reserved. |
| |
| Copyright 2007. Los Alamos National Security, LLC. |
| This software was produced under U.S. Government contract DE-AC52-06NA25396 |
| for Los Alamos National Laboratory (LANL), which is operated by |
| Los Alamos National Security, LLC for the U.S. Department of Energy. |
| The U.S. Government has rights to use, reproduce, and distribute this software. |
| NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, |
| EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. |
| If software is modified to produce derivative works, such modified software |
| should be clearly marked, so as not to confuse it with the version available |
| from LANL. |
| |
| Additionally, redistribution and use in source and binary forms, with or |
| without modification, are permitted provided that the following conditions |
| are met: |
| - Redistributions of source code must retain the above copyright notice, |
| this list of conditions and the following disclaimer. |
| - Redistributions in binary form must reproduce the above copyright notice, |
| this list of conditions and the following disclaimer in the documentation |
| and/or other materials provided with the distribution. |
| - Neither the name of Los Alamos National Security, LLC, Los Alamos National |
| Laboratory, LANL, the U.S. Government, nor the names of its contributors |
| may be used to endorse or promote products derived from this software |
| without specific prior written permission. |
| |
| THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS |
| "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, |
| THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR |
| CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; |
| OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, |
| WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR |
| OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF |
| ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| |
| =========================================================================*/ |
| |
| /*========================================================================= |
| |
| Copyright (c) 2011-2012 Argonne National Laboratory |
| All rights reserved. |
| |
| Redistribution and use in source and binary forms, with or without |
| modification, are permitted provided that the following conditions |
| are met: |
| 1. Redistributions of source code must retain the above copyright |
| notice, this list of conditions and the following disclaimer. |
| 2. Redistributions in binary form must reproduce the above copyright |
| notice, this list of conditions and the following disclaimer in the |
| documentation and/or other materials provided with the distribution. |
| |
| THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED |
| TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR |
| PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS |
| BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
| CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
| SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
| INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
| CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| POSSIBILITY OF SUCH DAMAGE. |
| |
| =========================================================================*/ |
| |
| |
| /* |
| BG/Q tuned version of HACC: 69.2% of peak performance on 96 racks of Sequoia |
| Argonne Leadership Computing Facility, Argonne, IL 60439 |
| Vitali Morozov (morozov@anl.gov) |
| Hal Finkel (hfinkel@anl.gov) |
| */ |
| |
| |
| #include "hip/hip_runtime.h" |
| #include "Timings.h" |
| #include "RCBForceTree.h" |
| #include "Partition.h" |
| |
| #include <cstring> |
| #include <cstdio> |
| #include <ctime> |
| #include <stdexcept> |
| #include <assert.h> |
| using namespace std; |
| |
| #ifdef __HIPCC__ |
| #include <cudaUtil.h> |
| |
| #define TILEX 4 //Unroll factor in the x dimension, best if 2 or 4, could also add 8 but that is too many registers |
| #define TILEY 4 //Unroll factor in the y dimension, best if 2 or 4, could add 8 but that is too many registers |
| #define BLOCKX 32 //Block size in the x dimension (should be 32) |
| #define BLOCKY 4 //Block size in the y dimension |
| #define MAXX 32 //Maximum blocks in the X dimension, smaller=more reuse but less parallelism |
| #define MAXY 256 //Maximum blocks in the Y dimension, there isn't much reason to make this smaller |
| |
| #define ALIGNX(n) ((n+TILEX-1)/TILEX*TILEX) //Rounds an integer to align with TILEX |
| #define ALIGNY(n) ((n+TILEY-1)/TILEY*TILEY) //Rounds an integer to align with TILEY |
| |
| cudaDeviceSelector __selector__; |
| #endif |
| |
| #ifdef __HIPCC__ |
| #define __HOST__ __host__ |
| #define __DEVICE__ __device__ |
| #else |
| #define __HOST__ |
| #define __DEVICE__ |
| #endif |
| |
| |
| // References: |
| // Emanuel Gafton and Stephan Rosswog. A fast recursive coordinate bisection tree for |
| // neighbour search and gravity. Mon. Not. R. Astron. Soc. to appear, 2011. |
| // http://arxiv.org/abs/1108.0028v1 |
| // |
| // Atsushi Kawai, Junichiro Makino and Toshikazu Ebisuzaki. |
| // Performance Analysis of High-Accuracy Tree Code Based on the Pseudoparticle |
| // Multipole Method. The Astrophysical Journal Supplement Series, 151:13-33, 2004. |
| // Related: http://arxiv.org/abs/astro-ph/0012041v1 |
| // |
| // R. H. Hardin and N. J. Sloane |
| // New Spherical 4-Designs. Discrete Math, 106/107 255-264, 1992. |
| // |
| // The library of spherical designs: |
| // http://www2.research.att.com/~njas/sphdesigns/ |
| namespace { |
| template <int TDPTS> |
| struct sphdesign {}; |
| |
| #define DECLARE_SPHDESIGN(TDPTS) \ |
| template <> \ |
| struct sphdesign<TDPTS> \ |
| { \ |
| static const POSVEL_T x[TDPTS]; \ |
| static const POSVEL_T y[TDPTS]; \ |
| static const POSVEL_T z[TDPTS]; \ |
| }; \ |
| /**/ |
| |
| DECLARE_SPHDESIGN(1) |
| DECLARE_SPHDESIGN(2) |
| DECLARE_SPHDESIGN(3) |
| DECLARE_SPHDESIGN(4) |
| DECLARE_SPHDESIGN(6) |
| DECLARE_SPHDESIGN(12) |
| DECLARE_SPHDESIGN(14) |
| |
| #undef DECLARE_SPHDESIGN |
| |
| /* this is not a t-design, but puts the monopole moment |
| at the center of mass. */ |
| const POSVEL_T sphdesign<1>::x[] = { |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<1>::y[] = { |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<1>::z[] = { |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<2>::x[] = { |
| 1.0, |
| -1.0 |
| }; |
| |
| const POSVEL_T sphdesign<2>::y[] = { |
| 0, |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<2>::z[] = { |
| 0, |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<3>::x[] = { |
| 1.0, |
| -.5, |
| -.5 |
| }; |
| |
| const POSVEL_T sphdesign<3>::y[] = { |
| 0, |
| .86602540378443864675, |
| -.86602540378443864675 |
| }; |
| |
| const POSVEL_T sphdesign<3>::z[] = { |
| 0, |
| 0, |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<4>::x[] = { |
| .577350269189625763, |
| .577350269189625763, |
| -.577350269189625763, |
| -.577350269189625763 |
| }; |
| |
| const POSVEL_T sphdesign<4>::y[] = { |
| .577350269189625763, |
| -.577350269189625763, |
| .577350269189625763, |
| -.577350269189625763 |
| }; |
| |
| const POSVEL_T sphdesign<4>::z[] = { |
| .577350269189625763, |
| -.577350269189625763, |
| -.577350269189625763, |
| .577350269189625763 |
| }; |
| |
| const POSVEL_T sphdesign<6>::x[] = { |
| 1.0, |
| -1.0, |
| 0, |
| 0, |
| 0, |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<6>::y[] = { |
| 0, |
| 0, |
| 1.0, |
| -1.0, |
| 0, |
| 0 |
| }; |
| |
| const POSVEL_T sphdesign<6>::z[] = { |
| 0, |
| 0, |
| 0, |
| 0, |
| 1.0, |
| -1.0 |
| }; |
| |
| // This is a 3-D 12-point spherical 4-design |
| // (the verticies of a icosahedron) from Hardin and Sloane. |
| const POSVEL_T sphdesign<12>::x[] = { |
| 0, |
| 0, |
| 0.525731112119134, |
| -0.525731112119134, |
| 0.85065080835204, |
| -0.85065080835204, |
| 0, |
| 0, |
| -0.525731112119134, |
| 0.525731112119134, |
| -0.85065080835204, |
| 0.85065080835204 |
| }; |
| |
| const POSVEL_T sphdesign<12>::y[] = { |
| 0.85065080835204, |
| 0.85065080835204, |
| 0, |
| 0, |
| 0.525731112119134, |
| 0.525731112119134, |
| -0.85065080835204, |
| -0.85065080835204, |
| 0, |
| 0, |
| -0.525731112119134, |
| -0.525731112119134 |
| }; |
| |
| const POSVEL_T sphdesign<12>::z[] = { |
| 0.525731112119134, |
| -0.525731112119134, |
| 0.85065080835204, |
| 0.85065080835204, |
| 0, |
| 0, |
| -0.525731112119134, |
| 0.525731112119134, |
| -0.85065080835204, |
| -0.85065080835204, |
| 0, |
| 0 |
| }; |
| |
| // This is a 3-D 14-point spherical 4-design by |
| // R. H. Hardin and N. J. A. Sloane. |
| const POSVEL_T sphdesign<14>::x[] = { |
| 1.0e0, |
| 5.947189772040725e-1, |
| 5.947189772040725e-1, |
| 5.947189772040725e-1, |
| -5.947189772040725e-1, |
| -5.947189772040725e-1, |
| -5.947189772040725e-1, |
| 3.012536847870683e-1, |
| 3.012536847870683e-1, |
| 3.012536847870683e-1, |
| -3.012536847870683e-1, |
| -3.012536847870683e-1, |
| -3.012536847870683e-1, |
| -1.0e0 |
| }; |
| |
| const POSVEL_T sphdesign<14>::y[] = { |
| 0.0e0, |
| 1.776539926025823e-1, |
| -7.678419429698292e-1, |
| 5.90187950367247e-1, |
| 1.776539926025823e-1, |
| 5.90187950367247e-1, |
| -7.678419429698292e-1, |
| 8.79474443923065e-1, |
| -7.588425179318781e-1, |
| -1.206319259911869e-1, |
| 8.79474443923065e-1, |
| -1.206319259911869e-1, |
| -7.588425179318781e-1, |
| 0.0e0 |
| }; |
| |
| const POSVEL_T sphdesign<14>::z[] = { |
| 0.0e0, |
| 7.840589244857197e-1, |
| -2.381765915652909e-1, |
| -5.458823329204288e-1, |
| -7.840589244857197e-1, |
| 5.458823329204288e-1, |
| 2.381765915652909e-1, |
| 3.684710570566285e-1, |
| 5.774116818882528e-1, |
| -9.458827389448813e-1, |
| -3.684710570566285e-1, |
| 9.458827389448813e-1, |
| -5.774116818882528e-1, |
| 0.0e0 |
| }; |
| } // anonymous namespace |
| |
| // Note: In Gafton and Rosswog the far-field force contribution is calculated |
| // per-cell (at the center of mass), and then a Taylor expansion about the center |
| // of mass is used to calculate the force on the individual particles. For this to |
| // work, the functional form of the force must be known (because the Jacobian |
| // and Hessian are required). Here, however, the functional form is not known, |
| // and so the pseudo-particle method of Makino is used instead. |
| |
| template <int TDPTS> |
| RCBForceTree<TDPTS>::RCBForceTree( |
| POSVEL_T* minLoc, |
| POSVEL_T* maxLoc, |
| POSVEL_T* minForceLoc, |
| POSVEL_T* maxForceLoc, |
| ID_T count, |
| POSVEL_T* xLoc, |
| POSVEL_T* yLoc, |
| POSVEL_T* zLoc, |
| POSVEL_T* xVel, |
| POSVEL_T* yVel, |
| POSVEL_T* zVel, |
| POSVEL_T* ms, |
| POSVEL_T* phiLoc, |
| ID_T *idLoc, |
| MASK_T *maskLoc, |
| POSVEL_T avgMass, |
| POSVEL_T fsm, |
| POSVEL_T r, |
| POSVEL_T oa, |
| ID_T nd, |
| ID_T ds, |
| ID_T tmin, |
| ForceLaw *fl, |
| float fcoeff, |
| POSVEL_T ppc) |
| { |
| #ifdef __HIPCC__ |
| //nvtxRangeId_t r0; |
| //r0=nvtxRangeStartA("RCBForceTree()"); |
| #endif |
| // Extract the contiguous data block from a vector pointer |
| particleCount = count; |
| |
| xx = xLoc; |
| yy = yLoc; |
| zz = zLoc; |
| vx = xVel; |
| vy = yVel; |
| vz = zVel; |
| mass = ms; |
| |
| numThreads=1; |
| |
| // static size for the interaction list |
| #define VMAX ALIGNY(16384) |
| nx_v=(POSVEL_T*)malloc(VMAX*sizeof(POSVEL_T)*numThreads); |
| ny_v=(POSVEL_T*)malloc(VMAX*sizeof(POSVEL_T)*numThreads); |
| nz_v=(POSVEL_T*)malloc(VMAX*sizeof(POSVEL_T)*numThreads); |
| nm_v=(POSVEL_T*)malloc(VMAX*sizeof(POSVEL_T)*numThreads); |
| |
| #ifdef __HIPCC__ |
| hipHostRegister(nx_v,VMAX*sizeof(POSVEL_T)*numThreads,0); |
| hipHostRegister(ny_v,VMAX*sizeof(POSVEL_T)*numThreads,0); |
| hipHostRegister(nz_v,VMAX*sizeof(POSVEL_T)*numThreads,0); |
| hipHostRegister(nm_v,VMAX*sizeof(POSVEL_T)*numThreads,0); |
| hipHostRegister(xx,count*sizeof(POSVEL_T),0); |
| hipHostRegister(yy,count*sizeof(POSVEL_T),0); |
| hipHostRegister(zz,count*sizeof(POSVEL_T),0); |
| hipHostRegister(mass,count*sizeof(POSVEL_T),0); |
| hipHostRegister(vx,count*sizeof(POSVEL_T),0); |
| hipHostRegister(vy,count*sizeof(POSVEL_T),0); |
| hipHostRegister(vz,count*sizeof(POSVEL_T),0); |
| |
| int size=ALIGNY(nd); |
| hipMalloc(&d_xx,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_yy,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_zz,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_vx,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_vy,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_vz,size*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_mass,size*sizeof(POSVEL_T)*numThreads); |
| |
| hipMalloc(&d_nx_v,VMAX*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_ny_v,VMAX*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_nz_v,VMAX*sizeof(POSVEL_T)*numThreads); |
| hipMalloc(&d_nm_v,VMAX*sizeof(POSVEL_T)*numThreads); |
| cudaCheckError(); |
| |
| |
| event_v=(hipEvent_t*)malloc(sizeof(hipEvent_t)*numThreads); |
| stream_v=(hipStream_t*)malloc(sizeof(hipStream_t)*numThreads); |
| for(int i=0;i<numThreads;i++) { |
| hipEventCreate(&event_v[i]); |
| hipStreamCreate(&stream_v[i]); |
| } |
| cudaCheckError(); |
| #endif |
| |
| phi = phiLoc; |
| id = idLoc; |
| mask = maskLoc; |
| |
| particleMass = avgMass; |
| fsrrmax = fsm; |
| rsm = r; |
| sinOpeningAngle = sinf(oa); |
| tanOpeningAngle = tanf(oa); |
| nDirect = nd; |
| depthSafety = ds; |
| taskPartMin = tmin; |
| ppContract = ppc; |
| |
| // Find the grid size of this chaining mesh |
| for (int dim = 0; dim < DIMENSION; dim++) { |
| minRange[dim] = minLoc[dim]; |
| maxRange[dim] = maxLoc[dim]; |
| minForceRange[dim] = minForceLoc[dim]; |
| maxForceRange[dim] = maxForceLoc[dim]; |
| } |
| |
| if (fl) { |
| m_own_fl = false; |
| m_fl = fl; |
| m_fcoeff = fcoeff; |
| } else { |
| //maybe change this to Newton's law or something |
| m_own_fl = true; |
| m_fl = new ForceLawNewton(); |
| m_fcoeff = 1.0; |
| } |
| |
| // Because the tree may be built in parallel, and no efficient way of locking |
| // the tree seems to be available in OpenMP (no reader/writer locks, etc.), |
| // we just estimate the number of tree nodes that will be needed. Hopefully, |
| // this will be an over estimate. If we need more than this, then tree nodes |
| // that really should be subdivided will not be. |
| // |
| // If the tree were perfectly balanced, then it would have a depth of |
| // log_2(particleCount/nDirect). The tree needs to have (2^depth)+1 entries. |
| // To that, a safety factor is added to the depth. |
| ID_T nds = (((ID_T)(particleCount/(POSVEL_T)nDirect)) << depthSafety) + 1; |
| tree.reserve(nds); |
| |
| int nthreads = 1; |
| |
| timespec b_start, b_end; |
| clock_gettime(CLOCK_THREAD_CPUTIME_ID, &b_start); |
| // Create the recursive RCB tree from the particle locations |
| #ifdef __HIPCC__ |
| //nvtxRangeId_t r1; |
| //r1=nvtxRangeStartA("createRCBForceTree"); |
| #endif |
| createRCBForceTree(); |
| |
| #ifdef __HIPCC__ |
| //nvtxRangeEnd(r1); |
| #endif |
| |
| clock_gettime(CLOCK_THREAD_CPUTIME_ID, &b_end); |
| double b_time = (b_end.tv_sec - b_start.tv_sec); |
| b_time += 1e-9*(b_end.tv_nsec - b_start.tv_nsec); |
| |
| printStats(b_time); |
| |
| // Interaction lists. |
| inx.resize(nthreads); |
| iny.resize(nthreads); |
| inz.resize(nthreads); |
| inm.resize(nthreads); |
| iq.resize(nthreads); |
| |
| #ifdef __HIPCC__ |
| //r1=nvtxRangeStartA("calcInternodeForces"); |
| #endif |
| calcInternodeForces(); |
| #ifdef __HIPCC__ |
| //nvtxRangeEnd(r1); |
| #endif |
| |
| #ifdef __HIPCC__ |
| //nvtxRangeEnd(r0); |
| #endif |
| } |
| |
| template <int TDPTS> |
| RCBForceTree<TDPTS>::~RCBForceTree() |
| { |
| #ifdef __HIPCC__ |
| //nvtxRangeId_t r0; |
| //r0=nvtxRangeStartA("~RCBForceTree"); |
| #endif |
| if (m_own_fl) { |
| delete m_fl; |
| } |
| #ifdef __HIPCC__ |
| hipHostUnregister(xx); |
| hipHostUnregister(yy); |
| hipHostUnregister(zz); |
| hipHostUnregister(mass); |
| hipHostUnregister(vx); |
| hipHostUnregister(vy); |
| hipHostUnregister(vz); |
| hipHostUnregister(nx_v); |
| hipHostUnregister(ny_v); |
| hipHostUnregister(nz_v); |
| hipHostUnregister(nm_v); |
| |
| hipFree(d_xx); |
| hipFree(d_yy); |
| hipFree(d_zz); |
| hipFree(d_vx); |
| hipFree(d_vy); |
| hipFree(d_vz); |
| hipFree(d_mass); |
| hipFree(d_nx_v); |
| hipFree(d_ny_v); |
| hipFree(d_nz_v); |
| hipFree(d_nm_v); |
| cudaCheckError(); |
| |
| for(int i=0;i<numThreads;i++) { |
| hipEventDestroy(event_v[i]); |
| hipStreamDestroy(stream_v[i]); |
| } |
| cudaCheckError(); |
| |
| free(event_v); |
| free(stream_v); |
| |
| #endif |
| free(nx_v); |
| free(ny_v); |
| free(nz_v); |
| free(nm_v); |
| #ifdef __HIPCC__ |
| //nvtxRangeEnd(r0); |
| #endif |
| } |
| |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::printStats(double buildTime) |
| { |
| size_t zeroLeafNodes = 0; |
| size_t nonzeroLeafNodes = 0; |
| size_t maxPPN = 0; |
| size_t leafParts = 0; |
| |
| for (ID_T tl = 1; tl < (ID_T) tree.size(); ++tl) { |
| if (tree[tl].cl == 0 && tree[tl].cr == 0) { |
| if (tree[tl].count > 0) { |
| ++nonzeroLeafNodes; |
| |
| leafParts += tree[tl].count; |
| maxPPN = std::max((size_t) tree[tl].count, maxPPN); |
| } else { |
| ++zeroLeafNodes; |
| } |
| } |
| } |
| |
| double localParticleCount = particleCount; |
| double localTreeSize = tree.size(); |
| double localTreeCapacity = tree.capacity(); |
| double localLeaves = zeroLeafNodes+nonzeroLeafNodes; |
| double localEmptyLeaves = zeroLeafNodes; |
| double localMeanPPN = leafParts/((double) nonzeroLeafNodes); |
| unsigned long localMaxPPN = maxPPN; |
| double localBuildTime = buildTime; |
| |
| /* |
| double globalParticleCount; |
| double globalTreeSize; |
| double globalTreeCapacity; |
| double globalLeaves; |
| double globalEmptyLeaves; |
| double globalMeanPPN; |
| unsigned long globalMaxPPN; |
| double globalBuildTime; |
| |
| bool printHere = true; |
| */ |
| |
| if ( Partition::getMyProc() == 0 ) { |
| printf("\ttree post-build statistics (local for rank 0):\n"); |
| printf("\t\tparticles: %.2f\n", localParticleCount); |
| printf("\t\tnodes: %.2f (allocated: %.2f)\n", localTreeSize, localTreeCapacity); |
| printf("\t\tleaves: %.2f (empty: %.2f)\n", localLeaves, localEmptyLeaves); |
| printf("\t\tmean ppn: %.2f (max ppn: %lu)\n", localMeanPPN, localMaxPPN); |
| printf("\t\tbuild time: %g s\n", localBuildTime); |
| } |
| } |
| |
| |
| extern "C" void cm(ID_T count, const POSVEL_T* xx, const POSVEL_T* yy, |
| const POSVEL_T* zz, const POSVEL_T* mass, |
| POSVEL_T* xmin, POSVEL_T* xmax, POSVEL_T* xc); |
| |
| static inline POSVEL_T pptdr(const POSVEL_T* xmin, const POSVEL_T* xmax, const POSVEL_T* xc) |
| { |
| return std::min(xmax[0] - xc[0], std::min(xmax[1] - xc[1], std::min(xmax[2] - xc[2], std::min(xc[0] - xmin[0], |
| std::min(xc[1] - xmin[1], xc[2] - xmin[2]))))); |
| } |
| |
| template <int TDPTS> |
| static inline void pppts(POSVEL_T tdr, const POSVEL_T* xc, |
| POSVEL_T* ppx, POSVEL_T* ppy, POSVEL_T* ppz) |
| { |
| for (int i = 0; i < TDPTS; ++i) { |
| ppx[i] = tdr*sphdesign<TDPTS>::x[i] + xc[0]; |
| ppy[i] = tdr*sphdesign<TDPTS>::y[i] + xc[1]; |
| ppz[i] = tdr*sphdesign<TDPTS>::z[i] + xc[2]; |
| } |
| } |
| |
| template <int TDPTS> |
| static inline void pp(ID_T count, const POSVEL_T* xx, const POSVEL_T* yy, |
| const POSVEL_T* zz, const POSVEL_T* mass, const POSVEL_T* xc, |
| const POSVEL_T* ppx, const POSVEL_T* ppy, const POSVEL_T* ppz, |
| POSVEL_T* ppm, POSVEL_T tdr) |
| { |
| POSVEL_T K = TDPTS; |
| POSVEL_T odr0 = 1/K; |
| |
| for (int i = 0; i < count; ++i) { |
| POSVEL_T xi = xx[i] - xc[0]; |
| POSVEL_T yi = yy[i] - xc[1]; |
| POSVEL_T zi = zz[i] - xc[2]; |
| POSVEL_T ri = sqrtf(xi*xi + yi*yi + zi*zi); |
| |
| for (int j = 0; j < TDPTS; ++j) { |
| POSVEL_T xj = ppx[j] - xc[0]; |
| POSVEL_T yj = ppy[j] - xc[1]; |
| POSVEL_T zj = ppz[j] - xc[2]; |
| POSVEL_T rj2 = xj*xj + yj*yj + zj*zj; |
| |
| POSVEL_T odr1 = 0, odr2 = 0; |
| if (rj2 != 0) { |
| POSVEL_T rj = sqrtf(rj2); |
| POSVEL_T aij = (xi*xj + yi*yj + zi*zj)/(ri*rj); |
| |
| odr1 = (3/K)*(ri/tdr)*aij; |
| odr2 = (5/K)*(ri/tdr)*(ri/tdr)*0.5*(3*aij*aij - 1); |
| } |
| |
| ppm[j] += mass[i]*(odr0 + odr1 + odr2); |
| } |
| } |
| } |
| |
| #ifdef __HIPCC__ |
| |
| typedef long long int int64; |
| |
| template<typename T> |
| __device__ __forceinline__ |
| T load(T *t) |
| { |
| return __ldg(t); //texture load |
| } |
| |
| template<int TILE_SIZE, typename T>__device__ void loadT(T * out, const T * in); |
| |
| //generic version (inefficient) |
| template<int TILE_SIZE, typename T> |
| __device__ __forceinline__ |
| void loadT(T * out, const T * in) { |
| #pragma unroll |
| for(int i=0;i<TILE_SIZE;i++) { |
| out[i]=__ldg(in+i); |
| } |
| } |
| |
| //Vector loads |
| template<> |
| __device__ __forceinline__ |
| void loadT<2,float>(float * out, const float * in) { |
| *reinterpret_cast<float2*>(out)=load(reinterpret_cast<const float2*>(in)); |
| } |
| template<> |
| __device__ __forceinline__ |
| void loadT<4,float>(float * out, const float * in) { |
| *reinterpret_cast<float4*>(out)=load(reinterpret_cast<const float4*>(in)); |
| } |
| |
| //static __device__ __forceinline__ float __internal_fast_rsqrtf(float a) |
| //{ |
| // float r; |
| // asm ("rsqrt.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a)); |
| // return r; |
| //} |
| |
| //computes the forces between tiles i and j, adds the change in force to xi,yi,zi |
| template<int TX, int TY> |
| __device__ __forceinline__ void computeForces(POSVEL_T xxi[], POSVEL_T yyi[], POSVEL_T zzi[], |
| POSVEL_T xxj[], POSVEL_T yyj[], POSVEL_T zzj[], POSVEL_T massj[], |
| POSVEL_T xi[], POSVEL_T yi[], POSVEL_T zi[], |
| POSVEL_T ma0, POSVEL_T ma1, POSVEL_T ma2, POSVEL_T ma3, POSVEL_T ma4, POSVEL_T ma5, |
| POSVEL_T mp_rsm2, POSVEL_T fsrrmax2) { |
| |
| #pragma unroll |
| for(int i=0;i<TY;i++) { |
| #pragma unroll |
| for(int j=0;j<TX;j++) { |
| POSVEL_T dxc = xxj[j] - xxi[i]; //1 FADD |
| POSVEL_T dyc = yyj[j] - yyi[i]; //1 FADD |
| POSVEL_T dzc = zzj[j] - zzi[i]; //1 FADD |
| |
| POSVEL_T r2 = dxc * dxc + dyc * dyc + dzc * dzc; //1 FMUL 2 FMA |
| POSVEL_T v=r2+mp_rsm2; //1 FADD |
| POSVEL_T v3=v*v*v; //2 FMUL |
| |
| POSVEL_T f = __frsqrt_rn(v3); //1 MUFU, |
| // MDS: Should ask someone why this line is dangling |
| // - ( ma0 + r2*(ma1 + r2*(ma2 + r2*(ma3 + r2*(ma4 + r2*ma5))))); //5 FMA, 1 FADD |
| //#define BUG |
| #ifndef BUG |
| f*=massj[j]*(r2<fsrrmax2 && r2>0.0f); //2 FMUL, 1 FSETP, 1 FCMP |
| #else |
| f*=massj[j]; //1 FMUL |
| f*=(r2<fsrrmax2 && r2>0.0f); //1 FMUL, 1 FSETP, 1 FCMP |
| #endif |
| |
| xi[i] = xi[i] + f * dxc; //1 FMA |
| yi[i] = yi[i] + f * dyc; //1 FMA |
| zi[i] = zi[i] + f * dzc; //1 FMA |
| } |
| } |
| } |
| |
| //loads a tile from memory. Use checkBounds and loadMass to disable bounds check or mass load at compile time |
| template <bool checkBounds, bool loadMass, int TILE_SIZE> |
| __device__ __forceinline__ |
| void loadTile(int i, int bounds, |
| const POSVEL_T* xx, const POSVEL_T* yy, const POSVEL_T* zz, const POSVEL_T* mass, |
| POSVEL_T xxi[], POSVEL_T yyi[], POSVEL_T zzi[], POSVEL_T massi[]) { |
| if(checkBounds) { |
| #pragma unroll |
| for(int64 u=0;u<TILE_SIZE;u++) { |
| int64 idx=TILE_SIZE*i+u; // 1 IMAD |
| |
| #if 1 |
| bool cond=idx<bounds; |
| xxi[u] = (cond) ? load(xx+idx) : 0.0f; // 1 ISETP, 1 LDG, 2 IMAD, 1 MOV |
| yyi[u] = (cond) ? load(yy+idx) : 0.0f; // 1 ISETP, 1 LDG, 2 IMAD, 1 MOV |
| zzi[u] = (cond) ? load(zz+idx) : 0.0f; // 1 ISETP, 1 LDG, 2 IMAD, 1 MOV |
| if(loadMass) massi[u] = (cond) ? load(mass+idx) : 0.0f; // 1 ISETP, 1 LDG, 2 IMAD, 1 MOV |
| #else |
| massi[u] = 0.0f; //1 MOV |
| if(idx<bounds) { //1 ISETP, 1 BRA |
| xxi[u] = load(xx+idx); //1 LDG, 2 IMAD |
| yyi[u] = load(yy+idx); //1 LDG, 2 IMAD |
| zzi[u] = load(zz+idx); //1 LDG, 2 IMAD |
| if(loadMass) massi[u] = load(mass+idx); //1 LDG, 2 IMAD |
| } |
| #endif |
| } |
| } else { |
| |
| int idx=TILE_SIZE*i; |
| loadT<TILE_SIZE>(xxi,xx+idx); //1 LDG, 2 IMAD |
| loadT<TILE_SIZE>(yyi,yy+idx); //1 LDG, 2 IMAD |
| loadT<TILE_SIZE>(zzi,zz+idx); //1 LDG, 2 IMAD |
| if(loadMass) loadT<TILE_SIZE>(massi,mass+idx); //1 LDG, 2 IMAD |
| } |
| } |
| |
| //applies the force in xi,yi,zi to update vx, vy, vz |
| //use checkBounds to disable bounds checking at compile time |
| template <bool checkBounds, int TILE_SIZE> |
| __device__ __forceinline__ |
| void applyForce(int i, int bounds,POSVEL_T fcoeff, |
| const POSVEL_T xi[], const POSVEL_T yi[], const POSVEL_T zi[], |
| POSVEL_T *vx, POSVEL_T *vy, POSVEL_T *vz) { |
| #pragma unroll |
| for(int u=0;u<TILE_SIZE;u++) { |
| int idx=TILE_SIZE*i+u; //1 IMAD |
| |
| if(!checkBounds || idx<bounds) |
| { //1 ISETP |
| atomicWarpReduceAndUpdate(vx+idx,fcoeff * xi[u]); //2 IMAD, 6 FADD |
| atomicWarpReduceAndUpdate(vy+idx,fcoeff * yi[u]); //2 IMAD, 6 FADD |
| atomicWarpReduceAndUpdate(vz+idx,fcoeff * zi[u]); //2 IMAD, 6 FADD |
| } |
| } |
| } |
| |
| //Tell the compiler how many blocks we expect to be active. |
| //This gives the compiler a better idea of how many registers to use. The second number is tunable. |
| __launch_bounds__(BLOCKX*BLOCKY,7) |
| __global__ |
| void Step10_cuda_kernel(int count, int count1, |
| const POSVEL_T* xx, const POSVEL_T* yy, |
| const POSVEL_T* zz, const POSVEL_T* mass, |
| const POSVEL_T* xx1, const POSVEL_T* yy1, |
| const POSVEL_T* zz1, const POSVEL_T* mass1, |
| POSVEL_T* vx, POSVEL_T* vy, |
| POSVEL_T* vz, POSVEL_T fsrrmax2, POSVEL_T mp_rsm2, POSVEL_T fcoeff) |
| { |
| const POSVEL_T ma0 = 0.269327, ma1 = -0.0750978, ma2 = 0.0114808, ma3 = -0.00109313, ma4 = 0.0000605491, ma5 = -0.00000147177; |
| |
| //Register arrays to hold tiles of data |
| POSVEL_T xxi[TILEY]; |
| POSVEL_T yyi[TILEY]; |
| POSVEL_T zzi[TILEY]; |
| POSVEL_T xxj[TILEX]; |
| POSVEL_T yyj[TILEX]; |
| POSVEL_T zzj[TILEX]; |
| POSVEL_T massj[TILEX]; |
| |
| //loop over interior region and calculate forces. |
| //for each tile i |
| for(int i=hipBlockIdx_y*hipBlockDim_y+hipThreadIdx_y;i<count/TILEY;i+=hipBlockDim_y*hipGridDim_y) //1 ISETP |
| { |
| POSVEL_T xi[TILEY]={0}; //TILEY MOV |
| POSVEL_T yi[TILEY]={0}; //TILEY MOV |
| POSVEL_T zi[TILEY]={0}; //TILEY MOV |
| |
| //load tile i,mass and bounds check are not needed |
| loadTile<false,false,TILEY>(i,count,xx,yy,zz,NULL,xxi,yyi,zzi,NULL); |
| |
| //for each tile j |
| for (int j=hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x;j<count1/TILEX;j+=hipBlockDim_x*hipGridDim_x) //1 ISETP |
| { |
| //load tile j, bounds check is not needed |
| loadTile<false,true,TILEX>(j,count1,xx1,yy1,zz1,mass1,xxj,yyj,zzj,massj); |
| |
| //compute forces between tile i and tile j |
| computeForces<TILEX,TILEY>(xxi,yyi,zzi,xxj,yyj,zzj,massj,xi,yi,zi,ma0,ma1,ma2,ma3,ma4,ma5,mp_rsm2,fsrrmax2); |
| } |
| |
| //process remaining elements at the end, use TILEX=1 |
| for (int j=count1/TILEX*TILEX+hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x;j<count1;j+=hipBlockDim_x*hipGridDim_x) //1 ISETP |
| { |
| //load tile j, bounds check is needed, mass is needed |
| loadTile<true,true,1>(j,count1,xx1,yy1,zz1,mass1,xxj,yyj,zzj,massj); |
| |
| //compute forces between tile i and tile j |
| computeForces<1,TILEY>(xxi,yyi,zzi,xxj,yyj,zzj,massj,xi,yi,zi,ma0,ma1,ma2,ma3,ma4,ma5,mp_rsm2,fsrrmax2); |
| } |
| |
| //apply the force we have calculated above, bounds check is not needed |
| applyForce<false,TILEY>(i,count,fcoeff,xi,yi,zi,vx,vy,vz); |
| } |
| |
| //At this point we have computed almost all interactions. |
| //However we still need to add contributions for particles at the end |
| |
| #if 1 |
| //process ramining elements in set TILEY=1 |
| //for each tile i |
| for(int i=count/TILEY*TILEY+hipBlockIdx_y*hipBlockDim_y+hipThreadIdx_y;i<count;i+=hipBlockDim_y*hipGridDim_y) //1 ISETP |
| { |
| POSVEL_T xi[1]={0}; //1 MOV |
| POSVEL_T yi[1]={0}; //1 MOV |
| POSVEL_T zi[1]={0}; //1 MOV |
| |
| //load xxi, yyi, zzi tiles, mass is not needed, bounds check is needed |
| loadTile<true,false,1>(i,count,xx,yy,zz,NULL,xxi,yyi,zzi,NULL); |
| |
| //for each tile j |
| for (int j=hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x;j<count1/TILEX;j+=hipBlockDim_x*hipGridDim_x) //1 ISETP |
| { |
| //load tile j, bounds check is not needed |
| loadTile<false,true,TILEX>(j,count1,xx1,yy1,zz1,mass1,xxj,yyj,zzj,massj); |
| |
| //compute forces between tile i and tile j |
| computeForces<TILEX,1>(xxi,yyi,zzi,xxj,yyj,zzj,massj,xi,yi,zi,ma0,ma1,ma2,ma3,ma4,ma5,mp_rsm2,fsrrmax2); |
| } |
| |
| //process remaining elements at the end, use TILEX=1 |
| for (int j=count1/TILEX*TILEX+hipBlockIdx_x*hipBlockDim_x+hipThreadIdx_x;j<count1;j+=hipBlockDim_x*hipGridDim_x) //1 ISETP |
| { |
| //load tile j, bounds check is needed, mass is needed |
| loadTile<true,true,1>(j,count1,xx1,yy1,zz1,mass1,xxj,yyj,zzj,massj); |
| |
| //compute forces between tile i and tile j |
| computeForces<1,1>(xxi,yyi,zzi,xxj,yyj,zzj,massj,xi,yi,zi,ma0,ma1,ma2,ma3,ma4,ma5,mp_rsm2,fsrrmax2); |
| } |
| |
| applyForce<true,1>(i,count,fcoeff,xi,yi,zi,vx,vy,vz); |
| } |
| #endif |
| |
| } |
| |
| |
| |
| #endif |
| |
| #ifdef __bgq__ |
| extern "C" Step16_int( int count1, float xxi, float yyi, float zzi, float fsrrmax2, float mp_rsm2, const float *xx1, const float *yy1, const float *zz1,const float *mass1, float *ax, float *ay, float *az ); |
| #endif |
| |
| static inline void nbody1(ID_T count, ID_T count1, const POSVEL_T* xx, const POSVEL_T* yy, |
| const POSVEL_T* zz, const POSVEL_T* mass, |
| const POSVEL_T* xx1, const POSVEL_T* yy1, |
| const POSVEL_T* zz1, const POSVEL_T* mass1, |
| POSVEL_T* vx, POSVEL_T* vy, POSVEL_T* vz, |
| ForceLaw *fl, float fcoeff, float fsrrmax, float rsm |
| #ifdef __HIPCC__ |
| , hipStream_t stream |
| #endif |
| ) |
| { |
| POSVEL_T fsrrmax2 = fsrrmax*fsrrmax; |
| POSVEL_T rsm2 = rsm*rsm; |
| |
| #ifdef __bgq__ |
| float ax = 0.0f, ay = 0.0f, az = 0.0f; |
| |
| for (int i = 0; i < count; ++i) |
| { |
| |
| Step16_int ( count1, xx[i],yy[i],zz[i], fsrrmax2,rsm2,xx1,yy1,zz1,mass1, &ax, &ay, &az ); |
| |
| vx[i] = vx[i] + ax * fcoeff; |
| vy[i] = vy[i] + ay * fcoeff; |
| vz[i] = vz[i] + az * fcoeff; |
| } |
| |
| #else |
| |
| #ifdef __HIPCC__ |
| //const int MAXX=64; |
| //const int MAXY=64; |
| |
| dim3 threads(BLOCKX,BLOCKY); |
| int blocksX=(count1+threads.x-1)/threads.x; |
| int blocksY=(count+threads.y-1)/threads.y; |
| dim3 blocks( min(blocksX,MAXX), min(blocksY,MAXY)); |
| |
| //call kernel |
| |
| cudaCheckError(); |
| //printf("count: %d, count1: %d\n", count,count1); |
| #if 0 |
| checkCudaPtr(xx,"xx"); |
| checkCudaPtr(yy,"yy"); |
| checkCudaPtr(zz,"zz"); |
| checkCudaPtr(mass,"mass"); |
| checkCudaPtr(xx1,"xx1"); |
| checkCudaPtr(yy1,"yy1"); |
| checkCudaPtr(zz1,"zz1"); |
| checkCudaPtr(mass1,"mass1"); |
| checkCudaPtr(vx,"vx"); |
| checkCudaPtr(vy,"vy"); |
| checkCudaPtr(vz,"vz"); |
| #endif |
| hipLaunchKernelGGL(Step10_cuda_kernel, dim3(blocks), dim3(threads), 0, stream, count,count1,xx,yy,zz,mass,xx1,yy1,zz1,mass1, vx, vy, vz, fsrrmax2, rsm2, fcoeff); |
| cudaCheckError(); |
| |
| //hipDeviceSynchronize(); |
| //exit(0); |
| #else |
| |
| for (int i = 0; i < count; ++i) |
| for (int j = 0; j < count1; ++j) { |
| POSVEL_T dx = xx1[j] - xx[i]; |
| POSVEL_T dy = yy1[j] - yy[i]; |
| POSVEL_T dz = zz1[j] - zz[i]; |
| POSVEL_T dist2 = dx*dx + dy*dy + dz*dz; |
| POSVEL_T f_over_r = mass[i]*mass1[j] * fl->f_over_r(dist2); |
| |
| POSVEL_T updateq = 1.0; |
| updateq *= (dist2 < fsrrmax2); |
| |
| vx[i] += updateq*fcoeff*f_over_r*dx; |
| vy[i] += updateq*fcoeff*f_over_r*dy; |
| vz[i] += updateq*fcoeff*f_over_r*dz; |
| } |
| #endif //end __HIPCC__ |
| |
| |
| #endif //end __bgq__ |
| } |
| |
| |
| static inline ID_T partition(ID_T n, |
| POSVEL_T* xx, POSVEL_T* yy, POSVEL_T* zz, |
| POSVEL_T* vx, POSVEL_T* vy, POSVEL_T* vz, |
| POSVEL_T* mass, POSVEL_T* phi, |
| ID_T* id, MASK_T* mask, POSVEL_T pv |
| ) |
| { |
| float t0, t1, t2, t3, t4, t5, t6, t7; |
| int32_t is, i, j; |
| long i0; |
| uint16_t i1; |
| |
| int idx[n]; |
| |
| is = 0; |
| for ( i = 0; i < n; i = i + 1 ) |
| { |
| if (xx[i] < pv) |
| { |
| idx[is] = i; |
| is = is + 1; |
| } |
| } |
| |
| #pragma unroll (4) |
| for ( j = 0; j < is; j++ ) |
| { |
| i = idx[j]; |
| |
| t6 = mass[i]; mass[i] = mass[j]; mass[j] = t6; |
| t7 = phi [i]; phi [i] = phi [j]; phi [j] = t7; |
| i1 = mask[i]; mask[i] = mask[j]; mask[j] = i1; |
| i0 = id [i]; id [i] = id [j]; id [j] = i0; |
| } |
| |
| #pragma unroll (4) |
| for ( j = 0; j < is; j++ ) |
| { |
| i = idx[j]; |
| |
| t0 = xx[i]; xx[i] = xx[j]; xx[j] = t0; |
| t1 = yy[i]; yy[i] = yy[j]; yy[j] = t1; |
| t2 = zz[i]; zz[i] = zz[j]; zz[j] = t2; |
| t3 = vx[i]; vx[i] = vx[j]; vx[j] = t3; |
| t4 = vy[i]; vy[i] = vy[j]; vy[j] = t4; |
| t5 = vz[i]; vz[i] = vz[j]; vz[j] = t5; |
| } |
| |
| return is; |
| } |
| |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::createRCBForceSubtree(int d, ID_T tl, ID_T tlcl, ID_T tlcr) |
| { |
| POSVEL_T *x1, *x2, *x3; |
| switch (d) { |
| case 0: |
| x1 = xx; |
| x2 = yy; |
| x3 = zz; |
| break; |
| case 1: |
| x1 = yy; |
| x2 = zz; |
| x3 = xx; |
| break; |
| /*case 2*/ default: |
| x1 = zz; |
| x2 = xx; |
| x3 = yy; |
| break; |
| } |
| |
| #ifdef __bgq__ |
| int tid = 0; |
| |
| #endif |
| const bool geoSplit = false; |
| POSVEL_T split = geoSplit ? (tree[tl].xmax[d]+tree[tl].xmin[d])/2 : tree[tl].xc[d]; |
| ID_T is = ::partition(tree[tl].count, x1 + tree[tl].offset, x2 + tree[tl].offset, x3 + tree[tl].offset, |
| vx + tree[tl].offset, vy + tree[tl].offset, vz + tree[tl].offset, |
| mass + tree[tl].offset, phi + tree[tl].offset, |
| id + tree[tl].offset, mask + tree[tl].offset, split |
| ); |
| |
| if (is == 0 || is == tree[tl].count) { |
| return; |
| } |
| |
| tree[tlcl].count = is; |
| tree[tlcr].count = tree[tl].count - tree[tlcl].count; |
| |
| if (tree[tlcl].count > 0) { |
| tree[tl].cl = tlcl; |
| tree[tlcl].offset = tree[tl].offset; |
| tree[tlcl].xmax[d] = split; |
| |
| createRCBForceTreeInParallel(tlcl); |
| } |
| |
| if (tree[tlcr].count > 0) { |
| tree[tl].cr = tlcr; |
| tree[tlcr].offset = tree[tl].offset + tree[tlcl].count; |
| tree[tlcr].xmin[d] = split; |
| |
| createRCBForceTreeInParallel(tlcr); |
| } |
| } |
| |
| // This is basically the algorithm from (Gafton and Rosswog, 2011). |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::createRCBForceTreeInParallel(ID_T tl) |
| { |
| ID_T cnt = tree[tl].count; |
| ID_T off = tree[tl].offset; |
| |
| // Compute the center-of-mass coordinates (and recompute the min/max) |
| ::cm(cnt, xx + off, yy + off, zz + off, mass + off, |
| tree[tl].xmin, tree[tl].xmax, tree[tl].xc); |
| |
| if (cnt <= nDirect) { |
| // The pseudoparticles |
| tree[tl].tdr = ppContract*::pptdr(tree[tl].xmin, tree[tl].xmax, tree[tl].xc); |
| memset(tree[tl].ppm, 0, sizeof(POSVEL_T)*TDPTS); |
| if (cnt > TDPTS) { // Otherwise, the pseudoparticles are never used |
| POSVEL_T ppx[TDPTS], ppy[TDPTS], ppz[TDPTS]; |
| pppts<TDPTS>(tree[tl].tdr, tree[tl].xc, ppx, ppy, ppz); |
| pp<TDPTS>(cnt, xx + off, yy + off, zz + off, mass + off, tree[tl].xc, |
| ppx, ppy, ppz, tree[tl].ppm, tree[tl].tdr); |
| } |
| |
| return; |
| } |
| |
| // Index of the right and left child levels |
| ID_T tlcl, tlcr; |
| { |
| tlcl = tree.size(); |
| tlcr = tlcl+1; |
| size_t newSize = tlcr+1; |
| tree.resize(newSize); |
| } |
| memset(&tree[tlcl], 0, sizeof(TreeNode)*2); |
| |
| // Both children have similar bounding boxes to the current node (the |
| // parent), so copy the bounding box here, and then overwrite the changed |
| // coordinate later. |
| for (int i = 0; i < DIMENSION; ++i) { |
| tree[tlcl].xmin[i] = tree[tl].xmin[i]; |
| tree[tlcr].xmin[i] = tree[tl].xmin[i]; |
| tree[tlcl].xmax[i] = tree[tl].xmax[i]; |
| tree[tlcr].xmax[i] = tree[tl].xmax[i]; |
| } |
| |
| // Split the longest edge at the center of mass. |
| POSVEL_T xlen[DIMENSION]; |
| for (int i = 0; i < DIMENSION; ++i) { |
| xlen[i] = tree[tl].xmax[i] - tree[tl].xmin[i]; |
| } |
| |
| int d; |
| if (xlen[0] > xlen[1] && xlen[0] > xlen[2]) { |
| d = 0; // Split in the x direction |
| } |
| else if (xlen[1] > xlen[2]) { |
| d = 1; // Split in the y direction |
| } |
| else { |
| d = 2; // Split in the z direction |
| } |
| |
| createRCBForceSubtree(d, tl, tlcl, tlcr); |
| |
| // Compute the pseudoparticles based on those of the children |
| POSVEL_T ppx[TDPTS], ppy[TDPTS], ppz[TDPTS]; |
| tree[tl].tdr = ppContract*::pptdr(tree[tl].xmin, tree[tl].xmax, tree[tl].xc); |
| pppts<TDPTS>(tree[tl].tdr, tree[tl].xc, ppx, ppy, ppz); |
| memset(tree[tl].ppm, 0, sizeof(POSVEL_T)*TDPTS); |
| |
| if (tree[tlcl].count > 0) { |
| if (tree[tlcl].count <= TDPTS) { |
| ID_T offc = tree[tlcl].offset; |
| pp<TDPTS>(tree[tlcl].count, xx + offc, yy + offc, zz + offc, mass + offc, |
| tree[tl].xc, ppx, ppy, ppz, tree[tl].ppm, tree[tl].tdr); |
| } else { |
| POSVEL_T ppxc[TDPTS], ppyc[TDPTS], ppzc[TDPTS]; |
| pppts<TDPTS>(tree[tlcl].tdr, tree[tlcl].xc, ppxc, ppyc, ppzc); |
| pp<TDPTS>(TDPTS, ppxc, ppyc, ppzc, tree[tlcl].ppm, tree[tl].xc, |
| ppx, ppy, ppz, tree[tl].ppm, tree[tl].tdr); |
| } |
| } |
| if (tree[tlcr].count > 0) { |
| if (tree[tlcr].count <= TDPTS) { |
| ID_T offc = tree[tlcr].offset; |
| pp<TDPTS>(tree[tlcr].count, xx + offc, yy + offc, zz + offc, mass + offc, |
| tree[tl].xc, ppx, ppy, ppz, tree[tl].ppm, tree[tl].tdr); |
| } else { |
| POSVEL_T ppxc[TDPTS], ppyc[TDPTS], ppzc[TDPTS]; |
| pppts<TDPTS>(tree[tlcr].tdr, tree[tlcr].xc, ppxc, ppyc, ppzc); |
| pp<TDPTS>(TDPTS, ppxc, ppyc, ppzc, tree[tlcr].ppm, tree[tl].xc, |
| ppx, ppy, ppz, tree[tl].ppm, tree[tl].tdr); |
| } |
| } |
| } |
| |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::createRCBForceTree() |
| { |
| // The top tree is the entire box |
| tree.resize(1); |
| memset(&tree[0], 0, sizeof(TreeNode)); |
| |
| tree[0].count = particleCount; |
| tree[0].offset = 0; |
| |
| for (int i = 0; i < DIMENSION; ++i) { |
| tree[0].xmin[i] = minRange[i]; |
| tree[0].xmax[i] = maxRange[i]; |
| } |
| |
| createRCBForceTreeInParallel(); |
| } |
| |
| |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::calcInternodeForce(ID_T tl, |
| const std::vector<ID_T> &parents) { |
| |
| POSVEL_T fsrrmax2 = fsrrmax*fsrrmax; |
| const TreeNode* tree_ = &tree[0]; |
| |
| int tid = 0; |
| |
| std::vector<ID_T> &q = iq[tid]; |
| q.clear(); |
| q.push_back(0); |
| |
| |
| POSVEL_T *nx=nx_v+tid*VMAX; |
| POSVEL_T *ny=ny_v+tid*VMAX; |
| POSVEL_T *nz=nz_v+tid*VMAX; |
| POSVEL_T *nm=nm_v+tid*VMAX; |
| |
| #ifdef __HIPCC__ |
| //Adjust pointers to this threads workspace |
| POSVEL_T *d_nx=d_nx_v+tid*VMAX; |
| POSVEL_T *d_ny=d_ny_v+tid*VMAX; |
| POSVEL_T *d_nz=d_nz_v+tid*VMAX; |
| POSVEL_T *d_nm=d_nm_v+tid*VMAX; |
| int size=ALIGNY(nDirect); |
| POSVEL_T *d_xxl=d_xx+tid*size; |
| POSVEL_T *d_yyl=d_yy+tid*size; |
| POSVEL_T *d_zzl=d_zz+tid*size; |
| POSVEL_T *d_massl=d_mass+tid*size; |
| POSVEL_T *d_vxl=d_vx+tid*size; |
| POSVEL_T *d_vyl=d_vy+tid*size; |
| POSVEL_T *d_vzl=d_vz+tid*size; |
| |
| hipEvent_t& event=event_v[tid]; |
| hipStream_t& stream=stream_v[tid]; |
| hipEventSynchronize(event); //wait for transfers from previous call to finish before overwriting nx,ny,nz,nm |
| cudaCheckError(); |
| #endif |
| |
| // The interaction list. |
| int SIZE = 0; // current size of these arrays |
| |
| while (!q.empty()) { |
| ID_T tln = q.back(); |
| q.pop_back(); |
| |
| // We should not interact with our own parents. |
| if (tln < tl) { |
| bool isParent = std::binary_search(parents.begin(), parents.end(), tln); |
| if (isParent) { |
| ID_T tlncr = tree_[tln].cr; |
| ID_T tlncl = tree_[tln].cl; |
| |
| if (tlncl != tl && tlncl > 0 && tree_[tlncl].count > 0) { |
| q.push_back(tlncl); |
| } |
| if (tlncr != tl && tlncr > 0 && tree_[tlncr].count > 0) { |
| q.push_back(tlncr); |
| } |
| |
| continue; |
| } |
| } |
| |
| // Is this node have a small enough opening angle to interact with? |
| POSVEL_T dx = tree_[tln].xc[0] - tree_[tl].xc[0]; |
| POSVEL_T dy = tree_[tln].xc[1] - tree_[tl].xc[1]; |
| POSVEL_T dz = tree_[tln].xc[2] - tree_[tl].xc[2]; |
| POSVEL_T dist2 = dx*dx + dy*dy + dz*dz; |
| |
| POSVEL_T sx = tree_[tln].xmax[0]-tree_[tln].xmin[0]; |
| POSVEL_T sy = tree_[tln].xmax[1]-tree_[tln].xmin[1]; |
| POSVEL_T sz = tree_[tln].xmax[2]-tree_[tln].xmin[2]; |
| POSVEL_T l2 = std::min(sx*sx, std::min(sy*sy, sz*sz)); // under-estimate |
| |
| POSVEL_T dtt2 = dist2*tanOpeningAngle*tanOpeningAngle; |
| bool looksBig; |
| // l2/dist2 is really tan^2 theta, for small theta, tan(theta) ~ theta |
| if (l2 > dtt2) { |
| // the under-estimate is too big, so this is definitely too big |
| looksBig = true; |
| } else { |
| // There are 8 corner points of the remote node, and the maximum angular |
| // size will be from one of those points to its opposite points. So there |
| // are 8 vector dot products to compute to determine the maximum angular |
| // size at any given reference point. (do we need to do this for each point |
| // in leaf node, or will the c.m. point be sufficient?). |
| looksBig = false; |
| for (int i = 0; i < 2; ++i) |
| for (int j = 0; j < 2; ++j) { |
| POSVEL_T x1 = (i == 0 ? tree_[tln].xmin : tree_[tln].xmax)[0] - tree_[tl].xc[0]; |
| POSVEL_T y1 = (j == 0 ? tree_[tln].xmin : tree_[tln].xmax)[1] - tree_[tl].xc[1]; |
| POSVEL_T z1 = tree_[tln].xmin[2] - tree_[tl].xc[2]; |
| |
| POSVEL_T x2 = (i == 0 ? tree_[tln].xmax : tree_[tln].xmin)[0] - tree_[tl].xc[0]; |
| POSVEL_T y2 = (j == 0 ? tree_[tln].xmax : tree_[tln].xmin)[1] - tree_[tl].xc[1]; |
| POSVEL_T z2 = tree_[tln].xmax[2] - tree_[tl].xc[2]; |
| |
| const bool useRealOA = false; |
| if (useRealOA) { |
| // |a x b| = a*b*sin(theta) |
| POSVEL_T cx = y1*z2 - z1*y2; |
| POSVEL_T cy = z1*x2 - x1*z2; |
| POSVEL_T cz = x1*y2 - y1*x2; |
| if ((cx*cx + cy*cy + cz*cz) > sinOpeningAngle*sinOpeningAngle* |
| (x1*x1 + y1*y1 + z1*z1)*(x2*x2 + y2*y2 + z2*z2) |
| ) { |
| looksBig = true; |
| break; |
| } |
| } else { |
| // Instead of using the real opening angle, use the tan approximation; this is |
| // better than the opening-angle b/c it incorporates depth information. |
| POSVEL_T ddx = x1 - x2, ddy = y1 - y2, ddz = z1 - z2; |
| POSVEL_T dh2 = ddx*ddx + ddy*ddy + ddz*ddz; |
| if (dh2 > dtt2) { |
| looksBig = true; |
| break; |
| } |
| } |
| } |
| } |
| |
| if (!looksBig) { |
| if (dist2 > fsrrmax2) { |
| // We could interact with this node, but it is too far away to make |
| // any difference, so it will be skipped, along with all of its |
| // children. |
| continue; |
| } |
| |
| // This node has fewer particles than pseudo particles, so just use the |
| // particles that are actually there. |
| if (tree_[tln].count <= TDPTS) { |
| ID_T offn = tree_[tln].offset; |
| ID_T cntn = tree_[tln].count; |
| |
| int start = SIZE; |
| SIZE = SIZE + cntn; |
| assert( SIZE < VMAX ); |
| |
| for ( int i = 0; i < cntn; ++i) { |
| nx[start + i] = xx[offn + i]; |
| ny[start + i] = yy[offn + i]; |
| nz[start + i] = zz[offn + i]; |
| nm[start + i] = mass[offn + i]; |
| } |
| |
| continue; |
| } |
| |
| // Interact the particles in this node with the pseudoparticles of the |
| // other node. |
| int start = SIZE; |
| SIZE = SIZE + TDPTS; |
| assert( SIZE < VMAX ); |
| |
| pppts<TDPTS>(tree_[tln].tdr, tree_[tln].xc, &nx[start], &ny[start], &nz[start]); |
| for ( int i = 0; i < TDPTS; ++i) { |
| nm[start + i] = tree_[tln].ppm[i]; |
| } |
| |
| continue; |
| } else if (tree_[tln].cr == 0 && tree_[tln].cl == 0) { |
| // This is a leaf node with which we must interact. |
| ID_T offn = tree_[tln].offset; |
| ID_T cntn = tree_[tln].count; |
| |
| int start = SIZE; |
| SIZE = SIZE + cntn; |
| assert( SIZE < VMAX ); |
| |
| for ( int i = 0; i < cntn; ++i) { |
| nx[start + i] = xx[offn + i]; |
| ny[start + i] = yy[offn + i]; |
| nz[start + i] = zz[offn + i]; |
| nm[start + i] = mass[offn + i]; |
| } |
| |
| continue; |
| } |
| |
| // This other node is not a leaf, but has too large an opening angle |
| // for an approx. interaction: queue its children. |
| |
| ID_T tlncr = tree_[tln].cr; |
| ID_T tlncl = tree_[tln].cl; |
| |
| if (tlncl > 0 && tree_[tlncl].count > 0) { |
| bool close = true; |
| for (int i = 0; i < DIMENSION; ++i) { |
| POSVEL_T dist = 0; |
| if (tree_[tl].xmax[i] < tree_[tlncl].xmin[i]) { |
| dist = tree_[tlncl].xmin[i] - tree_[tl].xmax[i]; |
| } else if (tree_[tl].xmin[i] > tree_[tlncl].xmax[i]) { |
| dist = tree_[tl].xmin[i] - tree_[tlncl].xmax[i]; |
| } |
| |
| if (dist > fsrrmax) { |
| close = false; |
| break; |
| } |
| } |
| |
| if (close) q.push_back(tlncl); |
| } |
| if (tlncr > 0 && tree_[tlncr].count > 0) { |
| bool close = true; |
| for (int i = 0; i < DIMENSION; ++i) { |
| POSVEL_T dist = 0; |
| if (tree_[tl].xmax[i] < tree_[tlncr].xmin[i]) { |
| dist = tree_[tlncr].xmin[i] - tree_[tl].xmax[i]; |
| } else if (tree_[tl].xmin[i] > tree_[tlncr].xmax[i]) { |
| dist = tree_[tl].xmin[i] - tree_[tlncr].xmax[i]; |
| } |
| |
| if (dist > fsrrmax) { |
| close = false; |
| break; |
| } |
| } |
| |
| if (close) q.push_back(tlncr); |
| } |
| } |
| |
| ID_T off = tree_[tl].offset; |
| ID_T cnt = tree_[tl].count; |
| |
| // Add self interactions... |
| int start = SIZE; |
| SIZE = SIZE + cnt; |
| assert( SIZE < VMAX ); |
| |
| for ( int i = 0; i < cnt; ++i) { |
| nx[start + i] = xx[off + i]; |
| ny[start + i] = yy[off + i]; |
| nz[start + i] = zz[off + i]; |
| nm[start + i] = mass[off + i]; |
| } |
| |
| #ifdef __HIPCC__ |
| hipMemcpyAsync(d_nx,nx,sizeof(POSVEL_T)*SIZE,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_ny,ny,sizeof(POSVEL_T)*SIZE,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_nz,nz,sizeof(POSVEL_T)*SIZE,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_nm,nm,sizeof(POSVEL_T)*SIZE,hipMemcpyHostToDevice,stream); |
| hipEventRecord(event,stream); //mark when transfers have finished |
| cudaCheckError(); |
| hipDeviceSynchronize(); |
| |
| |
| hipMemcpyAsync(d_xxl,xx+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_yyl,yy+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_zzl,zz+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_massl,mass+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_vxl,vx+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_vyl,vy+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| hipMemcpyAsync(d_vzl,vz+off,sizeof(POSVEL_T)*cnt,hipMemcpyHostToDevice,stream); |
| cudaCheckError(); |
| hipDeviceSynchronize(); |
| #endif |
| |
| // Process the interaction list... |
| #ifdef __HIPCC__ |
| ::nbody1(cnt, SIZE, d_xxl, d_yyl, d_zzl, d_massl, d_nx, d_ny, d_nz, d_nm, d_vxl, d_vyl, d_vzl, m_fl, m_fcoeff, fsrrmax, rsm, stream); |
| hipDeviceSynchronize(); |
| #else |
| ::nbody1(cnt, SIZE, xx + off, yy + off, zz + off, mass + off, nx, ny, nz, nm, vx + off, vy + off, vz + off, m_fl, m_fcoeff, fsrrmax, rsm); |
| #endif |
| |
| #ifdef __HIPCC__ |
| //transfer up vx vy vz |
| hipMemcpyAsync(vx+off,d_vxl,sizeof(POSVEL_T)*cnt,hipMemcpyDeviceToHost,stream); |
| hipMemcpyAsync(vy+off,d_vyl,sizeof(POSVEL_T)*cnt,hipMemcpyDeviceToHost,stream); |
| hipMemcpyAsync(vz+off,d_vzl,sizeof(POSVEL_T)*cnt,hipMemcpyDeviceToHost,stream); |
| cudaCheckError(); |
| hipDeviceSynchronize(); |
| #endif |
| |
| } |
| |
| // Iterate through the tree nodes, for each leaf node, start a task. |
| // That task iterates through the tree nodes, skipping any node (and all |
| // of its children) if all corners are too far away. Then it compares the |
| // opening angle. |
| template <int TDPTS> |
| void RCBForceTree<TDPTS>::calcInternodeForces() |
| { |
| |
| std::vector<ID_T> q(1, 0); |
| std::vector<ID_T> parents; |
| while (!q.empty()) { |
| ID_T tl = q.back(); |
| if (tree[tl].cr == 0 && tree[tl].cl == 0) { |
| // This is a leaf node. |
| q.pop_back(); |
| |
| bool inside = true; |
| for (int i = 0; i < DIMENSION; ++i) { |
| inside &= (tree[tl].xmax[i] < maxForceRange[i] && tree[tl].xmax[i] > minForceRange[i]) || |
| (tree[tl].xmin[i] < maxForceRange[i] && tree[tl].xmin[i] > minForceRange[i]); |
| } |
| |
| if (inside) { |
| calcInternodeForce(tl, parents); |
| } |
| } else if (parents.size() > 0 && parents.back() == tl) { |
| // This is second time here; we've done with all children. |
| parents.pop_back(); |
| q.pop_back(); |
| } else { |
| // This is the first time at this parent node, queue the children. |
| if (tree[tl].cl > 0) q.push_back(tree[tl].cl); |
| if (tree[tl].cr > 0) q.push_back(tree[tl].cr); |
| parents.push_back(tl); |
| } |
| } |
| } |
| |
| // Explicit template instantiation... |
| template class RCBForceTree<QUADRUPOLE_TDPTS>; |
| template class RCBForceTree<MONOPOLE_TDPTS>; |
| |
