tcorrected error in forcing function and velocity correction - sphere - GPU-bas… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit ced051a8ec4a7e3e2efee95b11c1e8b5aece3a73 | |
parent 1e3272182b4e8a3a6dbcb4988721d67853ef0281 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 22 Apr 2014 15:33:29 +0200 | |
corrected error in forcing function and velocity correction | |
Diffstat: | |
A doc/html/_images/math/177b8a21bb3d… | 0 | |
A doc/html/_images/math/22176d770731… | 0 | |
A doc/html/_images/math/470d00b97569… | 0 | |
A doc/html/_images/math/52f771ea3758… | 0 | |
A doc/html/_images/math/7b1a510ad1f3… | 0 | |
A doc/html/_images/math/964fa43eb7c2… | 0 | |
A doc/html/_images/math/e0ee7c8be01c… | 0 | |
M doc/html/_sources/cfd.txt | 27 ++++++++++++++------------- | |
M doc/html/cfd.html | 39 ++++++++++++++++-------------… | |
M doc/html/objects.inv | 0 | |
M doc/html/python_api.html | 4 ++-- | |
M doc/html/searchindex.js | 4 ++-- | |
M doc/pdf/sphere.pdf | 0 | |
M doc/sphinx/cfd.rst | 27 ++++++++++++++------------- | |
M src/device.cu | 1 + | |
M src/navierstokes.cuh | 21 +++++++++------------ | |
16 files changed, 62 insertions(+), 61 deletions(-) | |
--- | |
diff --git a/doc/html/_images/math/177b8a21bb3d0e2a33950cd8ceec4f96f5950bd6.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/22176d7707317572382e54a12e13095b61d90c26.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/470d00b97569e3a8976841debcd99aaf9c78d7a8.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/52f771ea3758e7653af6c471262bbc0ff0d9cf5a.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/7b1a510ad1f355daa1f756cfd47f75a3a29f5660.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/964fa43eb7c25fa23b82daf41b1b9e16e998ab55.png… | |
Binary files differ. | |
diff --git a/doc/html/_images/math/e0ee7c8be01c715a4e719394a4086ba228cbf51b.png… | |
Binary files differ. | |
diff --git a/doc/html/_sources/cfd.txt b/doc/html/_sources/cfd.txt | |
t@@ -416,7 +416,8 @@ The predicted velocity is corrected using the new pressure… | |
.. math:: | |
\boldsymbol{v}^{t+\Delta t} = \boldsymbol{v}^* | |
- - \frac{\Delta t}{\rho} \nabla \epsilon | |
+ %- \frac{\Delta t}{\rho} \nabla \epsilon | |
+ - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon | |
\quad \text{where} \quad | |
\epsilon = p^{t+\Delta t} - \beta p^t | |
t@@ -426,25 +427,25 @@ equation: | |
.. math:: | |
\Rightarrow | |
\phi^t \nabla \cdot | |
- \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+ \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \ri… | |
+ | |
- \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+ \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \ri… | |
\cdot \nabla \phi^t + \frac{\Delta \phi^t}{\Delta t} = 0 | |
.. math:: | |
\Rightarrow | |
\phi^t \nabla \cdot | |
- \boldsymbol{v}^* - \frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+ \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \phi^t \nabla^2 \epsilon | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
- - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+ - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t} = 0 | |
.. math:: | |
\Rightarrow | |
- \frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+ \frac{\Delta t}{\rho} \nabla^2 \epsilon | |
= \phi^t \nabla \cdot \boldsymbol{v}^* | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
- - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+ - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t} | |
The pressure difference in time becomes a `Poisson equation`_ with added terms: | |
t@@ -452,19 +453,19 @@ The pressure difference in time becomes a `Poisson equat… | |
.. math:: | |
\Rightarrow | |
\nabla^2 \epsilon | |
- = \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
+ = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
+ + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t} | |
- + \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t} | |
+ + \frac{\Delta \phi^t \rho}{\Delta t^2} | |
The right hand side of the above equation is termed the *forcing function* | |
:math:`f`, which is decomposed into two terms, :math:`f_1` and :math:`f_2`: | |
.. math:: | |
f_1 | |
- = \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
- + \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t} | |
+ = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
+ + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
+ + \frac{\Delta \phi^t \rho}{\Delta t^2} | |
f_2 = | |
\frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t} | |
diff --git a/doc/html/cfd.html b/doc/html/cfd.html | |
t@@ -406,48 +406,49 @@ equation:</p> | |
</div><p>The predicted velocity is corrected using the new pressure (Langtange… | |
2002):</p> | |
<div class="math"> | |
-<p><img src="_images/math/e3de9b8fce0b12c8459947bd22c727e379512fcc.png" alt="\… | |
-- \frac{\Delta t}{\rho} \nabla \epsilon | |
+<p><img src="_images/math/52f771ea3758e7653af6c471262bbc0ff0d9cf5a.png" alt="\… | |
+%- \frac{\Delta t}{\rho} \nabla \epsilon | |
+- \frac{\Delta t}{\rho \phi^t} \nabla \epsilon | |
\quad \text{where} \quad | |
\epsilon = p^{t+\Delta t} - \beta p^t"/></p> | |
</div><p>The above formulation of the future velocity is put into the continui… | |
equation:</p> | |
<div class="math"> | |
-<p><img src="_images/math/25961440701892ff5417a2474f7f52bb353e5d98.png" alt="\… | |
+<p><img src="_images/math/7b1a510ad1f355daa1f756cfd47f75a3a29f5660.png" alt="\… | |
\phi^t \nabla \cdot | |
-\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right) | |
+ | |
-\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right) | |
\cdot \nabla \phi^t + \frac{\Delta \phi^t}{\Delta t} = 0"/></p> | |
</div><div class="math"> | |
-<p><img src="_images/math/1c7a111b8952b4162797cce2cbd426d74083b834.png" alt="\… | |
+<p><img src="_images/math/964fa43eb7c25fa23b82daf41b1b9e16e998ab55.png" alt="\… | |
\phi^t \nabla \cdot | |
-\boldsymbol{v}^* - \frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+\boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \phi^t \nabla^2 \epsilon | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
-- \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+- \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t} = 0"/></p> | |
</div><div class="math"> | |
-<p><img src="_images/math/8a2e0d0a28c9ebec91bd12def7b3b36dd9e92da7.png" alt="\… | |
-\frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+<p><img src="_images/math/177b8a21bb3d0e2a33950cd8ceec4f96f5950bd6.png" alt="\… | |
+\frac{\Delta t}{\rho} \nabla^2 \epsilon | |
= \phi^t \nabla \cdot \boldsymbol{v}^* | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
-- \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+- \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t}"/></p> | |
</div><p>The pressure difference in time becomes a <a class="reference externa… | |
<div class="math"> | |
-<p><img src="_images/math/7572a9bb14a95da7f7b42fd4b5bec33c42800556.png" alt="\… | |
+<p><img src="_images/math/22176d7707317572382e54a12e13095b61d90c26.png" alt="\… | |
\nabla^2 \epsilon | |
-= \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
-+ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
+= \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
++ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t} | |
-+ \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t}"/></p> | |
++ \frac{\Delta \phi^t \rho}{\Delta t^2}"/></p> | |
</div><p>The right hand side of the above equation is termed the <em>forcing f… | |
<img class="math" src="_images/math/bb2c93730dbb48558bb3c4738c956c4e8f816437.p… | |
<div class="math"> | |
-<p><img src="_images/math/28d83ccd32811c0455355169ed8f6978d8d24dbe.png" alt="f… | |
-= \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
-+ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
-+ \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t} | |
+<p><img src="_images/math/470d00b97569e3a8976841debcd99aaf9c78d7a8.png" alt="f… | |
+= \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
++ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
++ \frac{\Delta \phi^t \rho}{\Delta t^2} | |
f_2 = | |
\frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t}"/></p> | |
diff --git a/doc/html/objects.inv b/doc/html/objects.inv | |
Binary files differ. | |
diff --git a/doc/html/python_api.html b/doc/html/python_api.html | |
t@@ -122,10 +122,10 @@ a <tt class="docutils literal"><span class="pre">sphere<… | |
57 | |
58 | |
59</pre></div></td><td class="code"><div class="highlight"><pre><span class="c… | |
-<span class="sd">"""</span> | |
+<span class="sd">'''</span> | |
<span class="sd">Example of two particles colliding.</span> | |
<span class="sd">Place script in sphere/python/ folder, and invoke with `pytho… | |
-<span class="sd">"""</span> | |
+<span class="sd">'''</span> | |
<span class="c"># Import the sphere module for setting up, running, and analyz… | |
<span class="c"># experiment. We also need the numpy module when setting array… | |
diff --git a/doc/html/searchindex.js b/doc/html/searchindex.js | |
t@@ -1 +1 @@ | |
-Search.setIndex({objects:{"":{sphere:[5,0,1,""]},sphere:{status:[5,2,1,""],con… | |
-\ No newline at end of file | |
+Search.setIndex({objects:{"":{sphere:[5,0,1,""]},sphere:{status:[5,2,1,""],con… | |
+\ No newline at end of file | |
diff --git a/doc/pdf/sphere.pdf b/doc/pdf/sphere.pdf | |
Binary files differ. | |
diff --git a/doc/sphinx/cfd.rst b/doc/sphinx/cfd.rst | |
t@@ -416,7 +416,8 @@ The predicted velocity is corrected using the new pressure… | |
.. math:: | |
\boldsymbol{v}^{t+\Delta t} = \boldsymbol{v}^* | |
- - \frac{\Delta t}{\rho} \nabla \epsilon | |
+ %- \frac{\Delta t}{\rho} \nabla \epsilon | |
+ - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon | |
\quad \text{where} \quad | |
\epsilon = p^{t+\Delta t} - \beta p^t | |
t@@ -426,25 +427,25 @@ equation: | |
.. math:: | |
\Rightarrow | |
\phi^t \nabla \cdot | |
- \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+ \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \ri… | |
+ | |
- \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho} \nabla \epsilon \right) | |
+ \left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \ri… | |
\cdot \nabla \phi^t + \frac{\Delta \phi^t}{\Delta t} = 0 | |
.. math:: | |
\Rightarrow | |
\phi^t \nabla \cdot | |
- \boldsymbol{v}^* - \frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+ \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \phi^t \nabla^2 \epsilon | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
- - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+ - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t} = 0 | |
.. math:: | |
\Rightarrow | |
- \frac{\Delta t}{\rho} \phi^t \nabla^2 \epsilon | |
+ \frac{\Delta t}{\rho} \nabla^2 \epsilon | |
= \phi^t \nabla \cdot \boldsymbol{v}^* | |
+ \nabla \phi^t \cdot \boldsymbol{v}^* | |
- - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho} | |
+ - \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t} | |
+ \frac{\Delta \phi^t}{\Delta t} | |
The pressure difference in time becomes a `Poisson equation`_ with added terms: | |
t@@ -452,19 +453,19 @@ The pressure difference in time becomes a `Poisson equat… | |
.. math:: | |
\Rightarrow | |
\nabla^2 \epsilon | |
- = \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
+ = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
+ + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t} | |
- + \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t} | |
+ + \frac{\Delta \phi^t \rho}{\Delta t^2} | |
The right hand side of the above equation is termed the *forcing function* | |
:math:`f`, which is decomposed into two terms, :math:`f_1` and :math:`f_2`: | |
.. math:: | |
f_1 | |
- = \frac{\nabla \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
- + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t \phi^t} | |
- + \frac{\Delta \phi^t \rho}{\Delta t^2 \phi^t} | |
+ = \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t} | |
+ + \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t} | |
+ + \frac{\Delta \phi^t \rho}{\Delta t^2} | |
f_2 = | |
\frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t} | |
diff --git a/src/device.cu b/src/device.cu | |
t@@ -1356,6 +1356,7 @@ __host__ void DEM::startTime() | |
dev_ns_p, | |
dev_ns_v, | |
dev_ns_v_p, | |
+ dev_ns_phi, | |
dev_ns_epsilon, | |
ns.beta, | |
ns.bc_bot, | |
diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh | |
t@@ -1746,16 +1746,6 @@ __global__ void findPredNSvelocities( | |
// Find the value of the forcing function. Only grad(epsilon) changes during | |
// the Jacobi iterations. The remaining, constant terms are only calculated | |
// during the first iteration. | |
-// The forcing function is: | |
-// f = (div(v_p)*rho)/dt | |
-// + (grad(phi) dot v_p*rho)/(dt*phi) | |
-// + (dphi*rho)/(dt*dt*phi) | |
-// - (grad(phi) dot grad(epsilon))/phi | |
-// The following is calculated in the first Jacobi iteration and saved: | |
-// f1 = (div(v_p)*rho)/dt | |
-// + (grad(phi) dot v_p*rho)/(dt*phi) | |
-// + (dphi*rho)/(dt*dt*phi) | |
-// f2 = grad(phi)/phi | |
// At each iteration, the value of the forcing function is found as: | |
// f = f1 - f2 dot grad(epsilon) | |
__global__ void findNSforcing( | |
t@@ -1811,9 +1801,13 @@ __global__ void findNSforcing( | |
// Find forcing function coefficients | |
//f1 = 0.0; | |
- f1 = div_v_p*devC_params.rho_f/devC_dt | |
+ /*f1 = div_v_p*devC_params.rho_f/devC_dt | |
+ 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/devC_dt | |
+ + dot(grad_phi, v_p)*devC_params.rho_f/devC_dt | |
+ + dphi*devC_params.rho_f/(devC_dt*devC_dt); | |
f2 = grad_phi/phi; | |
// Report values terms in the forcing function for debugging | |
t@@ -2135,6 +2129,7 @@ __global__ void updateNSvelocityPressure( | |
Float* dev_ns_p, | |
Float3* dev_ns_v, | |
Float3* dev_ns_v_p, | |
+ Float* dev_ns_phi, | |
Float* dev_ns_epsilon, | |
Float beta, | |
int bc_bot, | |
t@@ -2166,6 +2161,7 @@ __global__ void updateNSvelocityPressure( | |
const Float p_old = dev_ns_p[cellidx]; | |
const Float epsilon = dev_ns_epsilon[cellidx]; | |
const Float3 v_p = dev_ns_v_p[cellidx]; | |
+ const Float phi = dev_ns_phi[cellidx]; | |
// New pressure | |
Float p = beta*p_old + epsilon; | |
t@@ -2175,7 +2171,8 @@ __global__ void updateNSvelocityPressure( | |
= gradient(dev_ns_epsilon, x, y, z, dx, dy, dz); | |
// Find new velocity | |
- Float3 v = v_p - devC_dt/devC_params.rho_f*grad_epsilon; | |
+ //Float3 v = v_p - devC_dt/devC_params.rho_f*grad_epsilon; | |
+ Float3 v = v_p - devC_dt/(devC_params.rho_f*phi)*grad_epsilon; | |
// Print values for debugging | |
/* if (z == 0) { |