tDarcy flow seems correct, needs test against analytical solution - 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 b9a5fcbe8fbd1f8109ecb6a0c2a1e98aedc00bc4
(DIR) parent aa36e4b293f63604e6cd32a36ece1742a156fcc4
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Wed, 5 Jun 2013 12:02:09 +0200
Darcy flow seems correct, needs test against analytical solution
Diffstat:
M python/sphere.py | 5 +++++
M src/darcy.cpp | 94 +++++++++++++++++++++++++------
M src/device.cu | 1 -
M src/file_io.cpp | 8 ++++----
M src/porousflow.cpp | 12 +++++-------
M src/sphere.cpp | 2 --
M src/sphere.h | 19 +++++++++++++++----
7 files changed, 107 insertions(+), 34 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -333,10 +333,15 @@ class Spherebin:
for z in range(self.num[2]):
for y in range(self.num[1]):
for x in range(self.num[0]):
+ #print(str(x) + ',' + str(y) + ',' + str(z))
self.f_v[x,y,z,0] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ #print(numpy.fromfile(fh, dtype=numpy.float64, count=1))
self.f_v[x,y,z,1] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ #print(numpy.fromfile(fh, dtype=numpy.float64, count=1))
self.f_v[x,y,z,2] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ #print(numpy.fromfile(fh, dtype=numpy.float64, count=1))
self.f_rho[x,y,z] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ #print(numpy.fromfile(fh, dtype=numpy.float64, count=1))
finally:
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -2,14 +2,13 @@
#include <cstdio>
#include <cstdlib>
#include <string>
+#include <vector>
#include "typedefs.h"
#include "datatypes.h"
#include "constants.h"
#include "sphere.h"
-//#include "eigen-nvcc/Eigen/Core"
-
// Initialize memory
void DEM::initDarcyMem()
{
t@@ -20,6 +19,7 @@ void DEM::initDarcyMem()
d_K = new Float[ncells]; // hydraulic conductivity matrix
d_S = new Float[ncells]; // hydraulic storativity matrix
d_W = new Float[ncells]; // hydraulic recharge
+ d_n = new Float[ncells]; // cell porosity
}
// Free memory
t@@ -31,6 +31,7 @@ void DEM::freeDarcyMem()
free(d_K);
free(d_S);
free(d_W);
+ free(d_n);
}
// 3D index to 1D index
t@@ -187,6 +188,9 @@ void DEM::explDarcyStep()
}
}
}
+
+ // Find macroscopic cell fluid velocities
+ findDarcyVelocities();
}
// Print array values to file stream (stdout, stderr, other file)
t@@ -212,7 +216,7 @@ void DEM::printDarcyArray(FILE* stream, Float* arr, std::string desc)
}
// Print array values to file stream (stdout, stderr, other file)
-void DEM::printDarcyArray3(FILE* stream, Float3* arr)
+void DEM::printDarcyArray(FILE* stream, Float3* arr)
{
unsigned int x, y, z;
for (z=0; z<d_nz; z++) {
t@@ -230,10 +234,10 @@ void DEM::printDarcyArray3(FILE* stream, Float3* arr)
}
// Overload printDarcyArray to add optional description
-void DEM::printDarcyArray3(FILE* stream, Float3* arr, std::string desc)
+void DEM::printDarcyArray(FILE* stream, Float3* arr, std::string desc)
{
std::cout << "\n" << desc << ":\n";
- printDarcyArray3(stream, arr);
+ printDarcyArray(stream, arr);
}
// Find cell velocity
t@@ -247,8 +251,7 @@ void DEM::findDarcyVelocities()
Float nu = params.nu;
// Porosity [-]: n
- Float n = 0.5; // CHANGE THIS
-
+
unsigned int ix, iy, iz, cellidx;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
t@@ -257,22 +260,77 @@ void DEM::findDarcyVelocities()
cellidx = idx(ix,iy,iz);
dH = d_dH[cellidx];
- // Calculate flux
+ // Approximate cell porosity
+ Float n = cellPorosity(ix, iy, iz);
+
+ // Calculate flux (if 0,0,0 is lower left corner)
q.x = -d_K[cellidx]/nu * dH.x;
q.y = -d_K[cellidx]/nu * dH.y;
q.z = -d_K[cellidx]/nu * dH.z;
-
+
// Calculate velocity
v.x = q.x/n;
v.y = q.y/n;
v.z = q.z/n;
-
d_V[cellidx] = v;
}
}
}
}
+// Find particles with centres inside a spatial interval
+// NOTE: This function is untested and unused
+std::vector<unsigned int> DEM::particlesInCell(
+ const Float3 min, const Float3 max)
+{
+ // Particles radii inside cell will be stored in this vector
+ std::vector<unsigned int> pidx;
+
+ unsigned int i;
+ Float4 x;
+ for (i=0; i<np; ++i) {
+
+ // Read the position
+ x = k.x[i];
+
+ if (x.x >= min.x && x.y >= min.y && x.z >= min.z
+ && x.x < max.x && x.y < max.y && x.z < max.z) {
+ pidx.push_back(i);
+ }
+ }
+}
+
+// Find the porosity of a target cell
+Float DEM::cellPorosity(
+ const unsigned int x,
+ const unsigned int y,
+ const unsigned int z)
+{
+ const Float3 x_min = {x*d_dx, y*d_dy, z*d_dz};
+ const Float3 x_max = {x_min.x+d_dx, x_min.y+d_dy, x_min.z+d_dz};
+ const Float cell_volume = d_dx*d_dy*d_dz;
+
+ Float void_volume = cell_volume;
+
+ unsigned int i;
+ Float4 xr;
+ for (i=0; i<np; ++i) {
+
+ // Read the position and radius
+ xr = k.x[i];
+
+ if (xr.x >= x_min.x && xr.y >= x_min.y && xr.z >= x_min.z
+ && xr.x < x_max.x && xr.y < x_max.y && xr.z < x_max.z) {
+ void_volume -= 4.0/3.0*M_PI*xr.w*xr.w*xr.w;
+ }
+ }
+
+ // Return the porosity, which should always be between 0.0 and 1.0
+ Float n = fmin(1.0, fmax(0.0, (void_volume)/cell_volume));
+ d_n[idx(x,y,z)] = n;
+ return n;
+}
+
// Solve Darcy flow on a regular, cubic grid
void DEM::initDarcy(const Float cellsizemultiplier)
{
t@@ -294,12 +352,12 @@ void DEM::initDarcy(const Float cellsizemultiplier)
if (verbose == 1) {
std::cout << " - Fluid grid dimensions: "
- << d_nx << " * "
- << d_ny << " * "
+ << d_nx << "*"
+ << d_ny << "*"
<< d_nz << std::endl;
std::cout << " - Fluid grid cell size: "
- << d_dx << " * "
- << d_dy << " * "
+ << d_dx << "*"
+ << d_dy << "*"
<< d_dz << std::endl;
}
t@@ -307,9 +365,13 @@ void DEM::initDarcy(const Float cellsizemultiplier)
initDarcyVals();
}
-
+// Print final heads and free memory
void DEM::endDarcy()
{
- printDarcyArray(stdout, d_H, "d_H");
+ //printDarcyArray(stdout, d_H, "d_H");
+ //printDarcyArray(stdout, d_V, "d_V");
+ //printDarcyArray(stdout, d_n, "d_n");
freeDarcyMem();
}
+
+// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -860,7 +860,6 @@ __host__ void DEM::startTime()
if (darcy == 1)
explDarcyStep();
-
// Update particle kinematics
if (PROFILING == 1)
startTimer(&kernel_tic);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -246,7 +246,7 @@ void DEM::readbin(const char *target)
if (params.nu > 0.0 && darcy == 0) { // Lattice-Boltzmann flow
if (verbose == 1)
- cout << " - Reading LBM values:\t\t\t";
+ cout << " - Reading LBM values:\t\t\t\t ";
f = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
f_new = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
t@@ -461,9 +461,9 @@ void DEM::writebin(const char *target)
}
}
} else if (params.nu > 0.0 && darcy == 1) { // Darcy flow
- for (z = 0; z<d_dz; ++z) {
- for (y = 0; y<d_dy; ++y) {
- for (x = 0; x<d_dx; ++x) {
+ for (z = 0; z<d_nz; ++z) {
+ for (y = 0; y<d_ny; ++y) {
+ for (x = 0; x<d_nx; ++x) {
i = idx(x,y,z);
ofs.write(as_bytes(d_V[i].x), sizeof(Float));
ofs.write(as_bytes(d_V[i].y), sizeof(Float));
(DIR) diff --git a/src/porousflow.cpp b/src/porousflow.cpp
t@@ -30,7 +30,7 @@ int main(const int argc, const char *argv[])
int verbose = 1;
int dry = 0;
int nfiles = 0;
-
+ int darcy = 1;
// Process input parameters
int i;
t@@ -61,10 +61,9 @@ int main(const int argc, const char *argv[])
// Create DEM class, read data from input binary,
// check values, init cuda, transfer const mem, solve Darcy flow
- DEM dem(argvi, verbose, 1, dry, 1, 1, 1);
-
- dem.startTime();
-
+ DEM dem(argvi, verbose, 1, dry, 1, 1, darcy);
+ if (dry == 0)
+ dem.startTime();
}
}
t@@ -77,6 +76,5 @@ int main(const int argc, const char *argv[])
}
return 0; // Return successfull exit status
-}
-// END OF FILE
+}
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -738,8 +738,6 @@ void DEM::forcechains(const std::string format, const int threedim,
cout << '[' << x_min.y << ':' << x_max.y << "] ";
cout << '[' << x_min.z << ':' << x_max.z << "] " << "NaN notitle" << endl;
}
-
-
}
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -151,12 +151,13 @@ class DEM {
int d_nx, d_ny, d_nz; // Number of cells in each dim
Float d_dx, d_dy, d_dz; // Cell length in each dim
Float* d_H; // Cell hydraulic heads
+ Float3* d_V; // Cell fluid velocity
Float3* d_dH; // Cell spatial gradient in heads
Float* d_K; // Cell hydraulic conductivities (anisotropic)
Float* d_S; // Cell hydraulic storativity
Float* d_W; // Cell hydraulic recharge
- Float3* d_V; // Cell fluid velocity
-
+ Float* d_n; // Cell porosity
+
// Darcy functions
// Memory allocation
t@@ -172,6 +173,16 @@ class DEM {
// Set gradient to zero at grid edges
void setDarcyBCNeumannZero();
+ // Find particles in cell
+ std::vector<unsigned int> particlesInCell(
+ const Float3 min, const Float3 max);
+
+ // Find porosity of cell
+ Float cellPorosity(
+ const unsigned int x,
+ const unsigned int y,
+ const unsigned int z);
+
// Find darcy flow velocities from specific flux (q)
void findDarcyVelocities();
t@@ -266,8 +277,8 @@ class DEM {
// Print Darcy arrays to file stream
void printDarcyArray(FILE* stream, Float* arr);
void printDarcyArray(FILE* stream, Float* arr, std::string desc);
- void printDarcyArray3(FILE* stream, Float3* arr);
- void printDarcyArray3(FILE* stream, Float3* arr, std::string desc);
+ void printDarcyArray(FILE* stream, Float3* arr);
+ void printDarcyArray(FILE* stream, Float3* arr, std::string desc);
};