tDarcy porous flow (CPU) implemented, awaiting 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 4d546f1656814b970cd95e250c2771a5beee42f5
(DIR) parent 9f99cdfce0698a251618671dbb964b98bfe813f4
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Sat, 29 Jun 2013 12:11:33 +0200
Darcy porous flow (CPU) implemented, awaiting tests
Diffstat:
A python/darcy.py | 118 +++++++++++++++++++++++++++++++
M src/darcy.cpp | 83 +++++++++++++++++++++++--------
M src/sphere.h | 48 ++++++++++++++++++++++++++-----
M src/utility.cpp | 59 ++++++++++++++++++++++++++++++-
M src/utility.h | 6 ++++++
5 files changed, 286 insertions(+), 28 deletions(-)
---
(DIR) diff --git a/python/darcy.py b/python/darcy.py
t@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+
+# Import sphere functionality
+from sphere import *
+import sys
+
+### EXPERIMENT SETUP ###
+initialization = True
+consolidation = True
+#shearing = True
+rendering = False
+#plots = False
+
+
+
+# Number of particles
+#np = 1e2
+np = 1e4
+
+# Common simulation id
+sim_id = "darcy"
+
+# Deviatoric stress [Pa]
+#devs = 10e3
+devslist = [10.0e3]
+
+### INITIALIZATION ###
+
+# New class
+init = Spherebin(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
+
+# Save radii
+init.generateRadii(radius_mean = 0.05)
+
+# Use default params
+init.defaultParams(mu_s = 0.4, mu_d = 0.4, nu = 8.9e-4)
+
+# Initialize positions in random grid (also sets world size)
+#init.initRandomGridPos(gridnum = numpy.array([9, 9, 1000]), periodic = 1, contactmodel = 2)
+#init.initRandomGridPos(gridnum = numpy.array([32, 32, 1000]), periodic = 1, contactmodel = 2)
+init.initRandomGridPos(gridnum = numpy.array([32, 32, 1000]), periodic = 1, contactmodel = 1)
+
+# Bond ~30% of the particles
+#init.random2bonds(spacing=0.1)
+
+# Set duration of simulation
+init.initTemporal(total = 7.0)
+init.time_file_dt[0] = 0.05
+#init.time_file_dt[0] = init.time_dt[0]*0.99
+#init.time_total[0] = init.time_dt[0]*2.0
+#init.initTemporal(total = 0.5)
+#init.time_file_dt[0] = init.time_total[0]/5.0
+
+#init.f_rho[2,2,4] = 5.0
+#init.f_rho[6,6,10] = 1.1
+#init.f_rho[:,:,-1] = 1.0001
+
+if (initialization == True):
+
+ # Write input file for sphere
+ init.writebin()
+
+ # Run sphere
+ init.run(dry=True)
+ init.run(darcyflow=True)
+
+
+### CONSOLIDATION ###
+
+for devs in devslist:
+ # New class
+ cons = Spherebin(np = init.np, nw = 1, sid = sim_id + "-cons-devs{}".format(devs))
+
+ # Read last output file of initialization step
+ lastf = status(sim_id + "-init")
+ cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
+
+ # Setup consolidation experiment
+ cons.consolidate(deviatoric_stress = devs, periodic = init.periodic)
+
+
+ # Set duration of simulation
+ cons.initTemporal(total = 1.5)
+ #cons.initTemporal(total = 0.0019, file_dt = 0.00009)
+ #cons.initTemporal(total = 0.0019, file_dt = 1e-6)
+ #cons.initTemporal(total = 0.19, file_dt = 0.019)
+
+ cons.w_m[0] *= 0.001
+
+
+
+ if (consolidation == True):
+ # Write input file for sphere
+ cons.writebin()
+
+ # Run sphere
+ cons.run(dry=True) # show values, don't run
+ cons.run(darcyflow=True) # run
+
+ if (plots == True):
+ # Make a graph of energies
+ visualize(cons.sid, "energy", savefig=True, outformat='png')
+ visualize(cons.sid, "walls", savefig=True, outformat='png')
+
+ if (rendering == True):
+ # Render images with raytracer
+ cons.render(method = "pres", max_val = 2.0*devs, verbose = False)
+
+ project = cons.sid
+ lastfile = status(cons.sid)
+ sb = Spherebin()
+ for i in range(lastfile+1):
+ fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
+ sb.sid = project + ".output{:0=5}".format(i)
+ sb.readbin(fn, verbose = False)
+ for y in range(0,sb.num[1]):
+ sb.plotFluidDensities(y = y)
+ sb.plotFluidVelocities(y = y)
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -23,7 +23,7 @@ void DEM::initDarcyMem()
d_T = new Float3[ncells]; // hydraulic transmissivity matrix
d_Ss = new Float[ncells]; // hydraulic storativity matrix
d_W = new Float[ncells]; // hydraulic recharge
- d_n = new Float[ncells]; // cell porosity
+ d_phi = new Float[ncells]; // cell porosity
}
// Free memory
t@@ -37,7 +37,7 @@ void DEM::freeDarcyMem()
free(d_T);
free(d_Ss);
free(d_W);
- free(d_n);
+ free(d_phi);
}
// 3D index to 1D index
t@@ -57,7 +57,7 @@ void DEM::initDarcyVals()
const Float k = 1.0e-10;
// Density of the fluid [kg/m^3]
- const Float rho = 3600.0;
+ const Float rho = 1000.0;
unsigned int ix, iy, iz, cellidx;
for (ix=0; ix<d_nx; ++ix) {
t@@ -95,7 +95,7 @@ void DEM::copyDarcyVals(unsigned int read, unsigned int write)
d_T[write] = MAKE_FLOAT3(d_T[read].x, d_T[read].y, d_T[read].z);
d_Ss[write] = d_Ss[read];
d_W[write] = d_W[read];
- d_n[write] = d_n[read];
+ d_phi[write] = d_phi[read];
}
// Update ghost nodes from their parent cell values
t@@ -156,14 +156,32 @@ void DEM::setDarcyGhostNodes()
// Find cell transmissivities from hydraulic conductivities and cell dimensions
void DEM::findDarcyTransmissivities()
{
+ // Find porosities from cell particle content
+ findPorosities();
+
+ // Density of the fluid [kg/m^3]
+ const Float rho = 1000.0;
+
+ // Kozeny-Carman parameter
+ Float a = 1.0e-8;
+
unsigned int ix, iy, iz, cellidx;
- Float K;
+ Float K, k;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
cellidx = idx(ix,iy,iz);
- K = d_K[cellidx];
+
+ // Read cell porosity
+ Float phi = d_phi[cellidx];
+
+ // Calculate permeability from the Kozeny-Carman relationship
+ k = a*phi*phi*phi/(1.0 - phi*phi);
+
+ //K = d_K[cellidx];
+ // Save hydraulic conductivity [m/s]
+ d_K[cellidx] = k*rho*-params.g[2]/params.nu;
// Hydraulic transmissivity [m2/s]
Float3 T = {K*d_dx, K*d_dy, K*d_dz};
t@@ -293,6 +311,12 @@ Float hmean(Float a, Float b) {
// Boundary conditions are fixed values (Dirichlet)
void DEM::explDarcyStep()
{
+ // Find transmissivities from cell particle content
+ findDarcyTransmissivities();
+
+ // Check the time step length
+ checkDarcyTimestep();
+
// Cell dims squared
const Float dx2 = d_dx*d_dx;
const Float dy2 = d_dy*d_dy;
t@@ -404,14 +428,13 @@ void DEM::explDarcyStep()
}
}
- // Find macroscopic cell fluid velocities
- findDarcyVelocities();
-
// Swap d_H and d_H_new
Float* tmp = d_H;
d_H = d_H_new;
d_H_new = tmp;
+ // Find macroscopic cell fluid velocities
+ findDarcyVelocities();
}
// Print array values to file stream (stdout, stderr, other file)
t@@ -485,7 +508,7 @@ void DEM::findDarcyVelocities()
dH = d_dH[cellidx];
// Approximate cell porosity
- Float n = cellPorosity(ix, iy, iz);
+ Float phi = d_phi[cellidx];
// Calculate flux
// The sign might need to be reversed, depending on the
t@@ -495,9 +518,9 @@ void DEM::findDarcyVelocities()
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;
+ v.x = q.x/phi;
+ v.y = q.y/phi;
+ v.z = q.z/phi;
d_V[cellidx] = v;
}
}
t@@ -556,9 +579,21 @@ Float DEM::cellPorosity(
}
// 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;
+ Float phi = fmin(1.0, fmax(0.0, void_volume/cell_volume));
+ //Float phi = 0.1;
+ return phi;
+}
+
+void DEM::findPorosities()
+{
+ unsigned int ix, iy, iz, cellidx;
+ for (ix=0; ix<d_nx; ++ix) {
+ for (iy=0; iy<d_ny; ++iy) {
+ for (iz=0; iz<d_nz; ++iz) {
+ d_phi[idx(ix,iy,iz)] = cellPorosity(ix,iy,iz);
+ }
+ }
+ }
}
// Find particles with centres inside a spatial interval
t@@ -619,7 +654,7 @@ Float DEM::getTmax()
return max;
}
// Get maximum value in 1d array with ghost nodes
-Float DEM::getSmin()
+Float DEM::getSsmin()
{
Float min = 1.0e13; // initialize with a small number
unsigned int ix,iy,iz;
t@@ -704,9 +739,17 @@ void DEM::initDarcy(const Float cellsizemultiplier)
// Print final heads and free memory
void DEM::endDarcy()
{
- //printDarcyArray(stdout, d_H, "d_H");
- //printDarcyArray(stdout, d_V, "d_V");
- //printDarcyArray(stdout, d_n, "d_n");
+ FILE* Kfile;
+ if ((Kfile = fopen("d_K.txt","w"))) {
+ printDarcyArray(Kfile, d_K);
+ fclose(Kfile);
+ } else {
+ fprintf(stderr, "Error, could not open d_K.txt\n");
+ }
+ printDarcyArray(stdout, d_phi, "d_phi");
+ printDarcyArray(stdout, d_H, "d_H");
+ printDarcyArray(stdout, d_K, "d_K");
+ printDarcyArray(stdout, d_V, "d_V");
freeDarcyMem();
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -147,7 +147,7 @@ class DEM {
//// Darcy-flow
int darcy; // 0: no, 1: yes
- // Darcy values
+ // Darcy values, host
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
t@@ -155,11 +155,23 @@ class DEM {
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
+ Float3* d_T; // Cell hydraulic transmissivity
+ Float* d_Ss; // Cell hydraulic storativity
Float* d_W; // Cell hydraulic recharge
- Float* d_n; // Cell porosity
+ Float* d_phi; // Cell porosity
- // Darcy functions
+ // Darcy values, device
+ Float* dev_d_H; // Cell hydraulic heads
+ Float* dev_d_H_new; // Cell hydraulic heads
+ Float3* dev_d_V; // Cell fluid velocity
+ Float3* dev_d_dH; // Cell spatial gradient in heads
+ Float* dev_d_K; // Cell hydraulic conductivities
+ Float3* dev_d_T; // Cell hydraulic transmissivity
+ Float* dev_d_Ss; // Cell hydraulic storativity
+ Float* dev_d_W; // Cell hydraulic recharge
+ Float* dev_d_phi; // Cell porosity
+
+ //// Darcy functions
// Memory allocation
void initDarcyMem();
t@@ -168,6 +180,16 @@ class DEM {
// Set some values for the Darcy parameters
void initDarcyVals();
+ // Copy Darcy values from cell to cell (by index)
+ void copyDarcyVals(unsigned int read, unsigned int write);
+
+ // Update ghost nodes from their parent cell values
+ void setDarcyGhostNodes();
+
+ // Find cell transmissivities from hydraulic conductivities and cell
+ // dimensions
+ void findDarcyTransmissivities();
+
// Finds central difference gradients
void findDarcyGradients();
t@@ -202,6 +224,9 @@ class DEM {
const unsigned int y,
const unsigned int z);
+ // Find and save all cell porosities
+ void findPorosities();
+
// Find darcy flow velocities from specific flux (q)
void findDarcyVelocities();
t@@ -211,18 +236,27 @@ class DEM {
const unsigned int y,
const unsigned int z);
- // Get minimum value in 1D array
- Float minVal3dArr(Float* arr);
-
// Initialize Darcy values and arrays
void initDarcy(const Float cellsizemultiplier = 1.0);
// Clean up Darcy arrays
void endDarcy();
+ // Check whether the explicit integration is going to meet the
+ // stability criteria
+ Float getTmax();
+ Float getSsmin();
+ void checkDarcyTimestep();
+
// Perform a single time step, explicit integration
void explDarcyStep();
+ //// Darcy functions, device
+ void initDarcyMemDev();
+ void freeDarcyMemDev();
+ void transferDarcyToGlobalDeviceMemory(int statusmsg);
+ void transferDarcyFromGlobalDeviceMemory(int statusmsg);
+
public:
// Values and functions accessible from the outside
(DIR) diff --git a/src/utility.cpp b/src/utility.cpp
t@@ -4,7 +4,8 @@
//Round a / b to nearest higher integer value
-unsigned int iDivUp(unsigned int a, unsigned int b) {
+unsigned int iDivUp(unsigned int a, unsigned int b)
+{
return (a % b != 0) ? (a / b + 1) : (a / b);
}
t@@ -15,3 +16,59 @@ void swapFloatArrays(Float* arr1, Float* arr2)
arr1 = arr2;
arr2 = tmp;
}
+
+// Get minimum value in 1D array
+Float minVal(Float* arr, int length)
+{
+ Float min = 1.0e13; // initialize with a large number
+ unsigned int i;
+ Float val;
+ for (i=0; i<length; ++i) {
+ val = arr[i];
+ if (val < min) min = val;
+ }
+ return min;
+}
+
+// Get maximum value in 1d array
+Float maxVal(Float* arr, int length)
+{
+ Float max = -1.0e13; // initialize with a small number
+ unsigned int i;
+ Float val;
+ for (i=0; i<length; ++i) {
+ val = arr[i];
+ if (val > max) max = val;
+ }
+ return max;
+}
+
+// Get minimum value in 3d array
+Float minVal(Float3* arr, int length)
+{
+ Float min = 1.0e13; // initialize with a large number
+ unsigned int i;
+ Float3 val;
+ for (i=0; i<length; ++i) {
+ val = arr[i];
+ if (val.x < min) min = val.x;
+ if (val.y < min) min = val.y;
+ if (val.z < min) min = val.z;
+ }
+ return min;
+}
+
+// Get maximum value in 3d array
+Float maxVal(Float3* arr, int length)
+{
+ Float max = -1.0e13; // initialize with a small number
+ unsigned int i;
+ Float3 val;
+ for (i=0; i<length; ++i) {
+ val = arr[i];
+ if (val.x > max) max = val.x;
+ if (val.y > max) max = val.y;
+ if (val.z > max) max = val.z;
+ }
+ return max;
+}
(DIR) diff --git a/src/utility.h b/src/utility.h
t@@ -10,5 +10,11 @@ unsigned int iDivUp(unsigned int a, unsigned int b);
// Swap two arrays pointers
void swapFloatArrays(Float* arr1, Float* arr2);
+// Get minimum/maximum value in 1D or 3D array
+Float minVal(Float* arr, int length);
+Float minVal(Float3* arr, int length);
+Float maxVal(Float* arr, int length);
+Float maxVal(Float3* arr, int length);
+
#endif
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4