tFixed porosity calculator and added porosity tests - 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 d6e187230caddd28d3aeac6ebd09175c1067fa40
(DIR) parent b3a14695511b1c825ed7bd11f7ea0e0f07a3da98
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Wed, 2 Jan 2013 09:57:32 +0100
Fixed porosity calculator and added porosity tests
Diffstat:
M python/tests.py | 83 ++++++++++++++++++++++++++++---
M src/main.cpp | 15 +++++++++++----
M src/sphere.cpp | 45 +++++++++++++++++++++++--------
M src/sphere.h | 3 +++
4 files changed, 123 insertions(+), 23 deletions(-)
---
(DIR) diff --git a/python/tests.py b/python/tests.py
t@@ -2,18 +2,36 @@
from sphere import *
import subprocess
+def passed():
+ return "\tPassed"
+
+def failed():
+ return "\tFailed"
+
def compare(first, second, string):
if (first == second):
- print(string + "\tPassed")
+ print(string + passed())
else:
- print(string + "\tFailed")
+ print(string + failed())
+
+def compareFloats(first, second, string, criterion=1e-5):
+ if abs(first-second) < criterion:
+ print(string + passed())
+ else :
+ print(string + failed())
+
+def cleanup(spherebin):
+ 'Remove temporary files'
+ subprocess.call("rm -f ../input/" + spherebin.sid + ".bin", shell=True)
+ subprocess.call("rm -f ../output/" + spherebin.sid + ".*.bin", shell=True)
+ print("")
#### Input/output tests ####
print("### Input/output tests ###")
# Generate data in python
-orig = Spherebin(np = 100, nw = 1, sid = "test")
+orig = Spherebin(np = 100, nw = 1, sid = "test-initgrid")
orig.generateRadii(histogram = False)
orig.defaultParams()
orig.initRandomGridPos(g = numpy.zeros(orig.nd))
t@@ -24,24 +42,73 @@ orig.writebin(verbose=False)
# Test Python IO routines
py = Spherebin()
-py.readbin("../input/test.bin", verbose=False)
+py.readbin("../input/" + orig.sid + ".bin", verbose=False)
compare(orig, py, "Python IO:")
# Test C++ IO routines
orig.run(verbose=False, hideinputfile=True)
#orig.run()
cpp = Spherebin()
-cpp.readbin("../output/test.output00000.bin", verbose=False)
+cpp.readbin("../output/" + orig.sid + ".output00000.bin", verbose=False)
compare(orig, cpp, "C++ IO: ")
# Test CUDA IO routines
cuda = Spherebin()
-cuda.readbin("../output/test.output00001.bin", verbose=False)
+cuda.readbin("../output/" + orig.sid + ".output00001.bin", verbose=False)
cuda.time_current = orig.time_current
cuda.time_step_count = orig.time_step_count
compare(orig, cuda, "CUDA IO: ")
# Remove temporary files
-subprocess.call("rm ../input/" + orig.sid + ".bin", shell=True)
-subprocess.call("rm ../output/" + orig.sid + ".*.bin", shell=True)
+cleanup(orig)
+
+
+#### Porosity tests ####
+print("### porosity tests ###")
+
+def testPorosities(spherebin):
+
+ # Number of vertical slices
+ slicevals = [1, 2, 4]
+ i = 1 # iterator var
+ for slices in slicevals:
+
+ # Find correct value of bulk porosity
+ n_bulk = spherebin.bulkPorosity()
+ #print("Bulk: " + str(n_bulk))
+
+ porosity = spherebin.porosity(slices = slices)[0]
+ #print("Avg: " + str(numpy.average(porosity)))
+ #print(porosity)
+
+ # Check if average of porosity function values matches the bulk porosity
+ compareFloats(n_bulk, numpy.average(porosity), \
+ spherebin.sid + ": Porosity average to bulk porosity ("\
+ + str(i) + "/" + str(len(slicevals)) + "):")
+ i += 1
+
+# Test data from previous test
+testPorosities(orig)
+
+# Simple cubic packing of uniform spheres
+# The theoretical porosity is (4/3*pi*r^3)/(2r)^3 = 0.476
+sidelen = 10
+cubic = Spherebin(np = sidelen**3, sid='cubic')
+radius = 1.0
+cubic.generateRadii(psd='uni', radius_mean=radius, radius_variance=0.0, histogram=False)
+for ix in range(sidelen):
+ for iy in range(sidelen):
+ for iz in range(sidelen):
+ i = ix + sidelen * (iy + sidelen * iz) # linear index
+ cubic.x[i,0] = ix*radius*2.0 + radius
+ cubic.x[i,1] = iy*radius*2.0 + radius
+ cubic.x[i,2] = iz*radius*2.0 + radius
+cubic.L[:] = 2.0 * radius * sidelen
+
+cubic.initTemporal(0.2)
+cubic.initGrid()
+
+testPorosities(cubic)
+
+cleanup(cubic)
(DIR) diff --git a/src/main.cpp b/src/main.cpp
t@@ -53,7 +53,8 @@ int main(const int argc, const char *argv[])
<< "-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"
- << "\tpres, vel, angvel, xdisp, angpos\n"
+ << "\tnormal, pres, vel, angvel, xdisp, angpos\n"
+ << "\t'normal' is the default mode\n"
<< std::endl;
return 0; // Exit with success
}
t@@ -92,7 +93,9 @@ int main(const int argc, const char *argv[])
else if (argvi == "-m" || argvi == "--method") {
// Find out which
- if (std::string(argv[i+1]) == "pres")
+ if (std::string(argv[i+1]) == "normal")
+ method = 0;
+ else if (std::string(argv[i+1]) == "pres")
method = 1;
else if (std::string(argv[i+1]) == "vel")
method = 2;
t@@ -109,8 +112,12 @@ int main(const int argc, const char *argv[])
}
// Read max. value of colorbar as next argument
- max_val = atof(argv[i+2]);
- i += 2; // skip ahead
+ if (method != 0) {
+ max_val = atof(argv[i+2]);
+ i += 2; // skip ahead
+ } else {
+ i += 1;
+ }
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -306,6 +306,19 @@ Float sphericalCap(const Float h, const Float r)
return M_PI * h * h / 3.0 * (3.0 * r - h);
}
+// Returns the max. radius of any particle
+Float DEM::r_max()
+{
+ Float r_max = 0.0;
+ Float r;
+ for (unsigned int i=0; i<np; ++i) {
+ r = k.x[i].w;
+ if (r > r_max)
+ r_max = r;
+ }
+ return r;
+}
+
// Calculate the porosity with depth, and write to file in output directory
void DEM::porosity(const int z_slices)
{
t@@ -321,6 +334,15 @@ void DEM::porosity(const int z_slices)
// Calculate depth slice thickness
Float h_slice = (top - grid.origo[2]) / (Float)z_slices;
+ // Check that the vertical slice height does not exceed the
+ // max particle diameter, since this function doesn't work
+ // if the sphere intersects more than 1 boundary
+ if (h_slice <= r_max()*2.0) {
+ std::cerr << "Error! The number of z-slices is too high."
+ << std::endl;
+ exit(1);
+ }
+
// Calculate slice volume
Float V_slice = h_slice * grid.L[0] * grid.L[1];
t@@ -331,7 +353,7 @@ void DEM::porosity(const int z_slices)
Float porosity[z_slices];
// Loop over vertical slices
-#pragma omp parallel for if(np > 100)
+//#pragma omp parallel for if(np > 100)
for (int iz = 0; iz<z_slices; ++iz) {
// The void volume equals the slice volume, with the
t@@ -358,14 +380,14 @@ void DEM::porosity(const int z_slices)
// If the sphere is inside the slice and not intersecting the
// boundaries, subtract the entire sphere volume
- if (z_slice_low < z_sphere_low && z_sphere_high < z_slice_high) {
+ if (z_slice_low <= z_sphere_low && z_sphere_high <= z_slice_high) {
V_void -= V_sphere;
} else {
// If the sphere intersects with the lower boundary,
// and the centre is below the boundary
- if (z_slice_low > z_sphere_centre && z_slice_low < z_sphere_high) {
+ if (z_slice_low >= z_sphere_centre && z_slice_low <= z_sphere_high) {
// Subtract the volume of a spherical cap
V_void -= sphericalCap(z_sphere_high - z_slice_low, radius);
t@@ -381,24 +403,25 @@ void DEM::porosity(const int z_slices)
}
// If the sphere intersects with the upper boundary,
+ // and the centre is above the boundary
+ if (z_slice_high <= z_sphere_centre && z_slice_high >= z_sphere_low) {
+
+ // Subtract the volume of the spherical cap below
+ V_void -= sphericalCap(z_slice_high - z_sphere_low, radius);
+ }
+ // If the sphere intersects with the upper boundary,
// and the centre is below the boundary
- if (z_slice_high > z_sphere_centre && z_slice_high < z_sphere_high) {
+ else if (z_slice_high > z_sphere_centre && z_slice_high < z_sphere_high) {
// Subtract the volume of the sphere,
// then add the volume of the spherical cap above
V_void -= V_sphere + sphericalCap(z_sphere_high - z_slice_high, radius);
}
- // If the sphere intersects with the upper boundary,
- // and the centre is above the boundary
- else if (z_slice_high < z_sphere_centre && z_slice_high > z_sphere_low) {
-
- // Subtract the volume of the spherical cap below
- V_void -= sphericalCap(z_slice_high - z_sphere_low, radius);
- }
}
+
}
// Save the mid z-point
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -122,6 +122,9 @@ class DEM {
// Find and return the max. position of any particle in each dimension
float3 maxPos(void);
+ // Find and return the max. radius
+ Float r_max(void);
+
// Write porosities found in porosity() to text file
void writePorosities(
const char *target,