Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA
This commit is contained in:
commit
ef6d275161
38
docs/source/userGuide/models/color/analysis/basic.rst
Normal file
38
docs/source/userGuide/models/color/analysis/basic.rst
Normal file
|
@ -0,0 +1,38 @@
|
|||
======================================
|
||||
Color model -- Basic Analysis
|
||||
======================================
|
||||
|
||||
The basic analysis routine for the LBPM color model logs a time series
|
||||
of averaged information to the space-delimited CSV file ``timelog.csv``.
|
||||
The ``analysis_interval`` key should be specified to control the interval at
|
||||
to perform basic analysis, which corresponds to the number of simulation timesteps
|
||||
to perform between analyses. The basic analysis routine is designed to
|
||||
be lightweight so that it can be performed frequently, and it is almost always
|
||||
useful to enable it. Results will be immediately logged
|
||||
to ``timelog.csv`` which can be used to assess whether or not the simulation is
|
||||
behaving as expected. Furthermore, the analysis framework will check
|
||||
simulation measures to verfiy that no ``NaN`` are detected for the fluid
|
||||
pressure and flow rate.
|
||||
|
||||
The list of measures logged to ``timelog.csv`` are defined as follows.
|
||||
The region of space occupied by the wetting fluid is determined from the
|
||||
phase indicator field, :math:`\Omega_w:\phi<0`
|
||||
|
||||
|
||||
* ``sw`` -- water saturation (fluid component 2)
|
||||
* ``krw`` -- water effective permeability
|
||||
* ``krw`` -- effective permeability for fluid w
|
||||
* ``krn`` -- effective permeability for fluid n
|
||||
* ``krwf`` -- effective permeability for fluid w (with film uncertainty estimate)
|
||||
* ``krnf`` -- effective permeability for fluid n (with film uncertainty estimate)
|
||||
* ``vw`` -- speed for the non-wetting fluid
|
||||
* ``vn`` -- speed for the wetting fluid
|
||||
* ``force`` -- magnitude for effective driving force
|
||||
* ``pw`` -- average pressure for fluid w
|
||||
* ``pn`` -- average pressure for fluid n
|
||||
* ``wet`` -- total solid wetting energy
|
||||
|
||||
|
||||
More comprehensive analysis is performed in the ``subphase`` analysis module.
|
||||
|
||||
|
157
docs/source/userGuide/models/color/analysis/subphase.rst
Normal file
157
docs/source/userGuide/models/color/analysis/subphase.rst
Normal file
|
@ -0,0 +1,157 @@
|
|||
======================================
|
||||
Color model -- Subphase Analysis
|
||||
======================================
|
||||
|
||||
The subphase analysis routine for the LBPM color model logs a time series
|
||||
of averaged information to the space-delimited CSV file ``subphase.csv``.
|
||||
The ``subphase_analysis_interval`` key should be specified to control the interval at
|
||||
to perform basic analysis, which corresponds to the number of simulation timesteps
|
||||
to perform between analyses. The subphase analys routine performs the following functions
|
||||
|
||||
* analyzes the connectivity of fluid phases using a connected components algorithm
|
||||
* constructs iso-surfaces to represent interfaces within the system
|
||||
* computes averages of physical quantities based on the defined entities
|
||||
|
||||
Since it is more computationally expensive to carry out these operations compared to the
|
||||
basic analysis, it may be useful to choose ``subphase_analysis_interval`` to be larger than
|
||||
``analysis_interval``. Nevertheless, since analysis is performed in memory, it is orders of
|
||||
magnitude faster than to analyze data *in situ* rather than writing fields to disc. Monitoring
|
||||
the performance in MLUPS provides an easy way to tune the analysis interval so that the
|
||||
overall simulation performance is not subject to significant penalties.
|
||||
|
||||
There are three main entities that are constructed for subphase analysis
|
||||
|
||||
* :math:`\Omega_w:\phi<0` : region of space filled by fluid w
|
||||
* :math:`\Omega_n:\phi<0` : region of space filled by fluid n
|
||||
* :math:`\Omega_i: |\nabla \phi|<0 > \epsilon` : region of space for the interface region
|
||||
|
||||
The phase regions are defined identical to what is reported in ``timelog.csv``.
|
||||
The interface region is defined explicitly as the portion of space where
|
||||
significant composition gradients are present. For each entity, sub-entities are
|
||||
constructed by running a connected components algorithm on the set. This is done to
|
||||
separate the connected part of the phase from the disconnected part. This subset operation
|
||||
is performed by identifying the largest connected component (based on volume)
|
||||
and denoting this as the connected part of that region. The remaining portion of the
|
||||
phase is lumped into a disconnected sub-entity. Subphase analysis is therefore performed
|
||||
for the following six entities
|
||||
|
||||
* ``wc`` -- connected part of phase w
|
||||
* ``wd`` -- disconnected part of phase w
|
||||
* ``nc`` -- connected part of phase n
|
||||
* ``nd`` -- disconnected part of phase n
|
||||
* ``ic`` -- connected part of interface region
|
||||
* ``id`` -- disconnected part of interface region
|
||||
|
||||
For each entity :math:`\Omega_k` with :math:`k\in\{wc,wd,nc,nd,ic,id\}`
|
||||
an isosurface is constructed to approximate the boundary of the region,
|
||||
:math:`\Gamma_k`. Once each region is identified, the following measures are determined
|
||||
|
||||
**Geometric invariants**
|
||||
|
||||
* Volume -- :math:`V_k=\int_{\Omega_k} dV`
|
||||
* Surface area -- :math:`A_k=\int_{\Gamma_k} dS`
|
||||
* Integral mean curvature -- :math:`H_k=\int_{\Gamma_k} \frac 12 (\kappa_1 + \kappa_2) dS`
|
||||
* Euler characteristic -- :math:`\chi_k= \frac{1}{4\pi} \int_{\Gamma_k} \kappa_1 \kappa_2 dS`
|
||||
|
||||
**Conserved quantities**
|
||||
|
||||
* Total mass within the region :math:`M_k=\int_{\Omega_k} \rho dV`
|
||||
* Total momentum within the region :math:`\mathbf{P}_k=\int_{\Omega_k} \rho \mathbf{u} dV`
|
||||
* Total kinetic energy within the region :math:`K_k=\frac 12 \int_{\Omega_k} \rho \mathbf{u} \cdot \mathbf{u} dV`
|
||||
|
||||
**Thermodynamic quantities**
|
||||
|
||||
* Pressure -- :math:`p_k=\frac{1}{V_k}\int_{\Omega_k} p dV`
|
||||
* Solid wetting energy -- :math:`\gamma_s=\int_{\Gamma_s}\gamma dS`
|
||||
* Viscous dissipation -- :math:`\Phi_k=\int_{\Omega_k}\bm{\varepsilon} : \nabla \mathbf{u} dV`
|
||||
|
||||
The total solid wetting energy is determined by integrating the interfacial stresses in the
|
||||
immediate vicinity of the solid surface :math:`\Gamma_s`. The integral of the
|
||||
dissipation function is determined based on the viscous stress tensor, denoted by :math:`\bm{\varepsilon}`.
|
||||
|
||||
|
||||
The full list of measures are associated with the labels in ``subphase.csv``
|
||||
|
||||
* ``time`` -- timestep
|
||||
* ``rn`` -- density for phase n (input parameter)
|
||||
* ``rw`` -- density for phase w (input parameter)
|
||||
* ``nun`` -- kinematic viscosity for phase n (input parameter)
|
||||
* ``nuw`` -- kinematic viscosity for phase w (input parameter)
|
||||
* ``Fx`` -- external force in x direction (input parameter)
|
||||
* ``Fy`` -- external force in y direction (input parameter)
|
||||
* ``Fz`` -- external force in z direction (input parameter)
|
||||
* ``iftwn`` -- interfacial tension (input parameter)
|
||||
* ``wet`` -- total solid wetting energy
|
||||
* ``pwc`` -- average pressure for connected part of phase w
|
||||
* ``pwd`` -- average pressure for disconnected part of phase w
|
||||
* ``pnc`` -- average pressure for connected part of phase n
|
||||
* ``pnd`` -- average pressure for disconnected part of phase n
|
||||
* ``Mwc`` -- mass for connected part of phase w
|
||||
* ``Mwd`` --mass for disconnected part of phase w
|
||||
* ``Mwi`` -- mass for phase within diffuse interface region
|
||||
* ``Mnc`` -- mass for connected part of phase n
|
||||
* ``Mnd`` -- mass for disconnected part of phase n
|
||||
* ``Mni`` -- mass for phase n within diffuse interface region
|
||||
* ``Msw`` -- mass for component w within 2 voxels of solid
|
||||
* ``Msn`` -- mass for component n within 2 voxels of solid
|
||||
* ``Pwc_x`` -- x- momentum for connected part of phase w
|
||||
* ``Pwd_x`` -- x- momentum for disconnected part of phase w
|
||||
* ``Pwi_x`` -- x- momentum for phase w within diffuse interface
|
||||
* ``Pnc_x`` -- x- momentum for connected part of phase n
|
||||
* ``Pnd_x`` -- x- momentum for disconnected part of phase n
|
||||
* ``Pni_x`` -- x- momentum for phase n within diffuse interface
|
||||
* ``Psw_x`` -- x- momentum for phase w within 2 voxels of solid
|
||||
* ``Psn_x`` -- x- momentum for phase n within 2 voxels of solid
|
||||
* ``Pwc_y`` -- y- momentum for connected part of phase w
|
||||
* ``Pwd_y`` -- y- momentum for disconnected part of phase w
|
||||
* ``Pwi_y`` -- y- momentum for phase w within diffuse interface
|
||||
* ``Pnc_y`` -- y- momentum for connected part of phase n
|
||||
* ``Pnd_y`` -- y- momentum for connected part of phase n
|
||||
* ``Pni_y`` -- y- momentum for phase n within diffuse interface
|
||||
* ``Psw_y`` -- y- momentum for phase w within 2 voxels of solid
|
||||
* ``Psn_y`` -- y- momentum for phase n within 2 voxels of solid
|
||||
* ``Pwc_z`` -- z- momentum for connected part of phase w
|
||||
* ``Pwd_z`` -- z- momentum for disconnected part of phase w
|
||||
* ``Pwi_z`` -- z- momentum for phase w within diffuse interface
|
||||
* ``Pnc_z`` -- z- momentum for connected part of phase n
|
||||
* ``Pnd_z`` -- z- momentum for disconnected part of phase n
|
||||
* ``Pni_z`` -- z- momentum for phase n within diffuse interface
|
||||
* ``Psw_z`` -- z- momentum for phase w within 2 voxels of solid
|
||||
* ``Psn_z`` -- z- momentum for phase n within 2 voxels of solid
|
||||
* ``Kwc`` -- Kinetic energy for transport within connected part of phase w
|
||||
* ``Kwd`` -- Kinetic energy for transport within disconnected part of phase w
|
||||
* ``Kwi`` -- Kinetic energy for transport of phase w within diffuse interface region
|
||||
* ``Knc`` -- Kinetic energy for transport in connected part of phase n
|
||||
* ``Knd`` -- Kinetic energy for transport within disconnected part of phase n
|
||||
* ``Kni`` -- Kinetic energy for transport of phase n within diffuse interface region
|
||||
* ``Dwc`` -- Viscous dissipation for conneced pathway for phase w
|
||||
* ``Dwd`` -- Viscous dissipation for disconnected part of phase w
|
||||
* ``Dnc`` -- Viscous dissipation for connected pathway for phase n
|
||||
* ``Dnd`` -- Viscous dissipation for disconnected part of phase n
|
||||
* ``Vwc`` -- Volume for connected pathway for phase w
|
||||
* ``Awc`` -- Surface area for connected pathway for phase w
|
||||
* ``Hwc`` -- Integral mean curvature for connected pathway for phase w
|
||||
* ``Xwc`` -- Euler characteristic for connected pathway for phase w
|
||||
* ``Vwd`` -- Volume for disconnected phase w
|
||||
* ``Awd`` -- Surface area for disconnected phase w
|
||||
* ``Hwd`` -- Integral mean curvature for disconnected phase w
|
||||
* ``Xwd`` -- Euler characteristic for disconnected phase w
|
||||
* ``Nwd`` -- Number of connected components in disconnected phase w
|
||||
* ``Vnc`` -- Volume for connected pathway for phase n
|
||||
* ``Anc`` -- Surface area for connected pathway for phase n
|
||||
* ``Hnc`` -- Integral mean curvature for connected pathway for phase n
|
||||
* ``Xnc`` -- Euler characteristic for connected pathway for phase n
|
||||
* ``Vnd`` -- Volume for disconnected phase n
|
||||
* ``And`` -- Surface area for disconnected phase n
|
||||
* ``Hnd`` -- Integral mean curvature for disconnected phase n
|
||||
* ``Xnd`` -- Euler characteristic for disconnected phase n
|
||||
* ``Nnd`` -- number of connected components within disconnected phase n
|
||||
* ``Vi`` -- volume for diffuse interface region
|
||||
* ``Ai`` -- surface area for boundary of diffuse interface region
|
||||
* ``Hi`` -- integral mean curvature for boundary of diffuse interface region
|
||||
* ``Xi`` -- Euler characteristic for diffuse interface region
|
||||
* ``Vic`` -- volume for connected interface region
|
||||
* ``Aic`` -- surface area for boundary of connected interface region
|
||||
* ``Hic`` -- Integral mean curvature for connected interface region
|
||||
* ``Xic`` -- Euler characteristic for connected interface region
|
||||
* ``Nic`` -- number of connected components in connected interface region
|
|
@ -37,6 +37,16 @@ that can be over-ridden to develop customized simulations.
|
|||
|
||||
protocols/*
|
||||
|
||||
****************************
|
||||
Analysis capabilities
|
||||
****************************
|
||||
|
||||
.. toctree::
|
||||
:glob:
|
||||
:maxdepth: 2
|
||||
|
||||
analysis/*
|
||||
|
||||
|
||||
***************************
|
||||
Model parameters
|
||||
|
@ -44,10 +54,12 @@ Model parameters
|
|||
|
||||
The essential model parameters for the color model are
|
||||
|
||||
- :math:`\alpha` -- control the interfacial tension between fluids with key ``alpha``
|
||||
- :math:`\beta` -- control the width of the interface with key ``beta``
|
||||
- :math:`\tau_A` -- control the viscosity of fluid A with key ``tauA``
|
||||
- :math:`\tau_B` -- control the viscosity of fluid B with key ``tauB``
|
||||
- ``alpha`` -- control the interfacial tension between fluids with :math:`0 < \alpha < 0.01`
|
||||
- ``beta`` -- control the width of the interface with key :math:`\beta < 1`
|
||||
- ``tauA`` -- control the viscosity of fluid A with :math:`0.7 < \tau_A < 1.5`
|
||||
- ``tauB`` -- control the viscosity of fluid B with :math:`0.7 < \tau_B < 1.5`
|
||||
- ``rhoA`` -- control the viscosity of fluid A with :math:`0.05 < \rho_A < 1.0`
|
||||
- ``rhoB`` -- control the viscosity of fluid B with :math:`0.05 < \rho_B < 1.0`
|
||||
|
||||
****************************
|
||||
Model Formulation
|
||||
|
|
35
example/DropletCoalescence/Droplets.py
Normal file
35
example/DropletCoalescence/Droplets.py
Normal file
|
@ -0,0 +1,35 @@
|
|||
import numpy
|
||||
import math
|
||||
|
||||
nx=400
|
||||
ny=200
|
||||
nz=200
|
||||
N=nx*ny*nz
|
||||
Radius=64
|
||||
|
||||
mesh=(nx,ny,nz)
|
||||
data=numpy.ones(mesh,dtype=numpy.int8)
|
||||
|
||||
#print(data)
|
||||
print("Create two droplets")
|
||||
print("Mesh size: "+repr(mesh))
|
||||
print("Droplet radius: "+repr(Radius))
|
||||
|
||||
gap = 6
|
||||
c1x = nx/2 - gap/2 - Radius
|
||||
c2x = nx/2 + gap/2 + Radius
|
||||
|
||||
# assign a bubble on each side
|
||||
for x in range(0,200):
|
||||
for y in range(0,ny):
|
||||
for z in range(0,nz):
|
||||
if math.sqrt((x-c1x)*(x-c1x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
|
||||
data[x,y,z]=2
|
||||
|
||||
for x in range(200,nx):
|
||||
for y in range(0,ny):
|
||||
for z in range(0,nz):
|
||||
if math.sqrt((x-c2x)*(x-c2x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
|
||||
data[x,y,z]=2
|
||||
|
||||
data.tofile("Droplets.raw")
|
51
example/DropletCoalescence/input.db
Normal file
51
example/DropletCoalescence/input.db
Normal file
|
@ -0,0 +1,51 @@
|
|||
MRT {
|
||||
timestepMax = 10000
|
||||
tau = 0.7
|
||||
F = 1e-05, 0, 0
|
||||
Restart = false
|
||||
din = 1.0
|
||||
dout = 1.0
|
||||
flux = 0.0
|
||||
}
|
||||
|
||||
Color {
|
||||
tauA = 0.7;
|
||||
tauB = 1.0;
|
||||
rhoA = 1.0;
|
||||
rhoB = 1.0;
|
||||
alpha = 5e-3;
|
||||
beta = 0.95;
|
||||
F = 0, 0, 0
|
||||
Restart = false
|
||||
timestepMax = 40000
|
||||
ComponentLabels = -2
|
||||
ComponentAffinity = -0.5
|
||||
}
|
||||
|
||||
Domain {
|
||||
Filename = "Droplets.raw"
|
||||
nproc = 1, 1, 2 // Number of processors (Npx,Npy,Npz)
|
||||
n = 200, 200, 200 // Size of local domain (Nx,Ny,Nz)
|
||||
N = 200, 200, 400 // size of the input image
|
||||
voxel_length = 1.0
|
||||
BC = 0 // Boundary condition type
|
||||
Sw = 0.15
|
||||
ReadType = "8bit"
|
||||
ReadValues = -2, 1, 2
|
||||
WriteValues = -2, 1, 2
|
||||
ComponentLabels = -2
|
||||
HistoryLabels = -2
|
||||
}
|
||||
|
||||
Analysis {
|
||||
analysis_interval = 1000 // Frequency to perform analysis
|
||||
subphase_analysis_interval = 5000 // Frequency to perform analysis
|
||||
restart_interval = 60000 // Frequency to write restart data
|
||||
visualization_interval = 100000 // Frequency to write visualization 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 {
|
||||
}
|
|
@ -34,15 +34,15 @@ int main( int argc, char **argv )
|
|||
// Load the input database
|
||||
auto db = std::make_shared<Database>( argv[1] );
|
||||
if (argc > 2) {
|
||||
SimulationMode = "development";
|
||||
SimulationMode = "legacy";
|
||||
}
|
||||
|
||||
if ( rank == 0 ) {
|
||||
printf( "********************************************************\n" );
|
||||
printf( "Running Color LBM \n" );
|
||||
printf( "********************************************************\n" );
|
||||
if (SimulationMode == "development")
|
||||
printf("**** DEVELOPMENT MODE ENABLED *************\n");
|
||||
if (SimulationMode == "legacy")
|
||||
printf("**** LEGACY MODE ENABLED *************\n");
|
||||
}
|
||||
// Initialize compute device
|
||||
int device = ScaLBL_SetDevice( rank );
|
||||
|
@ -66,13 +66,16 @@ int main( int argc, char **argv )
|
|||
// structure and allocate variables
|
||||
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
|
||||
|
||||
if (SimulationMode == "development"){
|
||||
if (SimulationMode == "legacy"){
|
||||
ColorModel.Run();
|
||||
}
|
||||
else {
|
||||
double MLUPS=0.0;
|
||||
int timestep = 0;
|
||||
bool ContinueSimulation = true;
|
||||
|
||||
/* Variables for simulation protocols */
|
||||
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "none" );
|
||||
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "default" );
|
||||
/* image sequence protocol */
|
||||
int IMAGE_INDEX = 0;
|
||||
int IMAGE_COUNT = 0;
|
||||
|
@ -123,6 +126,7 @@ int main( int argc, char **argv )
|
|||
else{
|
||||
if (rank==0) printf("Finished simulating image sequence \n");
|
||||
ColorModel.timestep = ColorModel.timestepMax;
|
||||
ContinueSimulation = false;
|
||||
}
|
||||
}
|
||||
/*********************************************************/
|
||||
|
@ -144,7 +148,7 @@ int main( int argc, char **argv )
|
|||
double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
|
||||
/* stop simulation if previous point was sufficiently close to the endpoint*/
|
||||
if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
|
||||
if (ContinueSimulation){
|
||||
if (ContinueSimulation && SKIP_TIMESTEPS > 0 ){
|
||||
while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
|
||||
timestep += ANALYSIS_INTERVAL;
|
||||
if (PROTOCOL == "fractional flow") {
|
||||
|
@ -173,9 +177,8 @@ int main( int argc, char **argv )
|
|||
/*********************************************************/
|
||||
}
|
||||
}
|
||||
else
|
||||
ColorModel.Run();
|
||||
|
||||
|
||||
|
||||
PROFILE_STOP( "Main" );
|
||||
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
|
||||
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
|
||||
|
|
Loading…
Reference in New Issue
Block a user