tadd advection stability check - 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 9c151584116fc9b0a5b998903775373d1b7db9b0
(DIR) parent a911141e8df42ba8a6b2be040a989889cbd93a96
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 6 Nov 2014 12:36:06 +0100
add advection stability check
Diffstat:
M python/sphere.py | 10 ++++++----
M src/darcy.cpp | 40 ++++++++++++++++++++++++++++---
M src/sphere.h | 1 +
3 files changed, 44 insertions(+), 7 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -2785,9 +2785,10 @@ class sim:
self.cellSize()
#self.hydraulicPermeability()
#alpha_max = numpy.max(self.k/(self.beta_f*0.9*self.mu))
- k_max = 2.7e-8 # hardcoded
- phi_min = 0.1 # hardcoded
+ k_max = 2.7e-8 # hardcoded in darcy.cuh
+ phi_min = 0.1 # hardcoded in darcy.cuh
alpha_max = k_max/(self.beta_f*phi_min*self.mu)
+ print(alpha_max)
return safety * 1.0/(2.0*alpha_max)*1.0/(
1.0/(self.dx[0]**2) + \
1.0/(self.dx[1]**2) + \
t@@ -2907,8 +2908,9 @@ class sim:
# Check numerical stability of the fluid phase, by criteria derived by
# von Neumann stability analysis of the diffusion and advection terms
- elif self.fluid:
- self.time_dt[0] = self.largestFluidTimeStep(safety = 0.5)
+ if self.fluid:
+ fluid_time_dt = self.largestFluidTimeStep(safety = 0.5)
+ self.time_dt[0] = numpy.min([fluid_time_dt, self.time_dt[0]])
else:
raise Exception('Error: Could not automatically set a time step.')
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -109,6 +109,25 @@ Float DEM::smallestDarcyPorosity()
return phi_min;
}
+// Component-wise max absolute velocities
+Float3 DEM::largestDarcyVelocities()
+{
+ Float3 v_max_abs = MAKE_FLOAT3(0.0, 0.0, 0.0);
+ Float3 v;
+ for (unsigned int z=0; z<grid.num[2]; z++)
+ for (unsigned int y=0; y<grid.num[1]; y++)
+ for (unsigned int x=0; x<grid.num[0]; x++) {
+ v = darcy.v[d_idx(x,y,z)];
+ if (v.x > v_max_abs.x)
+ v_max_abs.x = fabs(v.x);
+ if (v.y > v_max_abs.y)
+ v_max_abs.y = fabs(v.y);
+ if (v.z > v_max_abs.z)
+ v_max_abs.z = fabs(v.z);
+ }
+ return v_max_abs;
+}
+
// Determine if the FTCS (forward time, central space) solution of the Navier
// Stokes equations is unstable
void DEM::checkDarcyStability()
t@@ -121,17 +140,32 @@ void DEM::checkDarcyStability()
const Float alpha_max = largestDarcyPermeability()
/(darcy.beta_f*smallestDarcyPorosity()*darcy.mu);
- // Check the diffusion term using von Neumann stability analysis
- //if (darcy.mu*time.dt/(dmin*dmin) > 0.5) {
+ // von Neumann stability analysis
if (time.dt >= 1.0/(2.0*alpha_max) *
1.0/(1.0/(dx*dx) + 1.0/(dy*dy) + 1.0/(dz*dz))) {
std::cerr
<< "\nError: The time step is too large to ensure stability in "
"the diffusive term of the fluid momentum equation.\n"
- "Increase the viscosity, decrease the time step, and/or increase "
+ "Decrease the time step, increase the fluid viscosity, increase "
+ "the fluid compressibility and/or increase "
"the fluid grid cell size." << std::endl;
//exit(1);
}
+
+ // Courant-Friedrichs-Lewy criteria
+ Float3 v_max_abs = largestDarcyVelocities();
+ if (v_max_abs.x*time.dt > dx ||
+ v_max_abs.y*time.dt > dy ||
+ v_max_abs.z*time.dt > dz) {
+ std::cerr
+ << "\nError: The time step is too large to ensure stability due to "
+ "large fluid velocities.\n v_max_abs = "
+ << v_max_abs.x << ", "
+ << v_max_abs.y << ", "
+ << v_max_abs.z <<
+ " m/s.\nDecrease the time step "
+ "and/or increase the fluid grid cell size." << std::endl;
+ }
}
// Print array values to file stream (stdout, stderr, other file)
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -310,6 +310,7 @@ class DEM {
void initDarcyMem();
Float largestDarcyPermeability();
Float smallestDarcyPorosity();
+ Float3 largestDarcyVelocities();
void initDarcyMemDev();
unsigned int darcyCells();
unsigned int darcyCellsVelocity();