tdoc contains sphere output, porosity implementation begun - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit b2bcb1d1fd37276c044d071197d7388cf7e1d392
(DIR) parent b9800d3a2ef6aa3ec0a0bbf31643f7968d03e01f
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Thu, 6 Dec 2012 21:39:41 +0100
doc contains sphere output, porosity implementation begun
Diffstat:
M doc/sphinx/conf.py | 2 +-
M doc/sphinx/introduction.rst | 17 +++++++++++++++++
M doc/sphinx/sphere_internals.rst | 6 ++++++
M src/Makefile | 16 ++++++++++++----
M src/main.cpp | 10 ++++++----
M src/sphere.cpp | 123 ++++++++++++++++++++++++++++---
M src/sphere.h | 6 +++++-
7 files changed, 160 insertions(+), 20 deletions(-)
---
(DIR) diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py
t@@ -32,7 +32,7 @@ sys.path.insert(0, os.path.abspath('../../python/'))
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
#extensions = []
-extensions = ['sphinx.ext.autodoc','breathe','sphinx.ext.pngmath']
+extensions = ['sphinx.ext.autodoc','breathe','sphinx.ext.pngmath','sphinxcontrib.programoutput']
breathe_projects = { "sphere": "../doxygen/xml/" }
breathe_default_project = "sphere"
(DIR) diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst
t@@ -30,8 +30,12 @@ Optional tools, required for simulation setup and data processing:
Optional tools, required for building the documentation:
* `Sphinx <http://sphinx-doc.org>`_
+
+ * `sphinxcontrib-programoutput <http://packages.python.org/sphinxcontrib-programoutput/>`_
+
* `Doxygen <http://www.stack.nl/~dimitri/doxygen/>`_
* `Breathe <http://michaeljones.github.com/breathe/>`_
+ * `dvipng <http://www.nongnu.org/dvipng/>`_
`Git <http://git-scm.com>`_ is used as the distributed version control system platform, and the source code is maintained at `Github <https://github.com/anders-dc/sphere/>`_. *sphere* is licensed under the `GNU Public License, v.3 <https://www.gnu.org/licenses/gpl.html>`_.
t@@ -47,6 +51,10 @@ If successfull, the GNU Makefile will create the required data folders, object f
$ ./sphere_* --version
+The output should look similar to this:
+
+.. program-output:: ../../sphere_linux_X86_64 --version
+
The documentation can be read in the `reStructuredText <http://docutils.sourceforge.net/docs/ref/rst/restructuredtext.html>`_-format in the ``doc/sphinx/`` folder, or build into e.g. HTML or PDF format with the following commands::
$ cd doc/sphinx
t@@ -57,3 +65,12 @@ To see all available output formats, execute::
$ make help
+
+Work flow
+---------
+After compiling the *sphere* binary, the procedure of a creating and handling a simulation is typically arranged in the following order:
+ * Setup of particle assemblage, physical properties and conditions using the Python API.
+ * Execution of *sphere* software, which simulates the particle behavior as a function of time, as a result of the conditions initially specified in the input file.
+ * Inspection, analysis, interpretation and visualization of *sphere* output in Python, and/or scene rendering using the built-in ray tracer.
+
+
(DIR) diff --git a/doc/sphinx/sphere_internals.rst b/doc/sphinx/sphere_internals.rst
t@@ -1,5 +1,11 @@
sphere internals
================
+
+The *sphere* executable has the following options:
+
+.. command-output:: ../../sphere_linux_X86_64 --help
+
+
After compiling the *sphere* binary (see sub-section \ref{subsec:compilation}), the procedure of a creating and handling a simulation is typically arranged in the following order:
\item Setup of particle assemblage, physical properties and conditions using the Python API, described in section \ref{sec:ModelSetup}, page \pageref{sec:ModelSetup}.
\item Execution of *sphere* software, which simulates the particle behavior as a function of time, as a result of the conditions initially specified in the input file. Described in section \ref{sec:Simulation}, page \pageref{sec:Simulation}.
(DIR) diff --git a/src/Makefile b/src/Makefile
t@@ -37,7 +37,7 @@ CCFILES=main.cpp file_io.cpp sphere.cpp
CUFILES=device.cu utility.cu
CCOBJECTS=$(CCFILES:.cpp=.o)
CUOBJECTS=$(CUFILES:.cu=.o)
- OBJECTS=$(CCOBJECTS) $(CUOBJECTS)
+OBJECTS=$(CCOBJECTS) $(CUOBJECTS)
# If a header-file changes, update all objects
CDEPS=*.h
t@@ -53,10 +53,10 @@ ifneq ($(DARWIN),) # If OS X
#GL_LIBS+=-framework glut -framework OpenGL
SDKPATH=/Developer/GPU\ Computing/C/common
EXECUTABLE=../sphere_darwin_$(ARCH)
- CCFLAGS+=-arch x86_64
+ CCFLAGS+=-arch x86_64 -lopenmp
LDFLAGS=-L$(SDKPATH)/lib/darwin
else # If Linux
- CCFLAGS+=-DUNIX
+ CCFLAGS+=-DUNIX -fopenmp
#GL_LIBS+=-L/usr/X11/lib -lglut -lGL -lGLU
SDKPATH=$(HOME)/NVIDIA_GPU_Computing_SDK/C/common
EXECUTABLE=../sphere_linux_$(ARCH)
t@@ -70,12 +70,18 @@ LDFLAGS+=-L$(SDKPATH)/../../shared/lib -L$(SDKPATH)/../lib
LDFLAGS+=-lcutil_x86_64 -lcuda -lcudart
LDFLAGS+=-fopenmp
-all: $(CCFILES) $(CUFILES) $(EXECUTABLE) folders
+all: $(CCFILES) $(CUFILES) $(EXECUTABLE) folders ../porosity
$(EXECUTABLE): $(OBJECTS)
$(LINKER) $(OBJECTS) $(LDFLAGS) -o $@
+../porosity: porosity.o sphere.o device.o file_io.o utility.o
+ $(LINKER) $^ $(LDFLAGS) -o $@
+
+porosity.o: porosity.cpp $(CDEPS)
+ $(CC) $(CCFLAGS) $(INCLUDES) -c $< -o $@
+
utility.o: utility.cu
$(NVCC) $(NVCCFLAGS) $(INCLUDES) -c $< -o $@
t@@ -105,6 +111,8 @@ folders: ../output ../img_out ../input
../img_out:
mkdir $@
+doc:
+ $(MAKE) -C ../doc
edit:
$(EDITOR) Makefile $(CCFILES) $(CUFILES) *.h *.cuh
(DIR) diff --git a/src/main.cpp b/src/main.cpp
t@@ -50,10 +50,10 @@ int main(const int argc, const char *argv[])
<< "-q, --quiet\t\tsuppress status messages to stdout\n"
<< "-n, --dry\t\tshow key experiment parameters and quit\n"
<< "-r, --render\t\trender input files instead of simulating temporal evolution\n"
- << "-dc, --dont-check\t\tdon't check values before running\n"
- << "Raytracer (-r) specific options:\n"
+ << "-dc, --dont-check\tdon't check values before running\n"
+ << "\nRaytracer (-r) specific options:\n"
<< "-m <method> <maxval>, --method <method> <maxval>\n\tcolor visualization method, possible values:\n"
- << "\t\t\t\tpres, vel, angvel, xdisp, angpos\n"
+ << "\tpres, vel, angvel, xdisp, angpos\n"
<< std::endl;
return 0; // Exit with success
}
t@@ -69,7 +69,9 @@ int main(const int argc, const char *argv[])
<< "| |___/ .__/|_| |_|\\___|_| \\___| |\n"
<< "| | | |\n"
<< "| |_| Version: " << VERS << " |\n"
- << "`-------------------------------------´\n";
+ << "`-------------------------------------´\n"
+ << " A discrete element method particle dynamics simulator.\n"
+ << " Written by Anders Damsgaard Christensen, license GPLv3+.\n";
return 0;
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -14,7 +14,8 @@
DEM::DEM(const std::string inputbin,
const int verbosity,
const int checkVals,
- const int dry)
+ const int dry,
+ const int initCuda)
: verbose(verbosity)
{
using std::cout;
t@@ -44,18 +45,20 @@ DEM::DEM(const std::string inputbin,
if (dry == 1)
exit(1);
- // Initialize CUDA
- initializeGPU();
+ if (initCuda == 1) {
+ // Initialize CUDA
+ initializeGPU();
- // Copy constant data to constant device memory
- transferToConstantDeviceMemory();
+ // Copy constant data to constant device memory
+ transferToConstantDeviceMemory();
- // Allocate device memory for particle variables,
- // tied to previously declared pointers in structures
- allocateGlobalDeviceMemory();
+ // Allocate device memory for particle variables,
+ // tied to previously declared pointers in structures
+ allocateGlobalDeviceMemory();
- // Transfer data from host to gpu device memory
- transferToGlobalDeviceMemory();
+ // Transfer data from host to gpu device memory
+ transferToGlobalDeviceMemory();
+ }
}
t@@ -264,4 +267,104 @@ void DEM::reportValues()
cout << " cells\n";
}
+// Returns the volume of a spherical cap
+Float sphericalCap(const Float h, const Float r)
+{
+ return M_PI * h * h / 3.0 * (3.0 * r - h);
+}
+
+// Calculate the porosity with depth, and write to file in output directory
+void DEM::porosity(const int z_slices)
+{
+
+ // Calculate depth slice thickness
+ Float h_slice = (walls.nx->w - grid.origo[2]) / (Float)z_slices;
+
+ // Calculate slice volume
+ Float V_slice = h_slice * grid.L[0] * grid.L[1];
+
+ // Array of porosity values
+ Float porosity[z_slices];
+
+ // Loop over vertical slices
+//#pragma omp parallel for if(np > 100)
+ for (int iz = 0; iz<z_slices; ++iz) {
+
+ // The void volume equals the slice volume, with the
+ // grain volumes subtracted
+ Float V_void = V_slice;
+
+ // Bottom and top position of depth slice
+ Float z_low = iz * h_slice;
+ Float z_high = z_low + h_slice;
+
+ // Loop over particles to see whether they are inside of the slice
+ for (unsigned int i = 0; i<np; ++i) {
+
+ Float z_pos = k.x[i].z;
+ Float radius = k.x[i].w;
+
+ // Sphere volume
+ Float V_sphere = 4.0/3.0 * M_PI * radius * radius * radius;
+
+
+ // If the sphere is inside the slice and not intersecting the
+ // boundaries, subtract the entire sphere volume
+ if (z_low < (z_pos - radius) && (z_pos + radius) < z_high) {
+ V_void -= V_sphere;
+ std::cout << "sphere " << i << " entirely inside volume" << std::endl;
+ } else {
+
+ // If the sphere intersects with the lower boundary,
+ // and the centre is below the boundary
+ if (z_low > z_pos && z_low < (z_pos + radius)) {
+
+ // Subtract the volume of a spherical cap
+ V_void -= sphericalCap(z_pos + radius - z_low, radius);
+
+ std::cout << "sphere " << i << " intersects lower boundary on upper hemisphere" << std::endl;
+ }
+
+ // If the sphere intersects with the lower boundary,
+ // and the centre is above the boundary
+ else if (z_low < z_pos && z_low > (z_pos - radius)) {
+
+ // Subtract the volume of the sphere,
+ // then add the volume of the spherical cap below
+ V_void -= V_sphere + sphericalCap(z_low - z_pos - radius, radius);
+ std::cout << "sphere " << i << " intersects lower boundary on lower hemisphere" << std::endl;
+ }
+
+ // If the sphere intersects with the upper boundary,
+ // and the centre is below the boundary
+ if (z_high > z_pos && z_high < (z_pos + radius)) {
+
+ // Subtract the volume of the sphere,
+ // then add the volume of the spherical cap above
+ V_void -= V_sphere + sphericalCap(z_pos + radius - z_high, radius);
+ std::cout << "sphere " << i << " intersects upper boundary on upper hemisphere" << std::endl;
+ }
+
+ // If the sphere intersects with the upper boundary,
+ // and the centre is above the boundary
+ else if (z_high < z_pos && z_high > (z_pos - radius)) {
+
+ // Subtract the volume of the spherical cap below
+ V_void -= sphericalCap(z_high - z_pos - radius, radius);
+ std::cout << "sphere " << i << " intersects upper boundary on lower hemisphere" << std::endl;
+ }
+
+ }
+ }
+
+ // Save the porosity
+ porosity[iz] = V_void / V_slice;
+ }
+
+ // Display results from the top down
+ for (int i=z_slices-1; i>=0; --i)
+ std::cout << porosity[i] << std::endl;
+
+}
+
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -130,7 +130,8 @@ class DEM {
DEM(std::string inputbin,
const int verbosity = 1,
const int checkVals = 1,
- const int dry = 0);
+ const int dry = 0,
+ const int initCuda = 1);
// Destructor
~DEM(void);
t@@ -162,6 +163,9 @@ class DEM {
// Write image data to PPM file
void writePPM(const char *target);
+ // Calculate porosity with depth and save as text file
+ void porosity(const int z_slices = 10);
+
};
#endif