diff --git a/README.md b/README.md index 09588d229..54d26c1fc 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,14 @@ uDeviceX comes with the GPLv2 LICENSE. === The uDeviceX folders is organized into: -balaprep: preprocessing tool for assigning ranks to compute nodes. -cell-placement: preprocessing tool for generating the initial RBC/CTC displacement -cuda-ctc: code for the CTC model -cuda-dpd: code for the DPD interactions -cuda-rbc: code for the RBC model -device-gen: preprocessing tool to generate new device geometries -halo-bench: OSU-like benchmark to measure latency and bandwidth across the MPI ranks. -mpi-dpd: the simulation code -proof-of-concept: tests and hacks. +* balaprep: preprocessing tool for assigning ranks to compute nodes. +* cell-placement: preprocessing tool for generating the initial RBC/CTC displacement +* cuda-ctc: code for the CTC model +* cuda-dpd: code for the DPD interactions +* cuda-rbc: code for the RBC model +* device-gen: preprocessing tool to generate new device geometries +* mpi-dpd: the simulation code +* postprocessing: auxiliary tools to extract quantities from simulations the output +* proof-of-concept: tests and hacks +* tests: accuracy and performance tests diff --git a/balaprep/cuda-dpdmassimo.cu b/balaprep/cuda-dpdmassimo.cu index 2e84c6e35..82182b127 100644 --- a/balaprep/cuda-dpdmassimo.cu +++ b/balaprep/cuda-dpdmassimo.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-02-26. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/balaprep/cudatimer.h b/balaprep/cudatimer.h index b7c27024e..d289d584e 100644 --- a/balaprep/cudatimer.h +++ b/balaprep/cudatimer.h @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-02-26. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ class CudaEventTimer { diff --git a/balaprep/mainmassimo.cu b/balaprep/mainmassimo.cu index f76cb1c93..4458a7703 100644 --- a/balaprep/mainmassimo.cu +++ b/balaprep/mainmassimo.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-02-25. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/balaprep/minmax.cu b/balaprep/minmax.cu index 441e22ed7..65e6b123b 100644 --- a/balaprep/minmax.cu +++ b/balaprep/minmax.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-13. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/balaprep/minmax.h b/balaprep/minmax.h index ff7592e68..59ac410a2 100644 --- a/balaprep/minmax.h +++ b/balaprep/minmax.h @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ typedef struct { diff --git a/balaprep/minmaxwrapper.cu b/balaprep/minmaxwrapper.cu index e63539a26..a3f1c915a 100644 --- a/balaprep/minmaxwrapper.cu +++ b/balaprep/minmaxwrapper.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/balaprep/prenompi.cc b/balaprep/prenompi.cc index b52eeb8c6..97fb7af31 100644 --- a/balaprep/prenompi.cc +++ b/balaprep/prenompi.cc @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-02-20. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/balaprep/scan_massimo.cu b/balaprep/scan_massimo.cu index 339cc8fba..74b574860 100644 --- a/balaprep/scan_massimo.cu +++ b/balaprep/scan_massimo.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "scanxdpd.h" diff --git a/balaprep/scanxdpd.h b/balaprep/scanxdpd.h index 0bf96b8aa..2c6d77139 100644 --- a/balaprep/scanxdpd.h +++ b/balaprep/scanxdpd.h @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ typedef struct { diff --git a/cell-placement/main.cpp b/cell-placement/main.cpp index bc3e9aeb8..a73a5a130 100644 --- a/cell-placement/main.cpp +++ b/cell-placement/main.cpp @@ -6,9 +6,6 @@ * Further edited by Dmitry Alexeev on 2014-03-25. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-ctc/misc.h b/cuda-ctc/misc.h index 18626d8d9..44897120a 100644 --- a/cuda-ctc/misc.h +++ b/cuda-ctc/misc.h @@ -8,4 +8,4 @@ #pragma once -typedef float real; \ No newline at end of file +typedef float real; diff --git a/cuda-dpd/cell-lists-faster.cu b/cuda-dpd/cell-lists-faster.cu index e78ca86c6..2b2c9a817 100644 --- a/cuda-dpd/cell-lists-faster.cu +++ b/cuda-dpd/cell-lists-faster.cu @@ -6,9 +6,6 @@ * Edited by Massimo Bernaschi on 2014-03-30. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/cell-lists.cu b/cuda-dpd/cell-lists.cu index 9a320eafc..e6e0b8410 100644 --- a/cuda-dpd/cell-lists.cu +++ b/cuda-dpd/cell-lists.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-21. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/cell-lists.h b/cuda-dpd/cell-lists.h index c9dc9b914..e4a236407 100644 --- a/cuda-dpd/cell-lists.h +++ b/cuda-dpd/cell-lists.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-21. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/dpd-rng.h b/cuda-dpd/dpd-rng.h index 9854d31dd..67e821a81 100644 --- a/cuda-dpd/dpd-rng.h +++ b/cuda-dpd/dpd-rng.h @@ -6,9 +6,6 @@ * Major editing (+Logistic RNG) from Yu-Hang Tang on 2015-03-19. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/dpd/Makefile b/cuda-dpd/dpd/Makefile index 81c92431f..3a8d0cc06 100644 --- a/cuda-dpd/dpd/Makefile +++ b/cuda-dpd/dpd/Makefile @@ -35,8 +35,8 @@ NVCCFLAGS += -DVISCOSITY_S_LEVEL=$(slevel) -lineinfo -Xptxas -v test-dpd: main.cpp libcuda-dpd.so $(CXX) $(CXXFLAGS) $^ -lcudart -lcurand -o test-dpd -libcuda-dpd.a: cuda-dpd.o cuda-dpd-bipartite.o ../profiler-dpd.o celllists - ar rcs $@ cuda-dpd.o cuda-dpd-bipartite.o ../profiler-dpd.o ../cell-lists.o ../cell-lists-faster.o +libcuda-dpd.a: cuda-dpd.o cuda-dpd-bipartite.o stress.o ../profiler-dpd.o celllists + ar rcs $@ cuda-dpd.o cuda-dpd-bipartite.o ../profiler-dpd.o ../cell-lists.o ../cell-lists-faster.o stress.o libcuda-dpd.so: cuda-dpd.o cuda-dpd-bipartite.o ../profiler-dpd.o celllists $(CXX) $(CXXFLAGS) -shared cuda-dpd.o cuda-dpd-bipartite.o ../profiler-dpd.o ../cell-lists.o ../cell-lists-faster.o -o libcuda-dpd.so -lcudart -lcurand @@ -48,6 +48,9 @@ cuda-dpd.o: $(CUDADPD) cuda-dpd.h cuda-dpd-bipartite.o: $(CUDADPDBIP) cuda-dpd.h $(NVCC) $(NVCCFLAGS) -c $(CUDADPDBIP) -o $@ +stress.o: stress.cu cuda-dpd.h + $(NVCC) $(NVCCFLAGS) -c $< -o $@ + ../%.o: make -C ../ $(@:../%=%) CXX="$(CXX)" NVCC="$(NVCC)" diff --git a/cuda-dpd/dpd/cuda-dpd-bipartite-floatized.cu b/cuda-dpd/dpd/cuda-dpd-bipartite-floatized.cu index 3f56424fe..fa57fd07e 100644 --- a/cuda-dpd/dpd/cuda-dpd-bipartite-floatized.cu +++ b/cuda-dpd/dpd/cuda-dpd-bipartite-floatized.cu @@ -5,9 +5,6 @@ * Created and authored by Yu-Hang Tang on 2015-03-18. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/cuda-dpd-bipartite.cu b/cuda-dpd/dpd/cuda-dpd-bipartite.cu index bcd3353e1..532fa0470 100644 --- a/cuda-dpd/dpd/cuda-dpd-bipartite.cu +++ b/cuda-dpd/dpd/cuda-dpd-bipartite.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-28. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/cuda-dpd-floatized-1.cu b/cuda-dpd/dpd/cuda-dpd-floatized-1.cu index 77824b65c..e18edfad7 100644 --- a/cuda-dpd/dpd/cuda-dpd-floatized-1.cu +++ b/cuda-dpd/dpd/cuda-dpd-floatized-1.cu @@ -6,9 +6,6 @@ * Created and authored by Yu-Hang Tang and Mauro Bisson on 2015-04-01. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/cuda-dpd-floatized-2.cu b/cuda-dpd/dpd/cuda-dpd-floatized-2.cu index 94f5af29f..ddc3751e7 100644 --- a/cuda-dpd/dpd/cuda-dpd-floatized-2.cu +++ b/cuda-dpd/dpd/cuda-dpd-floatized-2.cu @@ -6,9 +6,6 @@ * Created and authored by Yu-Hang Tang on 2015-03-18. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/cuda-dpd.cu b/cuda-dpd/dpd/cuda-dpd.cu index 34c77f4c5..4cc1a4b0e 100644 --- a/cuda-dpd/dpd/cuda-dpd.cu +++ b/cuda-dpd/dpd/cuda-dpd.cu @@ -6,9 +6,6 @@ * Major editing by Mauro Bisson on 2015-04-01. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/cuda-dpd.h b/cuda-dpd/dpd/cuda-dpd.h index dfeca681e..cee553e68 100644 --- a/cuda-dpd/dpd/cuda-dpd.h +++ b/cuda-dpd/dpd/cuda-dpd.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-04. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -29,7 +26,7 @@ template<> inline __device__ float viscosity_function<0>(float x){ return x; } void forces_dpd_cuda_nohost(const float * const xyzuvw, const float4 * const xyzouvwo, const ushort4 * const xyzo_half, float * const _axayaz, const int np, - const int * const cellsstart, const int * const cellscount, + const int * const cellsstart, const int * const cellscount, const float rc, const float XL, const float YL, const float ZL, const float aij, @@ -37,7 +34,7 @@ void forces_dpd_cuda_nohost(const float * const xyzuvw, const float4 * const xyz const float sigma, const float invsqrtdt, const float seed1, - cudaStream_t stream); + cudaStream_t stream); void forces_dpd_cuda(const float * const xp, const float * const yp, const float * const zp, const float * const xv, const float * const yv, const float * const zv, @@ -62,3 +59,13 @@ void forces_dpd_cuda_bipartite_nohost(cudaStream_t stream, const float2 * const const int3 halo_ncells, const float aij, const float gamma, const float sigmaf, const float seed, const int mask, float * const axayaz); + +void compute_stress(const float * const xyzuvw, + const int np, + const int * const cellsstart, const int * const cellscount, + const int XL, const int YL, const int ZL, + const float aij, const float gamma, const float sigmaf, const float seed, + float * const sigma_xx, float * const sigma_xy, float * const sigma_xz, + float * const sigma_yy, float * const sigma_yz, float * const sigma_zz, + float * const axayaz, + cudaStream_t stream); diff --git a/cuda-dpd/dpd/main.cpp b/cuda-dpd/dpd/main.cpp index 344d9f7be..934f9ddb2 100644 --- a/cuda-dpd/dpd/main.cpp +++ b/cuda-dpd/dpd/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-10. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/dpd/stress.cu b/cuda-dpd/dpd/stress.cu new file mode 100644 index 000000000..32c8aa51b --- /dev/null +++ b/cuda-dpd/dpd/stress.cu @@ -0,0 +1,261 @@ +/* + * stress.cu + * Part of uDeviceX/cuda-dpd-sem/dpd/ + * + * Created and authored by Diego Rossinelli on 2015-09-29. + * Copyright 2015. All rights reserved. + * + */ + +#include +#include + +#include "cuda-dpd.h" +#include "../dpd-rng.h" +#include "../hacks.h" + +namespace StressKernels +{ + struct InfoStress + { + int3 ncells; + float aij, gamma, sigmaf; + float *sigma_xx, *sigma_xy, *sigma_xz, *sigma_yy, *sigma_yz, *sigma_zz, *axayaz; + float seed; + }; + + __constant__ InfoStress info; + + texture texParticles; + texture texStart, texCount; + +#define _XCPB_ 2 +#define _YCPB_ 2 +#define _ZCPB_ 1 +#define CPB (_XCPB_ * _YCPB_ * _ZCPB_) + + __device__ float3 _dpd_interaction(const int dpid, const float3 xdest, const float3 udest, const int spid, const float2 stmp0, const float2 stmp1) + { + const int sentry = 3 * spid; + const float2 stmp2 = tex1Dfetch(texParticles, sentry + 2); + + const float _xr = xdest.x - stmp0.x; + const float _yr = xdest.y - stmp0.y; + const float _zr = xdest.z - stmp1.x; + const float rij2 = _xr * _xr + _yr * _yr + _zr * _zr; + assert(rij2 < 1); + + const float invrij = rsqrtf(rij2); + const float rij = rij2 * invrij; + const float argwr = 1 - rij; + const float wr = viscosity_function<-VISCOSITY_S_LEVEL>(argwr); + + const float xr = _xr * invrij; + const float yr = _yr * invrij; + const float zr = _zr * invrij; + + const float rdotv = + xr * (udest.x - stmp1.y) + + yr * (udest.y - stmp2.x) + + zr * (udest.z - stmp2.y); + + const float myrandnr = Logistic::mean0var1(info.seed, min(spid, dpid), max(spid, dpid)); + + const float strength = info.aij * argwr - (info.gamma * wr * rdotv + info.sigmaf * myrandnr) * wr; + + return make_float3(strength * xr, strength * yr, strength * zr); + } + +#define __IMOD(x,y) ((x)-((x)/(y))*(y)) + + template + __global__ void stress_kernel() + { + int mycount = 0, myscan = 0; + + __shared__ int volatile starts[CPB][16], scan[CPB][16]; + + if (threadIdx.x < 14) + { + const int cbase = blockIdx.x * blockDim.y + threadIdx.y; + + int dx, dy, dz; + dx = dy = dz = threadIdx.x / 3; + dx = threadIdx.x - dx * 3 - 1; + dy = __IMOD(dy, 3) - 1; + dz = __IMOD(dz / 3, 3) - 1; + + int cid = cbase + dz * info.ncells.x * info.ncells.y + dy * info.ncells.x + dx; + + const bool valid_cid = (cid >= 0) && (cid < info.ncells.x * info.ncells.y * info.ncells.z); + + starts[threadIdx.y][threadIdx.x] = (valid_cid) ? tex1Dfetch(texStart, cid) : 0; + + myscan = mycount = (valid_cid) ? tex1Dfetch(texCount, cid) : 0; + } + +#pragma unroll + for(int L = 1; L < 16; L <<= 1) + myscan += (threadIdx.x >= L) * __shfl_up(myscan, L); + + if (threadIdx.x < 15) + scan[threadIdx.y][threadIdx.x] = myscan - mycount; + + const int subtid = threadIdx.x % COLS; + const int slot = threadIdx.x / COLS; + + const int dststart = starts[threadIdx.y][13]; + const int lastdst = dststart + scan[threadIdx.y][14] - scan[threadIdx.y][13]; + + const int nsrc = scan[threadIdx.y][14]; + const int nsrcext = scan[threadIdx.y][13]; + + for(int pid = subtid; pid < nsrc; pid += COLS) + { + const int key9 = 9 * (pid >= scan[threadIdx.y][9]); + + int key3 = 3 * (pid >= scan[threadIdx.y][key9 + 3]); + key3 += (key9 < 9) ? 3 * (pid >= scan[threadIdx.y][key9 + 6]) : 0; + + int spid = pid - scan[threadIdx.y][key3 + key9] + starts[threadIdx.y][key3 + key9]; + + const int sentry = 3 * spid; + const float2 stmp0 = tex1Dfetch(texParticles, sentry); + const float2 stmp1 = tex1Dfetch(texParticles, sentry + 1); + + for(int dpid = dststart + slot; dpid < lastdst; dpid += ROWS) + { + float3 xdest, udest; + + float2 dtmp0 = tex1Dfetch(texParticles, 3 * dpid); + xdest.x = dtmp0.x; + xdest.y = dtmp0.y; + + dtmp0 = tex1Dfetch(texParticles, 3 * dpid + 1); + xdest.z = dtmp0.x; + udest.x = dtmp0.y; + + dtmp0 = tex1Dfetch(texParticles, 3 * dpid + 2); + udest.y = dtmp0.x; + udest.z = dtmp0.y; + + const float rx = xdest.x - stmp0.x; + const float ry = xdest.y - stmp0.y; + const float rz = xdest.z - stmp1.x; + + const float d2 = rx * rx + ry * ry + rz * rz; + + if ((dpid != spid) && (d2 < 1.0f)) + { + const float3 f = _dpd_interaction(dpid, xdest, udest, spid, stmp0, stmp1); + + atomicAdd(info.sigma_xx + dpid, f.x * rx); + atomicAdd(info.sigma_xy + dpid, f.x * ry); + atomicAdd(info.sigma_xz + dpid, f.x * rz); + atomicAdd(info.sigma_yy + dpid, f.y * ry); + atomicAdd(info.sigma_yz + dpid, f.y * rz); + atomicAdd(info.sigma_zz + dpid, f.z * rz); + + if (info.axayaz) + { + atomicAdd(info.axayaz + 3 * dpid , f.x); + atomicAdd(info.axayaz + 3 * dpid + 1, f.y); + atomicAdd(info.axayaz + 3 * dpid + 2, f.z); + } + + if (pid < nsrcext) + { + atomicAdd(info.sigma_xx + spid, f.x * rx); + atomicAdd(info.sigma_xy + spid, f.x * ry); + atomicAdd(info.sigma_xz + spid, f.x * rz); + atomicAdd(info.sigma_yy + spid, f.y * ry); + atomicAdd(info.sigma_yz + spid, f.y * rz); + atomicAdd(info.sigma_zz + spid, f.z * rz); + + if (info.axayaz) + { + atomicAdd(info.axayaz + 3*spid , -f.x); + atomicAdd(info.axayaz + 3*spid + 1, -f.y); + atomicAdd(info.axayaz + 3*spid + 2, -f.z); + } + } + } + } + } + } + + bool computestress_init = false; +} + +using namespace StressKernels; + +void compute_stress(const float * const xyzuvw, + const int np, + const int * const cellsstart, const int * const cellscount, + const int XL, const int YL, const int ZL, + const float aij, const float gamma, const float sigmaf, const float seed, + float * const sigma_xx, float * const sigma_xy, float * const sigma_xz, + float * const sigma_yy, float * const sigma_yz, float * const sigma_zz, + float * const axayaz, + cudaStream_t stream) +{ + if (np == 0) + { + printf("WARNING: stress_nohost called with np = %d\n", np); + return; + } + + if (!computestress_init) + { + texStart.channelDesc = cudaCreateChannelDesc(); + texStart.filterMode = cudaFilterModePoint; + texStart.mipmapFilterMode = cudaFilterModePoint; + texStart.normalized = 0; + + texCount.channelDesc = cudaCreateChannelDesc(); + texCount.filterMode = cudaFilterModePoint; + texCount.mipmapFilterMode = cudaFilterModePoint; + texCount.normalized = 0; + + texParticles.channelDesc = cudaCreateChannelDesc(); + texParticles.filterMode = cudaFilterModePoint; + texParticles.mipmapFilterMode = cudaFilterModePoint; + texParticles.normalized = 0; + + CUDA_CHECK(cudaFuncSetCacheConfig(stress_kernel<32, 1>, cudaFuncCachePreferL1)); + + computestress_init = true; + } + + size_t textureoffset; + CUDA_CHECK(cudaBindTexture(&textureoffset, &texParticles, xyzuvw, &texParticles.channelDesc, sizeof(float) * 6 * np)); + assert(textureoffset == 0); + + const int ncells = XL * YL * ZL; + + CUDA_CHECK(cudaBindTexture(&textureoffset, &texStart, cellsstart, &texStart.channelDesc, sizeof(int) * ncells)); + assert(textureoffset == 0); + CUDA_CHECK(cudaBindTexture(&textureoffset, &texCount, cellscount, &texCount.channelDesc, sizeof(int) * ncells)); + assert(textureoffset == 0); + + { + static InfoStress c = { make_int3(XL, YL, ZL), aij, gamma, sigmaf, + sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz, + axayaz, seed }; + + CUDA_CHECK(cudaMemcpyToSymbolAsync(info, &c, sizeof(c), 0, cudaMemcpyHostToDevice, stream)); + } + + if (axayaz) + CUDA_CHECK(cudaMemsetAsync(axayaz, 0, sizeof(float) * 3 * np, stream)); + + float * const ptrs[] = { sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz }; + + for(int c = 0; c < 6; ++c) + CUDA_CHECK(cudaMemsetAsync(ptrs[c], 0, sizeof(float) * np, stream)); + + stress_kernel<32, 1><<<(ncells + CPB - 1) / CPB, dim3(32, CPB), 0, stream>>>(); + + CUDA_CHECK(cudaPeekAtLastError()); +} + diff --git a/cuda-dpd/hacks.h b/cuda-dpd/hacks.h index 8d82a26b8..f3b4def27 100644 --- a/cuda-dpd/hacks.h +++ b/cuda-dpd/hacks.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-29. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/profiler-dpd.cpp b/cuda-dpd/profiler-dpd.cpp index a519845ac..7038baec8 100644 --- a/cuda-dpd/profiler-dpd.cpp +++ b/cuda-dpd/profiler-dpd.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-17. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/profiler-dpd.h b/cuda-dpd/profiler-dpd.h index fa1f6a7b3..8d532a81f 100644 --- a/cuda-dpd/profiler-dpd.h +++ b/cuda-dpd/profiler-dpd.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-23. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/sem/cell-factory.cu b/cuda-dpd/sem/cell-factory.cu index cabd7323c..e87899590 100644 --- a/cuda-dpd/sem/cell-factory.cu +++ b/cuda-dpd/sem/cell-factory.cu @@ -5,9 +5,6 @@ * Created and authored by Dmitry 2014-08-07 on 12:56:03. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/sem/cell-factory.h b/cuda-dpd/sem/cell-factory.h index 6aaae7082..19fb192bd 100644 --- a/cuda-dpd/sem/cell-factory.h +++ b/cuda-dpd/sem/cell-factory.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-08-08. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/sem/cuda-sem.cu b/cuda-dpd/sem/cuda-sem.cu index 3674f90f7..1771e3b39 100644 --- a/cuda-dpd/sem/cuda-sem.cu +++ b/cuda-dpd/sem/cuda-sem.cu @@ -5,9 +5,6 @@ * Created and authored by Dmitry 2014-08-05 on 17:25:23. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/sem/cuda-sem.h b/cuda-dpd/sem/cuda-sem.h index f80b30fd7..f2c544b67 100644 --- a/cuda-dpd/sem/cuda-sem.h +++ b/cuda-dpd/sem/cuda-sem.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-29. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/cuda-dpd/sem/main.cpp b/cuda-dpd/sem/main.cpp index d610429bd..a12da0b61 100644 --- a/cuda-dpd/sem/main.cpp +++ b/cuda-dpd/sem/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-10. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/sem/main.cu b/cuda-dpd/sem/main.cu index cba86b0cf..096f11e70 100644 --- a/cuda-dpd/sem/main.cu +++ b/cuda-dpd/sem/main.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-29. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/cuda-dpd/test-cell-lists.cu b/cuda-dpd/test-cell-lists.cu index 27d66f031..758e184f3 100644 --- a/cuda-dpd/test-cell-lists.cu +++ b/cuda-dpd/test-cell-lists.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-08-08. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -288,4 +285,4 @@ int main() } return 0; -} \ No newline at end of file +} diff --git a/cuda-dpd/tiny-float.h b/cuda-dpd/tiny-float.h index 7ff887306..9ffd7b916 100644 --- a/cuda-dpd/tiny-float.h +++ b/cuda-dpd/tiny-float.h @@ -5,9 +5,6 @@ * Created and authored by Yu-Hang Tang on 2015-03-14. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #ifndef _TINY_FLOAT_ diff --git a/cuda-rbc/misc.h b/cuda-rbc/misc.h index 18626d8d9..44897120a 100644 --- a/cuda-rbc/misc.h +++ b/cuda-rbc/misc.h @@ -8,4 +8,4 @@ #pragma once -typedef float real; \ No newline at end of file +typedef float real; diff --git a/device-gen/2Dto3D/Makefile b/device-gen/2Dto3D/Makefile deleted file mode 100644 index b6d431f68..000000000 --- a/device-gen/2Dto3D/Makefile +++ /dev/null @@ -1,7 +0,0 @@ -2Dto3D: main.cpp - g++ main.cpp -O0 -g3 -fopenmp -o 2Dto3D - -clean: - rm 2Dto3D - -.PHONY = clean diff --git a/device-gen/2Dto3D/main.cpp b/device-gen/2Dto3D/main.cpp deleted file mode 100644 index 6d8b48908..000000000 --- a/device-gen/2Dto3D/main.cpp +++ /dev/null @@ -1,90 +0,0 @@ -/* - * main.cpp - * Part of uDeviceX/device-gen/2Dto3D/ - * - * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. - * Copyright 2015. All rights reserved. - * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. - */ - -#include -#include -#include -#include -#include -#include - -using namespace std; - -int main(int argc, char ** argv) -{ - if (argc != 6) - { - printf("usage: ./2to3 \n"); - return 1; - } - - const float zextent = atof(argv[2]); - const float zmargin = atof(argv[3]); - const int NZ = atoi(argv[4]); - - int NX, NY; - float xextent, yextent; - - vector slice; - - { - printf("Reading file %s...\n", argv[1]); - FILE * f = fopen(argv[1], "r"); - assert(f != 0); - float zextentOld; - int NZOld; - fscanf(f, "%f %f %f\n", &xextent, &yextent, &zextentOld); - fscanf(f, "%d %d %d\n", &NX, &NY, &NZOld); - printf("Extent: [%f, %f, %f]. Grid size: [%d, %d, %d]\n", xextent, yextent, zextentOld, NX, NY,NZOld); - slice.resize(NX * NY, 0.0f); - fread(&slice[0], sizeof(float), slice.size(), f); - fclose(f); - } - - printf("Generating data with extent [%f, %f, %f], dimensions [%d, %d, %d], zmargin %f\n", - xextent, yextent, zextent + 2 * zmargin, NX, NY, NZ, zmargin); - vector volume(NX * NY * NZ, 0.0f); - - const float z0 = -zextent * 0.5 - zmargin; - const float dz = (zextent + 2 * zmargin) / (NZ - 1); - -//#pragma omp parallel for - for(int iz = 0; iz < NZ; ++iz) - { - const float z = z0 + iz * dz; - for(int iy = 0; iy < NY; ++iy) - for(int ix = 0; ix < NX; ++ix) - { - const float xysdf = slice[ix + NX * (NY - 1 - iy)]; // NY -1 to change Y-axis direction - const float zsdf = fabs(z) - zextent * 0.5; - float val; - if (xysdf < 0) - val = max(zsdf, xysdf); - else - val = (zsdf < 0) ? xysdf : sqrt(zsdf * zsdf + xysdf * xysdf); - - assert(iy + NY * (ix + NX * iz) < volume.size()); - assert(volume[iy + NY * (ix + NX * iz)] == 0.0f); - volume[iy + NY * (ix + NX * iz)] = val; - } - } - - { - FILE * f = fopen(argv[5], "w"); - assert(f != 0); - fprintf(f, "%f %f %f\n", yextent, xextent, zextent + 2.0f * zmargin); //exchange X and Y - fprintf(f, "%d %d %d\n", NY, NX, NZ); - fwrite(&volume[0], sizeof(float), volume.size(), f); - fclose(f); - } -} - diff --git a/device-gen/README.md b/device-gen/README.md new file mode 100644 index 000000000..105e34647 --- /dev/null +++ b/device-gen/README.md @@ -0,0 +1,30 @@ +# Generate microfluidic geometry + +Set of scripts to generate device geometries as Signed Distance Function in \*.dat format. +The dat format consists of header and the binary float data: +``` + + + +``` +Where size is the geometry length units (typically microns), grid size defines how many grid points are there. + +## Parabolic funnels +Geometry mimicing the microfluidic device by McFaul et al [Cell separation based on size and deformability using microfluidic funnel ratchets](http://www.ncbi.nlm.nih.gov/pubmed/22517056) +To generate a this geometry with 10 rows and 20 columns and with walls in z direction of width 4: +``` +cd funnels +make +./funnel -nColumns=20 -nRows=10 -zMargin=4 -out=geom.dat +``` + +## Later displacement device +Geometry reproducing CTC-iChip1 module by Karabacak et al [Microfluidic, marker-free isolation of circulating tumor cells from blood samples](http://www.nature.com/nprot/journal/v9/n3/full/nprot.2014.044.html) + +To build geometry constisting of 13 columns and 59 rows repeated twice, with wall widht 2 and such grid resolution that 0.5 grid points correspond to 1 unit of length: +``` +cd ctc-ichip +make +./ctc-ichip -nColumns=13 -nRows=59 -nRepeat=2 -zMargin=2.0 -out=13x59x2-05.dat -zResolution=0.5 +``` + diff --git a/device-gen/common/2Dto3D.cpp b/device-gen/common/2Dto3D.cpp new file mode 100644 index 000000000..b1daaf9b3 --- /dev/null +++ b/device-gen/common/2Dto3D.cpp @@ -0,0 +1,80 @@ +/* + * 2Dto3D.cpp + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ +#include "2Dto3D.h" +#include +#include +#include +#include +#include +#include +#include "common.h" + +using namespace std; + +void conver2Dto3D(const int NX, const int NY, const float xextent, const float yextent, const std::vector& slice, + const int NZ, const float zextent, const float zmargin, const std::string& fileName) +{ + printf("Generating data with extent [%f, %f, %f], dimensions [%d, %d, %d], zmargin %f\n", + xextent, yextent, zextent + 2 * zmargin, NX, NY, NZ, zmargin); + + vector outputslice(NX * NY, 0.0f); + + const float z0 = -zextent * 0.5 - zmargin; + const float dz = (zextent + 2 * zmargin) / (NZ - 1); + + FILE * f = fopen(fileName.c_str(), "w"); + assert(f != 0); + fprintf(f, "%f %f %f\n", yextent, xextent, zextent + 2.0f * zmargin); + fprintf(f, "%d %d %d\n", NY, NX, NZ); + + for(int iz = 0; iz < NZ; ++iz) + { + const float z = z0 + iz * dz; + +#pragma omp parallel for + for(int iy = 0; iy < NY; ++iy) + for(int ix = 0; ix < NX; ++ix) + { + const float xysdf = slice[ix + NX * (NY - 1 - iy)]; // NY -1 to change Y-axis direction + float val = xysdf; + if (zmargin != 0.0f) { + const float zsdf = fabs(z) - zextent * 0.5; + if (xysdf < 0) + val = max(zsdf, xysdf); + else + val = (zsdf < 0) ? xysdf : sqrt(zsdf * zsdf + xysdf * xysdf); + } + assert(iy + NY * (ix) < outputslice.size()); + + assert(fabs(val) < 1e3); // to check that the value has reasonable range + outputslice[iy + NY * ix] = val; + } + + if (iz == 0) + { + unsigned char * ptr = (unsigned char *)&outputslice[0]; + if ((ptr[0] >= 9 && ptr[0] <= 13) || ptr[0] == 32 ) + { + ptr[0] = (ptr[0] == 32) ? 33 : (ptr[0] < 11) ? 8 : 14; + printf("INFO: some symbols were changed while writing\n"); + } + } + + int result = fwrite(&outputslice.front(), sizeof(float), NX * NY, f); + + if (result != NX * NY) { + printf("ERROR: written less than expected"); + exit(3); + } + + } + + fclose(f); +} + diff --git a/device-gen/common/2Dto3D.h b/device-gen/common/2Dto3D.h new file mode 100644 index 000000000..e1de9a506 --- /dev/null +++ b/device-gen/common/2Dto3D.h @@ -0,0 +1,14 @@ +/* + * 2Dto3D.h + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ +#pragma once +#include +#include + +void conver2Dto3D(const int NX, const int NY, const float xextent, const float yextent, const std::vector& slice, + const int NZ, const float zextent, const float zmargin, const std::string& fileName); diff --git a/device-gen/common/collage.cpp b/device-gen/common/collage.cpp new file mode 100644 index 000000000..372efd973 --- /dev/null +++ b/device-gen/common/collage.cpp @@ -0,0 +1,127 @@ +/* + * collage.cpp + * Part of CTC/device-gen/sdf-collage/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ +#include "collage.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "common.h" +using namespace std; + +void collageSDF(int NX, int NY, const vector< vector >& sampleSDF, vector& outputSDF) +{ + outputSDF.resize(sampleSDF.size() * sampleSDF[0].size()); + printf("SIZE: %d\n", outputSDF.size()); + const int stride = NX; + for(int iy = 0; iy < sampleSDF.size() * NY; ++iy) + for(int ix = 0; ix < NX; ++ix) + { + const int dst = ix + stride * iy; + const int iobst = iy / NY; + assert(ix + NX * (iy - iobst*NY) < sampleSDF[iobst].size()); + outputSDF[dst] = sampleSDF[iobst][ix + NX * (iy - iobst*NY)]; + } +} + +void populateSDF(const int NX, const int NY, const float xextent, const float yextent, const std::vector& sampleSDF, + const int xtimes, const int ytimes, std::vector& outputSDF) +{ + printf("Populate %d * %d times\n", xtimes, ytimes); + const int stride = xtimes * NX; + outputSDF.resize(xtimes * ytimes * NX * NY); + + for(int ty = 0; ty < ytimes; ++ty) + for(int tx = 0; tx < xtimes; ++tx) { + for(int iy = 0; iy < NY; ++iy) + for(int ix = 0; ix < NX; ++ix) + { + const int gx = ix + NX * tx; + const int gy = iy + NY * ty; + const int dst = gx + stride * gy; + + assert(dst < outputSDF.size()); + assert(ix + NX * iy < sampleSDF.size()); + outputSDF[dst] = sampleSDF[ix + NX * iy]; + } + } +} + +void collageSDFWithWall(const int NX, const int NY, const float xextent, const float yextent, + const vector< vector >& sampleSDF, const int ytimes, + const float paddingAdd, vector& outputSDF) +{ + collageSDF(NX, NY, sampleSDF, outputSDF); + const int xtimes = 1; + int outputSDFNX = xtimes * NX; + int outputSDFNY = ytimes * NY; + const float x0 = -xtimes * xextent * 0.5; + const float dx = xtimes * xextent / (outputSDFNX - 1); + + const float y0 = -ytimes * yextent * 0.5; + const float dy = ytimes * yextent / (outputSDFNY - 1); + + const float angle = (1.8/180.)*M_PI; + const float normal[] = {-cos(angle), sin(angle)}; + const float wallWidth = -2*y0*tan(angle); + + float ypick = 25.0f; //15 + float widthOfBufferZone = paddingAdd - wallWidth; + float xpick = (wallWidth) * (-2.0f*y0 - ypick) / (-2.0f*y0); + const float angle2 = atan(xpick/ypick); + std::cout << "YY = " << xpick << ", " << y0 + ypick << " ANGLE = " << angle2/M_PI*180 << std::endl; + const float normal2[] = {-cos(angle2), -sin(angle2)}; + + const float linePoint[] = {-x0 - wallWidth, y0}; + const float linePoint2[] = {-x0 - xpick -wallWidth, -y0 - ypick}; + + for (int iy = 0; iy < outputSDFNY; ++iy) + for (int ix = 0; ix < outputSDFNX; ++ix) + { + const float signX = sign(dx*ix + x0); + float p[] = {dx*ix + x0, dy*iy + y0}; + float padding = signbit(-p[0])*widthOfBufferZone; + float xsdf = -1e6; + + if ((signX == -1 && p[1] > (y0 + ypick)) || (signX == 1 && p[1] < (-y0 - ypick))) { + xsdf = -(normal[0]*(fabs(p[0]) - linePoint[0] + padding) - signX*normal[1]*(p[1] - linePoint[1])); + } else { + xsdf = -(normal2[0]*(fabs(p[0]) - linePoint2[0] - signbit(p[0])*wallWidth + padding) - normal2[1]*(fabs(p[1]) - linePoint2[1])); + } + + outputSDF[ix + outputSDFNX*iy] = std::max(outputSDF[ix + outputSDFNX*iy], xsdf); + } +} + +void shiftSDF(const int NX, const int NY, const float xextent, const float yextent, const std::vector& inputGrid, + const float xshift, const float xpadding, int& newNX, float& newXextent, std::vector& outGrid) +{ + float h = xextent / (NX - 1); + int ixshift = xshift / h; + int ipadding = xpadding / h + 1; + newNX = NX + ipadding; + newXextent = xextent + xpadding; + + float minVal = *std::min(inputGrid.begin(), inputGrid.end()); + outGrid.resize(newNX * NY, -1e6); + + for(int iy = 0; iy < NY; ++iy) + for(int ix = 0; ix < NX ; ++ix) + { + int newIx = (ix + ixshift) % newNX; + assert(fabs(inputGrid[ix + NX * iy]) < 1e3); + outGrid[newIx + newNX * iy] = inputGrid[ix + NX * iy]; + } +} + diff --git a/device-gen/common/collage.h b/device-gen/common/collage.h new file mode 100644 index 000000000..78064c672 --- /dev/null +++ b/device-gen/common/collage.h @@ -0,0 +1,23 @@ +/* + * collage.h + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ +#pragma once + +#include + +void collageSDF(int NX, int NY, const std::vector< std::vector >& sampleSDF, std::vector& outputSDF); + +void collageSDFWithWall(const int NX, const int NY, const float xextent, const float yextent, + const std::vector< std::vector >& sampleSDF, const int ytimes, + const float paddingAdd, std::vector& outputSDF); + +void populateSDF(const int NX, const int NY, const float xextent, const float yextent, const std::vector& sampleSDF, + const int xtimes, const int ytimes, std::vector& outputSDF); + +void shiftSDF(const int NX, const int NY, const float xextent, const float yextent, const std::vector& inputGrid, + const float xshift, const float xpadding, int& newNX, float& newXextent, std::vector& outGrid); diff --git a/device-gen/common/common.h b/device-gen/common/common.h new file mode 100644 index 000000000..694a570c6 --- /dev/null +++ b/device-gen/common/common.h @@ -0,0 +1,64 @@ +/* + * common.h + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ + +#pragma once +#include +#include +#include +#include +#include + +inline float sign(float x) { + return 1.0f - 2.0f*std::signbit(x); +} + +inline void readDAT(const std::string& fileName, std::vector& data, + int& NX, int& NY, int& NZ, float& xextent, float& yextent, float& zextent) +{ + FILE * f = fopen(fileName.c_str(), "r"); + assert(f != 0); + int result = fscanf(f, "%f %f %f\n", &xextent, &yextent, &zextent); + assert(result == 3); + result = fscanf(f, "%d %d %d\n", &NX, &NY, &NZ); + assert(result == 3); + printf("Read file %s. Extent: [%f, %f, %f]. Grid size: [%d, %d, %d]\n", + fileName.c_str(), xextent, yextent, zextent, NX, NY, NZ); + data.resize(NX * NY * NZ, 0.0f); + result = fread(&data[0], sizeof(float), NX * NY * NZ, f); + if (result != data.size()) { + printf("ERROR: read less than expected"); + exit(3); + } + fclose(f); +} + +inline void writeDAT(const std::string& fileName, std::vector& data, + const int NX, const int NY, const int NZ, float xextent, float yextent, float zextent) +{ + FILE * f = fopen(fileName.c_str(), "w"); + assert(f != 0); + fprintf(f, "%f %f %f\n", xextent, yextent, zextent); + fprintf(f, "%d %d %d\n", NX, NY, NZ); + + unsigned char * ptr = (unsigned char *)&data[0]; + if ((ptr[0] >= 9 && ptr[0] <= 13) || ptr[0] == 32 ) + { + ptr[0] = (ptr[0] == 32) ? 33 : (ptr[0] < 11) ? 8 : 14; + printf("INFO: some symbols were changed while writing\n"); + } + + int result = fwrite(&data[0], sizeof(float), (int)data.size(), f); + if (result != data.size()) { + printf("ERROR: written less than expected"); + exit(3); + } + + fclose(f); +} + diff --git a/device-gen/common/device-builder.h b/device-gen/common/device-builder.h new file mode 100644 index 000000000..a0be05654 --- /dev/null +++ b/device-gen/common/device-builder.h @@ -0,0 +1,33 @@ + +/* + * device-builder.h + * Part of CTC/device-gen/ctc-ichip/ + * + * Created and authored by Kirill Lykov on 2015-09-7. + * Copyright 2015. All rights reserved. + * + */ + +#pragma once + +class DeviceBuilder +{ +protected: + typedef vector SDF; + int m_ncolumns, m_nrows; + float m_resolution, m_zmargin; + float m_unitSizeX, m_unitSizeY, m_unitSizeZ; // size of the egg with the empty space aroung it + std::string m_outFileName2D, m_outFileName3D; + int m_niterRedistance; + int m_unitNX, m_unitNY, m_unitNZ; +public: + DeviceBuilder(float unitSizeX, float unitSizeY, float unitSizeZ) + : m_ncolumns(0), m_nrows(0), m_resolution(0), m_zmargin(0), + m_unitSizeX(unitSizeX), m_unitSizeY(unitSizeY), m_unitSizeZ(unitSizeZ), + m_niterRedistance(1e3), m_unitNX(0), m_unitNY(0), m_unitNZ(0) + {} + + virtual void build() = 0; + + virtual ~DeviceBuilder() {} +}; diff --git a/device-gen/common/redistance.cpp b/device-gen/common/redistance.cpp new file mode 100644 index 000000000..1bd374a90 --- /dev/null +++ b/device-gen/common/redistance.cpp @@ -0,0 +1,131 @@ +/* + * resistance.h + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ +#include "redistance.h" +#include +#include + + float Redistance::sussman_scheme(int ix, int iy, float sgn0) + { + const float phicenter = _ACCESS(m_phi, ix, iy); + + const float dphidxm = phicenter - _ACCESS(m_phi, ix - 1, iy); + const float dphidxp = _ACCESS(m_phi, ix + 1, iy) - phicenter; + const float dphidym = phicenter - _ACCESS(m_phi, ix, iy - 1); + const float dphidyp = _ACCESS(m_phi, ix, iy + 1) - phicenter; + + if (sgn0 == 1) + { + const float xgrad0 = std::max( max(0.0f, dphidxm), -std::min(0.0f, dphidxp)) * m_invdx; + const float ygrad0 = std::max( max(0.0f, dphidym), -min(0.0f, dphidyp)) * m_invdy; + + const float G0 = sqrtf(xgrad0 * xgrad0 + ygrad0 * ygrad0) - 1.0f; + + return phicenter - m_dt * sgn0 * G0; + } + else + { + const float xgrad1 = std::max( -min(0.0f, dphidxm), std::max(0.0f, dphidxp)) * m_invdx; + const float ygrad1 = std::max( -min(0.0f, dphidym), std::max(0.0f, dphidyp)) * m_invdy; + + const float G1 = sqrtf(xgrad1 * xgrad1 + ygrad1 * ygrad1) - 1.0f; + + return phicenter - m_dt * sgn0 * G1; + } + } + +Redistance::Redistance(const float dt, const float dx, const float dy, + const int xsize, const int ysize) +: m_xsize(xsize), m_ysize(ysize), m_dt(dt), m_dx(dx), m_dy(dy), m_invdx(1.0f/dx), m_invdy(1.0f/dy), m_phi0(nullptr), m_phi(nullptr) +{} + + void Redistance::run(const int iterations, float * field) + { + for(int code = 0; code < 3 * 3; ++code) + { + if (code == 1 + 3) continue; + + const float deltax = m_dx * ((code % 3) - 1); + const float deltay = m_dy * ((code % 9) / 3 - 1); + + const float dl = sqrtf(deltax * deltax + deltay * deltay); + + m_dls[code] = dl; + } + + m_phi0 = new float[m_xsize * m_ysize]; + memcpy(m_phi0, field, sizeof(float) * m_xsize * m_ysize); + m_phi = field; + + float * tmp = new float[m_xsize * m_ysize]; + for(int t = 0; t < iterations; ++t) + { + if (t % 100 == 0) + printf("t: %d, size: %d %d\n", t, m_xsize, m_ysize); + +#pragma omp parallel for + for(int iy = 0; iy < m_ysize; ++iy) + for(int ix = 0; ix < m_xsize; ++ix) + { + const float myval0 = _ACCESS(m_phi0, ix, iy); + const float sgn0 = myval0 > 0 ? 1 : (myval0 < 0 ? -1 : 0); + + const bool boundary = ( + ix == 0 || ix == m_xsize - 1 || + iy == 0 || iy == m_ysize - 1); + if (boundary) + tmp[ix + m_xsize * iy] = simple_scheme(ix, iy, sgn0, myval0); + else + { + if (anycrossing(ix, iy, sgn0)) + tmp[ix + m_xsize * iy] = myval0; + else + tmp[ix + m_xsize * iy] = sussman_scheme(ix, iy, sgn0); + } + assert(fabs(tmp[ix + m_xsize * iy]) < 1e7); + } + + memcpy(field, tmp, sizeof(float) * m_xsize * m_ysize); + } + + delete [] tmp; + delete [] m_phi0; + m_phi0 = nullptr; + } + + float Redistance::simple_scheme(int ix, int iy, float sgn0, float myphi0) + { + float mindistance = 1e6f; + for(int code = 0; code < 3 * 3; ++code) + { + if (code == 1 + 3) continue; + + const int xneighbor = ix + (code % 3) - 1; + const int yneighbor = iy + (code % 9) / 3 - 1; + + if (xneighbor < 0 || xneighbor >= m_xsize) continue; + if (yneighbor < 0 || yneighbor >= m_ysize) continue; + + const float phi0_neighbor = _ACCESS(m_phi0, xneighbor, yneighbor); + const float phi_neighbor = _ACCESS(m_phi, xneighbor, yneighbor); + + const float dl = m_dls[code]; + + float distance = 0; + + if (sgn0 * phi0_neighbor < 0) + distance = - myphi0 * dl / (phi0_neighbor - myphi0); + else + distance = dl + abs(phi_neighbor); + + mindistance = std::min(mindistance, distance); + } + + return sgn0 * mindistance; + } + diff --git a/device-gen/common/redistance.h b/device-gen/common/redistance.h new file mode 100644 index 000000000..70f63698b --- /dev/null +++ b/device-gen/common/redistance.h @@ -0,0 +1,49 @@ +/* + * resistance.h + * Part of CTC/device-gen/common/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ + +#include +#include + +using namespace std; + +#define _ACCESS(f, x, y) f[(x) + m_xsize * (y)] + +class Redistance +{ + int m_xsize, m_ysize; + float * m_phi0, * m_phi; + float m_dt, m_dx, m_dy, m_invdx, m_invdy; + float m_dls[9]; + + template + inline bool anycrossing_dir(int ix, int iy, const float sgn0) + { + const int dx = d == 0, dy = d == 1, dz = d == 2; + + const float fm1 = _ACCESS(m_phi0, ix - dx, iy - dy); + const float fp1 = _ACCESS(m_phi0, ix + dx, iy + dy); + + return (fm1 * sgn0 < 0 || fp1 * sgn0 < 0); + } + + inline bool anycrossing(int ix, int iy, const float sgn0) + { + return + anycrossing_dir<0>(ix, iy, sgn0) || + anycrossing_dir<1>(ix, iy, sgn0); + } + + float simple_scheme(int ix, int iy, float sgn0, float myphi0); + + float sussman_scheme(int ix, int iy, float sgn0); +public: + Redistance(const float dt, const float dx, const float dy, + const int xsize, const int ysize); + void run(const int iterations, float * field); +}; diff --git a/device-gen/ctc-ichip/Makefile b/device-gen/ctc-ichip/Makefile new file mode 100644 index 000000000..93dccdc3e --- /dev/null +++ b/device-gen/ctc-ichip/Makefile @@ -0,0 +1,19 @@ +CXX = g++-4.9 +CXXFLAGS += -O0 -g3 -std=c++11 -fopenmp + +ctc-ichip: *.cpp collage.o redistance.o 2Dto3D.o + $(CXX) $(CXXFLAGS) -I../../mpi-dpd/ collage.o redistance.o 2Dto3D.o main.cpp -o ctc-ichip + +collage.o: ../common/collage.h ../common/collage.cpp + $(CXX) $(CXXFLAGS) -c $^ + +redistance.o: ../common/redistance.h ../common/redistance.cpp + $(CXX) $(CXXFLAGS) -c $^ + +2Dto3D.o: ../common/2Dto3D.h ../common/2Dto3D.cpp + $(CXX) $(CXXFLAGS) -c $^ + +clean: + rm -f ctc-ichip *.o *.d *.h.gch + + diff --git a/device-gen/ctc-ichip/main.cpp b/device-gen/ctc-ichip/main.cpp new file mode 100644 index 000000000..1f1c6c69b --- /dev/null +++ b/device-gen/ctc-ichip/main.cpp @@ -0,0 +1,294 @@ +/* + * main.cpp + * Part of CTC/device-gen/ctc-ichip/ + * + * Created and authored by Kirill Lykov on 2015-09-7. + * Copyright 2015. All rights reserved. + * + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include "../common/device-builder.h" +#include "../common/common.h" +#include "../common/collage.h" +#include "../common/redistance.h" +#include "../common/2Dto3D.h" + +using namespace std; + +struct Egg +{ + float r1, r2, alpha; + + Egg() + : r1(12.0f), r2(8.5f), alpha(0.03f) + { + } + + float x2y(float x) const { + return sqrt(r2*r2 * exp(-alpha * x) * (1.0f - x*x/r1/r1)); + } + + void run(vector& vx, vector& vy) { + int N = 500; + float dx = 2.0f * r1 / (N - 1); + for (int i = 0; i < N; ++i) { + float x = i * dx - r1; + float y = x2y(x); + vx.push_back(x); + vy.push_back(y); + } + + auto vxRev = vx; + vx.insert(vx.end(), vxRev.rbegin(), vxRev.rend()); + + auto vyRev = vy; + for_each(vyRev.begin(), vyRev.end(), [](float& i) { i *= -1.0f; }); + vy.insert(vy.end(), vyRev.rbegin(), vyRev.rend()); + } +}; + +class CTCiChip1Builder : public DeviceBuilder +{ + int m_nrepeat; + const float m_angle; + float m_desiredSubdomainSzX; +public: + CTCiChip1Builder() + : DeviceBuilder(56.0f, 32.0f, 128.0f), + m_nrepeat(0), m_angle(1.7f * M_PI / 180.0f) + {} + + CTCiChip1Builder& setNColumns(int ncolumns) + { + m_ncolumns = ncolumns; + return *this; + } + + CTCiChip1Builder& setNRows(int nrows) + { + m_nrows = nrows; + return *this; + } + + CTCiChip1Builder& setRepeat(float nrepeat) + { + m_nrepeat = nrepeat; + return *this; + } + + CTCiChip1Builder& setResolution(float resolution) + { + m_resolution = resolution; + return *this; + } + + CTCiChip1Builder& setZWallWidth(float zmargin) + { + m_zmargin = zmargin; + return *this; + } + + CTCiChip1Builder& setDiseredSubdomainX(float x) + { + m_desiredSubdomainSzX = x; + return *this; + } + + CTCiChip1Builder& setFileNameFor2D(const std::string& outFileName2D) + { + m_outFileName2D = outFileName2D; + return *this; + } + + CTCiChip1Builder& setFileNameFor3D(const std::string& outFileName3D) + { + m_outFileName3D = outFileName3D; + return *this; + } + + void build(); + +private: + void generateUnitSDF(vector& sdf) const; + + void shiftRows(int rowNX, int rowNY, float rowSizeX, float rowSizeY, const SDF& rowObstacles, + float& padding, float& addPadding, int& shiftedRowNX, float& shiftedRowSizeX,std::vector& shiftedRows) const; +}; + +void CTCiChip1Builder::build() +{ + if (m_ncolumns * m_nrows * m_nrepeat * m_resolution * m_zmargin == 0.0f || m_outFileName3D.length() == 0) + throw std::runtime_error("Invalid parameters"); + + // 1 Create 1 obstacle + m_unitNX = static_cast(m_unitSizeX * m_resolution); + m_unitNY = static_cast(m_unitSizeY * m_resolution); + m_unitNZ = static_cast(m_unitSizeZ * m_resolution); + + SDF eggSdf; + generateUnitSDF(eggSdf); + + // 2 Create 1 row of obstacles + int rowNX = m_ncolumns*m_unitNX; + int rowNY = m_unitNY; + int rowSizeX = m_ncolumns * m_unitSizeX; + int rowSizeY = m_unitSizeY; + SDF rowObstacles; + populateSDF(m_unitNX, m_unitNY, m_unitSizeX, m_unitSizeY, eggSdf, m_ncolumns, 1, rowObstacles); + + // 3 Shift rows + float padding = 0.0f; + float addPadding = 0.0f; + int shiftedRowNX = 0; // they are all the same length + float shiftedRowSizeX = 0.0f; + + std::vector shiftedRows; + shiftRows(rowNX, rowNY, rowSizeX, rowSizeY, rowObstacles, padding, addPadding, shiftedRowNX, shiftedRowSizeX, shiftedRows); + + // 4 Collage rows + SDF finalSDF; + collageSDFWithWall(shiftedRowNX, rowNY, shiftedRowSizeX, rowSizeY, shiftedRows, m_nrows, addPadding, finalSDF); + + // 5 Apply redistancing for the result + float finalExtent[] = {shiftedRowSizeX, static_cast(m_nrows * rowSizeY)}; + int finalN[] = {shiftedRowNX, m_nrows*rowNY}; + const float dx = finalExtent[0] / (finalN[0] - 1); + const float dy = finalExtent[1] / (finalN[1] - 1); + Redistance redistancer(0.25f * min(dx, dy), dx, dy, finalN[0], finalN[1]); + redistancer.run(m_niterRedistance, &finalSDF[0]); + + // 6 Repeat this pattern + SDF finalSDF2; + populateSDF(finalN[0], finalN[1], finalExtent[0], finalExtent[1], finalSDF, 1, m_nrepeat, finalSDF2); + std::swap(finalSDF, finalSDF2); + + // 6 Write result to the file + if (m_outFileName2D.length() != 0) + writeDAT(m_outFileName2D, finalSDF, finalN[0], m_nrepeat * finalN[1], 1, finalExtent[0], m_nrepeat * finalExtent[1], 1.0f); + + conver2Dto3D(finalN[0], m_nrepeat * finalN[1], finalExtent[0], m_nrepeat*finalExtent[1], finalSDF, + m_unitNZ, m_unitSizeZ - 2.0f*m_zmargin, m_zmargin, m_outFileName3D); +} + +void CTCiChip1Builder::generateUnitSDF(vector& sdf) const +{ + vector xs, ys; + Egg egg; + egg.run(xs, ys); + + const float xlb = -m_unitSizeX/2.0f; + const float ylb = -m_unitSizeY/2.0f; + + sdf.resize(m_unitNX * m_unitNY, 0.0f); + const float dx = m_unitSizeX / (m_unitNX - 1); + const float dy = m_unitSizeY / (m_unitNY - 1); + const int nsamples = xs.size(); + + for(int iy = 0; iy < m_unitNY; ++iy) + for(int ix = 0; ix < m_unitNX; ++ix) + { + const float x = xlb + ix * dx; + const float y = ylb + iy * dy; + + float distance2 = 1e6; + int iclosest = 0; + for(int i = 0; i < nsamples ; ++i) + { + const float xd = xs[i] - x; + const float yd = ys[i] - y; + const float candidate = xd * xd + yd * yd; + + if (candidate < distance2) + { + iclosest = i; + distance2 = candidate; + } + } + + float s = -1; + + { + const float ycurve = egg.x2y(x); + if (x >= -egg.r1 && x <= egg.r1 && fabs(y) <= ycurve) + s = +1; + } + + + sdf[ix + m_unitNX * iy] = s * sqrt(distance2); + } +} + +void CTCiChip1Builder::shiftRows(int rowNX, int rowNY, float rowSizeX, float rowSizeY, const SDF& rowObstacles, + float& padding, float& addPadding, int& shiftedRowNX, float& shiftedRowSizeX,std::vector& shiftedRows) const +{ + const int nRowsPerShift = static_cast(ceil(m_unitSizeX / (m_unitSizeY * tan(m_angle)))); + if (fabs(m_unitSizeX / (m_unitSizeY * tan(m_angle)) - nRowsPerShift) > 1e-1) { + throw std::runtime_error("Suggest changing the angle"); + } + + padding = float(ceil(m_nrows * m_unitSizeY * tan(m_angle))); + // TODO Do I need this nUniqueRows? + int nUniqueRows = m_nrows; + if (m_nrows > nRowsPerShift) { + nUniqueRows = nRowsPerShift; + padding = float(round(nRowsPerShift * m_unitSizeY * tan(m_angle))); + } + + // TODO fix this stupid workaround + if (padding < 32.0f) + padding = 0.0f; + if (padding == 57.0f) + padding = m_unitSizeX; + + // additional hack to have domain size in X direction to be devisible by desiredSubdomainSzX + { + float origSzX = m_ncolumns*m_unitSizeX + padding; + addPadding = (int(origSzX/m_desiredSubdomainSzX) + 1)*m_desiredSubdomainSzX - origSzX; + padding = padding + addPadding; // adjust padding to have desired size + } + + std::cout << "Launching rows generation. New size = "<< m_ncolumns*m_unitSizeX + padding << std::endl; + shiftedRows.resize(nUniqueRows); + for (int i = 0; i < nUniqueRows; ++i) { + float xshift = (nUniqueRows - i - 1) * 32.0f * tan(m_angle); + shiftSDF(rowNX, rowNY, rowSizeX, rowSizeY, rowObstacles, xshift, padding, shiftedRowNX, shiftedRowSizeX, shiftedRows[i]); + } +} + + +int main(int argc, char ** argv) +{ + ArgumentParser argp(vector(argv, argv + argc)); + + int nColumns = argp("-nColumns").asInt(1); + int nRows = argp("-nRows").asInt(1); + int nRepeat = argp("-nRepeat").asInt(1); + float zMargin = static_cast(argp("-zMargin").asDouble(5.0)); + float resolution = static_cast(argp("-zResolution").asDouble(1.0)); + std::string outFileName = argp("-out").asString("3d"); + + CTCiChip1Builder builder; + try { + builder.setNColumns(nColumns) + .setNRows(nRows) + .setRepeat(nRepeat) + .setResolution(resolution) + .setZWallWidth(zMargin) + .setFileNameFor2D("2d") + .setFileNameFor3D(outFileName) + .setDiseredSubdomainX(64.0f) + .build(); + } catch(const std::exception& ex) { + std::cout << "ERROR: " << ex.what() << std::endl; + } + return 0; +} + diff --git a/device-gen/funnels/Makefile b/device-gen/funnels/Makefile new file mode 100644 index 000000000..5c7799168 --- /dev/null +++ b/device-gen/funnels/Makefile @@ -0,0 +1,19 @@ +CXX = g++-4.9 +CXXFLAGS += -O0 -g3 -std=c++11 -fopenmp + +funnel: *.cpp collage.o redistance.o 2Dto3D.o + $(CXX) $(CXXFLAGS) -I../../mpi-dpd/ collage.o redistance.o 2Dto3D.o main.cpp -o funnel + +collage.o: ../common/collage.h ../common/collage.cpp + $(CXX) $(CXXFLAGS) -c $^ + +redistance.o: ../common/redistance.h ../common/redistance.cpp + $(CXX) $(CXXFLAGS) -c $^ + +2Dto3D.o: ../common/2Dto3D.h ../common/2Dto3D.cpp + $(CXX) $(CXXFLAGS) -c $^ + +clean: + rm -f test *.o *.d *.h.gch + + diff --git a/device-gen/funnels/main.cpp b/device-gen/funnels/main.cpp new file mode 100644 index 000000000..f8ea25a2b --- /dev/null +++ b/device-gen/funnels/main.cpp @@ -0,0 +1,231 @@ +/* + * main.cpp + * Part of CTC/device-gen/sdf-unit-par/ + * + * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. + * Copyright 2015. All rights reserved. + * + */ + +#include +#include +#include +#include +#include +#include +#include "../common/device-builder.h" +#include "../common/common.h" +#include "../common/collage.h" +#include "../common/redistance.h" +#include "../common/2Dto3D.h" + +using namespace std; + +struct Parabola +{ + float x0; + float ymax; // length of obstacle, for cutting the pick + float y0; + + Parabola(float gap, float xextent) : ymax(48.0) + { + x0 = xextent/2.0f - gap/2.0f; + if (gap > 9.0) + y0 = ymax; + else + y0 = ymax + 6.5; + } + + void line1(vector& vx, vector& vy) { + int N = 500; + float dx = 2.0f * fabs(x0) / (N - 1); + for (int i = 0; i < N; ++i) { + float x = i * dx - x0; + float y = 0.0f; + vx.push_back(x); + vy.push_back(y); + } + } + + void line2(vector& vx, vector& vy) { + int N = 1500; + float dx = 2.0f * fabs(x0) / (N - 1); + float alpha = -y0 / (x0 * x0); + for (int i = 0; i < N; ++i) { + float x = i * dx - x0; + float y = min(ymax, alpha * x*x + y0); + vx.push_back(x); + vy.push_back(y); + } + } +}; + +class FunnelsBuilder : public DeviceBuilder +{ + float m_gapSpace; // unit gap between obstacles +public: + FunnelsBuilder() + : DeviceBuilder(24.0f, 96.0f, 58.0f), m_gapSpace(1.0f) + {} + + FunnelsBuilder& setNColumns(int ncolumns) + { + m_ncolumns = ncolumns; + return *this; + } + + FunnelsBuilder& setNRows(int nrows) + { + m_nrows = nrows; + return *this; + } + + FunnelsBuilder& setResolution(float resolution) + { + m_resolution = resolution; + return *this; + } + + FunnelsBuilder& setZWallWidth(float zmargin) + { + m_zmargin = zmargin; + return *this; + } + + FunnelsBuilder& setFileNameFor2D(const std::string& outFileName2D) + { + m_outFileName2D = outFileName2D; + return *this; + } + + FunnelsBuilder& setFileNameFor3D(const std::string& outFileName3D) + { + m_outFileName3D = outFileName3D; + return *this; + } + + void build(); + +private: + void generateUnitSDF(float gap, vector& sdf) const; +}; + +void FunnelsBuilder::generateUnitSDF(float gap, vector& sdf) const +{ + assert(m_unitNX * m_unitNY * m_unitNZ != 0); + vector xs, ys; + Parabola par(gap, m_unitSizeX); + par.line1(xs, ys); + par.line2(xs, ys); + + const float xlb = -m_unitSizeX/2.0f; + const float ylb = -(m_unitSizeY - par.ymax)/2.0f; + + sdf.resize(m_unitNX * m_unitNY, 0.0f); + const float dx = m_unitSizeX / (m_unitNX - 1); //TODO NX-1 + const float dy = m_unitSizeY / (m_unitNY - 1); + const int nsamples = xs.size(); + + for(int iy = 0; iy < m_unitNY; ++iy) + for(int ix = 0; ix < m_unitNX; ++ix) + { + const float x = xlb + ix * dx; + const float y = ylb + iy * dy; + + float distance2 = 1e6; + int iclosest = 0; + for(int i = 0; i < nsamples ; ++i) + { + const float xd = xs[i] - x; + const float yd = ys[i] - y; + const float candidate = xd * xd + yd * yd; + + if (candidate < distance2) + { + iclosest = i; + distance2 = candidate; + } + } + + float s = -1; + + { + const float alpha = -par.y0 / (par.x0 * par.x0); + const float ycurve = min(par.ymax, alpha * x*x + par.y0); + + if (x >= -par.x0 && x <= par.x0 && y >= 0 && y <= ycurve) + s = +1; + } + + + sdf[ix + m_unitNX * iy] = s * sqrt(distance2); + } +} + +void FunnelsBuilder::build() +{ + if (m_ncolumns * m_nrows * m_resolution * m_zmargin == 0.0f || m_outFileName3D.length() == 0) + throw std::runtime_error("Invalid parameters"); + + // 1 Create obstacles with different gaps + m_unitNX = static_cast(m_unitSizeX * m_resolution); + m_unitNY = static_cast(m_unitSizeY * m_resolution); + m_unitNZ = static_cast(m_unitSizeZ * m_resolution); + + std::vector unitSDF(m_nrows); + for (int i = 0; i < m_nrows; ++i) { + float gap = m_gapSpace * (i + 3); + generateUnitSDF(gap, unitSDF[i]); + } + + // 2 Create rows of obstacles + std::vector rows(m_nrows); + for (int i = 0; i < m_nrows; ++i) { + populateSDF(m_unitNX, m_unitNY, m_unitSizeX, m_unitSizeY, unitSDF[i], m_ncolumns, 1, rows[i]); + } + + // 3 Collage rows + SDF finalSDF; + collageSDF(m_unitNX * m_ncolumns, m_unitNY, rows, finalSDF); + //collageSDF(m_unitNX * m_ncolumns, m_unitNY, m_unitSizeX * m_ncolumns, m_unitSizeY, rows, m_nrows, false, finalSDF); + + // 4 Apply redistancing for the result + float finalExtent[] = {m_unitSizeX * m_ncolumns, m_unitSizeY * m_nrows}; + int finalN[] = {m_unitNX * m_ncolumns, m_unitNY * m_nrows}; + const float dx = finalExtent[0] / (finalN[0] - 1); + const float dy = finalExtent[1] / (finalN[1] - 1); + Redistance redistancer(0.25f * min(dx, dy), dx, dy, finalN[0], finalN[1]); + redistancer.run(m_niterRedistance, &finalSDF[0]); + + if (m_outFileName2D.length() != 0) + writeDAT(m_outFileName2D.c_str(), finalSDF, finalExtent[0], finalExtent[1], 1.0f, finalN[0], finalN[1], 1); + + conver2Dto3D(finalN[0], finalN[1], finalExtent[0], finalExtent[1], finalSDF, m_unitNZ, m_unitSizeZ - 2.0f*m_zmargin, m_zmargin, m_outFileName3D); +} + +int main(int argc, char ** argv) +{ + ArgumentParser argp(vector(argv, argv + argc)); + + int nColumns = argp("-nColumns").asInt(1); + int nRows = argp("-nRows").asInt(1); + float zMargin = static_cast(argp("-zMargin").asDouble(5.0)); + float resolution = static_cast(argp("-zResolution").asDouble(1.0)); + + std::string outFileName = argp("-out").asString("3d"); + + FunnelsBuilder builder; + try { + builder.setNColumns(nColumns) + .setNRows(nRows) + .setResolution(1.0f) + .setZWallWidth(zMargin) + .setFileNameFor3D(outFileName) + .setResolution(resolution) + .build(); + } catch(const std::exception& ex) { + std::cout << "ERROR: " << ex.what() << std::endl; + return 1; + } + return 0; +} diff --git a/device-gen/pipe/Makefile b/device-gen/pipe/Makefile new file mode 100644 index 000000000..2d07290c1 --- /dev/null +++ b/device-gen/pipe/Makefile @@ -0,0 +1,7 @@ +sdf-unit: main.cpp + g++-4.9 main.cpp -O0 -g3 -o sdf-unit + +clean: + rm sdf-unit + +.PHONY = clean diff --git a/device-gen/pipe/main.cpp b/device-gen/pipe/main.cpp new file mode 100644 index 000000000..898c6031c --- /dev/null +++ b/device-gen/pipe/main.cpp @@ -0,0 +1,58 @@ +/* + * main.cpp + * Part of CTC/device-gen/sdf-unit-par/ + * + * Created and authored by Kirill Lykov on 2015-08-28. + * Copyright 2015. All rights reserved. + * + */ + +#include +#include +#include +#include +#include +#include "../common/common.h" +using namespace std; + +#define REAL float + +REAL distToSide(REAL x, REAL y, REAL radius) { + return sqrt(x*x + y*y) - radius; +} + +int main(int argc, char ** argv) +{ + if (argc != 4) + { + printf("usage: ./sdf-cylinder \n"); + return 1; + } + + const int N = atoi(argv[1]); + const REAL radius = atof(argv[2]);; + const REAL extent = 2.0f*radius + 4.0f; + + std:cout << "Will generate SDF with extent " << extent << " " << extent + << ". Grid size " << N << " x " << N << std::endl; + + const REAL xlb = -extent/2.0f; + const REAL ylb = -extent/2.0f; + + vector sdf(N * N, 0.0f); + const REAL dx = extent / (N-1); + const REAL dy = extent / (N-1); + + for(int iy = 0; iy < N; ++iy) + for(int ix = 0; ix < N; ++ix) + { + const REAL x = xlb + ix * dx; + const REAL y = ylb + iy * dy; + + sdf[ix + N * iy] = distToSide(x, y, radius); + } + + writeDAT(argv[3], sdf, extent, extent, REAL(1.0), N, N, 1); + + return 0; +} diff --git a/device-gen/plates/Makefile b/device-gen/plates/Makefile new file mode 100644 index 000000000..4e67bbfb1 --- /dev/null +++ b/device-gen/plates/Makefile @@ -0,0 +1,9 @@ +CXXFLAGS += -O3 -DNDEBUG + +plates: plates.cpp + $(CXX) $(CXXFLAGS) $< -o plates + +clean: + rm -f plates + +.PHONY = clean diff --git a/device-gen/plates/plates.cpp b/device-gen/plates/plates.cpp new file mode 100644 index 000000000..6e3a4546f --- /dev/null +++ b/device-gen/plates/plates.cpp @@ -0,0 +1,49 @@ +#include +#include +#include +#include + +int main(const int argc, const char * argv[]) +{ + if (argc != 5) + { + printf("usage: ./plates \n"); + return 1; + } + + const int NX = atoi(argv[1]); + const int NY = atoi(argv[2]); + const int NZ = atoi(argv[3]); + const int zmargin = atoi(argv[4]); + + float * data = new float[NX * NY * NZ]; + + for(int iz = 0; iz < NZ; ++iz) + { + const float zval = fabs(iz + 0.5 - NZ * 0.5) - (NZ / 2 - zmargin); + + for(int iy = 0; iy < NY; ++iy) + for(int ix = 0; ix < NX; ++ix) + data[ix + NX * (iy + NY * iz)] = zval; + } + + FILE * f = fopen("sdf.dat", "w"); + assert(f != 0); + fprintf(f, "%f %f %f\n", (float)NX, (float)NY, (float)NZ); + fprintf(f, "%d %d %d\n", NY, NX, NZ); + fwrite(data, sizeof(float), NX * NY * NZ, f); + fclose(f); + +#ifndef NDEBUG + { + FILE * f = fopen("sdf.raw", "w"); + assert(f != 0); + fwrite(data, sizeof(float), NX * NY * NZ, f); + fclose(f); + } +#endif + + delete [] data; + + return 0; +} diff --git a/device-gen/post-process/h52ply.py b/device-gen/post-process/h52ply.py new file mode 100755 index 000000000..c2a2dd31e --- /dev/null +++ b/device-gen/post-process/h52ply.py @@ -0,0 +1,37 @@ +#!/usr/bin/env /Applications/paraview.app/Contents/bin/pvpython + +''' + * Part of CTC/device-gen/post-processing/h52ply.py + * + * Created and authored by Kirill Lykov on 2015-08-28. + * Copyright 2015. All rights reserved. + * + * Users are NOT authorized + * to employ the present software for their own publications + * before getting a written permission from the author of this file. +''' + +import argparse +import os +from paraview.simple import * + +print("h52ply started") +parser = argparse.ArgumentParser(description='Transforms h5 which has xmf to ply using Paraview python lib.', + usage= './h52ply.py -i -o ') +parser.add_argument('-i','--inputFile', help='XMF', required=True) +parser.add_argument('-o','--outputFile', help='PLY', required=True) +args = vars(parser.parse_args()) + +# paraview wants to have absolute path +fullPath = os.path.dirname(os.path.abspath(args['inputFile'])) + '/' +print fullPath + +a13x59xmf = XDMFReader(FileNames=[fullPath + args['inputFile']]) +#a13x59xmf.GridStatus = ['Grid_26'] + +contour1 = Contour(Input=a13x59xmf) +contour1.Isosurfaces = [0.0] + +# save data +SaveData(fullPath + args['outputFile'], proxy=contour1) + diff --git a/device-gen/post-process/plyScale.py b/device-gen/post-process/plyScale.py new file mode 100755 index 000000000..82a06ea04 --- /dev/null +++ b/device-gen/post-process/plyScale.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python + +''' + * Part of CTC/device-gen/post-processing/plyScale.py + * + * Created and authored by Kirill Lykov on 2015-08-28. + * Copyright 2015. All rights reserved. + * + * Users are NOT authorized + * to employ the present software for their own publications + * before getting a written permission from the author of this file. +''' + +from plyfile import PlyData, PlyElement +import argparse +import copy +import numpy + +def computeExtent(vertices): + extentMax = [-10e6] * 3 + extentMin = [10e6] * 3 + for i in range(0, len(vertices)): + vi = vertices[i] + for dim in range(0, 3): + extentMax[dim] = max(extentMax[dim], vi[dim]) + extentMin[dim] = min(extentMin[dim], vi[dim]) + + origOrigin = [(extentMax[i] + extentMin[i])/2.0 for i in range(0, 3)] + origExtent = [extentMax[i] - extentMin[i] for i in range(0, 3)] + return (origOrigin, origExtent) + +parser = argparse.ArgumentParser(description='Modifies ply file to be used for rendering.\n Example: ./plyScale.py -f input.ply -o out.ply -r 210 -cutX 10.0') +parser.add_argument('-f','--inputFile', help='Input file name', required=True) +parser.add_argument('-o','--outputFile', help='Output file name', required=True) +parser.add_argument('--lx', help='Desired size of bounding box (X axis)', required=False, default="0") +parser.add_argument('--ly', help='Desired size of bounding box (Y axis)', required=False, default="0") +parser.add_argument('--lz', help='Desired size of bounding box (Z axis)', required=False, default="0") +parser.add_argument('-r','--order', help='Reorder axis. By default 012, to swap x and z use 210', required=False, default="012") +helpStringForCut = 'Remove all the faces which are above specified value for %s. The origin is in the center of mass. If axis reordering was applied, axis are in new coordinates.' +parser.add_argument('--cutX', help=helpStringForCut%('X'), required=False, default="none") +parser.add_argument('--cutY', help=helpStringForCut%('Y'), required=False, default="none") +parser.add_argument('--cutZ', help=helpStringForCut%('Z'), required=False, default="none") +args = vars(parser.parse_args()) + +desiredBox = [float(args['lx']), float(args['ly']), float(args['lz'])] + +plydata = PlyData.read(args['inputFile']) +vertices = plydata['vertex'].data + +# swap coords +order = args['order'] +if (order != "012"): + idx = [int(order[i]) for i in range(0, len(order))] + assert(len(idx) == 3) + print "Swapping axis!" + for i in range(0, len(vertices)): + v = copy.deepcopy(vertices[i]) + for dim in range(0, 3): + vertices[i][dim] = v[ idx[dim] ] + + +# Current box +(origOrigin, origExtent) = computeExtent(vertices) +for dim in range(0, 3): + if desiredBox[dim] == 0: + desiredBox[dim] = origExtent[dim] +print ("Extent is (%f, %f, %f). Center is (%f, %f, %f)."%(origExtent[0], origExtent[1], origExtent[2], + origOrigin[0], origOrigin[1], origOrigin[2])) +for i in range(0, len(vertices)): + for dim in range(0, 3): + vertices[i][dim] -= origOrigin[dim] + vertices[i][dim] *= desiredBox[dim]/origExtent[dim] + +if (args['cutX'] != "none" or args['cutY'] != "none" or args['cutZ'] != "none"): + print "Cut it!" + cut = [1e6]*3 + if args['cutX'] != "none": + cut[0] = float(args['cutX']) + if args['cutY'] != "none": + cut[1] = float(args['cutY']) + if args['cutZ'] != "none": + cut[2] = float(args['cutZ']) + + toDelete = list() + newInx = [None]*len(vertices) + j = 0 + for i in range(0, len(vertices)): + if (vertices[i][0] > cut[0] or vertices[i][1] > cut[1] or vertices[i][2] > cut[2]): + toDelete.append(i) + else: + newInx[i] = j + j += 1 + + plydata['vertex'].data = numpy.delete(plydata['vertex'].data, toDelete, axis=0) + # remove polygons containing these vertices + setVertToDel = set(toDelete) + faces = plydata['face'].data + facesToDelete = list() + for i in range(0, len(faces)): + curr = set(faces[i][0]) + common = curr & setVertToDel + if (common): + facesToDelete.append(i) + plydata['face'].data = numpy.delete(plydata['face'].data, facesToDelete, axis=0) + + #update vertices in polygons + faces = plydata['face'].data + for i in range(0, len(faces)): + f = faces[i][0] + for i in range(0, len(f)): + f[i] = newInx[ f[i] ] + +(finalOrigin, finalExtent) = computeExtent(vertices) +print ("Extent is (%f, %f, %f). Center is (%f, %f, %f)."%(finalExtent[0], finalExtent[1], finalExtent[2], + finalOrigin[0], finalOrigin[1], finalOrigin[2])) + +plydata.write(args['outputFile']) + diff --git a/device-gen/scripts/README b/device-gen/scripts/README deleted file mode 100644 index 6d30e57db..000000000 --- a/device-gen/scripts/README +++ /dev/null @@ -1,5 +0,0 @@ -Generates parabolic funnels. -./makeall.sh -./run.sh - -The result is in the file sdf.dat diff --git a/device-gen/scripts/cleanall.sh b/device-gen/scripts/cleanall.sh deleted file mode 100755 index b9e3a1b26..000000000 --- a/device-gen/scripts/cleanall.sh +++ /dev/null @@ -1,7 +0,0 @@ -#! /bin/bash -cd ../ -for d in */ ; do - pushd $d - make clean - popd -done diff --git a/device-gen/scripts/files.txt b/device-gen/scripts/files.txt deleted file mode 100644 index 217f12814..000000000 --- a/device-gen/scripts/files.txt +++ /dev/null @@ -1,12 +0,0 @@ -r4.dat -r5.dat -r6.dat -r7.dat -r8.dat -r9.dat -r10.dat -r11.dat -r12.dat -r13.dat -r14.dat -r15.dat \ No newline at end of file diff --git a/device-gen/scripts/makeall.sh b/device-gen/scripts/makeall.sh deleted file mode 100755 index 8c8723c2a..000000000 --- a/device-gen/scripts/makeall.sh +++ /dev/null @@ -1,7 +0,0 @@ -#! /bin/bash -cd ../ -for d in */ ; do - pushd $d - make - popd -done diff --git a/device-gen/scripts/run.sh b/device-gen/scripts/run.sh deleted file mode 100755 index ea81e33be..000000000 --- a/device-gen/scripts/run.sh +++ /dev/null @@ -1,19 +0,0 @@ -#! /usr/local/bin/bash - -unitXRes=32 -unitYRes=128 -unitZRes=64 - -nColumns=2 -# nRows is defined by files.txt - -for i in `seq 3 15`; do - ../sdf-unit-par/sdf-unit $unitXRes $unitYRes 24 96 $i gap$i.dat -done - -for i in `seq 3 15`; do - ../sdf-collage/sdf-collage gap$i.dat $nColumns 1 r$i.dat -done - -../sdf-collage/sdf-collage files.txt 1 1 collage.dat -../2Dto3D/2Dto3D collage.dat 40.0 4.0 $unitZRes sdf.dat diff --git a/device-gen/sdf-collage/Makefile b/device-gen/sdf-collage/Makefile deleted file mode 100644 index 306311353..000000000 --- a/device-gen/sdf-collage/Makefile +++ /dev/null @@ -1,7 +0,0 @@ -sdf-collage: main.cpp - g++ main.cpp -O0 -g3 -fopenmp -o sdf-collage - -clean: - rm sdf-collage - -.PHONY = clean diff --git a/device-gen/sdf-collage/main.cpp b/device-gen/sdf-collage/main.cpp deleted file mode 100644 index 2d2ac7dbc..000000000 --- a/device-gen/sdf-collage/main.cpp +++ /dev/null @@ -1,227 +0,0 @@ -/* - * main.cpp - * Part of uDeviceX/device-gen/sdf-collage/ - * - * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. - * Copyright 2015. All rights reserved. - * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. - */ - -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace std; - -#define _ACCESS(f, x, y) f[(x) + xsize * (y)] - -namespace Redistancing -{ - int xsize; - float * phi0, * phi; - float dt, invdx, invdy; - - template - inline bool anycrossing_dir(int ix, int iy, const float sgn0) - { - const int dx = d == 0, dy = d == 1, dz = d == 2; - - const float fm1 = _ACCESS(phi0, ix - dx, iy - dy); - const float fp1 = _ACCESS(phi0, ix + dx, iy + dy); - - return (fm1 * sgn0 < 0 || fp1 * sgn0 < 0); - } - - inline bool anycrossing(int ix, int iy, const float sgn0) - { - return - anycrossing_dir<0>(ix, iy, sgn0) || - anycrossing_dir<1>(ix, iy, sgn0) ; - } - - float sussman_scheme(int ix, int iy, float sgn0) - { - const float phicenter = _ACCESS(phi, ix, iy); - - const float dphidxm = phicenter - _ACCESS(phi, ix - 1, iy); - const float dphidxp = _ACCESS(phi, ix + 1, iy) - phicenter; - const float dphidym = phicenter - _ACCESS(phi, ix, iy - 1); - const float dphidyp = _ACCESS(phi, ix, iy + 1) - phicenter; - - if (sgn0 == 1) - { - const float xgrad0 = max( max((float)0, dphidxm), -min((float)0, dphidxp)) * invdx; - const float ygrad0 = max( max((float)0, dphidym), -min((float)0, dphidyp)) * invdy; - - const float G0 = sqrtf(xgrad0 * xgrad0 + ygrad0 * ygrad0) - 1; - - return phicenter - dt * sgn0 * G0; - } - else - { - const float xgrad1 = max( -min((float)0, dphidxm), max((float)0, dphidxp)) * invdx; - const float ygrad1 = max( -min((float)0, dphidym), max((float)0, dphidyp)) * invdy; - - const float G1 = sqrtf(xgrad1 * xgrad1 + ygrad1 * ygrad1) - 1; - - return phicenter - dt * sgn0 * G1; - } - } - - void redistancing(const int iterations, const float dt, const float dx, const float dy, - const int xsize, const int ysize, - float * field) - { - Redistancing::xsize = xsize; - Redistancing::dt = dt; - Redistancing::invdx = 1. / dx; - Redistancing::invdy = 1. / dy; - - Redistancing::phi0 = new float[xsize * ysize]; - memcpy(phi0, field, sizeof(float) * xsize * ysize); - Redistancing::phi = field; - - float * tmp = new float[xsize * ysize]; - for(int t = 0; t < iterations; ++t) - { - if (t % 30 == 0) - printf("t: %d\n", t); - -#pragma omp parallel for - for(int iy = 0; iy < ysize; ++iy) - for(int ix = 0; ix < xsize; ++ix) - { - const float myval0 = _ACCESS(phi0, ix, iy); - const float sgn0 = myval0 > 0 ? 1 : (myval0 < 0 ? -1 : 0); - - if (anycrossing(ix, iy, sgn0) || ix == 0 || ix == xsize - 1 || iy == 0 || iy == ysize - 1) - tmp[ix + xsize * iy] = myval0; - else - tmp[ix + xsize * iy] = sussman_scheme(ix, iy, sgn0); - } - - memcpy(field, tmp, sizeof(float) * xsize * ysize); - } - - delete [] tmp; - delete [] phi0; - phi0 = NULL; - } -} - -void mergeSDF(int NX, int NY, vector< vector >& cookie, vector& cake) -{ - cake.resize(cookie.size() * cookie[0].size()); - printf("SIZE: %d\n", cake.size()); - const int stride = NX; - for(int iy = 0; iy < cookie.size() * NY; ++iy) - for(int ix = 0; ix < NX; ++ix) - { - const int dst = ix + stride * iy; - const int iobst = iy / NY; - cake[dst] = cookie[iobst][ix + NX * (iy - iobst*NY)]; - } -} - -int main(int argc, char ** argv) -{ - if (argc != 5) - { - printf("usage: ./sdf-collage \n"); - return -1; - } - - const int xtimes = atoi(argv[2]); - int ytimes = atoi(argv[3]); - - float xextent, yextent, zextent; - int NX, NY,NZ; - vector< vector > cookie; - vector cake; - if (string(argv[1]) != "files.txt") - { - cookie.resize(1); - // for one file - FILE * f = fopen(argv[1], "r"); - assert(f != 0); - fscanf(f, "%f %f %f\n", &xextent, &yextent, &zextent); - fscanf(f, "%d %d %d\n", &NX, &NY, &NZ); - printf("Extent: [%f, %f, %f]. Grid size: [%d, %d, %d]\n", xextent, yextent, zextent, NX, NY,NZ); - assert(NZ == 1); - cookie[0].resize(NX * NY * NZ, 0.0f); - fread(&cookie[0][0], sizeof(float), NX * NY * NZ, f); - fclose(f); - - printf("Populate %d * %d times\n", xtimes, ytimes); - const int stride = xtimes * NX; - cake.resize(xtimes * ytimes * NX * NY); - - for(int ty = 0; ty < ytimes; ++ty) - for(int tx = 0; tx < xtimes; ++tx) { - for(int iy = 0; iy < NY; ++iy) - for(int ix = 0; ix < NX; ++ix) - { - const int gx = ix + NX * tx; - const int gy = iy + NY * ty; - const int dst = gx + stride * gy; - - assert(dst < cake.size()); - assert(ix + NX * iy < cookie[0].size()); - cake[dst] = cookie[0][ix + NX * iy]; - } - } - } else { - vector files; - - FILE* fs = fopen(argv[1], "r"); - assert(fs != 0); - string buf(127, ' '); - while(fscanf(fs, "%s\n", &buf[0]) == 1) { - files.push_back(buf); - } - fclose(fs); - - ytimes = files.size(); - cookie.resize(ytimes); - assert(xtimes == 1); - - for (int i = files.size() - 1; i >= 0; --i) - { - printf("Reading file %s ...\n", files[i].c_str()); - FILE * f = fopen(files[i].c_str(), "r"); - assert(f != 0); - fscanf(f, "%f %f %f\n", &xextent, &yextent, &zextent); - fscanf(f, "%d %d %d\n", &NX, &NY, &NZ); - printf("Extent: [%g, %g, %g]. Grid size: [%d, %d, %d]\n", xextent, yextent, zextent, NX, NY,NZ); - assert(NZ == 1); - cookie[i].resize(NX * NY * NZ, 0.0f); - fread(&cookie[i][0], sizeof(float), NX * NY * NZ, f); - fclose(f); - } - - mergeSDF(NX, NY, cookie, cake); - } - - const float dx = xextent / NX; - const float dy = yextent / NY; - Redistancing::redistancing(240, 0.25 * min(dx, dy), dx, dy, xtimes * NX, ytimes * NY, &cake[0]); - - { - FILE * f = fopen(argv[4], "w"); - assert(f != 0); - fprintf(f, "%f %f %f\n", xtimes * xextent, ytimes * yextent, 1.0f); - fprintf(f, "%d %d %d\n", xtimes * NX, ytimes * NY, 1); - fwrite(&cake[0], sizeof(float), cake.size(), f); - fclose(f); - } - - return 0; -} diff --git a/device-gen/sdf-unit-par/Makefile b/device-gen/sdf-unit-par/Makefile deleted file mode 100644 index 8d0bef89c..000000000 --- a/device-gen/sdf-unit-par/Makefile +++ /dev/null @@ -1,7 +0,0 @@ -sdf-unit: main.cpp - g++ main.cpp -o sdf-unit - -clean: - rm sdf-unit - -.PHONY = clean \ No newline at end of file diff --git a/device-gen/sdf-unit-par/main.cpp b/device-gen/sdf-unit-par/main.cpp deleted file mode 100644 index 8c5ceb639..000000000 --- a/device-gen/sdf-unit-par/main.cpp +++ /dev/null @@ -1,132 +0,0 @@ -/* - * main.cpp - * Part of uDeviceX/device-gen/sdf-unit-par/ - * - * Created and authored by Diego Rossinelli and Kirill Lykov on 2015-03-20. - * Copyright 2015. All rights reserved. - * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. - */ - -#include -#include -#include -#include -#include -using namespace std; - -struct Parabola -{ - float x0; - float ymax; // length of obstacle, for cutting the pick - float y0; - - Parabola(float gap, float xextent) : ymax(48.0) - { - x0 = xextent/2.0f - gap/2.0f; - if (gap > 9.0) - y0 = ymax; - else - y0 = ymax + 6.5; - } - - void line1(vector& vx, vector& vy) { - int N = 500; - float dx = 2.0f * fabs(x0) / (N - 1); - for (int i = 0; i < N; ++i) { - float x = i * dx - x0; - float y = 0.0f; - vx.push_back(x); - vy.push_back(y); - } - } - - void line2(vector& vx, vector& vy) { - int N = 1500; - float dx = 2.0f * fabs(x0) / (N - 1); - float alpha = -y0 / (x0 * x0); - for (int i = 0; i < N; ++i) { - float x = i * dx - x0; - float y = min(ymax, alpha * x*x + y0); - vx.push_back(x); - vy.push_back(y); - } - } -}; - -int main(int argc, char ** argv) -{ - if (argc != 7) - { - printf("usage: ./sdf-unit \n"); - return 1; - } - - const int NX = atoi(argv[1]); - const int NY = atoi(argv[2]); - const float xextent = atof(argv[3]); - const float yextent = atof(argv[4]); - const float gap = atof(argv[5]); - - vector xs, ys; - Parabola par(gap, xextent); - par.line1(xs, ys); - par.line2(xs, ys); - - const float xlb = -xextent/2.0f; - const float ylb = -(yextent - par.ymax)/2.0f; - printf("starting brute force sdf with %d x %d starting from %f %f to %f %f\n", - NX, NY, xlb, ylb, xlb + xextent, ylb + yextent); - - float * sdf = new float[NX * NY]; - const float dx = xextent / NX; - const float dy = yextent / NY; - const int nsamples = xs.size(); - - for(int iy = 0; iy < NY; ++iy) - for(int ix = 0; ix < NX; ++ix) - { - const float x = xlb + ix * dx; - const float y = ylb + iy * dy; - - float distance2 = 1e6; - int iclosest = 0; - for(int i = 0; i < nsamples ; ++i) - { - const float xd = xs[i] - x; - const float yd = ys[i] - y; - const float candidate = xd * xd + yd * yd; - - if (candidate < distance2) - { - iclosest = i; - distance2 = candidate; - } - } - - float s = -1; - - { - const float alpha = -par.y0 / (par.x0 * par.x0); - const float ycurve = min(par.ymax, alpha * x*x + par.y0); - - if (x >= -par.x0 && x <= par.x0 && y >= 0 && y <= ycurve) - s = +1; - } - - - sdf[ix + NX * iy] = s * sqrt(distance2); - } - - FILE * f = fopen(argv[6], "w"); - fprintf(f, "%f %f %f\n", xextent, yextent, 1.0f); - fprintf(f, "%d %d %d\n", NX, NY, 1); - fwrite(sdf, sizeof(float), NX * NY, f); - fclose(f); - - delete [] sdf; - - return 0; -} diff --git a/mpi-dpd/README.md b/mpi-dpd/README.md new file mode 100644 index 000000000..146275372 --- /dev/null +++ b/mpi-dpd/README.md @@ -0,0 +1,51 @@ +The main folder: mpi-dpd +======= + +This is the "main" folder containing the object files orchestrating the various kernels. +The generated object files are the following ones: + +* common.o: global simulation parameters, common datastructures like device arrays, cell lists etc. +* contact.o: computation of the contact/lubrication force across "touching" solute particles +* containers.o: particle arrays, collections encapsulating the data of RBCs and CTCs +* dpd.o: code coordinating the computation of cuda-dpd +* fsi.o: computation of the "flow-structure interaction" force between solvent and solute particles +* io.o: data dumps in XYZ, PLY (for cells), H5Part, HDF5 structured grids +* main.o: home sweet home +* minmax.o: computation of the extent of an array of RBCs or CTCs +* redistancing.o: computation of the distance transform for the implicit description of the wall geometry +* redistribute-particles.o: redistribution of the solvent across the MPI ranks once particles have moved +* redistribute-rbcs.o: redistribution of RBCs across MPI ranks once they have moved +* scan.o: computation of the prefix sum for the solvent cell lists count +* simulation.o: simulation "driver" coordinating the other object files, except for main.cu +* solute-exchange.o: exchange the "halo" solute particles across the MPI ranks close by to compute FSI and contact forces. +* solvent-exchange.o: exchange the "halo" solvent particles across the MPI ranks to compute the DPD interactions +* wall.o: computation of the particles interacting with the no-slip boundary conditions of the wall + +Compiling uDeviceX +------------- +The makefile will check for a .cache.Makefile, that can be optionally put in this folder. +For example my .cache.Makefile on Piz Daint is + +`h5part = 0 +NVCC = nvcc -I$(CRAY_MPICH2_DIR)/include -L$(CRAY_MPICH2_DIR)/lib -I/scratch/daint/diegor/h5part/include -I$(HDF5_DIR)/include -I/users/diegor/vtk/install/include/vtk-6.2 -I/users/diegor/h5part/include/ +CXX = CC $(CRAY_CUDATOOLKIT_POST_LINK_OPTS) $(CRAY_CUDATOOLKIT_INCLUDE_OPTS) -L/users/diegor/h5part/lib -L$(HDF5_DIR)/lib -L/users/diegor/vtk/install/lib` + +To clean uDeviceX entirely (mpi-dpd, cuda-dpd, cuda-rbc, cuda-ctc): +`make cleanall` + +To just cleanup the mpi-dpd folder: +`make clean` + +To compile: +`make -j` + +Running uDeviceX +----------- +Running uDeviceX consists of these steps: + +1. Generation of the geometry file (optional), the file should always be named `sdf.dat` (as Signed Distance Function). +See the folder `device-gen`. +2. Generation of the initial positioning of the RBCs and CTCs (optional). The IC files should be called +`rbcs-ic.txt and ctcs-ic.txt` See the folder `cell-placement`. +3. Execution of uDeviceX, for example `mpirun ./test 4 4 2 -walls -couette=1 -tend=5e4 -steps_per_dump=1000 -rbcs -contactforces` +4. Post processing of the simulation data (optional), for example `ls ./stress/* -rt1 | tail -n 50 | mpirun -n 32 -N 1 ../postprocessing/stress/stress -origin=0,0,5 -extent=192,192,85 -project=1,1,0 > stress-profile.txt` diff --git a/mpi-dpd/common-kernels.h b/mpi-dpd/common-kernels.h index a2b8a1572..59c165675 100644 --- a/mpi-dpd/common-kernels.h +++ b/mpi-dpd/common-kernels.h @@ -170,8 +170,7 @@ void write_AOS3f(float * const data, const int nparticles, float& s0, float& s1, data[laneid + 64] = s2; } -template -__global__ void subindex_local(const int nparticles, const float2 * particles, int * const partials, +__global__ static void subindex_local(const int nparticles, const float2 * particles, int * const partials, uchar4 * const subindices) { assert(blockDim.x == 128 && blockDim.x * gridDim.x >= nparticles); @@ -192,29 +191,18 @@ __global__ void subindex_local(const int nparticles, const float2 * particles, read_AOS6f(particles + 3 * base, nsrc, data0, data1, data2); - const bool inside = project || + const bool inside = (data0.x >= -XSIZE_SUBDOMAIN / 2 && data0.x < XSIZE_SUBDOMAIN / 2 && data0.y >= -YSIZE_SUBDOMAIN / 2 && data0.y < YSIZE_SUBDOMAIN / 2 && data1.x >= -ZSIZE_SUBDOMAIN / 2 && data1.x < ZSIZE_SUBDOMAIN / 2 ); if (lane < nsrc && inside) { - if (project) - { - const int xcid = min(XSIZE_SUBDOMAIN - 1, max(0, (int)floor((double)data0.x + XSIZE_SUBDOMAIN / 2))); - const int ycid = min(YSIZE_SUBDOMAIN - 1, max(0, (int)floor((double)data0.y + YSIZE_SUBDOMAIN / 2))); - const int zcid = min(ZSIZE_SUBDOMAIN - 1, max(0, (int)floor((double)data1.x + ZSIZE_SUBDOMAIN / 2))); + const int xcid = (int)floor((double)data0.x + XSIZE_SUBDOMAIN / 2); + const int ycid = (int)floor((double)data0.y + YSIZE_SUBDOMAIN / 2); + const int zcid = (int)floor((double)data1.x + ZSIZE_SUBDOMAIN / 2); - cid = xcid + XSIZE_SUBDOMAIN * (ycid + YSIZE_SUBDOMAIN * zcid); - } - else - { - const int xcid = (int)floor((double)data0.x + XSIZE_SUBDOMAIN / 2); - const int ycid = (int)floor((double)data0.y + YSIZE_SUBDOMAIN / 2); - const int zcid = (int)floor((double)data1.x + ZSIZE_SUBDOMAIN / 2); - - cid = xcid + XSIZE_SUBDOMAIN * (ycid + YSIZE_SUBDOMAIN * zcid); - } + cid = xcid + XSIZE_SUBDOMAIN * (ycid + YSIZE_SUBDOMAIN * zcid); } } diff --git a/mpi-dpd/common.cu b/mpi-dpd/common.cu index d6155f56c..aa86dce28 100644 --- a/mpi-dpd/common.cu +++ b/mpi-dpd/common.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-01-30. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/mpi-dpd/common.h b/mpi-dpd/common.h index 481592ee6..53f7ed2a6 100644 --- a/mpi-dpd/common.h +++ b/mpi-dpd/common.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli, on 2014-12-05. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -37,8 +34,8 @@ const float sigmaf = sigma / sqrt(dt); const float aij = 25; const float hydrostatic_a = 0.05; -extern float tend; -extern bool walls, pushtheflow, doublepoiseuille, rbcs, ctcs, xyz_dumps, hdf5field_dumps, hdf5part_dumps, is_mps_enabled, contactforces; +extern float tend, couette; +extern bool walls, pushtheflow, doublepoiseuille, rbcs, ctcs, xyz_dumps, hdf5field_dumps, hdf5part_dumps, is_mps_enabled, contactforces, stress; extern int steps_per_report, steps_per_dump, wall_creation_stepid, nvtxstart, nvtxstop; #include diff --git a/mpi-dpd/contact.cu b/mpi-dpd/contact.cu index 42059759f..8f4dbf67d 100644 --- a/mpi-dpd/contact.cu +++ b/mpi-dpd/contact.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ static const int maxsolutes = 32; @@ -42,11 +39,6 @@ namespace KernelsContact texture texCellsStart, texCellEntries; - __global__ void bulk_3tpp(const float2 * const particles, const int np, const int ncellentries, const int nsolutes, - float * const acc, const float seed, const int mysoluteid); - - __global__ void halo(const int nparticles_padded, const int ncellentries, const int nsolutes, const float seed); - void setup() { texCellsStart.channelDesc = cudaCreateChannelDesc(); @@ -58,10 +50,9 @@ namespace KernelsContact texCellEntries.filterMode = cudaFilterModePoint; texCellEntries.mipmapFilterMode = cudaFilterModePoint; texCellEntries.normalized = 0; - - CUDA_CHECK(cudaFuncSetCacheConfig(bulk_3tpp, cudaFuncCachePreferL1)); - CUDA_CHECK(cudaFuncSetCacheConfig(halo, cudaFuncCachePreferL1)); } + + __global__ void bulk_3tpp(const int nsolutes, const float seed); } ComputeContact::ComputeContact(MPI_Comm comm): @@ -77,11 +68,13 @@ cellsstart(KernelsContact::NCELLS + 16), cellscount(KernelsContact::NCELLS + 16) CUDA_CHECK(cudaMemcpyToSymbol(KernelsContact::params, ¶ms, sizeof(params))); CUDA_CHECK(cudaPeekAtLastError()); + + CUDA_CHECK(cudaFuncSetCacheConfig(KernelsContact::bulk_3tpp , cudaFuncCachePreferL1)); } namespace KernelsContact { - __global__ void populate(const uchar4 * const subindices, const int * const cellstart, + __global__ void populate(const uchar4 * const subindices, const int * const cellstart, const int nparticles, const int soluteid, const int ntotalparticles, CellEntry * const entrycells) { @@ -141,104 +134,66 @@ namespace KernelsContact const int ncells = XSIZE_SUBDOMAIN * YSIZE_SUBDOMAIN * ZSIZE_SUBDOMAIN; - CUDA_CHECK(cudaBindTexture(&textureoffset, &texCellsStart, cellsstart, &texCellsStart.channelDesc, sizeof(int) * ncells)); + CUDA_CHECK(cudaBindTexture(&textureoffset, &texCellsStart, cellsstart, &texCellsStart.channelDesc, sizeof(int) * (1 + ncells))); assert(textureoffset == 0); const int n = wsolutes.size(); - int ns[n]; - float2 * ps[n]; - float * as[n]; - - for(int i = 0; i < n; ++i) + if (n) { - ns[i] = wsolutes[i].n; - ps[i] = (float2 *)wsolutes[i].p; - as[i] = (float * )wsolutes[i].a; - } - - CUDA_CHECK(cudaMemcpyToSymbolAsync(cnsolutes, ns, sizeof(int) * n, 0, cudaMemcpyHostToDevice, stream)); - CUDA_CHECK(cudaMemcpyToSymbolAsync(csolutes, ps, sizeof(float2 *) * n, 0, cudaMemcpyHostToDevice, stream)); - CUDA_CHECK(cudaMemcpyToSymbolAsync(csolutesacc, as, sizeof(float *) * n, 0, cudaMemcpyHostToDevice, stream)); - } -} - -void ComputeContact::build_cells(std::vector wsolutes, cudaStream_t stream) -{ - this->nsolutes = wsolutes.size(); - - int ntotal = 0; - - for(int i = 0; i < wsolutes.size(); ++i) - ntotal += wsolutes[i].n; - - subindices.resize(ntotal); - cellsentries.resize(ntotal); + int ns[n]; + float2 * ps[n]; + float * as[n]; - CUDA_CHECK(cudaMemsetAsync(cellscount.data, 0, sizeof(int) * cellscount.size, stream)); - -#ifndef NDEBUG - CUDA_CHECK(cudaMemsetAsync(cellsentries.data, 0xff, sizeof(int) * cellsentries.capacity, stream)); - CUDA_CHECK(cudaMemsetAsync(subindices.data, 0xff, sizeof(int) * subindices.capacity, stream)); - CUDA_CHECK(cudaMemsetAsync(compressed_cellscount.data, 0xff, sizeof(unsigned char) * compressed_cellscount.capacity, stream)); - CUDA_CHECK(cudaMemsetAsync(cellsstart.data, 0xff, sizeof(int) * cellsstart.capacity, stream)); -#endif - - CUDA_CHECK(cudaPeekAtLastError()); - - int ctr = 0; - for(int i = 0; i < wsolutes.size(); ++i) - { - const ParticlesWrap it = wsolutes[i]; - - if (it.n) - subindex_local<<< (it.n + 127) / 128, 128, 0, stream >>> - (it.n, (float2 *)it.p, cellscount.data, subindices.data + ctr); + for(int i = 0; i < n; ++i) + { + ns[i] = wsolutes[i].n; + ps[i] = (float2 *)wsolutes[i].p; + as[i] = (float * )wsolutes[i].a; + } - ctr += it.n; + CUDA_CHECK(cudaMemcpyToSymbolAsync(cnsolutes, ns, sizeof(int) * n, 0, cudaMemcpyHostToDevice, stream)); + CUDA_CHECK(cudaMemcpyToSymbolAsync(csolutes, ps, sizeof(float2 *) * n, 0, cudaMemcpyHostToDevice, stream)); + CUDA_CHECK(cudaMemcpyToSymbolAsync(csolutesacc, as, sizeof(float *) * n, 0, cudaMemcpyHostToDevice, stream)); + } } - compress_counts<<< (compressed_cellscount.size + 127) / 128, 128, 0, stream >>> - (compressed_cellscount.size, (int4 *)cellscount.data, (uchar4 *)compressed_cellscount.data); + __global__ __launch_bounds__(128, 10) void bulk_3tpp(const int nsolutes, const float seed) + { + const int np = tex1Dfetch(texCellsStart, XCELLS * YCELLS * ZCELLS); - scan(compressed_cellscount.data, compressed_cellscount.size, stream, (uint *)cellsstart.data); + assert(blockDim.x * gridDim.x >= np * 3); - ctr = 0; - for(int i = 0; i < wsolutes.size(); ++i) - { - const ParticlesWrap it = wsolutes[i]; + const int gid = threadIdx.x + blockDim.x * blockIdx.x; + const int myslot = gid / 3; + const int zplane = gid % 3; - if (it.n) - KernelsContact::populate<<< (it.n + 127) / 128, 128, 0, stream >>> - (subindices.data + ctr, cellsstart.data, it.n, i, ntotal, (KernelsContact::CellEntry *)cellsentries.data); + if (myslot >= np) + return; - ctr += it.n; - } + float2 dst0, dst1, dst2; + int mysoluteid, actualpid; - CUDA_CHECK(cudaPeekAtLastError()); + { + CellEntry ce; + ce.pid = tex1Dfetch(texCellEntries, myslot); - KernelsContact::bind(cellsstart.data, cellsentries.data, ntotal, wsolutes, stream, cellscount.data); -} + mysoluteid = ce.code.w; -namespace KernelsContact -{ - __global__ __launch_bounds__(128, 10) - void bulk_3tpp(const float2 * const particles, - const int np, const int ncellentries, const int nsolutes, - float * const acc, const float seed, const int mysoluteid) - { - assert(blockDim.x * gridDim.x >= np * 3); + ce.code.w = 0; + actualpid = ce.pid; - const int gid = threadIdx.x + blockDim.x * blockIdx.x; - const int pid = gid / 3; - const int zplane = gid % 3; + assert(mysoluteid < nsolutes); + assert(actualpid >= 0 && actualpid < cnsolutes[mysoluteid]); - if (pid >= np) - return; + dst0 = _ACCESS(csolutes[mysoluteid] + 3 * actualpid + 0); + dst1 = _ACCESS(csolutes[mysoluteid] + 3 * actualpid + 1); + dst2 = _ACCESS(csolutes[mysoluteid] + 3 * actualpid + 2); - const float2 dst0 = _ACCESS(particles + 3 * pid + 0); - const float2 dst1 = _ACCESS(particles + 3 * pid + 1); - const float2 dst2 = _ACCESS(particles + 3 * pid + 2); + assert(dst0.x >= -XOFFSET && dst0.x < XOFFSET); + assert(dst0.y >= -YOFFSET && dst0.y < YOFFSET); + assert(dst1.x >= -ZOFFSET && dst1.x < ZOFFSET); + } int scan1, scan2, ncandidates, spidbase; int deltaspid1, deltaspid2; @@ -295,13 +250,15 @@ namespace KernelsContact float xforce = 0, yforce = 0, zforce = 0; -#pragma unroll 3 for(int i = 0; i < ncandidates; ++i) { const int m1 = (int)(i >= scan1); const int m2 = (int)(i >= scan2); const int slot = i + (m2 ? deltaspid2 : m1 ? deltaspid1 : spidbase); - assert(slot >= 0 && slot < ncellentries); + assert(slot >= 0 && slot < np); + + if (slot >= myslot) + continue; CellEntry ce; ce.pid = tex1Dfetch(texCellEntries, slot); @@ -313,9 +270,6 @@ namespace KernelsContact const int spid = ce.pid; assert(spid >= 0 && spid < cnsolutes[soluteid]); - if (mysoluteid < soluteid || mysoluteid == soluteid && pid <= spid) - continue; - const int sentry = 3 * spid; const float2 stmp0 = _ACCESS(csolutes[soluteid] + sentry ); const float2 stmp1 = _ACCESS(csolutes[soluteid] + sentry + 1); @@ -339,7 +293,7 @@ namespace KernelsContact const float t2 = ljsigma2 * invr2; const float t4 = t2 * t2; const float t6 = t4 * t2; - const float lj = min(1e3f, max(0.f, 24.f * invrij * t6 * (2.f * t6 - 1.f))); + const float lj = min(1e4f, max(0.f, 24.f * invrij * t6 * (2.f * t6 - 1.f))); const float wr = viscosity_function<-VISCOSITY_S_LEVEL>(1.f - rij); @@ -352,7 +306,7 @@ namespace KernelsContact yr * (dst2.x - stmp2.x) + zr * (dst2.y - stmp2.y); - const float myrandnr = Logistic::mean0var1(seed, pid, spid); + const float myrandnr = Logistic::mean0var1(seed, myslot, slot); const float strength = lj + (- params.gamma * wr * rdotv + params.sigmaf * myrandnr) * wr; @@ -368,93 +322,57 @@ namespace KernelsContact assert(!isnan(yinteraction)); assert(!isnan(zinteraction)); - assert(fabs(xinteraction) < 1e4); - assert(fabs(yinteraction) < 1e4); - assert(fabs(zinteraction) < 1e4); + assert(fabs(xinteraction) < 1e5); + assert(fabs(yinteraction) < 1e5); + assert(fabs(zinteraction) < 1e5); atomicAdd(csolutesacc[soluteid] + sentry , -xinteraction); atomicAdd(csolutesacc[soluteid] + sentry + 1, -yinteraction); atomicAdd(csolutesacc[soluteid] + sentry + 2, -zinteraction); } - atomicAdd(acc + 3 * pid + 0, xforce); - atomicAdd(acc + 3 * pid + 1, yforce); - atomicAdd(acc + 3 * pid + 2, zforce); + const float xacc = atomicAdd(csolutesacc[mysoluteid] + 3 * actualpid + 0, xforce); + const float yacc = atomicAdd(csolutesacc[mysoluteid] + 3 * actualpid + 1, yforce); + const float zacc = atomicAdd(csolutesacc[mysoluteid] + 3 * actualpid + 2, zforce); - for(int c = 0; c < 3; ++c) - assert(!isnan(acc[3 * pid + c])); + assert(!isnan(xacc)); + assert(!isnan(yacc)); + assert(!isnan(zacc)); } -} - -void ComputeContact::bulk(std::vector wsolutes, cudaStream_t stream) -{ - NVTX_RANGE("Contact/bulk", NVTX_C6); - if (wsolutes.size() == 0) - return; - - for(int i = 0; i < wsolutes.size(); ++i) + __global__ void halo(const float2 * halo, const int nhalo, const int nsolutes, const float seed, float * const acc) { - ParticlesWrap it = wsolutes[i]; - - if (it.n) - KernelsContact::bulk_3tpp<<< (3 * it.n + 127) / 128, 128, 0, stream >>> - ((float2 *)it.p, it.n, cellsentries.size, wsolutes.size(), (float *)it.a, local_trunk.get_float(), i); - - CUDA_CHECK(cudaPeekAtLastError()); - } -} - -namespace KernelsContact -{ - __constant__ int packstarts_padded[27], packcount[26]; - __constant__ Particle * packstates[26]; - __constant__ Acceleration * packresults[26]; + const int nbulk = tex1Dfetch(texCellsStart, XCELLS * YCELLS * ZCELLS); - __global__ void halo(const int nparticles_padded, const int ncellentries, const int nsolutes, const float seed) - { - assert(blockDim.x * gridDim.x >= nparticles_padded); + assert(blockDim.x * gridDim.x >= nhalo); const int laneid = threadIdx.x & 0x1f; const int warpid = threadIdx.x >> 5; - const int localbase = 32 * (warpid + 4 * blockIdx.x); - const int pid = localbase + laneid; - - if (localbase >= nparticles_padded) - return; + const int unpackbase = 32 * (warpid + 4 * blockIdx.x); + const int nunpack = min(32, nhalo - unpackbase); - int nunpack; float2 dst0, dst1, dst2; - float * dst = NULL; + read_AOS6f((float2 *)(halo + 3 * unpackbase), nunpack, dst0, dst1, dst2); - { - const uint key9 = 9 * (localbase >= packstarts_padded[9]) + 9 * (localbase >= packstarts_padded[18]); - const uint key3 = 3 * (localbase >= packstarts_padded[key9 + 3]) + 3 * (localbase >= packstarts_padded[key9 + 6]); - const uint key1 = (localbase >= packstarts_padded[key9 + key3 + 1]) + (localbase >= packstarts_padded[key9 + key3 + 2]); - const int code = key9 + key3 + key1; - assert(code >= 0 && code < 26); - assert(localbase >= packstarts_padded[code] && localbase < packstarts_padded[code + 1]); - - const int unpackbase = localbase - packstarts_padded[code]; - assert (unpackbase >= 0); - assert(unpackbase < packcount[code]); - - nunpack = min(32, packcount[code] - unpackbase); - - if (nunpack == 0) - return; + float xforce, yforce, zforce; + read_AOS3f(acc + 3 * unpackbase, nunpack, xforce, yforce, zforce); - read_AOS6f((float2 *)(packstates[code] + unpackbase), nunpack, dst0, dst1, dst2); + const bool outside_plus = + dst0.x >= XOFFSET || + dst0.x >= -XOFFSET && dst0.y >= YOFFSET || + dst0.x >= -XOFFSET && dst0.y >= -YOFFSET && dst1.x >= ZOFFSET; - dst = (float*)(packresults[code] + unpackbase); - } + const bool inside_outerhalo = + dst0.x < XOFFSET + 1 && + dst0.y < YOFFSET + 1 && + dst1.x < ZOFFSET + 1 ; - float xforce, yforce, zforce; - read_AOS3f(dst, nunpack, xforce, yforce, zforce); + const bool valid = laneid < nunpack && outside_plus && inside_outerhalo; - const int nzplanes = laneid < nunpack ? 3 : 0; + if (!valid) + return; - for(int zplane = 0; zplane < nzplanes; ++zplane) + for(int zplane = 0; zplane < 3; ++zplane) { int scan1, scan2, ncandidates, spidbase; int deltaspid1, deltaspid2; @@ -470,8 +388,8 @@ namespace KernelsContact assert(xcount >= 0); const int ycenter = YOFFSET + (int)floorf(dst0.y); - const int zcenter = ZOFFSET + (int)floorf(dst1.x); + const int zmy = zcenter - 1 + zplane; const bool zvalid = zmy >= 0 && zmy < ZCELLS; @@ -515,7 +433,7 @@ namespace KernelsContact const int m2 = (int)(i >= scan2); const int slot = i + (m2 ? deltaspid2 : m1 ? deltaspid1 : spidbase); - assert(slot >= 0 && slot < ncellentries); + assert(slot >= 0 && slot < nbulk); CellEntry ce; ce.pid = tex1Dfetch(texCellEntries, slot); const int soluteid = ce.code.w; @@ -548,7 +466,7 @@ namespace KernelsContact const float t2 = ljsigma2 * invr2; const float t4 = t2 * t2; const float t6 = t4 * t2; - const float lj = min(1e3f, max(0.f, 24.f * invrij * t6 * (2.f * t6 - 1.f))); + const float lj = min(1e4f, max(0.f, 24.f * invrij * t6 * (2.f * t6 - 1.f))); const float wr = viscosity_function<-VISCOSITY_S_LEVEL>(1.f - rij); @@ -561,7 +479,7 @@ namespace KernelsContact yr * (dst2.x - stmp2.x) + zr * (dst2.y - stmp2.y); - const float myrandnr = Logistic::mean0var1(seed, pid, spid); + const float myrandnr = Logistic::mean0var1(seed, unpackbase + laneid, spid); const float strength = lj + (- params.gamma * wr * rdotv + params.sigmaf * myrandnr) * wr; @@ -577,9 +495,9 @@ namespace KernelsContact assert(!isnan(yinteraction)); assert(!isnan(zinteraction)); - assert(fabs(xinteraction) < 1e4); - assert(fabs(yinteraction) < 1e4); - assert(fabs(zinteraction) < 1e4); + assert(fabs(xinteraction) < 1e5); + assert(fabs(yinteraction) < 1e5); + assert(fabs(zinteraction) < 1e5); atomicAdd(csolutesacc[soluteid] + sentry , -xinteraction); atomicAdd(csolutesacc[soluteid] + sentry + 1, -yinteraction); @@ -587,54 +505,87 @@ namespace KernelsContact } } - write_AOS3f(dst, nunpack, xforce, yforce, zforce); + acc[3 * (unpackbase + laneid) + 0] = xforce; + acc[3 * (unpackbase + laneid) + 1] = yforce; + acc[3 * (unpackbase + laneid) + 2] = zforce; } } -void ComputeContact::halo(ParticlesWrap halos[26], cudaStream_t stream) +void ComputeContact::halo(ParticlesWrap halowrap, cudaStream_t stream) { NVTX_RANGE("Contact/halo", NVTX_C7); - int nremote_padded = 0; + CUDA_CHECK(cudaPeekAtLastError()); - { - int recvpackcount[26], recvpackstarts_padded[27]; + wsolutes.push_back(halowrap); - for(int i = 0; i < 26; ++i) - recvpackcount[i] = halos[i].n; + int ntotal = 0; - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsContact::packcount, recvpackcount, - sizeof(recvpackcount), 0, cudaMemcpyHostToDevice, stream)); + for(int i = 0; i < wsolutes.size(); ++i) + ntotal += wsolutes[i].n; - recvpackstarts_padded[0] = 0; - for(int i = 0, s = 0; i < 26; ++i) - recvpackstarts_padded[i + 1] = (s += 32 * ((halos[i].n + 31) / 32)); + subindices.resize(ntotal); + cellsentries.resize(ntotal); - nremote_padded = recvpackstarts_padded[26]; + CUDA_CHECK(cudaMemsetAsync(cellscount.data, 0, sizeof(int) * cellscount.size, stream)); - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsContact::packstarts_padded, recvpackstarts_padded, - sizeof(recvpackstarts_padded), 0, cudaMemcpyHostToDevice, stream)); +#ifndef NDEBUG + CUDA_CHECK(cudaMemsetAsync(cellsentries.data, 0xff, sizeof(int) * cellsentries.capacity, stream)); + CUDA_CHECK(cudaMemsetAsync(subindices.data, 0xff, sizeof(int) * subindices.capacity, stream)); + CUDA_CHECK(cudaMemsetAsync(compressed_cellscount.data, 0xff, sizeof(unsigned char) * compressed_cellscount.capacity, stream)); + CUDA_CHECK(cudaMemsetAsync(cellsstart.data, 0xff, sizeof(int) * cellsstart.capacity, stream)); +#endif - const Particle * recvpackstates[26]; + CUDA_CHECK(cudaPeekAtLastError()); - for(int i = 0; i < 26; ++i) - recvpackstates[i] = halos[i].p; + int ctr = 0; + for(int i = 0; i < wsolutes.size(); ++i) + { + const ParticlesWrap it = wsolutes[i]; - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsContact::packstates, recvpackstates, - sizeof(recvpackstates), 0, cudaMemcpyHostToDevice, stream)); + if (it.n) + subindex_local<<< (it.n + 127) / 128, 128, 0, stream >>> + (it.n, (float2 *)it.p, cellscount.data, subindices.data + ctr); - Acceleration * packresults[26]; + ctr += it.n; + } + + compress_counts<<< (compressed_cellscount.size + 127) / 128, 128, 0, stream >>> + (compressed_cellscount.size, (int4 *)cellscount.data, (uchar4 *)compressed_cellscount.data); + + scan(compressed_cellscount.data, compressed_cellscount.size, stream, (uint *)cellsstart.data); + + ctr = 0; + for(int i = 0; i < wsolutes.size(); ++i) + { + const ParticlesWrap it = wsolutes[i]; - for(int i = 0; i < 26; ++i) - packresults[i] = halos[i].a; + if (it.n) + KernelsContact::populate<<< (it.n + 127) / 128, 128, 0, stream >>> + (subindices.data + ctr, cellsstart.data, it.n, i, ntotal, (KernelsContact::CellEntry *)cellsentries.data); - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsContact::packresults, packresults, - sizeof(packresults), 0, cudaMemcpyHostToDevice, stream)); + ctr += it.n; } - if(nremote_padded) - KernelsContact::halo<<< (nremote_padded + 127) / 128, 128, 0, stream>>> - (nremote_padded, cellsentries.size, nsolutes, local_trunk.get_float()); + CUDA_CHECK(cudaPeekAtLastError()); + + KernelsContact::bind(cellsstart.data, cellsentries.data, ntotal, wsolutes, stream, cellscount.data); + + if (cellsentries.size) + KernelsContact::bulk_3tpp<<< (3 * cellsentries.size + 127) / 128, 128, 0, stream >>> + (wsolutes.size(), local_trunk.get_float()); + + ctr = 0; + for(int i = 0; i < wsolutes.size(); ++i) + { + const ParticlesWrap it = wsolutes[i]; + + if (it.n) + KernelsContact::halo<<< (it.n + 127) / 128, 128, 0, stream>>> + ((float2 *)it.p, it.n, wsolutes.size(), local_trunk.get_float(), (float *)it.a); + + ctr += it.n; + } CUDA_CHECK(cudaPeekAtLastError()); } diff --git a/mpi-dpd/contact.h b/mpi-dpd/contact.h index d02aa692e..fb51b3c58 100644 --- a/mpi-dpd/contact.h +++ b/mpi-dpd/contact.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -21,24 +18,20 @@ class ComputeContact : public SoluteExchange::Visitor { - //cudaEvent_t evuploaded; - - int nsolutes; - + std::vector wsolutes; + SimpleDeviceBuffer subindices; SimpleDeviceBuffer compressed_cellscount; SimpleDeviceBuffer cellsentries, cellsstart, cellscount; Logistic::KISS local_trunk; - + public: ComputeContact(MPI_Comm comm); - void build_cells(std::vector wsolutes, cudaStream_t stream); - - void bulk(std::vector wsolutes, cudaStream_t stream); + void attach_bulk(std::vector wsolutes) { this->wsolutes = wsolutes; } /*override of SoluteExchange::Visitor::halo*/ - void halo(ParticlesWrap solutes[26], cudaStream_t stream); + void halo(ParticlesWrap allhalos, cudaStream_t stream); }; diff --git a/mpi-dpd/containers.cu b/mpi-dpd/containers.cu index dfd96c564..21a06f7d2 100644 --- a/mpi-dpd/containers.cu +++ b/mpi-dpd/containers.cu @@ -6,9 +6,6 @@ * Created and authored by Diego Rossinelli on 2014-12-05. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/mpi-dpd/containers.h b/mpi-dpd/containers.h index 0bd850d44..70224df51 100644 --- a/mpi-dpd/containers.h +++ b/mpi-dpd/containers.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-05. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/ctc.h b/mpi-dpd/ctc.h index b9f2d4293..09a4b567a 100644 --- a/mpi-dpd/ctc.h +++ b/mpi-dpd/ctc.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-18. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/dpd.cu b/mpi-dpd/dpd.cu index c573642b8..b71eace14 100644 --- a/mpi-dpd/dpd.cu +++ b/mpi-dpd/dpd.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-04. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -20,11 +17,15 @@ using namespace std; -ComputeDPD::ComputeDPD(MPI_Comm cartcomm): SolventExchange(cartcomm, 0), local_trunk(0, 0, 0, 0) +ComputeDPD::ComputeDPD(MPI_Comm cartcomm): +SolventExchange(cartcomm, 0), local_trunk(0, 0, 0, 0), +sigma_xx(NULL), sigma_xy(NULL), sigma_xz(NULL), sigma_yy(NULL), sigma_yz(NULL), sigma_zz(NULL) { int myrank; MPI_CHECK(MPI_Comm_rank(cartcomm, &myrank)); + local_trunk = Logistic::KISS(9078 - 2 * myrank, 321 - myrank, 552, 456); + for(int i = 0; i < 26; ++i) { int d[3] = { (i + 2) % 3 - 1, (i / 3 + 2) % 3 - 1, (i / 9 + 2) % 3 - 1 }; @@ -79,10 +80,17 @@ void ComputeDPD::local_interactions(const Particle * const xyzuvw, const float4 NVTX_RANGE("DPD/local", NVTX_C5); if (n > 0) + { forces_dpd_cuda_nohost((float*)xyzuvw, xyzouvwo, xyzo_half, (float *)a, n, cellsstart, cellscount, 1, XSIZE_SUBDOMAIN, YSIZE_SUBDOMAIN, ZSIZE_SUBDOMAIN, aij, gammadpd, - sigma, 1. / sqrt(dt), local_trunk.get_float(), stream); + sigma, 1. / sqrt(dt), current_lseed = local_trunk.get_float(), stream); + + if (sigma_xx) + compute_stress((float *)xyzuvw, n, cellsstart, cellscount, XSIZE_SUBDOMAIN, YSIZE_SUBDOMAIN, ZSIZE_SUBDOMAIN, + aij, gammadpd, sigmaf, current_lseed, + sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz, (float *)a, stream); + } } namespace BipsBatch @@ -102,9 +110,17 @@ namespace BipsBatch __constant__ BatchInfo batchinfos[26]; - __global__ void - interaction_kernel(const float aij, const float gamma, const float sigmaf, - const int ndstall, float * const adst, const int sizeadst) + struct StressInfo + { + float *sigma_xx, *sigma_xy, *sigma_xz, *sigma_yy, *sigma_yz, *sigma_zz; + }; + + __constant__ StressInfo stressinfo; + + + template < bool computestresses > __global__ + void interaction_kernel(const float aij, const float gamma, const float sigmaf, + const int ndstall, float * const adst, const int sizeadst) { #if !defined(__CUDA_ARCH__) #warning __CUDA_ARCH__ not defined! assuming 350 @@ -151,7 +167,8 @@ namespace BipsBatch const float vp = info.xdst[4 + dpid * 6]; const float wp = info.xdst[5 + dpid * 6]; - const int dstbase = 3 * info.scattered_entries[dpid]; + const int dstentry = info.scattered_entries[dpid]; + const int dstbase = 3 * dstentry; assert(dstbase < sizeadst * 3); uint scan1, scan2, ncandidates, spidbase; @@ -281,6 +298,16 @@ namespace BipsBatch xforce += strength * xr; yforce += strength * yr; zforce += strength * zr; + + if (computestresses) + { + atomicAdd(stressinfo.sigma_xx + dstentry, strength * xr * _xr); + atomicAdd(stressinfo.sigma_xy + dstentry, strength * xr * _yr); + atomicAdd(stressinfo.sigma_xz + dstentry, strength * xr * _zr); + atomicAdd(stressinfo.sigma_yy + dstentry, strength * yr * _yr); + atomicAdd(stressinfo.sigma_yz + dstentry, strength * yr * _zr); + atomicAdd(stressinfo.sigma_zz + dstentry, strength * zr * _zr); + } } atomicAdd(adst + dstbase + 0, xforce); @@ -295,12 +322,15 @@ namespace BipsBatch cudaEvent_t evhalodone; void interactions(const float aij, const float gamma, const float sigma, const float invsqrtdt, - const BatchInfo infos[20], cudaStream_t computestream, cudaStream_t uploadstream, float * const acc, const int n) + const BatchInfo infos[20], cudaStream_t computestream, cudaStream_t uploadstream, float * const acc, + float * const sigma_xx, float * const sigma_xy, float * const sigma_xz, float * const sigma_yy, + float * const sigma_yz, float * const sigma_zz, const int n) { if (firstcall) { CUDA_CHECK(cudaEventCreate(&evhalodone, cudaEventDisableTiming)); - CUDA_CHECK(cudaFuncSetCacheConfig(interaction_kernel, cudaFuncCachePreferL1)); + CUDA_CHECK(cudaFuncSetCacheConfig(interaction_kernel, cudaFuncCachePreferL1)); + CUDA_CHECK(cudaFuncSetCacheConfig(interaction_kernel, cudaFuncCachePreferL1)); firstcall = false; } @@ -316,12 +346,24 @@ namespace BipsBatch const int nthreads = 2 * hstart_padded[26]; + if (sigma_xx) + { + StressInfo strinfo = {sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz }; + + CUDA_CHECK(cudaMemcpyToSymbolAsync(stressinfo, &strinfo, sizeof(strinfo), 0, cudaMemcpyHostToDevice, uploadstream)); + } + CUDA_CHECK(cudaEventRecord(evhalodone, uploadstream)); CUDA_CHECK(cudaStreamWaitEvent(computestream, evhalodone, 0)); if (nthreads) - interaction_kernel<<< (nthreads + 127) / 128, 128, 0, computestream>>>(aij, gamma, sigma * invsqrtdt, nthreads, acc, n); + { + if (sigma_xx) + interaction_kernel<<< (nthreads + 127) / 128, 128, 0, computestream>>>(aij, gamma, sigma * invsqrtdt, nthreads, acc, n); + else + interaction_kernel<<< (nthreads + 127) / 128, 128, 0, computestream>>>(aij, gamma, sigma * invsqrtdt, nthreads, acc, n); + } CUDA_CHECK(cudaPeekAtLastError()); } @@ -346,7 +388,7 @@ void ComputeDPD::remote_interactions(const Particle * const p, const int n, Acce const int m2 = 0 == dz; BipsBatch::BatchInfo entry = { - (float *)sendhalos[i].dbuf.data, (float2 *)recvhalos[i].dbuf.data, interrank_trunks[i].get_float(), + (float *)sendhalos[i].dbuf.data, (float2 *)recvhalos[i].dbuf.data, current_rseeds[i] = interrank_trunks[i].get_float(), sendhalos[i].dbuf.size, recvhalos[i].dbuf.size, interrank_masks[i], recvhalos[i].dcellstarts.data, sendhalos[i].scattered_entries.data, dx, dy, dz, @@ -357,7 +399,8 @@ void ComputeDPD::remote_interactions(const Particle * const p, const int n, Acce infos[i] = entry; } - BipsBatch::interactions(aij, gammadpd, sigma, 1. / sqrt(dt), infos, stream, uploadstream, (float *)a, n); + BipsBatch::interactions(aij, gammadpd, sigma, 1. / sqrt(dt), infos, stream, uploadstream, (float *)a, + sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz, n); CUDA_CHECK(cudaPeekAtLastError()); } diff --git a/mpi-dpd/dpd.h b/mpi-dpd/dpd.h index 290fc7c49..4dfb9e4c3 100644 --- a/mpi-dpd/dpd.h +++ b/mpi-dpd/dpd.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-14. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -25,18 +22,43 @@ //see the vanilla version of this code for details about how this class operates class ComputeDPD : public SolventExchange -{ +{ Logistic::KISS local_trunk; Logistic::KISS interrank_trunks[26]; + float current_lseed, current_rseeds[26], + * sigma_xx, * sigma_xy, * sigma_xz, * sigma_yy, + * sigma_yz, * sigma_zz; + bool interrank_masks[26]; - + public: - + ComputeDPD(MPI_Comm cartcomm); + void set_stress_buffers(float * const stress_xx, float * const stress_xy, float * const stress_xz, float * const stress_yy, + float * const stress_yz, float * const stress_zz) + { + sigma_xx = stress_xx; + sigma_xy = stress_xy; + sigma_xz = stress_xz; + sigma_yy = stress_yy; + sigma_yz = stress_yz; + sigma_zz = stress_zz; + } + + void clr_stress_buffers() + { + sigma_xx = NULL; + sigma_xy = NULL; + sigma_xz = NULL; + sigma_yy = NULL; + sigma_yz = NULL; + sigma_zz = NULL; + } + void remote_interactions(const Particle * const p, const int n, Acceleration * const a, cudaStream_t stream, cudaStream_t uploadstream); - void local_interactions(const Particle * const xyzuvw, const float4 * const xyzouvwo, const ushort4 * const xyzo_half, const int n, Acceleration * const a, - const int * const cellsstart, const int * const cellscount, cudaStream_t stream); + void local_interactions(const Particle * const xyzuvw, const float4 * const xyzouvwo, const ushort4 * const xyzo_half, const int n, + Acceleration * const a, const int * const cellsstart, const int * const cellscount, cudaStream_t stream); }; diff --git a/mpi-dpd/fsi.cu b/mpi-dpd/fsi.cu index 4c8cdbec8..bef853471 100644 --- a/mpi-dpd/fsi.cu +++ b/mpi-dpd/fsi.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include <../dpd-rng.h> @@ -31,7 +28,7 @@ ComputeFSI::ComputeFSI(MPI_Comm comm) //TODO: use CUDA_CHECK(cudaEventCreateWithFlags(&evuploaded, cudaEventDisableTiming)); - KernelsFSI::Params params = {12.5 , gammadpd, sigmaf}; + KernelsFSI::Params params = {aij , gammadpd, sigmaf}; CUDA_CHECK(cudaMemcpyToSymbol(KernelsFSI::params, ¶ms, sizeof(params))); @@ -157,6 +154,7 @@ namespace KernelsFSI const float _zr = dst1.x - stmp1.x; const float rij2 = _xr * _xr + _yr * _yr + _zr * _zr; + assert(rij2 > 0); const float invrij = rsqrtf(rij2); @@ -271,188 +269,7 @@ void ComputeFSI::bulk(std::vector wsolutes, cudaStream_t stream) CUDA_CHECK(cudaPeekAtLastError()); } -namespace KernelsFSI -{ - __constant__ int packstarts_padded[27], packcount[26]; - __constant__ Particle * packstates[26]; - __constant__ Acceleration * packresults[26]; - - __global__ void interactions_halo(const int nparticles_padded, const int nsolvent, float * const accsolvent, const float seed) - { - assert(blockDim.x * gridDim.x >= nparticles_padded); - - const int laneid = threadIdx.x & 0x1f; - const int warpid = threadIdx.x >> 5; - const int localbase = 32 * (warpid + 4 * blockIdx.x); - const int pid = localbase + laneid; - - if (localbase >= nparticles_padded) - return; - - int nunpack; - float2 dst0, dst1, dst2; - float * dst = NULL; - - { - const uint key9 = 9 * (localbase >= packstarts_padded[9]) + 9 * (localbase >= packstarts_padded[18]); - const uint key3 = 3 * (localbase >= packstarts_padded[key9 + 3]) + 3 * (localbase >= packstarts_padded[key9 + 6]); - const uint key1 = (localbase >= packstarts_padded[key9 + key3 + 1]) + (localbase >= packstarts_padded[key9 + key3 + 2]); - const int code = key9 + key3 + key1; - assert(code >= 0 && code < 26); - assert(localbase >= packstarts_padded[code] && localbase < packstarts_padded[code + 1]); - - const int unpackbase = localbase - packstarts_padded[code]; - assert (unpackbase >= 0); - assert(unpackbase < packcount[code]); - - nunpack = min(32, packcount[code] - unpackbase); - - if (nunpack == 0) - return; - - read_AOS6f((float2 *)(packstates[code] + unpackbase), nunpack, dst0, dst1, dst2); - - dst = (float*)(packresults[code] + unpackbase); - } - - float xforce = 0, yforce = 0, zforce = 0; - - const int nzplanes = laneid < nunpack ? 3 : 0; - - for(int zplane = 0; zplane < nzplanes; ++zplane) - { - int scan1, scan2, ncandidates, spidbase; - int deltaspid1, deltaspid2; - - { - enum - { - XCELLS = XSIZE_SUBDOMAIN, - YCELLS = YSIZE_SUBDOMAIN, - ZCELLS = ZSIZE_SUBDOMAIN, - XOFFSET = XCELLS / 2, - YOFFSET = YCELLS / 2, - ZOFFSET = ZCELLS / 2 - }; - - const int NCELLS = XSIZE_SUBDOMAIN * YSIZE_SUBDOMAIN * ZSIZE_SUBDOMAIN; - const int xcenter = XOFFSET + (int)floorf(dst0.x); - const int xstart = max(0, xcenter - 1); - const int xcount = min(XCELLS, xcenter + 2) - xstart; - - if (xcenter - 1 >= XCELLS || xcenter + 2 <= 0) - continue; - - assert(xcount >= 0); - - const int ycenter = YOFFSET + (int)floorf(dst0.y); - - const int zcenter = ZOFFSET + (int)floorf(dst1.x); - const int zmy = zcenter - 1 + zplane; - const bool zvalid = zmy >= 0 && zmy < ZCELLS; - - int count0 = 0, count1 = 0, count2 = 0; - - if (zvalid && ycenter - 1 >= 0 && ycenter - 1 < YCELLS) - { - const int cid0 = xstart + XCELLS * (ycenter - 1 + YCELLS * zmy); - assert(cid0 >= 0 && cid0 + xcount <= NCELLS); - spidbase = tex1Dfetch(texCellsStart, cid0); - count0 = ((cid0 + xcount == NCELLS) ? nsolvent : tex1Dfetch(texCellsStart, cid0 + xcount)) - spidbase; - } - - if (zvalid && ycenter >= 0 && ycenter < YCELLS) - { - const int cid1 = xstart + XCELLS * (ycenter + YCELLS * zmy); - assert(cid1 >= 0 && cid1 + xcount <= NCELLS); - deltaspid1 = tex1Dfetch(texCellsStart, cid1); - count1 = ((cid1 + xcount == NCELLS) ? nsolvent : tex1Dfetch(texCellsStart, cid1 + xcount)) - deltaspid1; - } - - if (zvalid && ycenter + 1 >= 0 && ycenter + 1 < YCELLS) - { - const int cid2 = xstart + XCELLS * (ycenter + 1 + YCELLS * zmy); - deltaspid2 = tex1Dfetch(texCellsStart, cid2); - assert(cid2 >= 0 && cid2 + xcount <= NCELLS); - count2 = ((cid2 + xcount == NCELLS) ? nsolvent : tex1Dfetch(texCellsStart, cid2 + xcount)) - deltaspid2; - } - - scan1 = count0; - scan2 = count0 + count1; - ncandidates = scan2 + count2; - - deltaspid1 -= scan1; - deltaspid2 -= scan2; - } - - for(int i = 0; i < ncandidates; ++i) - { - const int m1 = (int)(i >= scan1); - const int m2 = (int)(i >= scan2); - const int spid = i + (m2 ? deltaspid2 : m1 ? deltaspid1 : spidbase); - - assert(spid >= 0 && spid < nsolvent); - - const int sentry = 3 * spid; - const float2 stmp0 = tex1Dfetch(texSolventParticles, sentry ); - const float2 stmp1 = tex1Dfetch(texSolventParticles, sentry + 1); - const float2 stmp2 = tex1Dfetch(texSolventParticles, sentry + 2); - - const float _xr = dst0.x - stmp0.x; - const float _yr = dst0.y - stmp0.y; - const float _zr = dst1.x - stmp1.x; - - const float rij2 = _xr * _xr + _yr * _yr + _zr * _zr; - - const float invrij = rsqrtf(rij2); - - const float rij = rij2 * invrij; - - if (rij2 >= 1) - continue; - - const float argwr = 1.f - rij; - const float wr = viscosity_function<-VISCOSITY_S_LEVEL>(argwr); - - const float xr = _xr * invrij; - const float yr = _yr * invrij; - const float zr = _zr * invrij; - - const float rdotv = - xr * (dst1.y - stmp1.y) + - yr * (dst2.x - stmp2.x) + - zr * (dst2.y - stmp2.y); - - const float myrandnr = Logistic::mean0var1(seed, pid, spid); - - const float strength = params.aij * argwr + (- params.gamma * wr * rdotv + params.sigmaf * myrandnr) * wr; - - const float xinteraction = strength * xr; - const float yinteraction = strength * yr; - const float zinteraction = strength * zr; - - xforce += xinteraction; - yforce += yinteraction; - zforce += zinteraction; - - assert(!isnan(xinteraction)); - assert(!isnan(yinteraction)); - assert(!isnan(zinteraction)); - assert(fabs(xinteraction) < 1e4); - assert(fabs(yinteraction) < 1e4); - assert(fabs(zinteraction) < 1e4); - - atomicAdd(accsolvent + sentry , -xinteraction); - atomicAdd(accsolvent + sentry + 1, -yinteraction); - atomicAdd(accsolvent + sentry + 2, -zinteraction); - } - } - - write_AOS3f(dst, nunpack, xforce, yforce, zforce); - } -} - -void ComputeFSI::halo(ParticlesWrap halos[26], cudaStream_t stream) +void ComputeFSI::halo(ParticlesWrap halowrap, cudaStream_t stream) { NVTX_RANGE("FSI/halo", NVTX_C7); @@ -460,50 +277,9 @@ void ComputeFSI::halo(ParticlesWrap halos[26], cudaStream_t stream) CUDA_CHECK(cudaPeekAtLastError()); - int nremote_padded = 0; - - { - int recvpackcount[26], recvpackstarts_padded[27]; - - for(int i = 0; i < 26; ++i) - recvpackcount[i] = halos[i].n; - - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsFSI::packcount, recvpackcount, - sizeof(recvpackcount), 0, cudaMemcpyHostToDevice, stream)); - - recvpackstarts_padded[0] = 0; - for(int i = 0, s = 0; i < 26; ++i) - recvpackstarts_padded[i + 1] = (s += 32 * ((halos[i].n + 31) / 32)); - - nremote_padded = recvpackstarts_padded[26]; - - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsFSI::packstarts_padded, recvpackstarts_padded, - sizeof(recvpackstarts_padded), 0, cudaMemcpyHostToDevice, stream)); - } - - { - const Particle * recvpackstates[26]; - - for(int i = 0; i < 26; ++i) - recvpackstates[i] = halos[i].p; - - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsFSI::packstates, recvpackstates, - sizeof(recvpackstates), 0, cudaMemcpyHostToDevice, stream)); - } - - { - Acceleration * packresults[26]; - - for(int i = 0; i < 26; ++i) - packresults[i] = halos[i].a; - - CUDA_CHECK(cudaMemcpyToSymbolAsync(KernelsFSI::packresults, packresults, - sizeof(packresults), 0, cudaMemcpyHostToDevice, stream)); - } - - if(nremote_padded) - KernelsFSI::interactions_halo<<< (nremote_padded + 127) / 128, 128, 0, stream>>> - (nremote_padded, wsolvent.n, (float *)wsolvent.a, local_trunk.get_float()); + if (halowrap.n) + KernelsFSI::interactions_3tpp<<< (3 * halowrap.n + 127) / 128, 128, 0, stream >>> + ((float2 *)halowrap.p, halowrap.n, wsolvent.n, (float *)halowrap.a, (float *)wsolvent.a, local_trunk.get_float()); CUDA_CHECK(cudaPeekAtLastError()); } diff --git a/mpi-dpd/fsi.h b/mpi-dpd/fsi.h index 130e753da..ba83f2e17 100644 --- a/mpi-dpd/fsi.h +++ b/mpi-dpd/fsi.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -36,5 +33,5 @@ class ComputeFSI : public SoluteExchange::Visitor void bulk(std::vector wsolutes, cudaStream_t stream); /*override of SoluteExchange::Visitor::halo*/ - void halo(ParticlesWrap solutes[26], cudaStream_t stream); + void halo(ParticlesWrap halowrap, cudaStream_t stream); }; diff --git a/mpi-dpd/io.cu b/mpi-dpd/io.cu index 11b963c29..e4c8b3288 100644 --- a/mpi-dpd/io.cu +++ b/mpi-dpd/io.cu @@ -6,9 +6,6 @@ * Major bug in H5 dump fixed by Panotelli on 2015-03-24. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -181,6 +178,54 @@ void ply_dump(MPI_Comm comm, MPI_Comm cartcomm, const char * filename, MPI_CHECK( MPI_File_close(&f)); } +void stress_dump(MPI_Comm cartcomm, const char * filename, const int nparticles, + const Particle * const particles, + const float * const stress_xx, const float * const stress_xy, const float * const stress_xz, + const float * const stress_yy, const float * const stress_yz, const float * const stress_zz) +{ + std::vector buf(nparticles * 12); + + int rank; + MPI_CHECK( MPI_Comm_rank(cartcomm, &rank) ); + + int dims[3], periods[3], coords[3]; + MPI_CHECK( MPI_Cart_get(cartcomm, 3, dims, periods, coords) ); + + int NALL = 0; + const int n = nparticles; + MPI_CHECK( MPI_Allreduce(&n, &NALL, 1, MPI_INT, MPI_SUM, cartcomm) ); + + MPI_File f; + MPI_CHECK( MPI_File_open(cartcomm, filename , MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &f) ); + + MPI_CHECK( MPI_File_set_size (f, sizeof(float) * 12 * NALL )); + + const int L[3] = { XSIZE_SUBDOMAIN, YSIZE_SUBDOMAIN, ZSIZE_SUBDOMAIN }; + + for(int i = 0; i < n; ++i) + { + const int base = 12 * i; + + for(int c = 0; c < 3; ++c) + buf[base + c] = particles[i].x[c] + L[c] / 2 + coords[c] * L[c]; + + for(int c = 0; c < 3; ++c) + buf[base + 3 + c] = particles[i].u[c]; + + buf[base + 6] = stress_xx[i] / 2; + buf[base + 7] = stress_xy[i] / 2; + buf[base + 8] = stress_xz[i] / 2; + buf[base + 9] = stress_yy[i] / 2; + buf[base + 10] = stress_yz[i] / 2; + buf[base + 11] = stress_zz[i] / 2; + } + + _write_bytes(buf.data(), sizeof(float) * 12 * n, f, cartcomm); + + MPI_CHECK( MPI_File_close(&f)); +} + + H5PartDump::H5PartDump(const string fname, MPI_Comm comm, MPI_Comm cartcomm): tstamp(0), disposed(false) { _initialize(fname, comm, cartcomm); diff --git a/mpi-dpd/io.h b/mpi-dpd/io.h index a970f5983..4d6fefd2a 100644 --- a/mpi-dpd/io.h +++ b/mpi-dpd/io.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-01-30. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -23,6 +20,11 @@ void ply_dump(MPI_Comm comm, MPI_Comm cartcomm, const char * filename, int (*mesh_indices)[3], const int ninstances, const int ntriangles_per_instance, Particle * _particles, int nvertices_per_instance, bool append); +void stress_dump(MPI_Comm comm, const char * filename, const int nparticles, + const Particle * const particles, + const float * const stress_xx, const float * const stress_xy, const float * const stress_xz, + const float * const stress_yy, const float * const stress_yz, const float * const stress_zz); + class H5PartDump { float origin[3]; diff --git a/mpi-dpd/main.cu b/mpi-dpd/main.cu index 3e5a7e5f4..57026523d 100644 --- a/mpi-dpd/main.cu +++ b/mpi-dpd/main.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-14. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -23,8 +20,9 @@ #include "simulation.h" bool currently_profiling = false; -float tend; -bool walls, pushtheflow, doublepoiseuille, rbcs, ctcs, xyz_dumps, hdf5field_dumps, hdf5part_dumps, is_mps_enabled, adjust_message_sizes, contactforces; +float tend, couette; +bool walls, pushtheflow, doublepoiseuille, rbcs, ctcs, xyz_dumps, hdf5field_dumps, + hdf5part_dumps, is_mps_enabled, adjust_message_sizes, contactforces, stress; int steps_per_report, steps_per_dump, wall_creation_stepid, nvtxstart, nvtxstop; LocalComm localcomm; @@ -84,6 +82,8 @@ int main(int argc, char ** argv) nvtxstop = argp("-nvtxstop").asInt(10500); adjust_message_sizes = argp("-adjust_message_sizes").asBool(false); contactforces = argp("-contactforces").asBool(false); + stress = argp("-stress").asBool(false); + couette = argp("-couette").asDouble(0); #ifndef _NO_DUMPS_ const bool mpi_thread_safe = argp("-mpi_thread_safe").asBool(true); diff --git a/mpi-dpd/minmax.cu b/mpi-dpd/minmax.cu index 57aba10b1..58520a9df 100644 --- a/mpi-dpd/minmax.cu +++ b/mpi-dpd/minmax.cu @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-23. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "minmax.h" @@ -258,4 +255,4 @@ void minmax(const Particle * const rbc, int size, int n, float3 *minrbc, float3 minmaxmba<<>>(rbc, minrbc, maxrbc, size, ptoblockds); } -} \ No newline at end of file +} diff --git a/mpi-dpd/minmax.h b/mpi-dpd/minmax.h index a5317be31..7eca7f67d 100644 --- a/mpi-dpd/minmax.h +++ b/mpi-dpd/minmax.h @@ -5,9 +5,6 @@ * Created and authored by Massimo Bernaschi on 2015-03-24. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/redistancing.cu b/mpi-dpd/redistancing.cu index 8747df39b..2e2878576 100644 --- a/mpi-dpd/redistancing.cu +++ b/mpi-dpd/redistancing.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-17. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "redistancing.h" diff --git a/mpi-dpd/redistancing.h b/mpi-dpd/redistancing.h index ede725383..0a0057d84 100644 --- a/mpi-dpd/redistancing.h +++ b/mpi-dpd/redistancing.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-17. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/redistribute-particles.cu b/mpi-dpd/redistribute-particles.cu index e31ec5d3b..624e84629 100644 --- a/mpi-dpd/redistribute-particles.cu +++ b/mpi-dpd/redistribute-particles.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-02-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -744,7 +741,7 @@ void RedistributeParticles::bulk(const int nparticles, int * const cellstarts, i subindices.resize(nparticles); if (nparticles) - subindex_local<<< (nparticles + 127) / 128, 128, 0, mystream>>> + subindex_local<<< (nparticles + 127) / 128, 128, 0, mystream>>> (nparticles, RedistributeParticlesKernels::texparticledata, cellcounts, subindices.data); /* #ifndef NDEBUG diff --git a/mpi-dpd/redistribute-particles.h b/mpi-dpd/redistribute-particles.h index f7e281b41..dd9e5c86d 100644 --- a/mpi-dpd/redistribute-particles.h +++ b/mpi-dpd/redistribute-particles.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-14. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/redistribute-rbcs.cu b/mpi-dpd/redistribute-rbcs.cu index cb3998844..b812e47e8 100644 --- a/mpi-dpd/redistribute-rbcs.cu +++ b/mpi-dpd/redistribute-rbcs.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-01. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/mpi-dpd/redistribute-rbcs.h b/mpi-dpd/redistribute-rbcs.h index 0917578d4..0f98d3f45 100644 --- a/mpi-dpd/redistribute-rbcs.h +++ b/mpi-dpd/redistribute-rbcs.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-01. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/scan.cu b/mpi-dpd/scan.cu index ef5245f4e..71ea7094a 100644 --- a/mpi-dpd/scan.cu +++ b/mpi-dpd/scan.cu @@ -5,9 +5,6 @@ * Created and authored by Mauro Bisson on 2015-07-28. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ template diff --git a/mpi-dpd/simulation.cu b/mpi-dpd/simulation.cu index 4ecf0f4a5..690094d67 100644 --- a/mpi-dpd/simulation.cu +++ b/mpi-dpd/simulation.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-24. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -257,7 +254,7 @@ void Simulation::_create_walls(const bool verbose, bool & termination_request) int nsurvived = 0; ExpectedMessageSizes new_sizes; - wall = new ComputeWall(cartcomm, particles->xyzuvw.data, particles->size, nsurvived, new_sizes, verbose); + wall = new ComputeWall(cartcomm, particles->xyzuvw.data, particles->size, nsurvived, new_sizes, couette); //adjust the message sizes if we're pushing the flow in x { @@ -374,9 +371,6 @@ void Simulation::_forces() CUDA_CHECK(cudaPeekAtLastError()); - if (contactforces) - contact.build_cells(wsolutes, mainstream); - dpd.local_interactions(particles->xyzuvw.data, xyzouvwo.data, xyzo_half.data, particles->size, particles->axayaz.data, cells.start, cells.count, mainstream); @@ -400,17 +394,17 @@ void Simulation::_forces() dpd.recv(mainstream, uploadstream); - solutex.recv_p(uploadstream); + solutex.recv_p(uploadstream, mainstream); - solutex.halo(uploadstream, mainstream); + if (contactforces) + contact.attach_bulk(wsolutes); + + solutex.halo(uploadstream, mainstream, downloadstream); dpd.remote_interactions(particles->xyzuvw.data, particles->size, particles->axayaz.data, mainstream, uploadstream); fsi.bulk(wsolutes, mainstream); - if (contactforces) - contact.bulk(wsolutes, mainstream); - CUDA_CHECK(cudaPeekAtLastError()); if (rbcscoll) @@ -441,6 +435,15 @@ void Simulation::_datadump(const int idtimestep) int n = particles->size; + if (stress) + { + for(int c = 0; c < 6; ++c) + stresses_datadump[c].resize(n); + + for(int c = 0; c < 6; ++c) + CUDA_CHECK(cudaMemcpyAsync(stresses_datadump[c].data, stresses[c].data, sizeof(float) * n, cudaMemcpyDeviceToHost,0)); + } + if (rbcscoll) n += rbcscoll->pcount(); @@ -568,6 +571,43 @@ void Simulation::_datadump_async() xyz_dump(myactivecomm, mycartcomm, "xyz/particles->xyz", "all-particles", p, n, datadump_idtimestep > 0); } + if (stress) + { + char filename[1024]; + sprintf(filename, "stress/stresses-%05d.data", iddatadump); + + if(rank == 0 && iddatadump == 0) + mkdir("stress", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + + const int nsolvent = stresses_datadump[0].size; + + stress_dump(mycartcomm, filename, stresses_datadump[0].size, p, + stresses_datadump[0].data, stresses_datadump[1].data, stresses_datadump[2].data, + stresses_datadump[3].data, stresses_datadump[4].data, stresses_datadump[5].data); + + //lets quickly compute the average stress in the system + const int v1[6] = {0, 0, 0, 1, 1, 2}; + const int v2[6] = {0, 1, 2, 1, 2, 2}; + + float avgstress[6] = {0, 0, 0, 0, 0, 0}; + for(int c = 0; c < 6; ++c) + for(int i = 0; i < nsolvent; ++i) + avgstress[c] += stresses_datadump[c].data[i] + p[i].u[v1[c]] * p[i].u[v2[c]]; + + int ntotsolvent; + MPI_CHECK( MPI_Reduce(&nsolvent, &ntotsolvent, 1, MPI_INT, MPI_SUM, 0, myactivecomm)); + + float totavgstress[6]; + MPI_CHECK( MPI_Reduce(avgstress, totavgstress, 6, MPI_FLOAT, MPI_SUM, 0, myactivecomm)); + + for(int c = 0; c < 6; ++c) + totavgstress[c] /= ntotsolvent; + + if (rank == 0) + printf("average stress: sxx:%.3e sxy:%.3e sxz:%.3e syy:%.3e syz:%.3e szz:%.3e\n", + totavgstress[0], totavgstress[1], totavgstress[2], totavgstress[3], totavgstress[4], totavgstress[5]); + } + if (hdf5part_dumps) { NVTX_RANGE("h5part dump", NVTX_C3); @@ -679,7 +719,6 @@ Simulation::Simulation(MPI_Comm cartcomm, MPI_Comm activecomm, bool (*check_term if (contactforces) solutex.attach_halocomputation(contact); - //localcomm.initialize(activecomm); int dims[3], periods[3], coords[3]; MPI_CHECK( MPI_Cart_get(cartcomm, 3, dims, periods, coords) ); @@ -792,9 +831,6 @@ void Simulation::_lockstep() dpd.local_interactions(particles->xyzuvw.data, xyzouvwo.data, xyzo_half.data, particles->size, particles->axayaz.data, cells.start, cells.count, mainstream); - if (contactforces) - contact.build_cells(wsolutes, mainstream); - solutex.post_p(mainstream, downloadstream); dpd.post(particles->xyzuvw.data, particles->size, mainstream, downloadstream); @@ -809,17 +845,17 @@ void Simulation::_lockstep() dpd.recv(mainstream, uploadstream); - solutex.recv_p(uploadstream); + solutex.recv_p(uploadstream, mainstream); - solutex.halo(uploadstream, mainstream); + if (contactforces) + contact.attach_bulk(wsolutes); + + solutex.halo(uploadstream, mainstream, downloadstream); dpd.remote_interactions(particles->xyzuvw.data, particles->size, particles->axayaz.data, mainstream, uploadstream); fsi.bulk(wsolutes, mainstream); - if (contactforces) - contact.bulk(wsolutes, mainstream); - CUDA_CHECK(cudaPeekAtLastError()); if (rbcscoll) @@ -946,7 +982,6 @@ void Simulation::run() int it; - for(it = 0; it < nsteps; ++it) { const bool verbose = it > 0 && rank == 0; @@ -1025,8 +1060,27 @@ void Simulation::run() printf("the simulation begins now and it consists of %.3e steps\n", (double)(nsteps - it)); } + if(stress && it % steps_per_dump == 0) + { + for(int c = 0; c < 6; ++c) + stresses[c].resize(particles->size); + + dpd.set_stress_buffers(stresses[0].data, stresses[1].data, stresses[2].data, stresses[3].data, stresses[4].data, stresses[5].data); + + if (wall) + wall->set_stress_buffers(stresses[0].data, stresses[1].data, stresses[2].data, stresses[3].data, stresses[4].data, stresses[5].data); + } + _forces(); + if (stress && it % steps_per_dump == 0 ) + { + dpd.clr_stress_buffers(); + + if (wall) + wall->clr_stress_buffers(); + } + #ifndef _NO_DUMPS_ if (it % steps_per_dump == 0) _datadump(it); diff --git a/mpi-dpd/simulation.h b/mpi-dpd/simulation.h index 5a097d522..f3c6e718e 100644 --- a/mpi-dpd/simulation.h +++ b/mpi-dpd/simulation.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-03-24. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -41,6 +38,7 @@ class Simulation ParticleArray * particles, * newparticles; SimpleDeviceBuffer xyzouvwo; SimpleDeviceBuffer xyzo_half; + SimpleDeviceBuffer stresses[6]; CellLists cells; CollectionRBC * rbcscoll; @@ -93,6 +91,7 @@ class Simulation PinnedHostBuffer particles_datadump; PinnedHostBuffer accelerations_datadump; + PinnedHostBuffer stresses_datadump[6]; cudaEvent_t evdownloaded; diff --git a/mpi-dpd/solute-exchange.cu b/mpi-dpd/solute-exchange.cu index e2d63c03d..c57c0bfd6 100644 --- a/mpi-dpd/solute-exchange.cu +++ b/mpi-dpd/solute-exchange.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ //#define _DUMBCRAY_ @@ -59,7 +56,8 @@ iterationcount(-1), packstotalstart(27), host_packstotalstart(27), host_packstot _adjust_packbuffers(); CUDA_CHECK(cudaEventCreateWithFlags(&evPpacked, cudaEventDisableTiming | cudaEventBlockingSync)); - CUDA_CHECK(cudaEventCreateWithFlags(&evAcomputed, cudaEventDisableTiming | cudaEventBlockingSync)); + CUDA_CHECK(cudaEventCreateWithFlags(&evAcomputed, cudaEventDisableTiming)); + CUDA_CHECK(cudaEventCreateWithFlags(&evAdownloaded, cudaEventDisableTiming | cudaEventBlockingSync)); CUDA_CHECK(cudaPeekAtLastError()); } @@ -351,7 +349,7 @@ void SoluteExchange::pack_p(cudaStream_t stream) if (wsolutes.size() == 0) return; - NVTX_RANGE("FSI/pack", NVTX_C4); + NVTX_RANGE("SOLUTEX/pack", NVTX_C4); ++iterationcount; @@ -371,7 +369,7 @@ void SoluteExchange::post_p(cudaStream_t stream, cudaStream_t downloadstream) //consolidate the packing { - NVTX_RANGE("FSI/consolidate", NVTX_C5); + NVTX_RANGE("SOLUTEX/consolidate", NVTX_C5); CUDA_CHECK(cudaEventSynchronize(evPpacked)); @@ -446,7 +444,7 @@ void SoluteExchange::post_p(cudaStream_t stream, cudaStream_t downloadstream) //post the sending of the packs { - NVTX_RANGE("FSI/send", NVTX_C6); + NVTX_RANGE("SOLUTEX/send", NVTX_C6); reqsendC.resize(26); @@ -485,13 +483,13 @@ void SoluteExchange::post_p(cudaStream_t stream, cudaStream_t downloadstream) } } -void SoluteExchange::recv_p(cudaStream_t uploadstream) +void SoluteExchange::recv_p(cudaStream_t uploadstream, cudaStream_t computestream) { if (wsolutes.size() == 0) return; - NVTX_RANGE("FSI/recv-p", NVTX_C7); - + NVTX_RANGE("SOLUTEX/recv-p", NVTX_C7); + _wait(reqrecvC); _wait(reqrecvP); @@ -504,7 +502,6 @@ void SoluteExchange::recv_p(cudaStream_t uploadstream) remote[i].preserve_resize(count); #ifndef NDEBUG - CUDA_CHECK(cudaMemsetAsync(remote[i].dstate.data, 0xff, sizeof(Particle) * remote[i].dstate.capacity, uploadstream)); CUDA_CHECK(cudaMemsetAsync(remote[i].result.data, 0xff, sizeof(Acceleration) * remote[i].result.capacity, uploadstream)); #endif @@ -524,39 +521,71 @@ void SoluteExchange::recv_p(cudaStream_t uploadstream) } _postrecvC(); - - for(int i = 0; i < 26; ++i) - CUDA_CHECK(cudaMemcpyAsync(remote[i].dstate.data, remote[i].hstate.data, sizeof(Particle) * remote[i].hstate.size, - cudaMemcpyHostToDevice, uploadstream)); + + //collate halos + { + int c = 0; + for(int i = 0; i < 26; ++i) + c += remote[i].hstate.size; + +#ifndef NDEBUG + CUDA_CHECK(cudaMemsetAsync(allremotehalosacc.data, 0xff, sizeof(Acceleration) * allremotehalosacc.capacity, computestream)); + CUDA_CHECK(cudaMemsetAsync(allremotehalos.data, 0xff, sizeof(Particle) * allremotehalos.capacity, uploadstream)); +#endif + + allremotehalos.resize(c); + allremotehalosacc.resize(c); + + CUDA_CHECK(cudaMemsetAsync(allremotehalosacc.data, 0, sizeof(Acceleration) * allremotehalosacc.size, computestream)); + + c = 0; + for(int i = 0; i < 26; ++i) + { + CUDA_CHECK(cudaMemcpyAsync(allremotehalos.data + c, remote[i].hstate.data, sizeof(Particle) * remote[i].hstate.size, cudaMemcpyHostToDevice, uploadstream)); + + c += remote[i].hstate.size; + } + } } -void SoluteExchange::halo(cudaStream_t uploadstream, cudaStream_t stream) +void SoluteExchange::halo(cudaStream_t uploadstream, cudaStream_t computestream, cudaStream_t downloadstream) { - NVTX_RANGE("FSI/halo", NVTX_C7); + NVTX_RANGE("SOLUTEX/halo", NVTX_C7); if (wsolutes.size() == 0) return; - + if (iterationcount) _wait(reqsendA); - - ParticlesWrap halos[26]; - - for(int i = 0; i < 26; ++i) - halos[i] = ParticlesWrap(remote[i].dstate.data, remote[i].dstate.size, remote[i].result.devptr); - + + ParticlesWrap halowrap(allremotehalos.data, allremotehalos.size, allremotehalosacc.data); + CUDA_CHECK(cudaStreamSynchronize(uploadstream)); - + for(int i = 0; i < visitors.size(); ++i) - visitors[i]->halo(halos, stream); - + visitors[i]->halo(halowrap, computestream); + CUDA_CHECK(cudaPeekAtLastError()); - - CUDA_CHECK(cudaEventRecord(evAcomputed, stream)); - + + CUDA_CHECK(cudaEventRecord(evAcomputed, computestream)); + + CUDA_CHECK(cudaStreamWaitEvent(downloadstream, evAcomputed, 0)); + + //split back halos + { + int c = 0; + for(int i = 0; i < 26; ++i) + { + CUDA_CHECK(cudaMemcpyAsync(remote[i].result.data, allremotehalosacc.data + c, sizeof(Acceleration) * remote[i].hstate.size, cudaMemcpyDeviceToHost, downloadstream)); + c += remote[i].hstate.size; + } + } + + CUDA_CHECK(cudaEventRecord(evAdownloaded, downloadstream)); + for(int i = 0; i < 26; ++i) local[i].update(); - + #ifndef _DUMBCRAY_ _postrecvP(); #endif @@ -567,9 +596,9 @@ void SoluteExchange::post_a() if (wsolutes.size() == 0) return; - NVTX_RANGE("FSI/send-a", NVTX_C1); + NVTX_RANGE("SOLUTEX/send-a", NVTX_C1); - CUDA_CHECK(cudaEventSynchronize(evAcomputed)); + CUDA_CHECK(cudaEventSynchronize(evAdownloaded)); reqsendA.resize(26); for(int i = 0; i < 26; ++i) @@ -632,7 +661,7 @@ void SoluteExchange::recv_a(cudaStream_t stream) if (wsolutes.size() == 0) return; - NVTX_RANGE("FSI/merge", NVTX_C2); + NVTX_RANGE("SOLUTEX/merge", NVTX_C2); { float * recvbags[26]; @@ -667,4 +696,5 @@ SoluteExchange::~SoluteExchange() CUDA_CHECK(cudaEventDestroy(evPpacked)); CUDA_CHECK(cudaEventDestroy(evAcomputed)); + CUDA_CHECK(cudaEventDestroy(evAdownloaded)); } diff --git a/mpi-dpd/solute-exchange.h b/mpi-dpd/solute-exchange.h index bf315c628..a1b5e50b4 100644 --- a/mpi-dpd/solute-exchange.h +++ b/mpi-dpd/solute-exchange.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-12-02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -19,11 +16,11 @@ class SoluteExchange { enum { TAGBASE_C = 113, TAGBASE_P = 365, TAGBASE_A = 668, TAGBASE_P2 = 1055, TAGBASE_A2 = 1501 }; - + public: - - struct Visitor { virtual void halo(ParticlesWrap solutehalos[26], cudaStream_t stream) = 0; }; - + + struct Visitor { virtual void halo(ParticlesWrap allhalos, cudaStream_t stream) = 0; }; + protected: MPI_Comm cartcomm; @@ -34,20 +31,20 @@ class SoluteExchange dims[3], periods[3], coords[3], myrank, recv_tags[26], recv_counts[26], send_counts[26]; - cudaEvent_t evPpacked, evAcomputed; + cudaEvent_t evPpacked, evAcomputed, evAdownloaded; SimpleDeviceBuffer packscount, packsstart, packsoffset, packstotalstart; PinnedHostBuffer host_packstotalstart, host_packstotalcount; - + SimpleDeviceBuffer packbuf; PinnedHostBuffer host_packbuf; - + std::vector wsolutes; - + std::vector reqsendC, reqrecvC, reqsendP, reqrecvP, reqsendA, reqrecvA; std::vector visitors; - + class TimeSeriesWindow { static const int N = 200; @@ -77,14 +74,12 @@ class SoluteExchange public: - SimpleDeviceBuffer dstate; PinnedHostBuffer hstate; PinnedHostBuffer result; std::vector pmessage; void preserve_resize(int n) { - dstate.resize(n); hstate.preserve_resize(n); result.resize(n); history.update(n); @@ -92,10 +87,13 @@ class SoluteExchange int expected() const { return (int)ceil(history.max() * 1.1); } - int capacity() const { assert(hstate.capacity == dstate.capacity); return dstate.capacity; } + int capacity() const { return hstate.capacity; } } remote[26]; + SimpleDeviceBuffer allremotehalos; + SimpleDeviceBuffer allremotehalosacc; + class LocalHalo { TimeSeriesWindow history; @@ -202,7 +200,7 @@ class SoluteExchange void _pack_attempt(cudaStream_t stream); public: - + SoluteExchange(MPI_Comm cartcomm); void bind_solutes(std::vector wsolutes) { this->wsolutes = wsolutes; } @@ -213,9 +211,9 @@ class SoluteExchange void post_p(cudaStream_t stream, cudaStream_t downloadstream); - void recv_p(cudaStream_t uploadstream); - - void halo(cudaStream_t uploadstream, cudaStream_t stream); + void recv_p(cudaStream_t uploadstream, cudaStream_t computestream); + + void halo(cudaStream_t uploadstream, cudaStream_t computestream, cudaStream_t downloadstream); void post_a(); diff --git a/mpi-dpd/solvent-exchange.cu b/mpi-dpd/solvent-exchange.cu index 3fdd19b5b..246901ada 100644 --- a/mpi-dpd/solvent-exchange.cu +++ b/mpi-dpd/solvent-exchange.cu @@ -6,9 +6,6 @@ * Major editing from Massimo-Bernaschi on 2015-03-20. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/mpi-dpd/solvent-exchange.h b/mpi-dpd/solvent-exchange.h index d4ab34a6d..e74620e5b 100644 --- a/mpi-dpd/solvent-exchange.h +++ b/mpi-dpd/solvent-exchange.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-18. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/mpi-dpd/wall.cu b/mpi-dpd/wall.cu index aea985207..882c992a8 100644 --- a/mpi-dpd/wall.cu +++ b/mpi-dpd/wall.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-19. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include @@ -57,8 +54,10 @@ namespace SolidWallsKernel texture texWallParticles; texture texWallCellStart, texWallCellCount; + template __global__ void interactions_3tpp(const float2 * const particles, const int np, const int nsolid, - float * const acc, const float seed, const float sigmaf); + float * const acc, const float seed, const float sigmaf, const float xvel, const float y0); + void setup() { texSDF.normalized = 0; @@ -83,7 +82,8 @@ namespace SolidWallsKernel texWallCellCount.mipmapFilterMode = cudaFilterModePoint; texWallCellCount.normalized = 0; - CUDA_CHECK(cudaFuncSetCacheConfig(interactions_3tpp, cudaFuncCachePreferL1)); + CUDA_CHECK(cudaFuncSetCacheConfig(interactions_3tpp, cudaFuncCachePreferL1)); + CUDA_CHECK(cudaFuncSetCacheConfig(interactions_3tpp, cudaFuncCachePreferL1)); } __device__ float sdf(float x, float y, float z) @@ -201,8 +201,6 @@ namespace SolidWallsKernel return make_float3(xmygrad, ymygrad, zmygrad); } - - __global__ void fill_keys(const Particle * const particles, const int n, int * const key) { assert(blockDim.x * gridDim.x >= n); @@ -377,8 +375,18 @@ namespace SolidWallsKernel } } + struct StressInfo + { + float *sigma_xx, *sigma_xy, *sigma_xz, *sigma_yy, *sigma_yz, *sigma_zz; + }; + + __constant__ StressInfo stressinfo; + + + template __global__ __launch_bounds__(128, 16) void interactions_3tpp(const float2 * const particles, const int np, const int nsolid, - float * const acc, const float seed, const float sigmaf) + float * const acc, const float seed, const float sigmaf, + const float xvelocity_wall, const float z0) { assert(blockDim.x * gridDim.x >= np * 3); @@ -487,9 +495,11 @@ namespace SolidWallsKernel const float xr = _xr * invrij; const float yr = _yr * invrij; const float zr = _zr * invrij; - + + const float xvel = zq > z0 ? xvelocity_wall : 0; + const float rdotv = - xr * (dst1.y - 0) + + xr * (dst1.y - xvel) + yr * (dst2.x - 0) + zr * (dst2.y - 0); @@ -500,6 +510,16 @@ namespace SolidWallsKernel xforce += strength * xr; yforce += strength * yr; zforce += strength * zr; + + if (computestresses) + { + atomicAdd(stressinfo.sigma_xx + pid, strength * xr * _xr); + atomicAdd(stressinfo.sigma_xy + pid, strength * xr * _yr); + atomicAdd(stressinfo.sigma_xz + pid, strength * xr * _zr); + atomicAdd(stressinfo.sigma_yy + pid, strength * yr * _yr); + atomicAdd(stressinfo.sigma_yz + pid, strength * yr * _zr); + atomicAdd(stressinfo.sigma_zz + pid, strength * zr * _zr); + } } atomicAdd(acc + 3 * pid + 0, xforce); @@ -756,8 +776,10 @@ struct FieldSampler }; ComputeWall::ComputeWall(MPI_Comm cartcomm, Particle* const p, const int n, int& nsurvived, - ExpectedMessageSizes& new_sizes, const bool verbose): - cartcomm(cartcomm), arrSDF(NULL), solid4(NULL), solid_size(0), + ExpectedMessageSizes& new_sizes, const float xvelocity): + sigma_xx(NULL), sigma_xy(NULL), sigma_xz(NULL), sigma_yy(NULL), sigma_yz(NULL), sigma_zz(NULL), + cartcomm(cartcomm), arrSDF(NULL), solid4(NULL), solid_size(0), xvelocity(xvelocity), + cells(XSIZE_SUBDOMAIN + 2 * XMARGIN_WALL, YSIZE_SUBDOMAIN + 2 * YMARGIN_WALL, ZSIZE_SUBDOMAIN + 2 * ZMARGIN_WALL) { MPI_CHECK( MPI_Comm_rank(cartcomm, &myrank)); @@ -766,6 +788,7 @@ ComputeWall::ComputeWall(MPI_Comm cartcomm, Particle* const p, const int n, int& float * field = new float[ XTEXTURESIZE * YTEXTURESIZE * ZTEXTURESIZE]; + static const bool verbose = false; FieldSampler sampler("sdf.dat", cartcomm, verbose); const int L[3] = { XSIZE_SUBDOMAIN, YSIZE_SUBDOMAIN, ZSIZE_SUBDOMAIN }; @@ -1104,7 +1127,6 @@ void ComputeWall::interactions(const Particle * const p, const int n, Accelerati const int * const cellsstart, const int * const cellscount, cudaStream_t stream) { NVTX_RANGE("WALL/interactions", NVTX_C3); - //cellsstart and cellscount IGNORED for now if (n > 0 && solid_size > 0) { @@ -1121,8 +1143,21 @@ void ComputeWall::interactions(const Particle * const p, const int n, Accelerati &SolidWallsKernel::texWallCellCount.channelDesc, sizeof(int) * cells.ncells)); assert(textureoffset == 0); - SolidWallsKernel::interactions_3tpp<<< (3 * n + 127) / 128, 128, 0, stream>>> - ((float2 *)p, n, solid_size, (float *)acc, trunk.get_float(), sigmaf); + const float z0 = (dims[2] - 1 - 2 * coords[2]) * ZSIZE_SUBDOMAIN / 2; + + if (sigma_xx) + { + SolidWallsKernel::StressInfo strinfo = { sigma_xx, sigma_xy, sigma_xz, sigma_yy, sigma_yz, sigma_zz }; + + CUDA_CHECK(cudaMemcpyToSymbolAsync(SolidWallsKernel::stressinfo, &strinfo, sizeof(strinfo), 0, cudaMemcpyHostToDevice, stream)); + + SolidWallsKernel::interactions_3tpp<<< (3 * n + 127) / 128, 128, 0, stream>>> + ((float2 *)p, n, solid_size, (float *)acc, trunk.get_float(), sigmaf, xvelocity, z0); + } + else + SolidWallsKernel::interactions_3tpp<<< (3 * n + 127) / 128, 128, 0, stream>>> + ((float2 *)p, n, solid_size, (float *)acc, trunk.get_float(), sigmaf, xvelocity, z0); + CUDA_CHECK(cudaUnbindTexture(SolidWallsKernel::texWallParticles)); CUDA_CHECK(cudaUnbindTexture(SolidWallsKernel::texWallCellStart)); diff --git a/mpi-dpd/wall.h b/mpi-dpd/wall.h index 466c54977..8e9db4228 100644 --- a/mpi-dpd/wall.h +++ b/mpi-dpd/wall.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-19. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once @@ -31,6 +28,8 @@ class ComputeWall int solid_size; float4 * solid4; + float * sigma_xx, * sigma_xy, * sigma_xz, * sigma_yy, + * sigma_yz, * sigma_zz, xvelocity; cudaArray * arrSDF; @@ -38,12 +37,33 @@ class ComputeWall public: - ComputeWall(MPI_Comm cartcomm, Particle* const p, const int n, int& nsurvived, ExpectedMessageSizes& new_sizes, const bool verbose); + ComputeWall(MPI_Comm cartcomm, Particle* const p, const int n, int& nsurvived, ExpectedMessageSizes& new_sizes, const float xvelocity); ~ComputeWall(); void bounce(Particle * const p, const int n, cudaStream_t stream); + void set_stress_buffers(float * const stress_xx, float * const stress_xy, float * const stress_xz, float * const stress_yy, + float * const stress_yz, float * const stress_zz) + { + sigma_xx = stress_xx; + sigma_xy = stress_xy; + sigma_xz = stress_xz; + sigma_yy = stress_yy; + sigma_yz = stress_yz; + sigma_zz = stress_zz; + } + + void clr_stress_buffers() + { + sigma_xx = NULL; + sigma_xy = NULL; + sigma_xz = NULL; + sigma_yy = NULL; + sigma_yz = NULL; + sigma_zz = NULL; + } + void interactions(const Particle * const p, const int n, Acceleration * const acc, const int * const cellsstart, const int * const cellscount, cudaStream_t stream); }; diff --git a/postprocessing/argument-parser.h b/postprocessing/argument-parser.h new file mode 100644 index 000000000..d01b9f87d --- /dev/null +++ b/postprocessing/argument-parser.h @@ -0,0 +1,243 @@ +/* + * ArgumentParser.h + * Cubism + * + *This argument parser assumes that all arguments are optional ie, each of the argument names is preceded by a '-' + *all arguments are however NOT optional to avoid a mess with default values and returned values when not found! + * + *More converter could be required: + *add as needed + *TypeName as{TypeName}() in Value + * + * Created by Christian Conti on 6/7/10. That is a long time ago. + * Modified by Diego Rossinelli several times after his dreadlocks hair cut. + * Copyright 2010 ETH Zurich. All rights reserved. + * + */ + +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +class Value +{ +private: + string content; + +public: + +Value() : content("") {} + +Value(string content_) : content(content_) { /*printf("%s\n",content.c_str());*/ } + + double asDouble(double def=0) const + { + if (content == "") return def; + return (double) atof(content.c_str()); + } + + int asInt(int def=0) const + { + if (content == "") return def; + return atoi(content.c_str()); + } + + bool asBool(bool def=false) const + { + if (content == "") return def; + if (content == "0") return false; + if (content == "false") return false; + + return true; + } + + string asString(string def="") const + { + if (content == "") return def; + + return content; + } + + vector asVecFloat(const int musthave_size = -1) const + { + //printf("mycontent is %s\n", content.c_str()); + std::stringstream ss(content); + //assert(ss.good()); + vector retval; + double e; + + while (ss >> e) + { + retval.push_back(e); + // printf("reading %f\n", e); + if (ss.peek() == ',') + ss.ignore(); + } + + if (musthave_size > 0) + assert(musthave_size == (int)retval.size()); + + return retval; + } +}; + +class ArgumentParser +{ +private: + + map mapArguments; + + const int iArgC; + const char** vArgV; + bool bStrictMode, bVerbose; + + const char delimiter; +public: + + Value operator()(const string arg) + { + map::const_iterator it = mapArguments.find(arg); + + if (bStrictMode) + { + if (it == mapArguments.end()) + { + printf("Runtime option NOT SPECIFIED! ABORTING! name: %s\n",arg.data()); + abort(); + } + } + + if (bVerbose) + printf("%s is %s\n", arg.data(), mapArguments[arg].asString().data()); + + if (it != mapArguments.end()) + return mapArguments[arg]; + else + return Value(); + } + + bool check(const string arg) const + { + return mapArguments.find(arg) != mapArguments.end(); + } + +ArgumentParser(const int argc, const char ** argv, bool bVerbose = false, const char delimiter = '=') : + mapArguments(), iArgC(argc), vArgV(argv), bStrictMode(false), bVerbose(bVerbose), delimiter(delimiter) + { + for (int i = 1; i args, bool bVerbose = false, const char delimiter = '='): + mapArguments(), iArgC(args.size()), vArgV(NULL), bStrictMode(false), bVerbose(bVerbose), delimiter(delimiter) + { + for(vector::iterator it = args.begin(); it != args.end(); ++it) + { + const char * arg = it->c_str(); + + int sep = 0; + while(arg[sep] != '\0' && arg[sep] != delimiter) + ++sep; + + string value; + + if (arg[sep] != '\0') + value = string(arg + sep + 1); + else + value = "1"; + + mapArguments[string(arg, sep)] = Value(value); + } + + mute(); + } + + int getargc() const { return iArgC; } + + const char** getargv() const { return vArgV; } + + void set_strict_mode() + { + bStrictMode = true; + } + + void unset_strict_mode() + { + bStrictMode = false; + } + + void mute() + { + bVerbose = false; + } + + void loud() + { + bVerbose = true; + } + + void print_arguments(FILE * f = stdout) + { + printf("PRINTOUT OF THE RUNTIME OPTIONS\n"); + + for(map::const_iterator it=mapArguments.begin(); it!=mapArguments.end(); it++) + fprintf(f, "%s: <%s>\n", it->first.c_str(), it->second.asString().c_str()); + + printf("END OF THE PRINTOUT.\n"); + } + + void print_arguments(string path2log) + { + FILE * f = fopen(path2log.c_str(), "w"); + + if (f == NULL) + { + printf("could not save the log to <%s>. Exiting now\n", path2log.c_str()); + exit(-1); + } + + print_arguments(f); + + fclose(f); + } + + vector find(string name) + { + map::iterator itb = mapArguments.lower_bound(name); + map::iterator ite = mapArguments.end(); + + vector retval; + for(map::iterator it = itb; it != ite; ++it) + { + if (it->first.find(name) == string::npos) + break; + + retval.push_back(it->first); + } + return retval; + } +}; diff --git a/postprocessing/mpi-check.h b/postprocessing/mpi-check.h new file mode 100644 index 000000000..31900eaab --- /dev/null +++ b/postprocessing/mpi-check.h @@ -0,0 +1,19 @@ +#include + +#include + +#define MPI_CHECK(ans) do { mpiAssert((ans), __FILE__, __LINE__); } while(0) + +inline void mpiAssert(int code, const char *file, int line, bool abort=true) +{ + if (code != MPI_SUCCESS) + { + char error_string[2048]; + int length_of_error_string = sizeof(error_string); + MPI_Error_string(code, error_string, &length_of_error_string); + + printf("mpiAssert: %s %d %s\n", file, line, error_string); + + MPI_Abort(MPI_COMM_WORLD, code); + } +} diff --git a/postprocessing/ply2vtkpts/Makefile b/postprocessing/ply2vtkpts/Makefile new file mode 100644 index 000000000..ec76ea0e7 --- /dev/null +++ b/postprocessing/ply2vtkpts/Makefile @@ -0,0 +1,17 @@ +CXX ?= CC + +ply2vtk: main.cpp + $(CXX) -Ofast -std=c++11 -fopenmp main.cpp \ + -I/apps/daint/VTK/6.2/gnu_491/include/vtk-6.2 -L/apps/daint/VTK/6.2/gnu_491/lib \ + -lvtkIOImage-6.2 -lvtkCommonDataModel-6.2 -lvtkpng-6.2 -lvtktiff-6.2 \ + -lvtkmetaio-6.2 -lvtkDICOMParser-6.2 -lvtkzlib-6.2 -lvtksys-6.2 \ + -lvtkIOXMLParser-6.2 -lvtkCommonExecutionModel-6.2 -lvtkCommonTransforms-6.2 \ + -lvtkCommonCore-6.2 -lvtkIOXML-6.2 -lvtkexpat-6.2 -lvtkjpeg-6.2 -lvtkIOCore-6.2 \ + -lvtkCommonSystem-6.2 -lvtkCommonTransforms-6.2 -lvtkCommonMath-6.2 \ + -lvtkCommonMisc-6.2 \ + -o ply2vtk + +clean: + rm -f ply2vtk + +.PHONY = clean diff --git a/postprocessing/ply2vtkpts/daint-ply2vtkpts.sh b/postprocessing/ply2vtkpts/daint-ply2vtkpts.sh new file mode 100644 index 000000000..b5c9c0ade --- /dev/null +++ b/postprocessing/ply2vtkpts/daint-ply2vtkpts.sh @@ -0,0 +1,45 @@ +module swap PrgEnv-cray PrgEnv-gnu +module load cray-hdf5-parallel +module unload cray-mpich/7.0.4 +module load cray-mpich/7.1.1 +module unload gcc +module load gcc/4.9.1 +module load vtk + +convert_some() +{ + MYFOLDER=$1 + SRCPATH=$2 + SRCPATTERN=$3 + NVERTPERCELLS=$4 + + mkdir -p $MYFOLDER + + echo `date` "convert_some: $*" >> ${MYFOLDER}/log.txt + + find "$SRCPATH" -name "$SRCPATTERN" > /tmp/asd.txt + + for F in $(cat /tmp/asd.txt) + do + SRC=`basename $F` + + DST=${MYFOLDER}/${SRC%.ply}.vtp + + aprun ./ply2vtk $NVERTPERCELLS $F $DST + done +} + +if (( $# != 3)) +then + echo "usage ./convert-all.sh " + + exit 1 +fi + + +MYFOLDER=$1 #for example "ichip31" +NVERTRBC=$2 #for example 498 +NVERTCTC=$3 #for example 5220 + +convert_some "$MYFOLDER" /scratch/daint/alexeedm/ctc/"$MYFOLDER"/ply/ "rbcs-*.ply" $NVERTRBC +convert_some "$MYFOLDER" /scratch/daint/alexeedm/ctc/"$MYFOLDER"/ply/ "ctcs-*.ply" $NVERTCTC \ No newline at end of file diff --git a/postprocessing/ply2vtkpts/main.cpp b/postprocessing/ply2vtkpts/main.cpp new file mode 100644 index 000000000..c1465f927 --- /dev/null +++ b/postprocessing/ply2vtkpts/main.cpp @@ -0,0 +1,223 @@ +#include +#include + +#include +#include +#include +#include + + +#define MPI_CHECK(ans) do { mpiAssert((ans), __FILE__, __LINE__); } while(0) + +inline void mpiAssert(int code, const char *file, int line, bool abort=true) +{ + if (code != MPI_SUCCESS) + { + char error_string[2048]; + int length_of_error_string = sizeof(error_string); + MPI_Error_string(code, error_string, &length_of_error_string); + + printf("mpiAssert: %s %d %s\n", file, line, error_string); + + MPI_Abort(MPI_COMM_WORLD, code); + } +} + +using namespace std; + +#include +#include +#include +#include +#include +#include + +int dump_vtk_points(const char * dstpath, const int nrbcs, const float * xs, const float * ys, const float * zs ) +{ + vtkSmartPointer points = + vtkSmartPointer::New(); + + for ( unsigned int i = 0; i < nrbcs; ++i ) + points->InsertNextPoint ( xs[i], ys[i], zs[i] ); + + // Create a polydata object and add the points to it. + vtkSmartPointer polydata = + vtkSmartPointer::New(); + polydata->SetPoints(points); + + // Write the file + vtkSmartPointer writer = + vtkSmartPointer::New(); + writer->SetFileName(dstpath); +#if VTK_MAJOR_VERSION <= 5 + writer->SetInput(polydata); +#else + writer->SetInputData(polydata); +#endif + + writer->Write(); + + return EXIT_SUCCESS; +} + +int main(int argc, char ** argv) +{ + MPI_CHECK(MPI_Init(&argc, &argv)); + + int nranks, rank; + MPI_CHECK(MPI_Comm_size(MPI_COMM_WORLD, &nranks)); + MPI_CHECK(MPI_Comm_rank(MPI_COMM_WORLD, &rank)); + + const bool verbose = false; + + if (argc != 4) + { + if (rank == 0) + printf("usage: test \n"); + + exit(EXIT_FAILURE); + } + + const int nvpc = atoi(argv[1]); + const char * path = argv[2]; + const char * dstpath = argv[3]; + + if (rank == 0) + printf("reading at location <%s>\n", path); + + int nallvertices, nallrbcs, headersize; + + const double tstart = omp_get_wtime(); + + if (rank == 0) + { + FILE * f = fopen(path, "r"); + assert(f); + char line[2048]; + + auto eat_line = [&] () + { + fgets(line, 2048, f); + + if (verbose) + printf("reading <%s>\n", line); + }; + + for(int i = 0; i < 3; ++i) + eat_line(); + + int retval = sscanf(line, "element vertex %d\n", &nallvertices); + assert(retval == 1); + + nallrbcs = nallvertices / nvpc; + + if (verbose) + printf("*** nvertices: %d\n", nallvertices); + + for(int i = 0; i < 7; ++i) + eat_line(); + + int nfaces = -1; + retval = sscanf(line, "element face %d\n", &nfaces); + assert(retval == 1); + + if (verbose) + printf("*** nfaces: %d\n", nfaces); + + for(int i = 0; i < 2; ++i) + eat_line(); + + headersize = ftell(f); + + fclose(f); + } + + MPI_CHECK(MPI_Bcast(&nallvertices, 1, MPI_INT, 0, MPI_COMM_WORLD)); + MPI_CHECK(MPI_Bcast(&nallrbcs, 1, MPI_INT, 0, MPI_COMM_WORLD)); + MPI_CHECK(MPI_Bcast(&headersize, 1, MPI_INT, 0, MPI_COMM_WORLD)); + + const double theader = omp_get_wtime(); + + const int myrbcs_size = nallrbcs / nranks + (int)(rank < (nallrbcs % nranks)); + const int myrbcs_start = nallrbcs / nranks * rank + min(rank, nallrbcs % nranks); + + float * data = new float[6 * myrbcs_size * nvpc]; + + { + MPI_File filehandle; + MPI_CHECK( MPI_File_open(MPI_COMM_WORLD, path, MPI_MODE_RDONLY, MPI_INFO_NULL, &filehandle) ); + + MPI_Status status; + MPI_CHECK( MPI_File_read_at(filehandle, headersize + myrbcs_start * nvpc * 6 * sizeof(float), + data, myrbcs_size * nvpc * 6, MPI_FLOAT, &status)); + + MPI_CHECK( MPI_File_close(&filehandle)); + } + + vector coords[3]; + + for(int i = 0; i < 3; ++i) + coords[i].resize(myrbcs_size); + +#pragma omp parallel for + for(int r = 0; r < myrbcs_size; ++r) + { + float com[3] = {0, 0, 0}; + + for(int v = 0; v < nvpc; ++v) + for(int c = 0; c < 3; ++c) + com[c] += data[c + 6 * (v + nvpc * r)]; + + for(int i = 0; i < 3; ++i) + coords[i][r] = com[i] /nvpc; + } + + delete [] data; + + vector allcoords[3]; + + if (rank == 0) + for(int i = 0; i < 3; ++i) + allcoords[i].resize(nranks * ((nallrbcs + nranks - 1) / nranks)); + + for(int i = 0; i < 3; ++i) + { + MPI_CHECK( MPI_Gather(&coords[i].front(), nallrbcs / nranks, MPI_FLOAT, + &allcoords[i].front(), nallrbcs / nranks, MPI_FLOAT, + 0, MPI_COMM_WORLD) ); + + MPI_CHECK( MPI_Gather(&coords[i].back(), 1, MPI_FLOAT, + (&allcoords[i].front()) + nranks * (nallrbcs / nranks), 1, MPI_FLOAT, + 0, MPI_COMM_WORLD) ); + } + + const double tthroughput = omp_get_wtime(); + + if (rank == 0) + dump_vtk_points(dstpath, nallrbcs, &allcoords[0].front(), &allcoords[1].front(), &allcoords[2].front()); + + const double tvtk = omp_get_wtime(); + + if (rank == 0) + { + const double ttotal = tvtk - tstart; + + printf("TOTAL TIME: %.2f\n", ttotal); + + printf("TDISTRIBUTION: HEADER:%.1f%%\tI/O+REDUCE:%.1f%%\tVTK:%.1f%%\t\n", + 100 / ttotal * (theader - tstart), + 100 / ttotal * (tthroughput - theader), + 100 / ttotal * (tvtk - tthroughput)); + + const double memfp_ply = 6. * nallvertices * sizeof(float) / pow(1024., 3); + const double memfp_vtk = 3. * nallrbcs * sizeof(float) / pow(1024., 3); + + printf("THROUGHPUT: %.1f GB/s\n", (memfp_ply + memfp_vtk) / (tthroughput - theader)); + printf("VTK DUMP: %.1f GB/s\n", memfp_vtk / (tvtk - tthroughput)); + } + + MPI_CHECK(MPI_Finalize()); + + return 0; +} + diff --git a/postprocessing/stress/Makefile b/postprocessing/stress/Makefile new file mode 100644 index 000000000..1707fd490 --- /dev/null +++ b/postprocessing/stress/Makefile @@ -0,0 +1,9 @@ +CXX = mpicxx + +stress: main.cpp + $(CXX) main.cpp -I../ -g -O3 -Wno-deprecated-declarations -Wno-unused-result -o stress + +clean: + rm -f stress + +.PHONY = clean diff --git a/postprocessing/stress/example.run b/postprocessing/stress/example.run new file mode 100644 index 000000000..8123d29ce --- /dev/null +++ b/postprocessing/stress/example.run @@ -0,0 +1,28 @@ +make + +echo EXAMPLE1: REDUCE ALL +find ../../mpi-dpd/stress/* | ./stress -origin=0,0,0 -extent=48,48,48 -project=1,1,1 + +echo EXAMPLE2: 1D PROFILE + GNUPLOT +find ../../mpi-dpd/stress/* | ./stress -origin=0,0,0 -extent=48,48,48 -project=1,0,1 > profile.txt +gnuplot -persist <<- END_GNUPLOT + plot "profile.txt" u 1:2 w lp title "sigma_xx", "profile.txt" u 1:3 w lp title "sigma_xy", "profile.txt" u 1:4 w lp title "sigma_xz", "profile.txt" u 1:5 w lp title "sigma_yy", "profile.txt" u 1:6 w lp title "sigma_yz", "profile.txt" u 1:7 w lp title "sigma_zz" +END_GNUPLOT + +echo EXAMPLE3: HEIGHT FIELD + MATPLOTLIB +find ../../mpi-dpd/stress/* | ./stress -origin=0,0,0 -extent=48,48,32 -project=0,1,0 | csplit --suppress-matched - '/^$/' {*} -f channel. +python <<- END_PYTHON +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.cm as cm + +for path in ["channel.00", "channel.01", "channel.02", "channel.03", "channel.04", "channel.05"]: + ncols, nrows = 48, 32 + temp = np.loadtxt(path).T + grid = temp.reshape((nrows, ncols)) + plt.figure(path) + plt.imshow(grid, interpolation='nearest', cmap=cm.gist_rainbow) + +plt.show() + +END_PYTHON \ No newline at end of file diff --git a/postprocessing/stress/main.cpp b/postprocessing/stress/main.cpp new file mode 100644 index 000000000..c90aeb02f --- /dev/null +++ b/postprocessing/stress/main.cpp @@ -0,0 +1,295 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +using namespace std; + +int main(int argc, const char ** argv) +{ + MPI_CHECK( MPI_Init(&argc, (char ***)&argv) ); + + int nranks, rank; + MPI_CHECK( MPI_Comm_rank(MPI_COMM_WORLD, &rank)); + MPI_CHECK( MPI_Comm_size(MPI_COMM_WORLD, &nranks)); + + ArgumentParser argp(argc, argv); + + const bool verbose = argp("-verbose").asBool(false); + const bool avg = argp("-average").asBool(true); + vector origin = argp("-origin").asVecFloat(3); + vector extent = argp("-extent").asVecFloat(3); + vector projectf = argp("-project").asVecFloat(3); + string contributions = argp("-contributions").asString("uf"); + + const double ufactor = contributions.find("u") != string::npos; + const double ffactor = contributions.find("f") != string::npos; + + bool project[3]; + for(int c = 0; c < 3; ++c) + project[c] = projectf[c] != 0; + + int nprojections = 0; + for(int c = 0; c < 3; ++c) + nprojections += project[c]; + + const int noutputchannels = 9; + const size_t chunksize = (1 << 29) / 12 / sizeof(float); + + float * const pbuf = new float[12 * chunksize]; + + float binsize[3]; + for(int c = 0; c < 3; ++c) + binsize[c] = project[c] ? extent[c] : 1; + + int nbins[3]; + for(int c = 0; c < 3; ++c) + nbins[c] = extent[c] / binsize[c]; + + const int ntotbins = nbins[0] * nbins[1] * nbins[2]; + + int * const bincount = new int[ntotbins]; + memset(bincount, 0, sizeof(int) * ntotbins); + + const int noutput = noutputchannels * ntotbins; + + double * const bindata = new double[noutput]; + memset(bindata, 0, sizeof(double) * noutput); + + vector paths; + + { + string myinput; + + if (rank == 0) + for (string line; getline(cin, line);) + myinput += line + "\n"; + + int inputsize = myinput.size(); + MPI_CHECK( MPI_Bcast(&inputsize, 1, MPI_INTEGER, 0, MPI_COMM_WORLD)); + + myinput.resize(inputsize); + + MPI_CHECK( MPI_Bcast(&myinput[0], inputsize, MPI_CHAR, 0, MPI_COMM_WORLD)); + + int c = 0; + istringstream iss(myinput); + + for (string line; getline(iss, line); ++c) + if (c % nranks == rank) + paths.push_back(line); + } + + int numfiles = paths.size(); + + size_t totalfootprint = 0; + double timeIO = 0; + + for(int ipath = 0; ipath < (int)paths.size(); ++ipath) + { + const char * const path = paths[ipath].c_str(); + + if (verbose) + fprintf(stderr, "working on <%s>\n", path); + + int fdin = open(path, O_RDONLY); + + if (!fdin) + { + fprintf(stderr, "can't access <%s> , exiting now.\n", path); + exit(-1); + } + + if (verbose) + perror("reading...\n"); + + const size_t filesize = lseek(fdin, 0, SEEK_END); + + totalfootprint += filesize; + + lseek(fdin, 0, SEEK_SET); + + const size_t nparticles = filesize / 12 / sizeof(float); + assert(filesize % (12 * sizeof(float)) == 0); + + if (verbose) + { + fprintf(stderr, "i have found %d particles\n", (int)nparticles); + fprintf(stderr, "particle chunk %d\n", (int)chunksize); + } + + for(size_t base = 0; base < nparticles; base += chunksize) + { + const int nhotparticles = min(nparticles - base, chunksize); + const size_t nhotbytes = nhotparticles * sizeof(float) * 12; + + size_t nreadbytes = 0; + int start = 0; + + while(start < nhotparticles) + { + const double tstart = MPI_Wtime(); + nreadbytes += read(fdin, pbuf, nhotbytes - nreadbytes); + timeIO += MPI_Wtime() - tstart; + + const int stop = nreadbytes / sizeof(float) / 12; + +#ifndef NDEBUG + if (verbose) + { + float avgs[12]; + for(int i = 0; i < 12; ++i) + avgs[i] = 0; + + for(int i = 0; i < nhotparticles; ++i) + for(int c = 0; c < 12; ++c) + avgs[c] += pbuf[12 * i + c]; + + for(int i = 0; i < 12; ++i) + printf("AVG %d: %.3e\n", i, avgs[i] / nhotparticles); + } +#endif + + for(int i = start; i < stop; ++i) + { + const int srcbase = 12 * i; + + int index[3]; + for(int c = 0; c < 3; ++c) + index[c] = (int)((pbuf[srcbase + c] - origin[c]) / binsize[c]); + + bool valid = true; + for(int c = 0; c < 3; ++c) + valid &= index[c] >= 0 && index[c] < nbins[c]; + + if (!valid) + continue; + + const int binid = index[0] + nbins[0] * (index[1] + nbins[1] * index[2]); + ++bincount[binid]; + + const int dstbase = noutputchannels * binid; + + const int v1[6] = {0, 0, 0, 1, 1, 2}; + const int v2[6] = {0, 1, 2, 1, 2, 2}; + + for(int c = 0; c < 6; ++c) + bindata[dstbase + c] += + ffactor * pbuf[srcbase + 6 + c] + + ufactor * pbuf[srcbase + 3 + v1[c]] * pbuf[srcbase + 3 + v2[c]]; + + for(int c = 0; c < 3; ++c) + bindata[dstbase + 6 + c] += pbuf[srcbase + 3 + c]; + } + + start = stop; + } + } + + close(fdin); + } + + if (rank == 0 && !numfiles) + { + perror("ooops zero files were read. Exiting now.\n"); + exit(-1); + } + + MPI_CHECK( MPI_Reduce(rank ? bincount : MPI_IN_PLACE, bincount, ntotbins, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD) ); + MPI_CHECK( MPI_Reduce(rank ? bindata : MPI_IN_PLACE, bindata, noutput, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) ); + MPI_CHECK( MPI_Reduce(rank ? &timeIO : MPI_IN_PLACE, &timeIO, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) ); + MPI_CHECK( MPI_Reduce(rank ? &totalfootprint : MPI_IN_PLACE, &totalfootprint, 1, MPI_OFFSET, MPI_SUM, 0, MPI_COMM_WORLD) ); + + if (rank) + goto finalize; + + if (avg) + for(int i = 0; i < ntotbins; ++i) + for(int c = 0; c < noutputchannels; ++c) + bindata[noutputchannels * i + c] /= bincount[i]; + + if (nprojections == 3) + { + assert(noutput == noutputchannels); + + for(int c = 0; c < noutputchannels; ++c) + printf("%+.3e\t", bindata[c]); + + printf("\n"); + } + else if (nprojections == 2) + { + int ctr = 0; + for(int iz = 0; iz < nbins[2]; ++iz) + for(int iy = 0; iy < nbins[1]; ++iy) + for(int ix = 0; ix < nbins[0]; ++ix) + { + printf("%03d ", ctr); + + for(int c = 0; c < noutputchannels; ++c) + printf("%+.4e ", bindata[noutputchannels * ctr + c]); + + printf("\n"); + + ++ctr; + } + } + else if (nprojections == 1) + { + int nx = 0; + for(int c = 0; c < 3; ++c) + if (nbins[c] > 1) + { + nx = nbins[c]; + break; + } + + for(int c = 0; c < noutputchannels; ++c) + { + int ctr = 0; + + for(int iz = 0; iz < nbins[2]; ++iz) + for(int iy = 0; iy < nbins[1]; ++iy) + for(int ix = 0; ix < nbins[0]; ++ix) + { + printf("%+.5e ", bindata[noutputchannels * ctr + c]); + + ++ctr; + + if (ctr % nx == 0) + printf("\n"); + } + + if (c < noutputchannels - 1) + printf("\n"); + } + } + else + { + perror("woops invalid number of projections. Exiting now...\n"); + exit(-1); + } + + if (verbose) + perror("all is done. ciao.\n"); + + fprintf(stderr, "total footprint: %.3f MB, I/O time: %.3f ms\n", totalfootprint * 1. / 1024 / 1024, timeIO * 1e3); + fprintf(stderr, "read throughput: %.3f GB/s\n", totalfootprint /( 1024 * 1024) / timeIO / 1024); + +finalize: + + delete [] pbuf; + delete [] bincount; + delete [] bindata; + + MPI_CHECK( MPI_Finalize() ); + + return 0; +} diff --git a/proof-of-concept/commets/README.md b/proof-of-concept/commets/README.md new file mode 100644 index 000000000..e26e4190e --- /dev/null +++ b/proof-of-concept/commets/README.md @@ -0,0 +1 @@ +Tools to transform comment blocks in uDeviceX diff --git a/proof-of-concept/commets/comments.awk b/proof-of-concept/commets/comments.awk new file mode 100755 index 000000000..2affd7c51 --- /dev/null +++ b/proof-of-concept/commets/comments.awk @@ -0,0 +1,113 @@ +#!/usr/bin/awk -f + +function err(s) { + printf "(comments.awk) %s\n", s + exit +} + +function logg(s, cmd) { + cmd = "cat 1>&2" + printf "(comments.awk) %s\n", s | cmd + close(cmd) +} + +function read_file(fn, line, sep) { # reads file content into string `a' + while (getline line < fn > 0) { + a = a sep line + sep = "\n" + } + if (line != sep) a = a sep +} + + +function write_file() { + printf "%s", a +} + +function ss(i, l) { # substring of `a' starting at `i' with length `l' + return substr(a, i, l) +} + +function ch(i) { # character `i' in `a' + return ss(i, 1) +} + +function spacep(c) { + return c == " " || c == "\n" || c == "\t" +} + +function find_comment_block( i, n) { # sets `lo', `hi' for the + # comments block + # find comment start `/*' + n = length(a); i = 1 + while (spacep(ch(i))) i++; # eat spaces + if (ss(i, 2) != "/*") return !HAS_COMMENT_BLOCK + + lo = i # start of the comment block + + while (1) { + i++ + if (i == n) return !HAS_COMMENT_BLOCK + if (ss(i, 2) == "*/") {hi = i + 1; return HAS_COMMENT_BLOCK} + } +} + +function remove(s, lo, hi, head, tail) { + head = substr(s, 1 , lo - 1) + tail = substr(s, hi + 1 ) + return head tail +} + +function replace(s, lo, hi, t, head, tail) { + head = substr(s, 1 , lo - 1) + tail = substr(s, hi + 1 ) + return head t tail +} + +function end_index(s, t, i) { + i = index(s, t) + if (i == 0) return 0 + return i + length(t) +} + +function transform_comment( lo, hi, be, end) { + beg = " * Users are NOT authorized" + end = "permission from the author of this file." + + lo = index(cm, beg) + if (lo == 0) return 0 + hi = end_index(cm, end) + if (hi == 0) err("I found beg but cannot find end in comments block in " fn) + + cm = remove(cm, lo, hi) + return 1 +} + +function transfrom( rc) { + if (!find_comment_block()) return + + cm = substr(a, lo, hi - lo + 1) # get comment block as a string + rc = transform_comment() # transfrom comment block + + if (!rc) return + + logg("trans: " fn) + a = replace(a, lo, hi, cm) +} + +function read_cpy(fn, line, sep) { # reads file content into string `a' + while (getline line < fn > 0) { + cpy = cpy sep line + sep = "\n" + } + if (line != sep) cpy = cpy sep +} + +BEGIN { + HAS_COMMENT_BLOCK = 1 + #read_cpy(fcpy = ARGV[1]) + + read_file(fn = ARGV[1]) + transfrom() + write_file() +} diff --git a/proof-of-concept/commets/comments.sh b/proof-of-concept/commets/comments.sh new file mode 100755 index 000000000..3ad55beba --- /dev/null +++ b/proof-of-concept/commets/comments.sh @@ -0,0 +1,14 @@ +#!/bin/bash + +f=mpi-dpd/common.cu +c=proof-of-concept/commets/comments.awk +t=/tmp/comments.$$.cu + +$c $f > d +diff $f d + +for f in `find . -type f -name '*.cu' -or -name '*.h' -or -name '*.cpp' -or -name '*.cc'` +do + $c $f > $t + cp $t $f +done > d diff --git a/proof-of-concept/create-hdf5-file/main.cpp b/proof-of-concept/create-hdf5-file/main.cpp index c0a624e88..c44f12781 100644 --- a/proof-of-concept/create-hdf5-file/main.cpp +++ b/proof-of-concept/create-hdf5-file/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2015-01-27. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/cuda-aware-mpi/main.cu b/proof-of-concept/cuda-aware-mpi/main.cu index 5dcf7e1f1..790b41d1f 100644 --- a/proof-of-concept/cuda-aware-mpi/main.cu +++ b/proof-of-concept/cuda-aware-mpi/main.cu @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-24. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/funnel-obstacle/funnel-obstacle.cpp b/proof-of-concept/funnel-obstacle/funnel-obstacle.cpp index 78ffa671c..256329941 100644 --- a/proof-of-concept/funnel-obstacle/funnel-obstacle.cpp +++ b/proof-of-concept/funnel-obstacle/funnel-obstacle.cpp @@ -5,9 +5,6 @@ * Created and authored by Kirill Lykov on 2014-07-31. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "funnel-obstacle.h" diff --git a/proof-of-concept/funnel-obstacle/funnel-obstacle.h b/proof-of-concept/funnel-obstacle/funnel-obstacle.h index 65f507373..79505dd83 100644 --- a/proof-of-concept/funnel-obstacle/funnel-obstacle.h +++ b/proof-of-concept/funnel-obstacle/funnel-obstacle.h @@ -5,9 +5,6 @@ * Created and authored by Kirill Lykov on 2014-07-31. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #ifndef LS_OBSTACLE_H_ #define LS_OBSTACLE_H_ diff --git a/proof-of-concept/funnel-obstacle/obst-check.cpp b/proof-of-concept/funnel-obstacle/obst-check.cpp index 9670368d9..335197aaa 100644 --- a/proof-of-concept/funnel-obstacle/obst-check.cpp +++ b/proof-of-concept/funnel-obstacle/obst-check.cpp @@ -5,9 +5,6 @@ * Created and authored by Kirill Lykov on 2014-07-30. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include #include diff --git a/proof-of-concept/logistic_rng/logistic.h b/proof-of-concept/logistic_rng/logistic.h index 49893064f..c3296f9c8 100644 --- a/proof-of-concept/logistic_rng/logistic.h +++ b/proof-of-concept/logistic_rng/logistic.h @@ -5,9 +5,6 @@ * Created and authored by Yu-Hang Tang on 2015-03-20. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #ifndef __LOGISTIC_RNG__ diff --git a/proof-of-concept/logistic_rng/test.cu b/proof-of-concept/logistic_rng/test.cu index f96c1fb7a..b89bbe9fc 100644 --- a/proof-of-concept/logistic_rng/test.cu +++ b/proof-of-concept/logistic_rng/test.cu @@ -5,9 +5,6 @@ * Created and authored by Yu-Hang Tang on 2015-03-20. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ diff --git a/proof-of-concept/poiseuille-dpd/main.cpp b/proof-of-concept/poiseuille-dpd/main.cpp index 3744ea017..46697007e 100644 --- a/proof-of-concept/poiseuille-dpd/main.cpp +++ b/proof-of-concept/poiseuille-dpd/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-dpd/main.cpp b/proof-of-concept/vanilla-dpd/main.cpp index 293cbae7a..7f9325391 100644 --- a/proof-of-concept/vanilla-dpd/main.cpp +++ b/proof-of-concept/vanilla-dpd/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-07-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-microfluidics/main.cpp b/proof-of-concept/vanilla-microfluidics/main.cpp index a6014fd4e..07ba2dff2 100644 --- a/proof-of-concept/vanilla-microfluidics/main.cpp +++ b/proof-of-concept/vanilla-microfluidics/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-08-08. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-mpi-dpd/common.cpp b/proof-of-concept/vanilla-mpi-dpd/common.cpp index eabc39919..fa41eb938 100644 --- a/proof-of-concept/vanilla-mpi-dpd/common.cpp +++ b/proof-of-concept/vanilla-mpi-dpd/common.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "common.h" diff --git a/proof-of-concept/vanilla-mpi-dpd/common.h b/proof-of-concept/vanilla-mpi-dpd/common.h index ffaedc26c..43451820c 100644 --- a/proof-of-concept/vanilla-mpi-dpd/common.h +++ b/proof-of-concept/vanilla-mpi-dpd/common.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.cpp b/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.cpp index f3dc18345..8a5ede318 100644 --- a/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.cpp +++ b/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.h b/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.h index bd2905405..58adc11ec 100644 --- a/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.h +++ b/proof-of-concept/vanilla-mpi-dpd/dpd-interactions.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/proof-of-concept/vanilla-mpi-dpd/main.cpp b/proof-of-concept/vanilla-mpi-dpd/main.cpp index 01bb8dc2c..c54672b89 100644 --- a/proof-of-concept/vanilla-mpi-dpd/main.cpp +++ b/proof-of-concept/vanilla-mpi-dpd/main.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-06. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.cpp b/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.cpp index 96653c941..25f1e33a0 100644 --- a/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.cpp +++ b/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.cpp @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.h b/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.h index 809bfe678..2467c924c 100644 --- a/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.h +++ b/proof-of-concept/vanilla-mpi-dpd/redistribute-particles.h @@ -5,9 +5,6 @@ * Created and authored by Diego Rossinelli on 2014-11-07. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/proof-of-concept/vanilla-walls/funnel-bouncer.cpp b/proof-of-concept/vanilla-walls/funnel-bouncer.cpp index 6183fa12d..4f9a39de2 100644 --- a/proof-of-concept/vanilla-walls/funnel-bouncer.cpp +++ b/proof-of-concept/vanilla-walls/funnel-bouncer.cpp @@ -6,9 +6,6 @@ * Authored by Diego Rossinelli on 2014-08-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "funnel-bouncer.h" diff --git a/proof-of-concept/vanilla-walls/funnel-bouncer.h b/proof-of-concept/vanilla-walls/funnel-bouncer.h index 7983e8a10..cc6f2dce3 100644 --- a/proof-of-concept/vanilla-walls/funnel-bouncer.h +++ b/proof-of-concept/vanilla-walls/funnel-bouncer.h @@ -6,9 +6,6 @@ * Major editing from Kirill Lykov on 2014-08-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/proof-of-concept/vanilla-walls/main.cpp b/proof-of-concept/vanilla-walls/main.cpp index 84c81d3b0..0080a04aa 100644 --- a/proof-of-concept/vanilla-walls/main.cpp +++ b/proof-of-concept/vanilla-walls/main.cpp @@ -6,9 +6,6 @@ * Major editing from Kirill Lykov on 2014-08-04 * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/proof-of-concept/vanilla-walls/particles.cpp b/proof-of-concept/vanilla-walls/particles.cpp index 68e654d02..4ce2757b9 100644 --- a/proof-of-concept/vanilla-walls/particles.cpp +++ b/proof-of-concept/vanilla-walls/particles.cpp @@ -6,9 +6,6 @@ * Authored by Diego Rossinelli on 2014-08-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include "particles.h" diff --git a/proof-of-concept/vanilla-walls/particles.h b/proof-of-concept/vanilla-walls/particles.h index acb0b63c4..167be31e1 100644 --- a/proof-of-concept/vanilla-walls/particles.h +++ b/proof-of-concept/vanilla-walls/particles.h @@ -6,9 +6,6 @@ * Authored by Diego Rossinelli on 2014-08-12. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #pragma once diff --git a/tests/contact/Makefile b/tests/contact/Makefile new file mode 100644 index 000000000..e14c7da00 --- /dev/null +++ b/tests/contact/Makefile @@ -0,0 +1,26 @@ +-include ../../mpi-dpd/.cache.Makefile + +NVCC ?= nvcc -ccbin $(CXX) +ARCH_VAL ?= compute_35 +CODE_VAL ?= sm_35 + + +NVCCFLAGS += -I$(HDF5_DIR)/include -Xcudafe "--diag_suppress=unrecognized_gcc_pragma" +NVCCFLAGS += -arch $(ARCH_VAL) -code $(CODE_VAL) -O3 -use_fast_math -g -DNDEBUG -Xcompiler "-fopenmp" +CXXFLAGS += -L../../cuda-dpd/dpd -L../../cuda-rbc/ -L../../cuda-ctc/ -O3 -g -std=c++11 -DNDEBUG -fopenmp +NVCCFLAGS += -I../../cuda-dpd/dpd -I../../cuda-dpd/ -I../../cuda-rbc/ -I../../cuda-ctc -I../../mpi-dpd +LIBS = -lcuda-dpd -lcuda-rbc -lcuda-ctc -lcudart -ldl -lz -fopenmp + +testcontact: ../../mpi-dpd/test testcontact.o + make -C ../../mpi-dpd + rm -f ../../mpi-dpd/main.o + $(CXX) $(CXXFLAGS) testcontact.o ../../mpi-dpd/*.o $(LIBS) -o testcontact + +../test: + make -C ../ + cp ../*.o . + +testcontact.o: testcontact.cu ../../mpi-dpd/argument-parser.h ../../mpi-dpd/common.h ../../mpi-dpd/containers.h ../../mpi-dpd/contact.h ../../cuda-dpd/dpd-rng.h + $(NVCC) $(NVCCFLAGS) testcontact.cu -c -o testcontact.o + +.PHONY: ../test diff --git a/tests/contact/testcontact.cu b/tests/contact/testcontact.cu new file mode 100644 index 000000000..3d799fc1d --- /dev/null +++ b/tests/contact/testcontact.cu @@ -0,0 +1,300 @@ +/* + * main.cu + * ctc PANDA + * + * Created by Dmitry Alexeev on Oct 20, 2015 + * Copyright 2015 ETH Zurich. All rights reserved. + * + */ + +/* + * main.cu + * Part of uDeviceX/mpi-dpd/ + * + * Created and authored by Diego Rossinelli on 2014-11-14. + * Copyright 2015. All rights reserved. + * + * Users are NOT authorized + * to employ the present software for their own publications + * before getting a written permission from the author of this file. + */ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + enum + { + XCELLS = XSIZE_SUBDOMAIN, + YCELLS = YSIZE_SUBDOMAIN, + ZCELLS = ZSIZE_SUBDOMAIN, + XOFFSET = XCELLS / 2, + YOFFSET = YCELLS / 2, + ZOFFSET = ZCELLS / 2 + }; +using namespace std; + +float tend, couette; +bool walls, pushtheflow, doublepoiseuille, rbcs, ctcs, xyz_dumps, hdf5field_dumps, hdf5part_dumps, is_mps_enabled, adjust_message_sizes, contactforces, stress; +int steps_per_report, steps_per_dump, wall_creation_stepid, nvtxstart, nvtxstop; + +LocalComm localcomm; + +static const float ljsigma = 0.5; +static const float ljsigma2 = ljsigma * ljsigma; + +template +inline float _viscosity_function(float x) +{ + return sqrtf(viscosity_function(x)); +} + +template<> inline float _viscosity_function<1>(float x) { return sqrtf(x); } +template<> inline float _viscosity_function<0>(float x){ return x; } + +int main(int argc, char ** argv) +{ + CUDA_CHECK(cudaSetDevice(0)); + CUDA_CHECK(cudaDeviceReset()); + + { + is_mps_enabled = false; + + const char * mps_variables[] = { + "CRAY_CUDA_MPS", + "CUDA_MPS", + "CRAY_CUDA_PROXY", + "CUDA_PROXY" + }; + + for(int i = 0; i < 4; ++i) + is_mps_enabled |= getenv(mps_variables[i])!= NULL && atoi(getenv(mps_variables[i])) != 0; + } + + int nranks, rank; + MPI_CHECK(MPI_Init(&argc, &argv)); + MPI_CHECK( MPI_Comm_size(MPI_COMM_WORLD, &nranks) ); + MPI_CHECK( MPI_Comm_rank(MPI_COMM_WORLD, &rank) ); + MPI_Comm activecomm = MPI_COMM_WORLD; + + bool reordering = true; + const char * env_reorder = getenv("MPICH_RANK_REORDER_METHOD"); + + MPI_Comm cartcomm; + int periods[] = {1, 1, 1}; + int ranks[] = {1, 1, 1}; + + + MPI_CHECK( MPI_Cart_create(activecomm, 3, ranks, periods, (int)reordering, &cartcomm) ); + activecomm = cartcomm; + + { + MPI_CHECK(MPI_Barrier(activecomm)); + localcomm.initialize(activecomm); + + MPI_CHECK(MPI_Barrier(activecomm)); + + // test here + const size_t myseed = 0x563d00cf;//time(NULL); + srand48(myseed); + printf("myseed: 0x%x\n", myseed); + + int n = 25e3; + vector ic(n); + vector acc(n); + for (int i=0; i gpuacc(n); + + Logistic::KISS local_trunk = Logistic::KISS(7119 - rank, 187 + rank, 18278, 15674); + + const double center[3] = { -XSIZE_SUBDOMAIN/2, -YSIZE_SUBDOMAIN/2, -ZSIZE_SUBDOMAIN/2} ;//YSIZE_SUBDOMAIN/2 -}; + const double halfwidth[3] = {XSIZE_SUBDOMAIN/10., YSIZE_SUBDOMAIN/10., ZSIZE_SUBDOMAIN / 10.}; + + for(int i = 0; i < n; ++i) + { + ic[i].x[0] = center[0] + halfwidth[0] * 2 * (drand48() - 0.5); + ic[i].x[1] = center[1] + halfwidth[1] * 2 * (drand48() - 0.5); + ic[i].x[2] = center[2] + halfwidth[2] * 2 * (drand48() - 0.5); + ic[i].u[0] = 0.5 - drand48(); + ic[i].u[1] = 0.5 - drand48(); + ic[i].u[2] = 0.5 - drand48(); + } + + if (false)//if (true) + { + ic.resize(2); + acc.resize(2); + gpuacc.resize(2); + n = 2; + //ic[0].x[0] = 24.413; ic[0].x[1] = +14.924; ic[0].x[2] = +7.326; + //ic[1].x[0] = +23.895; ic[1].x[1] = +14.887; ic[1].x[2] = +7.455 ; + + ic[0].x[0] = +23.670 ; ic[0].x[1] =+23.494; ic[0].x[2] =-18.980; + ic[1].x[0] = -23.851 ; ic[1].x[1] =+23.696; ic[1].x[2] =-18.790; + } + + float seed = local_trunk.get_float(); + +#pragma omp parallel for + for (int i=0; i= 1) + continue; + + const double invr2 = invrij * invrij; + const double t2 = ljsigma2 * invr2; + const double t4 = t2 * t2; + const double t6 = t4 * t2; + const double lj = min(1e4f, max(0.f, 24.f * invrij * t6 * (2.f * t6 - 1.f))); + + const double wr = _viscosity_function<0>(1.f - rij); + + const double xr = _xr * invrij; + const double yr = _yr * invrij; + const double zr = _zr * invrij; + + const double strength = lj; + + const double xinteraction = strength * xr; + const double yinteraction = strength * yr; + const double zinteraction = strength * zr; + + acc[i].a[0] += xinteraction; + acc[i].a[1] += yinteraction; + acc[i].a[2] += zinteraction; + } + + ParticleArray p; + p.resize(n); + + CUDA_CHECK( cudaMemcpy(p.xyzuvw.data, &ic[0], n * sizeof(Particle), cudaMemcpyHostToDevice) ); + CUDA_CHECK( cudaMemset(p.axayaz.data, 0, n * sizeof(Acceleration)) ); + std::vector wsolutes; + wsolutes.push_back(ParticlesWrap(p.xyzuvw.data, n, p.axayaz.data)); + + ComputeContact contact(cartcomm); + SoluteExchange solutex(cartcomm); + + solutex.attach_halocomputation(contact); + contact.attach_bulk(wsolutes); + + solutex.bind_solutes(wsolutes); + solutex.pack_p(0); + solutex.post_p(0, 0); + solutex.recv_p(0); + solutex.halo(0, 0); + solutex.post_a(); + solutex.recv_a(0); + + CUDA_CHECK( cudaMemcpy(&gpuacc[0], p.axayaz.data, n * sizeof(Acceleration), cudaMemcpyDeviceToHost) ); + + { + double fx = 0, fy = 0, fz = 0; + double hfx = 0, hfy = 0, hfz = 0; + + for (int i=0; i= tol && fabs(err) >= tol; + + if (failed) + printf("p %d c %d: %e ref: %e -> %e %e\n", i / 3, i % 3, res[i], ref[i], err, relerr); + + if (i % 3 == 2 && failed) + { + const int pid = i/3; + + const bool inside = + ic[i].x[0] >= -XOFFSET && ic[pid].x[0] < XOFFSET && + ic[i].x[1] >= -YOFFSET && ic[pid].x[1] < YOFFSET && + ic[i].x[2] >= -ZOFFSET && ic[pid].x[2] < ZOFFSET ; + + printf("%d: CPU [%+.3f %+.3f %+.3f] GPU [%+.3f %+.3f %+.3f] -> p %+.3f %+.3f %+.3f -> inside: %d\n", + i, acc[pid].a[0], acc[pid].a[1], acc[pid].a[2], + gpuacc[pid].a[0], gpuacc[pid].a[1], gpuacc[pid].a[2], + ic[pid].x[0], ic[pid].x[1], ic[pid].x[2], inside); + + failed = false; + } + + assert(fabs(relerr) < tol || fabs(err) < tol); + + l1 += fabs(err); + l1_rel += fabs(relerr); + + linf = std::max(linf, fabs(err)); + linf_rel = std::max(linf_rel, fabs(relerr)); + } + + printf("l-infinity errors: %.03e (absolute) %.03e (relative)\n", linf, linf_rel); + printf(" l-1 errors: %.03e (absolute) %.03e (relative)\n", l1, l1_rel); + } + } + + if (activecomm != cartcomm) + MPI_CHECK(MPI_Comm_free(&activecomm)); + + MPI_CHECK(MPI_Comm_free(&cartcomm)); + + MPI_CHECK(MPI_Finalize()); + + CUDA_CHECK(cudaDeviceSynchronize()); + + CUDA_CHECK(cudaDeviceReset()); + + return 0; +} diff --git a/halo-bench/Makefile b/tests/halo-bench/Makefile similarity index 100% rename from halo-bench/Makefile rename to tests/halo-bench/Makefile diff --git a/halo-bench/byte_latency.cpp b/tests/halo-bench/byte_latency.cpp similarity index 100% rename from halo-bench/byte_latency.cpp rename to tests/halo-bench/byte_latency.cpp diff --git a/halo-bench/halo-exchanger.cu b/tests/halo-bench/halo-exchanger.cu similarity index 99% rename from halo-bench/halo-exchanger.cu rename to tests/halo-bench/halo-exchanger.cu index 9b540b282..a31b0a237 100644 --- a/halo-bench/halo-exchanger.cu +++ b/tests/halo-bench/halo-exchanger.cu @@ -5,9 +5,6 @@ * Created and authored by Panagiotis Chatzidoukas on 2015-03-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/halo-bench/halo_bench.cpp b/tests/halo-bench/halo_bench.cpp similarity index 99% rename from halo-bench/halo_bench.cpp rename to tests/halo-bench/halo_bench.cpp index 63bdc01e2..9ef3905de 100644 --- a/halo-bench/halo_bench.cpp +++ b/tests/halo-bench/halo_bench.cpp @@ -5,9 +5,6 @@ * Created and authored by Panagiotis Chatzidoukas on 2015-03-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/halo-bench/hpm.cpp b/tests/halo-bench/hpm.cpp similarity index 94% rename from halo-bench/hpm.cpp rename to tests/halo-bench/hpm.cpp index 752004e06..fd780e5c6 100644 --- a/halo-bench/hpm.cpp +++ b/tests/halo-bench/hpm.cpp @@ -5,9 +5,6 @@ * Created and authored by Panagiotis Chatzidoukas on 2015-03-09. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #include diff --git a/halo-bench/mesh_distances.cpp b/tests/halo-bench/mesh_distances.cpp similarity index 100% rename from halo-bench/mesh_distances.cpp rename to tests/halo-bench/mesh_distances.cpp diff --git a/halo-bench/mesh_topo.cpp b/tests/halo-bench/mesh_topo.cpp similarity index 100% rename from halo-bench/mesh_topo.cpp rename to tests/halo-bench/mesh_topo.cpp diff --git a/halo-bench/osu_latency.c b/tests/halo-bench/osu_latency.c similarity index 100% rename from halo-bench/osu_latency.c rename to tests/halo-bench/osu_latency.c diff --git a/halo-bench/osu_latency_rdp.cpp b/tests/halo-bench/osu_latency_rdp.cpp similarity index 96% rename from halo-bench/osu_latency_rdp.cpp rename to tests/halo-bench/osu_latency_rdp.cpp index 46b165237..1c6c9a455 100644 --- a/halo-bench/osu_latency_rdp.cpp +++ b/tests/halo-bench/osu_latency_rdp.cpp @@ -5,9 +5,6 @@ * Created and authored by phadjido 2015-03-12 on 18:59:02. * Copyright 2015. All rights reserved. * - * Users are NOT authorized - * to employ the present software for their own publications - * before getting a written permission from the author of this file. */ #define BENCHMARK "OSU MPI%s Latency Test" diff --git a/halo-bench/scripts/env.sh b/tests/halo-bench/scripts/env.sh similarity index 100% rename from halo-bench/scripts/env.sh rename to tests/halo-bench/scripts/env.sh diff --git a/halo-bench/scripts/exp_12x12x2.sh b/tests/halo-bench/scripts/exp_12x12x2.sh similarity index 100% rename from halo-bench/scripts/exp_12x12x2.sh rename to tests/halo-bench/scripts/exp_12x12x2.sh diff --git a/halo-bench/scripts/exp_14x14x4.sh b/tests/halo-bench/scripts/exp_14x14x4.sh similarity index 100% rename from halo-bench/scripts/exp_14x14x4.sh rename to tests/halo-bench/scripts/exp_14x14x4.sh diff --git a/halo-bench/scripts/exp_3x3x3.sh b/tests/halo-bench/scripts/exp_3x3x3.sh similarity index 100% rename from halo-bench/scripts/exp_3x3x3.sh rename to tests/halo-bench/scripts/exp_3x3x3.sh diff --git a/halo-bench/scripts/exp_3x3x3_new.sh b/tests/halo-bench/scripts/exp_3x3x3_new.sh similarity index 100% rename from halo-bench/scripts/exp_3x3x3_new.sh rename to tests/halo-bench/scripts/exp_3x3x3_new.sh diff --git a/halo-bench/scripts/modules.sh b/tests/halo-bench/scripts/modules.sh similarity index 100% rename from halo-bench/scripts/modules.sh rename to tests/halo-bench/scripts/modules.sh diff --git a/halo-bench/scripts/unsetenv.sh b/tests/halo-bench/scripts/unsetenv.sh similarity index 100% rename from halo-bench/scripts/unsetenv.sh rename to tests/halo-bench/scripts/unsetenv.sh