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

fixing merge

Browse Source
This commit is contained in:
James McClure 2023-04-03 09:05:33 -04:00
commit 2989e890fa
14 changed files with 680 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
analysis/SubPhase.cpp
View File

@ -430,11 +430,11 @@ void SubPhase::Basic() {
//double total_flow_rate = water_flow_rate + not_water_flow_rate;
//double fractional_flow = water_flow_rate / total_flow_rate;
double h = Dm->voxel_length;
double krn = h * h * nu_n * Porosity * not_water_flow_rate / force_mag;
double krw = h * h * nu_w * Porosity* water_flow_rate / force_mag;
double krn = h * h * nu_n * Porosity* Porosity * not_water_flow_rate / force_mag;
double krw = h * h * nu_w * Porosity* Porosity* water_flow_rate / force_mag;
/* not counting films */
double krnf = krn - h * h * nu_n * Porosity * not_water_film_flow_rate / force_mag;
double krwf = krw - h * h * nu_w * Porosity * water_film_flow_rate / force_mag;
double krnf = krn - h * h * nu_n * Porosity* Porosity * not_water_film_flow_rate / force_mag;
double krwf = krw - h * h * nu_w * Porosity* Porosity * water_film_flow_rate / force_mag;
double eff_pressure = 1.0 / (krn + krw); // effective pressure drop
fprintf(TIMELOG,

213
docs/source/examples/membrane/cellModel.rst Normal file
View File

@ -0,0 +1,213 @@
********************************
Membrane Charging Dynamics
********************************
In this example, we consider membrane charging dynamics for a simple cell.
For the case considered in ``example/SingleCell`` an input membrane geometry is provided in the
file ``Bacterium.swc``, which specifies an oblong cell shape, relying on the ``.swc`` file format that
is commonly used to approximate neuron structures. The case considered is the four ion membrane transport
problem considered in Figure 4 from McClure & Li
The cell simulation is performed by the executable ``lbpm_nernst_planck_cell_simulator``, which is launched
in the same way as other parallel tools
.. code:: bash
mpirun -np 2 $LBPM_BIN/lbpm_nernst_planck_cell_simulator Bacterium.db
The input file ``Bacterium.db`` specifies the following
.. code:: c
MultiphysController {
timestepMax = 25000
num_iter_Ion_List = 4
analysis_interval = 100
tolerance = 1.0e-9
visualization_interval = 1000 // Frequency to write visualization data
}
Ions {
use_membrane = true
Restart = false
MembraneIonConcentrationList = 150.0e-3, 10.0e-3, 15.0e-3, 155.0e-3 //user-input unit: [mol/m^3]
temperature = 293.15 //unit [K]
number_ion_species = 4 //number of ions
tauList = 1.0, 1.0, 1.0, 1.0
IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 4.0e-3, 20.0e-3, 16.0e-3, 0.0e-3 //user-input unit: [mol/m^3]
BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
BC_InletList = 0, 0, 0, 0
BC_OutletList = 0, 0, 0, 0
}
Poisson {
lattice_scheme = "D3Q19"
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density
SolidLabels = 0 //solid labels for assigning solid boundary condition
SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2]
WriteLog = true //write convergence log for LB-Poisson solver
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "Bacterium.swc"
nproc = 2, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 128, 64, 64 // size of the input image
voxel_length = 0.01 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "swc"
ReadValues = 0, 1, 2
WriteValues = 0, 1, 2
}
Analysis {
analysis_interval = 100
subphase_analysis_interval = 50 // Frequency to perform analysis
restart_interval = 5000 // Frequency to write restart data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
save_electric_potential = true
save_concentration = true
save_velocity = false
}
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0, 0.0, 0.0
MassFractionIn = 1e-1, 1.0, 5e-3, 0.0
MassFractionOut = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionIn = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionOut = 1e-1, 1.0, 5e-3, 0.0
}
*******************
Example Output
*******************
Successful output looks like the following
.. code:: bash
********************************************************
Running LBPM Nernst-Planck Membrane solver
********************************************************
.... Read membrane permeability (MassFractionIn)
.... Read membrane permeability (MassFractionOut)
.... Read membrane permeability (ThresholdMassFractionIn)
.... Read membrane permeability (ThresholdMassFractionOut)
.... Read MembraneIonConcentrationList
voxel length = 0.010000 micron
voxel length = 0.010000 micron
Reading SWC file...
Number of lines in SWC file: 7
Number of lines extracted is: 7
shift swc data by 0.150000, 0.140000, 0.140000
Media porosity = 1.000000
LB Ion Solver: Initialized solid phase & converting to Signed Distance function
Domain set.
LB Ion Solver: Create ScaLBL_Communicator
LB Ion Solver: Set up memory efficient layout
LB Ion Solver: Allocating distributions
LB Ion Solver: Setting up device map and neighbor list
**** Creating membrane data structure ******
Number of active lattice sites (rank = 0): 262160
Membrane labels: 1
label=2, volume fraction = 0.133917
Creating membrane data structure...
Copy initial neighborlist...
Cut membrane links...
(cut 7105 links crossing membrane)
Construct membrane data structures...
Create device data structures...
Construct communication data structures...
Ion model setup complete
Analyze system with sub-domain size = 66 x 66 x 66
Set up analysis routines for 4 ions
LB Ion Solver: initializing D3Q7 distributions
...initializing based on membrane list
.... Set concentration(0): inside=0.15 [mol/m^3], outside=0.004 [mol/m^3]
.... Set concentration(1): inside=0.01 [mol/m^3], outside=0.02 [mol/m^3]
.... Set concentration(2): inside=0.015 [mol/m^3], outside=0.016 [mol/m^3]
.... Set concentration(3): inside=0.155 [mol/m^3], outside=0 [mol/m^3]
LB Ion Solver: initializing charge density
LB Ion Solver: solid boundary: non-flux boundary is assigned
LB Ion Solver: inlet boundary for Ion 1 is periodic
LB Ion Solver: outlet boundary for Ion 1 is periodic
LB Ion Solver: inlet boundary for Ion 2 is periodic
LB Ion Solver: outlet boundary for Ion 2 is periodic
LB Ion Solver: inlet boundary for Ion 3 is periodic
LB Ion Solver: outlet boundary for Ion 3 is periodic
LB Ion Solver: inlet boundary for Ion 4 is periodic
LB Ion Solver: outlet boundary for Ion 4 is periodic
*****************************************************
LB Ion Transport Solver:
Ion 1: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 2: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 3: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 4: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
*****************************************************
Ion model initialized
Main loop time_conv computed from ion 1: 2.5e-08[s/lt]
Main loop time_conv computed from ion 2: 2.5e-08[s/lt]
Main loop time_conv computed from ion 3: 2.5e-08[s/lt]
Main loop time_conv computed from ion 4: 2.5e-08[s/lt]
***********************************************************************************
LB-Poisson Solver: steady-state MaxTimeStep = 4000; steady-state tolerance = 1e-10
LB relaxation tau = 3.5
***********************************************************************************
LB-Poisson Solver: Use averaged MSE to check solution convergence.
LB-Poisson Solver: Use D3Q19 lattice structure.
voxel length = 0.010000 micron
voxel length = 0.010000 micron
Reading SWC file...
Number of lines in SWC file: 7
Number of lines extracted is: 7
shift swc data by 0.150000, 0.140000, 0.140000
Media porosity = 1.000000
LB-Poisson Solver: Initialized solid phase & converting to Signed Distance function
Domain set.
LB-Poisson Solver: Create ScaLBL_Communicator
LB-Poisson Solver: Set up memory efficient layout
LB-Poisson Solver: Allocating distributions
LB-Poisson Solver: Setting up device map and neighbor list
.... LB-Poisson Solver: check neighbor list
.... LB-Poisson Solver: copy neighbor list to GPU
Poisson solver created
LB-Poisson Solver: initializing D3Q19 distributions
LB-Poisson Solver: number of Poisson solid labels: 1
label=0, surface potential=0 [V], volume fraction=0
LB-Poisson Solver: number of Poisson initial-value labels: 2
label=1, initial potential=0 [V], volume fraction=0.96
label=2, initial potential=0 [V], volume fraction=0.13
POISSON MODEL: Reading restart file!
Poisson solver initialized
... getting Poisson solver error
-------------------------------------------------------------------
set coefficients
********************************************************
CPU time = 0.008526
Lattice update rate (per core)= 30.749833 MLUPS
Lattice update rate (total)= 61.499666 MLUPS
********************************************************

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

@ -20,3 +20,5 @@ a basic introduction to working with LBPM.
color/*
analysis/*
membrane/*

4
docs/source/publications/publications.rst
View File

@ -1,4 +1,4 @@
\============
============
Publications
============
@ -13,3 +13,5 @@ Publications
* J.E. McClure, M. Berrill, W. Gray, C.T. Miller, C.T. "Tracking interface and common curve dynamics for two-fluid flow in porous media. Journal of Fluid Mechanics, 796, 211-232 (2016) https://doi.org/10.1017/jfm.2016.212
* J.E.McClure, D.Adalsteinsson, C.Pan, W.G.Gray, C.T.Miller "Approximation of interfacial properties in multiphase porous medium systems" Advances in Water Resources, 30 (3): 354-365 (2007) https://doi.org/10.1016/j.advwatres.2006.06.010
* J.E. McClure, Z. Li "Capturing membrane structure and function in lattice Boltzmann models" arXiv preprint arXiv:2208.14122 https://arxiv.org/pdf/2208.14122.pdf

6
docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst
View File

@ -1,6 +0,0 @@
=============================================
Poisson-Boltzmann model
=============================================
The LBPM Poisson-Boltzmann solver is designed to solve the Poisson-Boltzmann equation
to solve for the electric field in an ionic fluid.

366
docs/source/userGuide/models/cell/cell.rst Normal file
View File

@ -0,0 +1,366 @@
=============================================
Cell model
=============================================
LBPM includes a whole-cell simulator based on a coupled solution of the Nernst-Planck equations with Gauss's law.
The resulting model is fully non-equilibrium, and can resolve the dynamics of how ions diffuse through the cellular
environment when subjected to complex membrane responses.
The lattice Boltzmann formulation is described below.
*********************
Nernst-Planck model
*********************
The Nernst-Planck model is designed to model ion transport based on the
Nernst-Planck equation.
.. math::
:nowrap:
$$
\frac{\partial C_k}{\partial t} + \nabla \cdot \mathbf{j}_k = 0
$$
where
.. math::
:nowrap:
$$
\mathbf{j}_k = C_k \mathbf{u} - D_k \Big( \nabla C_k + \frac{z_k C_k}{V_T} \nabla \psi\Big)
$$
A LBM solution is developed using a three-dimensional, seven velocity (D3Q7) lattice structure for each species. Each distribution is associated with a particular discrete velocity, such that the concentration is given by their sum,
.. math::
:nowrap:
$$
C_k = \sum_{q=0}^{6} f^k_q \;.
$$
Lattice Boltzmann equations (LBEs) are defined to determine the evolution of the distributions
.. math::
:nowrap:
$$
f^{k}_q (\mathbf{x}_n + \bm{\xi}_q \Delta t, t+ \Delta t)-
f^{k}_q (\mathbf{x}_n, t) = \frac{1}{\lambda_k}
\Big( f^{k}_q - f^{eq}_q \Big)\;,
$$
where the relaxation time :math:`\lambda_k` controls the bulk diffusion coefficient,
.. math::
:nowrap:
$$
D_k = c_s^2\Big( \lambda_k - \frac 12\Big)\;.
$$
The speed of sound for the D3Q7 lattice model is :math:`c_s^2 = \frac 14` and the weights are :math:`W_0 = 1/4` and :math:`W_1,\ldots, W_6 = 1/8`.
Equilibrium distributions are established from the fact that molecular velocity distribution follows a Gaussian distribution within the bulk fluids,
.. math::
:nowrap:
$$
f^{eq}_q = W_q C_k \Big[ 1 + \frac{\bm{\xi_q}\cdot \mathbf{u}^\prime}{c_s^2} \Big]\;.
$$
The velocity is given by
.. math::
:nowrap:
$$
\mathbf{u}^\prime = \mathbf{u} - \frac{z_k D_k}{V_T} \nabla \psi \;.
$$
Keys for the Nernst-Planck solver are provided in the ``Ion`` section of the input file database. Supported keys are
- ``use_membrane`` -- set up a membrane structure (defaults to ``true`` if not specified)
- ``Restart`` -- read concentrations from restart file (defaults to ``false`` if not specified)
- ``number_ion_species`` -- number of ions to use in the model
- ``temperature`` -- temperature to use for the thermal voltage (:math:`V_T=k_B T / e`, where the electron charge is :math:`e=1.6\times10^{-19}` Coulomb)
- ``FluidVelDummy`` -- vector providing a dummy fluid velocity field (for advection component)
- ``ElectricFieldDummy`` -- vectory providing a dummy electric field (for force component)
- ``tauList`` -- list of relaxation times to set the diffusion coefficient based on :math:`\lambda_k`.
- ``IonDiffusivityList`` -- list of physical ion diffusivities in units :math:`\mbox{m}^2/\mbox{second}`.
- ``IonValenceList`` -- list of ion valence charges for each ion in the model.
- ``IonConcentrationList`` -- list of concentrations to set for each ion.
- ``MembraneIonConcentrationList`` -- list of concentrations to set for each ion inside the membrane.
- ``BC_InletList`` -- boundary conditions for each ion at the z-inlet (``0`` for periodic, ``1`` to set concentration)
- ``BC_OutletList`` -- boundary conditions for each ion at the z-outlet
- ``InletValueList`` -- concentration value to set at the inlet (if not periodic)
- ``OutletValueList`` -- concentration value to set at the outlet (if not periodic)
*********************
Gauss's Law Model
*********************
The LBPM Gauss's law solver is designed to solve for the electric field in an ionic fluid.
.. math::
:nowrap:
$$
\nabla^2_{fe} \psi (\mathbf{x}_i) = \frac{1}{6 \Delta x^2}
\Bigg( 2 \sum_{q=1}^{6} \psi(\mathbf{x}_i + \bm{\xi}_q \Delta t)
+ \sum_{q=7}^{18} \psi(\mathbf{x}_i + \bm{\xi}_q \Delta t)
- 24 \psi (\mathbf{x}_i) \Bigg) \;,
$$
The equilibrium functions are defined as
.. math::
:nowrap:
$$
g_q^{eq} = w_q \psi\;,
$$
where :math:`w_0=1/2`, :math:`w_q=1/24` for :math:`q=1,\ldots,6` and :math:`w_q=1/48` for :math:`q=7,\ldots,18`
which implies that
.. math::
:nowrap:
$$
\psi = \sum_{q=0}^{Q} g_q^{eq}\;.
$$
Given a particular initial condition for :math:`\psi`, let us consider application of the standard D3Q19 streaming step based on the equilibrium distributions
.. math::
:nowrap:
$$
g_q^\prime(\mathbf{x}, t) = g_q^{eq}(\mathbf{x}-\bm{\xi}_q\Delta t, t+ \Delta t)\;.
$$
Relative to the solution of Gauss's law, the error is given by
.. math::
:nowrap:
$$
\varepsilon_{\psi} =
8 \Big[ -g_0 + \sum_{q=1}^Q g_q^\prime(\mathbf{x}, t) \Big]
+ \frac{\rho_e}{\epsilon_r \epsilon_0} \;.
$$
Using the fact that :math:`f_0 = W_0 \psi`, we can compute the value
:math:`\psi^\prime` that would kill the error. We set :math:`\varepsilon_{\psi}=0`
and rearrange terms to obtain
.. math::
:nowrap:
$$
\psi^\prime (\mathbf{x},t) = \frac{1}{W_0}\Big[ \sum_{q=1}^Q g_q^\prime(\mathbf{x}, t)
+ \frac{1}{8}\frac{\rho_e}{\epsilon_r \epsilon_0}\Big] \;.
$$
The local value of the potential is then updated based on a relaxation scheme, which is controlled by the relaxation time :math:`\tau_\psi`
.. math::
:nowrap:
$$
\psi(\mathbf{x},t+\Delta t) \leftarrow \Big(1 - \frac{1}{\tau_\psi} \Big )\psi (\mathbf{x},t)
+ \frac{1}{\tau_\psi} \psi^\prime (\mathbf{x},t)\;.
$$
The algorithm can then proceed to the next timestep.
Keys to control the Gauss's law solver are specified in the ``Poisson`` section of the input database.
Supported keys are:
- ``Restart`` -- read electric potential from a restart file (default ``false``)
- ``timestepMax`` -- maximum number of timesteps to run before exiting
- ``tau`` -- relaxation time
- ``analysis_interval`` -- how often to check solution for steady state
- ``tolerance`` -- controls the required accuracy
- ``epsilonR`` -- controls the electric permittivity
- ``WriteLog`` -- write a convergence log
***************************
Membrane Model
***************************
The LBPM membrane model provides the basis to model cellular dynamics.
There are currently two supported ways to specify the membrane location:
1. provide a segemented image that is labeled to differentiate the cell
interior and exterior. See the script ``NaCl-cell.py`` and input file ``NaCl.db`` as a reference for how to use labeled images.
- ``IonConcentrationFile`` -- list of files that specify the initial concentration for each ion
- ``Filename`` -- 8-bit binary file provided in the ``Domain`` section of the input database
- ``ReadType`` -- this should be ``"8bit"`` (this is the default)
2. provide a ``.swc`` file that specifies the geometry (see example input file below).
- ``Filename`` -- swc file name should be provided in the ``Domain`` section of the input database
- ``ReadType`` -- this should be ``"swc"`` (required since ``"8bit"`` is the internal default)
Example input files for both cases are stored within the LBPM repository, located at ``example/SingleCell/``
The membrane simply prevents the diffusion of ions. All lattice links crossing the membrane are stored in a dedicated data structure so that transport is decoupled from the bulk regions. Suppose that site :math:`\mathbf{x}_{q\ell}` is inside the membrane and :math:`\mathbf{x}_{p\ell}` is outside the membrane, with :math:`\mathbf{x}_{p \ell } = \mathbf{x}_{q\ell} + \bm{\xi}_q \Delta t`. For each species :math:`k`, transport across each link :math:`\ell` is controlled by a pair of coefficients, :math:`\alpha^k_{\ell p}` and :math:`\alpha^k_{\ell q}`. Ions transported from the outside to the inside are transported by the particular distribution that is associated with the direction :math:`\xi_q`
.. math::
:nowrap:
$$
{ f_{q}^{k \prime} (\mathbf{x}_{q \ell}) \gets (1-\alpha^k_{\ell q}) f_{q}^{k} (\mathbf{x}_{q\ell}) + \alpha^k_{\ell p } f_{ p}^{k} (\mathbf{x}_{p\ell})}
$$
Similarly, for ions transported from the inside to the outside
.. math::
:nowrap:
$$
{f_{p}^{k \prime} (\mathbf{x}_{p\ell}) \gets (1-\alpha^k_{\ell p}) f_{p}^{k} (\mathbf{x}_{p\ell}) + \alpha^k_{\ell q } f_{q}^{k} (\mathbf{x}_{q\ell})}
$$
The basic closure relationship that is implemented is for voltage-gated ion channels.
Let :math:`\Delta \psi_\ell = \psi(\mathbf{x}_{p\ell} ,t) - \psi(\mathbf{x}_{q\ell},t)` be the membrane potential across link :math:`\ell`. Since :math:`\psi` is determined based on the charge density, :math:`\Delta \psi_\ell` can vary with both space and time. The behavior of the gate is implmented as follows,
.. math::
:nowrap:
$$
\Delta \psi_\ell > \tilde{V}_m\; \Rightarrow \; \mbox{gate is open} \; \Rightarrow \; \alpha^{k}_{q \ell} = \alpha_{1} + \alpha_2\;,
$$
and
.. math::
:nowrap:
$$
\Delta \psi_\ell \le \tilde{V}_m\; \Rightarrow \; \mbox{gate is closed}\; \Rightarrow \; \alpha^{{k}}_{q \ell} = \alpha_1\;
$$
where :math:`\tilde{V}_m` is the membrane voltage threshold that controls gate. Mass conservation dictates that
.. math::
:nowrap:
$$
\alpha_1 \ge 0\;, \quad \alpha_2 \ge 0\;, \quad \alpha_1 + \alpha_2 \le 1\;.
$$
The rule is enforced based on the Heaviside function, as follows
.. math::
:nowrap:
$$
\alpha_{\ell q}^{k} (\Delta \psi_\ell) = \alpha_1 + \alpha_2 H\big(\Delta \psi_\ell - \tilde{V}_m \big)\;.
$$
Note that different coefficients are specified for each ion in the model.
Keys for the membrane model are set in the ``Membrane`` section of the input file database. Supported keys are
- ``VoltageThreshold`` -- voltage threshold (may be different for each ion)
- ``MassFractionIn`` -- value of :math:`\alpha^k_{\ell p}` when the voltage threshold is not met
- ``MassFractionOut`` -- value of :math:`\alpha^k_{\ell q}` when the voltage threshold is not met
- ``ThresholdMassFractionIn`` -- value of :math:`\alpha^k_{\ell p}` when the voltage threshold is met
- ``ThresholdMassFractionOut`` -- value of :math:`\alpha^k_{\ell q}` when the voltage threshold is met
****************************
Example Input File
****************************
.. code-block:: c
MultiphysController {
timestepMax = 25000
num_iter_Ion_List = 4
analysis_interval = 100
tolerance = 1.0e-9
visualization_interval = 1000 // Frequency to write visualization data
}
Ions {
use_membrane = true
Restart = false
MembraneIonConcentrationList = 150.0e-3, 10.0e-3, 15.0e-3, 155.0e-3 //user-input unit: [mol/m^3]
temperature = 293.15 //unit [K]
number_ion_species = 4 //number of ions
tauList = 1.0, 1.0, 1.0, 1.0
IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 4.0e-3, 20.0e-3, 16.0e-3, 0.0e-3 //user-input unit: [mol/m^3]
BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration
//SolidLabels = 0 //solid labels for assigning solid boundary condition; ONLY for BC_Solid=1
//SolidValues = 1.0e-5 // user-input surface ion concentration unit: [mol/m^2]; ONLY for BC_Solid=1
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
BC_InletList = 0, 0, 0, 0
BC_OutletList = 0, 0, 0, 0
}
Poisson {
lattice_scheme = "D3Q19"
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density
SolidLabels = 0 //solid labels for assigning solid boundary condition
SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2]
WriteLog = true //write convergence log for LB-Poisson solver
// ------------------------------- Testing Utilities ----------------------------------------
// ONLY for code debugging; the followings test sine/cosine voltage BCs; disabled by default
TestPeriodic = false
TestPeriodicTime = 1.0 //unit:[sec]
TestPeriodicTimeConv = 0.01 //unit:[sec]
TestPeriodicSaveInterval = 0.2 //unit:[sec]
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "Bacterium.swc"
nproc = 2, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 128, 64, 64 // size of the input image
voxel_length = 0.01 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "swc"
ReadValues = 0, 1, 2
WriteValues = 0, 1, 2
}
Analysis {
analysis_interval = 100
subphase_analysis_interval = 50 // Frequency to perform analysis
restart_interval = 5000 // Frequency to write restart data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
save_electric_potential = true
save_concentration = true
save_velocity = false
}
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0, 0.0, 0.0
MassFractionIn = 1e-1, 1.0, 5e-3, 0.0
MassFractionOut = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionIn = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionOut = 1e-1, 1.0, 5e-3, 0.0
}

4
docs/source/userGuide/models/index.rst
View File

@ -12,9 +12,7 @@ Currently supported lattice Boltzmann models
mrt/*
nernstPlanck/*
PoissonBoltzmann/*
cell/*
greyscale/*

6
docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst
View File

@ -1,6 +0,0 @@
=============================================
Nernst-Planck model
=============================================
The Nernst-Planck model is designed to model ion transport based on the
Nernst-Planck equation.

39
example/SingleCell/NaCl-cell.py
View File

@ -1,41 +1,38 @@
import numpy as np
import matplotlib.pylab as plt
Nx = 64
Ny = 64
Nz = 64
cx = Nx/2
cy = Ny/2
cz = Nz/2
radius = 12
D=np.ones((40,40,40),dtype="uint8")
D=np.ones((Nx,Ny,Nz),dtype="uint8")
cx = 20
cy = 20
cz = 20
radius = 8
for i in range(0, Nx):
for j in range (0, Ny):
for k in range (0,Nz):
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
if (dist < radius ) :
D[i,j,k] = 2
D.tofile("cell_64x64x64.raw")
D.tofile("cell_40x40x40.raw")
C1=np.zeros((Nx,Ny,Nz),dtype="double")
C2=np.zeros((Nx,Ny,Nz),dtype="double")
C1=np.zeros((40,40,40),dtype="double")
C2=np.zeros((40,40,40),dtype="double")
for i in range(0, Nx):
for j in range (0, Ny):
for k in range (0,Nz):
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
#outside the cell
C1[i,j,k] = 125.0e-6 # Na
C2[i,j,k] = 125.0e-6 # Cl
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
# inside the cell
if (dist < radius ) :
C1[i,j,k] = 110.0e-6
C2[i,j,k] = 110.0e-6
C1[i,j,k] = 5.0e-6
C2[i,j,k] = 5.0e-6
C1.tofile("cell_concentration_Na_64x64x64.raw")
C2.tofile("cell_concentration_Cl_64x64x64.raw")
C1.tofile("cell_concentration_Na_40x40x40.raw")
C2.tofile("cell_concentration_Cl_40x40x40.raw")

38
example/SingleCell/NaCl.db
View File

@ -1,9 +1,10 @@
MultiphysController {
timestepMax = 2000
timestepMax = 20
num_iter_Ion_List = 2
analysis_interval = 40
analysis_interval = 50
tolerance = 1.0e-9
visualization_interval = 40 // Frequency to write visualization data
visualization_interval = 100 // Frequency to write visualization data
analysis_interval = 50 // Frequency to perform analysis
}
Stokes {
tau = 1.0
@ -12,7 +13,9 @@ Stokes {
nu_phys = 0.889e-6 //fluid kinematic viscosity; user-input unit: [m^2/sec]
}
Ions {
IonConcentrationFile = "cell_concentration_Na_64x64x64.raw", "double", "cell_concentration_Cl_64x64x64.raw", "double"
use_membrane = true
Restart = false
IonConcentrationFile = "cell_concentration_Na_40x40x40.raw", "double", "cell_concentration_Cl_40x40x40.raw", "double"
temperature = 293.15 //unit [K]
number_ion_species = 2 //number of ions
tauList = 1.0, 1.0
@ -25,6 +28,8 @@ Ions {
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
}
Poisson {
lattice_scheme = "D3Q19"
Restart = false
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
@ -41,16 +46,19 @@ Poisson {
TestPeriodicTimeConv = 0.01 //unit:[sec]
TestPeriodicSaveInterval = 0.2 //unit:[sec]
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
timestepMax = 100000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 200 //timestep checking steady-state convergence
tolerance = 1.0e-6 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "cell_64x64x64.raw"
Filename = "cell_40x40x40.raw"
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 64, 64, 64 // size of the input image
voxel_length = 0.1 //resolution; user-input unit: [um]
n = 40, 40, 40 // Size of local domain (Nx,Ny,Nz)
N = 40, 40, 40 // size of the input image
voxel_length = 1.0 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "8bit"
ReadValues = 0, 1, 2
@ -72,9 +80,9 @@ Visualization {
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0
MassFractionIn = 1e-2, 1e-8
MassFractionOut = 1e-2, 1e-8
ThresholdMassFractionIn = 1e-2, 1e-8
ThresholdMassFractionOut = 1e-2, 1e-8
MassFractionIn = 0e-2, 5e-2
MassFractionOut = 0e-2, 5e-2
ThresholdMassFractionIn = 0e-2, 5e-2
ThresholdMassFractionOut = 0e-2, 5e-2
}

50
models/IonModel.cpp
View File

@ -300,7 +300,7 @@ void ScaLBL_IonModel::ReadParams(string filename, vector<int> &num_iter) {
void ScaLBL_IonModel::ReadParams(string filename) {
//NOTE: the maximum iteration timesteps for ions are left unspecified
// it relies on the multiphys controller to compute the max timestep
USE_MEMBRANE = false;
USE_MEMBRANE = true;
Restart = false;
// read the input database
db = std::make_shared<Database>(filename);
@ -1081,20 +1081,42 @@ void ScaLBL_IonModel::Initialize() {
printf("LB Ion Solver: initializing D3Q7 distributions\n");
//USE_MEMBRANE = true;
if (USE_MEMBRANE){
double *Ci_host;
if (rank == 0)
printf(" ...initializing based on membrane list \n");
Ci_host = new double[number_ion_species * Np];
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentrationMembrane( &Ci_host[ic * Np], ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
}
delete[] Ci_host;
double *Ci_host;
Ci_host = new double[number_ion_species * Np];
if (ion_db->keyExists("IonConcentrationFile")) {
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");
if (File_ion.size() == 2*number_ion_species) {
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentration_FromFile(&Ci_host[ic * Np], File_ion,
ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host,
number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np],
Np);
}
}
}
else{
if (rank == 0)
printf(" ...initializing based on membrane list \n");
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentrationMembrane( &Ci_host[ic * Np], ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
}
}
delete[] Ci_host;
}
else if (ion_db->keyExists("IonConcentrationFile")) {
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");

8
models/MRTModel.cpp
View File

@ -392,19 +392,17 @@ void ScaLBL_MRTModel::Run() {
double h = Dm->voxel_length;
double absperm =
h * h * mu * Mask->Porosity() * flow_rate / force_mag;
double absperm_adjusted =
h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag;
absperm_adjusted *= 1013.0; // Convert to mDarcy
absperm *= 1013.0; // Convert to mDarcy
if (rank == 0) {
printf(" %f\n", absperm);
FILE *log_file = fopen("Permeability.csv", "a");
fprintf(log_file,
"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g "
"%.8g %.8g %.8g\n",
"%.8g %.8g\n",
timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As,
h * Hs, Xs, vax, vay, vaz, absperm, absperm_adjusted);
h * Hs, Xs, vax, vay, vaz, absperm);
fclose(log_file);
}
}

17
models/PoissonSolver.cpp
View File

@ -995,7 +995,6 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo
delete [] host_Error;
//************************************************************************/
if(WriteLog==true){
getConvergenceLog(timestep,error);
}
@ -1051,22 +1050,22 @@ void ScaLBL_Poisson::SolveElectricPotentialAAodd(int timestep_from_Study){
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
@ -1121,22 +1120,22 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}

4
tests/lbpm_color_simulator.cpp
View File

@ -96,8 +96,10 @@ int main( int argc, char **argv )
SKIP_TIMESTEPS = 0;
if (PROTOCOL == "fractional flow")
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
if (PROTOCOL == "seed water")
if (PROTOCOL == "seed water"){
SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
}
}
/* analysis keys*/
int ANALYSIS_INTERVAL = ColorModel.timestepMax;