OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
JamesEMcclure
2021-09-03 09:55:03 -04:00
parent 2ceaa23afb d3a14c21ff
commit 3d64e1a8a1
17 changed files with 5925 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
common/Array.h
View File

@@ -658,7 +658,7 @@ public: // Math operations
void cat( const Array &x, int dim = 0 );
//! Initialize the array with random values (defined from the function table)
void rand();
//void rand();
//! Return true if NaNs are present
bool NaNs() const;

3
common/Array.hpp
View File

@@ -1324,11 +1324,12 @@ TYPE Array<TYPE, FUN, Allocator>::interp( const double *x ) const
/********************************************************
* Math operations (should call the Math class) *
********************************************************/
template<class TYPE, class FUN, class Allocator>
/*template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::rand()
{
FUN::rand( *this );
}
*/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &
Array<TYPE, FUN, Allocator>::operator+=( const Array<TYPE, FUN, Allocator> &rhs )

4
common/FunctionTable.cpp
View File

@@ -4,7 +4,7 @@
/********************************************************
* Random number generation *
********************************************************/
template<> char genRand<char>()
/*template<> char genRand<char>()
{
static std::random_device rd;
static std::mt19937 gen( rd() );
@@ -88,7 +88,7 @@ template<> long double genRand<long double>()
static std::uniform_real_distribution<double> dis;
return dis( gen );
}
*/
/********************************************************
* axpy *

7
common/FunctionTable.hpp
View File

@@ -39,21 +39,20 @@
#include <algorithm>
#include <cstring>
#include <limits>
#include <random>
//#include <random>
/********************************************************
* Random number initialization *
********************************************************/
template<class TYPE> TYPE genRand();
/*template<class TYPE> TYPE genRand();
template<class TYPE, class FUN>
inline void FunctionTable::rand( Array<TYPE, FUN> &x )
{
for ( size_t i = 0; i < x.length(); i++ )
x( i ) = genRand<TYPE>();
}
*/
/********************************************************
* Reduction *

7
docs/source/examples/running.rst
View File

@@ -1,9 +1,10 @@
============
Running LBPM
============
==============
LBPM examples
==============
There are two main components to running LBPM simulators.
First is understanding how to launch MPI tasks on your system,
which depends on the particular implementation of MPI that you are using,
as well as other details of the local configuration. The second component is
understanding the LBPM input file structure.

167
docs/source/userGuide/models/color/index.rst
View File

@@ -66,7 +66,172 @@ Model Formulation
****************************
Two LBEs are constructed to model the mass transport, incorporating the anti-diffusion
.. math::
:nowrap:
$$
A_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) = w_q N_a \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
+ \beta \frac{N_b}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big] \;
$$
.. math::
:nowrap:
$$
B_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) =
w_q N_b \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
- \beta \frac{N_a}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big]\;,
$$
The number density for each fluid is obtained from the sum of the mass transport distributions
.. math::
:nowrap:
$$
N_a = \sum_q A_q\;, \quad N_b = \sum_q B_q\;
$$
The phase indicator field is then defined as
.. math::
:nowrap:
$$
\phi = \frac{N_a-N_b}{N_a+N_b}
$$
The fluid density and kinematic viscosity are determined based on linear interpolation
.. math::
:nowrap:
$$
\rho_0 = \frac{(1+\phi) \rho_n}{2}+ \frac{(1-\phi) \rho_w}{2} \;,
$$
.. math::
:nowrap:
$$
s_\nu = \frac{(1+\phi)}{2\tau_n} +\frac{(1-\phi)}{2\tau_w} \;,
$$
where
.. math::
:nowrap:
$$
\nu_w = \frac{1}{3}\Big(\tau_w - \frac{1}{2} \Big) \;, \quad
\nu_n = \frac{1}{3}\Big(\tau_n - \frac{1}{2} \Big) \;.
$$
These values are then used to model the momentum transport.
The LBE governing momentum transport is defined based on a MRT relaxation process with additional
terms to account for the interfacial stresses
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = \sum^{Q-1}_{k=0} M^{-1}_{qk} \lambda_{k} (m_k^{eq}-m_k) + t_q \bm{\xi}_q \cdot \frac{\bm{F}}{c_s^2} \;,
$$
The moments are linearly indepdendent:
.. math::
:nowrap:
$$
m_k = \sum_{q=0}^{18} M_{qk} f_q\;.
$$
The relaxation parameters are determined from the relaxation time:
.. math::
:nowrap:
$$
\lambda_1 = \lambda_2= \lambda_9 = \lambda_{10}= \lambda_{11}= \lambda_{12}= \lambda_{13}= \lambda_{14}= \lambda_{15} = s_\nu \;,
$$
.. math::
:nowrap:
$$
\lambda_{4}= \lambda_{6}= \lambda_{8} = \lambda_{16} = \lambda_{17} = \lambda_{18}= \frac{8(2-s_\nu)}{8-s_\nu} \;,
$$
The non-zero equilibrium moments are defined as
.. math::
:nowrap:
$$
m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\
$$
.. math::
:nowrap:
$$
m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;,
$$
where the color gradient is determined from the phase indicator field
.. math::
:nowrap:
$$
\textbf{C}=\nabla \phi\;.
$$
and the unit normal vector is
.. math::
:nowrap:
$$
\bm{n} = \frac{\textbf{C}}{|\textbf{C}|}\;.
$$
****************************
Boundary Conditions
****************************

202
docs/source/userGuide/models/mrt/mrt.rst
View File

@@ -1,6 +1,200 @@
=============================================
###############################################################################
MRT model
=============================================
###############################################################################
The LBPM single fluid model is implemented by combining a multi-relaxation time (MRT) D3Q19
lattice Boltzmann equation (LBE) to solve for the momentum transport, recovering the Navier-Stokes
equations to second order based on the Chapman-Enskog expansion. The MRT model is used to assess the
permeability of digital rock images in either the Darcy or non-Darcy flow regimes.
A typical command to launch the LBPM color simulator is as follows
```
mpirun -np $NUMPROCS lbpm_permeability_simulator input.db
```
where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is
the name of the input database that provides the simulation parameters.
Note that the specific syntax to launch MPI tasks may vary depending on your system.
For additional details please refer to your local system documentation.
***************************
Model parameters
***************************
The essential model parameters for the color model are
- ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5`
****************************
Model Formulation
****************************
The LBE governing momentum transport is defined based on a MRT relaxation based on the D3Q19 discrete
velocity set, which determines the values :math:`\bm{\xi}_q`
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = \sum^{Q-1}_{k=0} M^{-1}_{qk} \lambda_{k} (m_k^{eq}-m_k) + t_q \bm{\xi}_q \cdot \frac{\bm{F}}{c_s^2} \;,
$$
Where :math:`\bm{F}` an external body force and :math:`c_s^2 = 1/3` is the speed of sound for the LB model.
The moments are linearly indepdendent functions of the distributions:
.. math::
:nowrap:
$$
m_k = \sum_{q=0}^{18} M_{qk} f_q\;.
$$
The non-zero equilibrium moments are
.. math::
:nowrap:
$$
m_1^{eq} = 19\frac{j_x^2+j_y^2+j_z^2}{\rho} - 11\rho \;,
$$
.. math::
:nowrap:
$$
m_2^{eq} = 3\rho - \frac{11}{2} \frac{j_x^2+j_y^2+j_z^2}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_4^{eq} = -\frac 2 3 j_x \;,
$$
.. math::
:nowrap:
$$
m_6^{eq} = -\frac 2 3 j_y \;,
$$
.. math::
:nowrap:
$$
m_8^{eq} = -\frac 2 3 j_z \;,
$$
.. math::
:nowrap:
$$
m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho}\;,
$$
.. math::
:nowrap:
$$
m_{10}^{eq} = -\frac{2j_x^2-j_y^2-j_z^2)}{2\rho} \;,
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{12}^{eq} = -\frac{j_y^2-j_z^2}{2\rho} \;,
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = \frac{j_x j_y}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = \frac{j_y j_z}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = \frac{j_x j_z}{\rho} \;,
$$
The relaxation parameters are determined based on the relaxation time :math:`\tau`
.. math::
:nowrap:
$$
\lambda_1 = \lambda_2= \lambda_9 = \lambda_{10}= \lambda_{11}= \lambda_{12}= \lambda_{13}= \lambda_{14}= \lambda_{15} = s_\nu = \frac{1}{\tau} \;,
$$
.. math::
:nowrap:
$$
\lambda_{4}= \lambda_{6}= \lambda_{8} = \lambda_{16} = \lambda_{17} = \lambda_{18}= \frac{8(2-s_\nu)}{8-s_\nu} \;,
$$
****************************
Example Input File
****************************
****************************
Boundary Conditions
****************************
The following external boundary conditions are supported by ``lbpm_permeability_simulator``
and can be set by setting the ``BC`` key values in the ``Domain`` section of the
input file database
- ``BC = 0`` -- fully periodic boundary conditions
- ``BC = 3`` -- constant pressure boundary condition
- ``BC = 4`` -- constant volumetric flux boundary condition
For ``BC = 0`` any mass that exits on one side of the domain will re-enter at the other
side. If the pore-structure for the image is tight, the mismatch between the inlet and
outlet can artificially reduce the permeability of the sample due to the blockage of
flow pathways at the boundary. LBPM includes an internal utility that will reduce the impact
of the boundary mismatch by eroding the solid labels within the inlet and outlet layers
(https://doi.org/10.1007/s10596-020-10028-9) to create a mixing layer.
The number mixing layers to use can be set using the key values in the ``Domain`` section
of the input database
- ``InletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the inlet
- ``OUtletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the outlet
For the other boundary conditions a thin reservoir of fluid (default ``3`` voxels)
is established at either side of the domain. The inlet is defined as the boundary face
where ``z = 0`` and the outlet is the boundary face where ``z = nprocz*nz``. By default a
reservoir of fluid A is established at the inlet and a reservoir of fluid B is established at
the outlet, each with a default thickness of three voxels. To over-ride the default label at
the inlet or outlet, the ``Domain`` section of the database may specify the following key values
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet
The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the
solution of the Navier-Stokes equations.

6
example/Bubble/input.db
View File

@@ -47,4 +47,8 @@ Analysis {
}
Visualization {
}
}
FlowAdaptor {
}

2
example/MicroModel/BlobSize.in
View File

@@ -1,2 +0,0 @@
1
32 32 32

6
example/MicroModel/Color.in
View File

@@ -1,6 +0,0 @@
0.7
1.0e-2 0.95 -1.0
0.7
0.0 0.0 0.0
0 1 1.04 0.96
900000 20000 1e-5

144
example/MicroModel/DiscPack.in
View File

@@ -1,144 +0,0 @@
33.13248649 29.93643243 15.41628378
75.07072973 26.66563068 16.99835135
118.9166622 26.51018919 13.54747297
155.9032432 26.54839365 18.37214865
194.4647297 25.64848784 13.23047838
236.4466216 27.04415726 18.86951351
276.9660811 27.18572081 18.66178378
316.6336486 27.16016042 17.04760811
360.6031081 29.63283514 14.62525676
400.4937838 26.33180541 15.91213514
436.6832432 29.38145946 15.30281081
479.5852703 27.09887028 18.76375676
15.17307297 69.60128378 14.82944595
53.0137973 64.87014865 16.21706757
94.52747297 68.84901351 17.6057027
135.5813919 65.70918919 16.92474324
173.2494595 69.96463514 14.05135135
215.6589189 68.16136486 17.94732432
259.9767568 70.38693243 14.80982432
297.0702703 67.83155405 14.68744595
336.5431081 67.57271622 17.27831081
378.3740541 67.72916216 18.13277027
419.2805405 66.79117568 17.42808108
455.1844595 72.06137838 13.27178784
33.72422973 108.8409189 17.24389189
74.43232432 108.0137838 18.75
114.1253784 106.4523919 17.1622027
153.882027 112.3773514 13.02102703
199.0522973 111.1353243 15.15556757
241.1697297 107.2128378 13.80640541
277.315 108.4926486 18.42555405
320.4716216 106.1688919 15.76947297
354.8572973 105.7547973 13.44669865
400.797973 110.3642297 16.28766216
439.5955405 107.9637838 18.33110811
482.8833784 107.3503784 13.76293243
10.30476351 149.3496351 15.64486486
55.28555405 147.8185946 16.92521622
95.24567568 149.0092432 17.58089189
131.8110405 150.626527 13.8865
179.3263514 146.4626757 14.22206757
216.0708108 152.3747568 14.82614865
261.7548649 147.6711757 13.01481622
299.3959459 147.4904324 13.76047297
339.8782432 149.7220541 16.66758108
373.3464865 152.5632297 13.45467297
417.9712162 147.855973 17.36840541
458.4025676 144.5771622 14.61705405
33.32618919 191.2175676 12.4216527
74.30910811 189.2091892 18.52312162
114.6439054 190.2966216 16.95763514
158.2144595 191.8867568 15.55116216
194.6536486 189.0091892 15.71814865
237.9637838 189.242027 17.27986486
271.4318919 186.4898649 13.0739527
317.7631081 188.4841892 16.13393243
358.1116216 191.2745946 14.83137838
398.7190541 189.4286486 18.13921622
439.5593243 189.8759459 16.8382973
482.1985135 194.4328378 12.19701081
13.90195446 229.4845946 18.36889189
55.69509459 226.9613514 16.02147297
96.98522973 232.4032432 14.46966216
137.4593243 227.1902703 13.4368527
174.2690541 228.4505405 16.73159459
216.2844595 229.237027 18.25581081
259.6358108 233.6717568 12.89051757
298.1172973 229.1837838 15.79105405
338.9263514 229.9989189 16.83474324
378.4617568 230.0036486 18.49422973
418.9406757 228.4566216 13.19394865
459.3066216 230.6732432 17.26921622
35.56972973 273.1122973 15.46828378
72.91724324 272.0410811 16.95751351
114.3243243 273.7995946 15.17266216
158.1804054 267.382973 13.8992027
194.7917568 269.7763514 17.70321622
236.57 271.7548649 12.18126081
272.412973 266.6922973 13.18252973
317.5945946 270.2414865 18.8062027
359.0810811 268.9272973 17.10451351
401.6145946 268.1164865 13.57114865
434.9897297 265.9721622 12.38897297
479.9452703 270.4847297 12.98088378
13.27824108 310.6754054 18.66364865
53.25386486 312.0385135 14.46228378
92.49813514 308.3451351 15.20474324
136.1003514 309.7871622 17.76025676
174.5641892 311.4295946 13.06244865
219.4805405 310.8893243 14.40940541
258.0839189 310.0159459 17.41775676
297.6322973 310.9914865 18.51806757
337.9448649 310.7555405 17.73358108
374.9710811 312.2205405 14.25995946
422.8451351 314.0112162 14.40621622
459.8067568 311.2154054 13.47336486
34.99885135 351.7381081 17.47081081
75.33945946 351.677027 15.52572973
114.8586486 351.3459459 18.89374324
155.7139189 350.8706757 18.07736486
191.6158108 352.282973 12.90032838
234.2532432 350.5644595 14.65478378
275.127027 350.5708108 16.79345946
317.652973 352.3060811 16.43624324
363.8067568 347.8472973 12.7394473
401.5713514 349.2537838 14.74306757
436.0778378 349.6571622 14.30860811
482.565 353.0252703 15.15274324
10.92297703 391.8589189 15.69475676
52.44290541 394.7612162 15.2132027
93.64185135 390.0597297 15.0062027
139.127527 389.3374324 14.73910811
175.5394595 391.7282432 18.06841892
216.4606757 392.4894595 17.91210811
251.7971622 390.4447297 12.52764324
298.7958108 393.1183784 16.58391892
336.7648649 392.7705405 16.06710811
379.9774324 392.8135135 16.87058108
417.1693243 392.1031081 16.41077027
459.4625676 391.9014865 18.88939189
32.60527027 431.387973 13.79572973
73.08481081 428.7002703 14.2487973
114.5770541 432.5933784 17.63783784
156.0283784 432.7117568 17.82533784
198.9813514 432.5832432 14.1537027
237.5352703 434.3167568 15.54622973
277.7677027 433.0071622 15.02437838
316.0105405 434.197973 16.34566216
359.23 437.6654054 13.07579189
401.887027 430.8952703 15.38594595
442.0587838 434.6848649 12.87081892
481.2160811 438.0641892 13.1349027
13.43983284 473.0802703 17.79305405
58.56352703 468.0102703 12.698
95.47004054 472.9767568 16.32575676
134.9401622 472.807027 18.70305405
173.3897297 469.9112162 15.53236486
216.2163514 473.1051351 18.61358108
256.3348649 474.0559459 14.83377027
294.8790541 471.1768919 14.53377027
339.6306757 472.1563514 16.27194595
377.4156757 471.6939189 17.32056757
423.6267568 474.1678378 12.76715541
459.5241892 471.4267568 16.65775676

3
example/MicroModel/Domain.in
View File

@@ -1,3 +0,0 @@
1 1 1
10 400 400
20.0 500.0 500.0

5466
example/MicroModel/FullMicromodel.discs Normal file
View File

File diff suppressed because it is too large Load Diff

47
example/MicroModel/GenerateCases.sh Executable file
View File

@@ -0,0 +1,47 @@
#!/bin/bash
tau1=1.18
tau2=0.7
alpha=0.005
Q="1.179757e-08 1.179757e-07 1.179757e-06"
#Cases for drainage
DrainWet="0.79 0.47 0.0"
# Cases for imbibition
ImbWet="0.92 0.47"
for q in $Q; do
echo $q;
flux=$(echo $q | sed 's/1.179757e-08/0.002/g')
flux=$(echo $flux | sed 's/1.179757e-07/0.02/g')
flux=$(echo $flux | sed 's/1.179757e-06/0.2/g')
for i in $DrainWet; do
NAME="Juanes_drain_Q"$flux"_wet"$i
echo $NAME
mkdir $NAME
echo "$tau1 $tau2" > $NAME/Color.in
echo "$alpha 0.95 $i" >> $NAME/Color.in
echo "0.0" >> $NAME/Color.in
echo "0.0 0.0 0.0" >> $NAME/Color.in
echo "0 4 $q 0.0" >> $NAME/Color.in
echo "5000000 25000 1e-5" >> $NAME/Color.in
echo "1000" >> $NAME/Color.in
done
for i in $ImbWet; do
NAME="Juanes_imb_Q"$flux"_wet"$i
echo $NAME
mkdir $NAME
echo "$tau1 $tau2" > $NAME/Color.in
echo "$alpha 0.95 $i" >> $NAME/Color.in
echo "0.0" >> $NAME/Color.in
echo "0.0 0.0 0.0" >> $NAME/Color.in
echo "0 4 $q 0.0" >> $NAME/Color.in
echo "5000000 25000 1e-5" >> $NAME/Color.in
echo "1000" >> $NAME/Color.in
done
done

21
example/MicroModel/ReadDiscPack.R Normal file
View File

@@ -0,0 +1,21 @@
require("ggplot2")
Discs<-read.csv("FullMicromodel.discs",head=FALSE,sep=" ")
colnames(Discs)<-c("cx","cy","radius")
L=0.45
# Extract some subset from the interior of the discs
SubDiscs<-subset(Discs,Discs$cy>0.9-L)
SubDiscs<-subset(SubDiscs,SubDiscs$cy<0.9+L)
SubDiscs<-subset(SubDiscs,SubDiscs$cx>0.9-L)
SubDiscs<-subset(SubDiscs,SubDiscs$cx<0.9+L)
SubDiscs$cx<-SubDiscs$cx-0.9+L
SubDiscs$cy<-SubDiscs$cy-0.9+L
write.table(SubDiscs,file="DiscPack.in",quote=FALSE,row.names=FALSE,col.names=FALSE,sep=" ")
ShowPlot<-ggplot(SubDiscs)+geom_circle(aes(x0=cx,y0=cy,r=radius))

6
example/Piston/input.db
View File

@@ -37,4 +37,8 @@ Analysis {
Visualization {
}
}
FlowAdaptor {
}

6
example/Plates/input.db
View File

@@ -36,4 +36,8 @@ Analysis {
}
Visualization {
}
}
FlowAdaptor {
}