tfixed minor error in forcing function, added optional output - 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 53bf17b35e77b163aa6ece5eaf32ca46541bb272
(DIR) parent ea7d8aa29788bcbca20557265e50c4247dfc4a26
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 6 May 2014 13:44:53 +0200
fixed minor error in forcing function, added optional output
Diffstat:
M src/device.cu | 8 ++++----
M src/navierstokes.cuh | 31 +++++++++++++++++++++++--------
M tests/stokes_law.py | 2 +-
3 files changed, 28 insertions(+), 13 deletions(-)
---
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -675,10 +675,8 @@ __host__ void DEM::startTime()
// MAIN CALCULATION TIME LOOP
while (time.current <= time.total) {
-
// Print current step number to terminal
- //printf("Step: %d\n", time.step_count);
-
+ //printf("\n\n@@@ DEM time step: %ld\n", iter);
// Routine check for errors
checkForCudaErrors("Start of main while loop");
t@@ -690,7 +688,7 @@ __host__ void DEM::startTime()
// in the fine, uniform and homogenous grid.
if (PROFILING == 1)
startTimer(&kernel_tic);
- calcParticleCellID<<<dimGrid, dimBlock>>>(dev_gridParticleCellID,
+ calcParticleCellID<<<dimGrid, dimBlock>>>(dev_gridParticleCellID,
dev_gridParticleIndex,
dev_x);
t@@ -1148,6 +1146,8 @@ __host__ void DEM::startTime()
for (unsigned int nijac = 0; nijac<ns.maxiter; ++nijac) {
+ //printf("### Jacobi iteration %d\n", nijac);
+
// Only grad(epsilon) changes during the Jacobi iterations.
// The remaining terms of the forcing function are only
// calculated during the first iteration.
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -6,13 +6,16 @@
//#include <cutil_math.h>
#include <helper_math.h>
-#include "vector_arithmetic.h" // for arbitrary prec. vectors
+#include "vector_arithmetic.h" // for arbitrary precision vectors
#include "sphere.h"
#include "datatypes.h"
#include "utility.h"
#include "constants.cuh"
#include "debug.h"
+// Enable reporting of forcing function terms to stdout
+//#define REPORT_FORCING_TERMS
+
// Arithmetic mean of two numbers
__inline__ __device__ Float amean(Float a, Float b) {
return (a+b)*0.5;
t@@ -2048,6 +2051,10 @@ __global__ void findNSforcing(
Float f1;
Float3 f2;
+#ifdef REPORT_FORCING_TERMS
+ Float t1, t2, t3, t4;
+#endif
+
// Check if this is the first Jacobi iteration. If it is, find f1 and f2
if (nijac == 0) {
t@@ -2069,16 +2076,19 @@ __global__ void findNSforcing(
+ dot(grad_phi, v_p)*devC_params.rho_f/(devC_dt*phi)
+ dphi*devC_params.rho_f/(devC_dt*devC_dt*phi);
f2 = grad_phi/phi;*/
- f1 = div_v_p*devC_params.rho_f*phi/(ndem*devC_dt)
- + dot(grad_phi, v_p)*devC_params.rho_f/(ndem*devC_dt)
- + dphi*devC_params.rho_f/(ndem*devC_dt*devC_dt);
+ const Float dt = devC_dt*ndem;
+ f1 = div_v_p*devC_params.rho_f*phi/dt
+ + dot(grad_phi, v_p)*devC_params.rho_f/dt
+ + dphi*devC_params.rho_f/(dt*dt);
f2 = grad_phi/phi;
+#ifdef REPORT_FORCING_TERMS
// Report values terms in the forcing function for debugging
+ t1 = div_v_p*phi*devC_params.rho_f/dt;
+ t2 = dot(grad_phi, v_p)*devC_params.rho_f/dt;
+ t4 = dphi*devC_params.rho_f/(dt*dt);
+#endif
/*
- const Float f1t1 = div_v_p*devC_params.rho_f/devC_dt;
- const Float f1t2 = dot(grad_phi, v_p)*devC_params.rho_f/(devC_dt*phi);
- const Float f1t3 = dphi*devC_params.rho_f/(devC_dt*devC_dt*phi);
printf("[%d,%d,%d] f1 = %f\t"
"f1t1 = %f\tf1t2 = %f\tf1t3 = %f\tf2 = %f\n",
x,y,z, f1, f1t1, f1t2, f1t3, f2);
t@@ -2122,7 +2132,12 @@ __global__ void findNSforcing(
// Forcing function value
const Float f = f1 - dot(f2, grad_epsilon);
- //printf("[%d,%d,%d]\tf1 = %f\tf2 = %f\tf = %f\n", x,y,z, f1, f2, f);
+
+#ifdef REPORT_FORCING_TERMS
+ t3 = -dot(f2, grad_epsilon);
+ printf("[%d,%d,%d]\tt1 = %f\tt2 = %f\tt3 = %f\tt4 = %f\n",
+ x,y,z, t1, t2, t3, t4);
+#endif
// Save forcing function value
__syncthreads();
(DIR) diff --git a/tests/stokes_law.py b/tests/stokes_law.py
t@@ -8,7 +8,7 @@ import matplotlib.pyplot as plt
print("### Stokes test - single sphere sedimentation ###")
## Stokes drag
-orig = sphere.sim(sid = "stokes_law_DEMCFD1", fluid = True)
+orig = sphere.sim(sid = "stokes_law_gradP", fluid = True)
cleanup(orig)
orig.defaultParams()
orig.addParticle([0.5,0.5,1.46], 0.05)