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

docs/Makefile

@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = source
+BUILDDIR      = $(HOME)/local/doc/build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

docs/README.md

@@ -0,0 +1,19 @@
+Dependencies for LBPM documentation
+
+# install sphinx
+pip install Sphinx
+
+# foamatting requires sphinx read-the-docs-theme
+pip install sphinx-rtd-theme
+
+# equation rendering requires latex and dvipng command
+sudo apt-get install dvipng
+sudo apt-get install texlive texstudio
+sudo apt-get install texlive-latex-recommended texlive-pictures texlive-latex-extra
+
+
+# To build the docs
+Step 1) install dependencies listed above
+Step 2) type 'make html' from the docs/ directory
+Step 3) point your browser at ~/local/doc/build/html/index.html
+# 

docs/source/conf.py

@@ -0,0 +1,70 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# This file only contains a selection of the most common options. For a full
+# list see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+import os
+import sys
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'LBPM'
+copyright = '2021, James E McClure'
+author = 'James E McClure'
+
+# The full version, including alpha/beta/rc tags
+release = '1.0'
+
+
+# -- General configuration ---------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.imgmath'
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = []
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'alabaster'
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+
+## Read the docs style:
+if os.environ.get('READTHEDOCS') != 'True':
+    try:
+        import sphinx_rtd_theme
+    except ImportError:
+        pass  # assume we have sphinx >= 1.3
+    else:
+        html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
+    html_theme = 'sphinx_rtd_theme'
+    
+#def setup(app):
+#    app.add_stylesheet("fix_rtd.css")

docs/source/developerGuide/buildingModels/overview.rst

@@ -0,0 +1,20 @@
+===========================
+Implementing a new LB model
+===========================
+
+While LBPM includes a range of fully-functioning lattice Boltzmann models, the commonly used
+Bhatnager-Gross-Krook (BGK) model has been deliberately excluded. While the physical limitations
+of this model are well-known, implementing the BGK model is an excellent way to understand
+how to implement new LB models within the more general framework of LBPM. In this excercise
+you will
+
+* learn "what goes where"
+
+
+
+* don't modify core data structures (unless you have a really good reason)
+
+* re-use existing components whenever possible
+
+
+

docs/source/developerGuide/designOverview/dataStructures.rst

@@ -0,0 +1,6 @@
+===============
+Data Structures
+===============
+
+LBPM includes a variety of generalized data structures to facilitate the implementation
+of different lattice Boltzmann models. 

docs/source/developerGuide/index.rst

@@ -0,0 +1,18 @@
+###############################################################################
+Developer guide
+###############################################################################
+
+The LBPM developer guide provides essential information on how to add new physics
+into the framework. 
+
+.. toctree::
+   :glob:
+   :maxdepth: 2
+
+   designOverview/*
+   
+   buildingModels/*
+
+   testingModels/*
+
+ 

docs/source/developerGuide/testingModels/unitTests.rst

@@ -0,0 +1,9 @@
+=================
+Adding unit tests
+=================
+
+Unit tests in LBPM are implemented using ctest
+
+* general overview
+
+* launching unit tests for GPU (MPI flags etc.)

docs/source/examples/running.rst

@@ -0,0 +1,9 @@
+============
+Running LBPM
+============
+
+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. 

docs/source/index.rst

@@ -0,0 +1,25 @@
+.. LBPM documentation master file, created by
+   sphinx-quickstart on Thu May 20 12:19:14 2021.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+LBPM -- Documentation
+===================================================
+
+.. toctree::
+   :glob:
+   :maxdepth: 2
+   :caption: Contents:
+
+   install
+   examples/*
+   userGuide/*
+   developerGuide/*
+   publications/*
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`

docs/source/install.rst

@@ -0,0 +1,10 @@
+============
+Installation
+============
+
+LBPM requires CMake, MPI and HDF5 as required dependencies. 
+
+.. code-block:: bash
+
+    ls $LBPM_SOURCE/sample_scripts
+    

docs/source/publications/publications.rst

@@ -0,0 +1,13 @@
+============
+Publications
+============
+
+* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
+
+
+* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
+
+
+* Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966
+
+ 

docs/source/userGuide/IO/fileformat.rst

@@ -0,0 +1,13 @@
+========================
+I/O conventions for LBPM
+========================
+
+There are three main kinds of output file that are supported by LBPM.
+
+
+* CSV files --
+
+* formatted binary files --
+  
+* unformatted binary files --
+

docs/source/userGuide/analysis/analysis.rst

@@ -0,0 +1,5 @@
+===========================
+Internal Analysis Framework
+===========================
+
+placeholder for analysis 

docs/source/userGuide/index.rst

@@ -0,0 +1,17 @@
+###############################################################################
+User Guide
+###############################################################################
+
+Welcome to the LBPM user guide.
+
+.. toctree::
+   :glob:
+   :maxdepth: 2
+
+   models/*
+
+   analysis/*
+
+   visualization/*
+
+   IO/*

docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst

@@ -0,0 +1,6 @@
+=============================================
+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. 

docs/source/userGuide/models/color/index.rst

@@ -0,0 +1,91 @@
+###############################################################################
+Color model
+###############################################################################
+
+The LBPM color model is implemented by combining a multi-relaxation time D3Q19
+lattice Boltzmann equation (LBE) to solve for the momentum transport with two D3Q7
+LBEs for the mass transport. The color model will obey strict mass and momentum
+conservation while minimizing diffusive fluxes across the interface between fluids.
+The color model is a good choice for modeling dense fluids that are strongly immiscible
+(e.g. water-oil systems). Due to the strong anti-diffusion in the interface region,
+the color model is not suitable for modeling processes such as Ostwald ripening that
+depend on diffusive fluxes between fluid phases.
+
+A typical command to launch the LBPM color simulator is as follows
+
+```
+mpirun -np $NUMPROCS lbpm_color_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.
+
+****************************
+Simulation protocols
+****************************
+
+Simulation protocols are designed to make it simpler to design and execute common
+computational experiments. Protocols will automatically determine boundary conditions
+needed to perform a particular simulation. LBPM will internall set default simulation paramaters
+that can be over-ridden to develop customized simulations.
+
+.. toctree::
+   :glob:
+   :maxdepth: 2
+
+   protocols/*
+
+
+***************************
+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``
+
+****************************
+Model Formulation
+****************************
+
+
+  
+****************************
+Boundary Conditions
+****************************
+
+The following external boundary conditions are supported by ``lbpm_color_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
+

docs/source/userGuide/models/color/protocols/centrifuge.rst

@@ -0,0 +1,12 @@
+======================================
+Color model -- Centrifuge Protocol
+======================================
+
+The centrifuge protocol is designed to mimic SCAL centrifuge experiments that
+are used to infer the capillary pressure. The centrifuge protocol 
+
+.. code-block:: bash
+
+    image_sequence = "image1.raw", "image2.raw"
+    
+

docs/source/userGuide/models/color/protocols/coreFlooding.rst

@@ -0,0 +1,63 @@
+======================================
+Color model -- Core Flooding
+======================================
+
+The core flooding protocol is designed to mimic SCAL experiments where one
+immiscible fluid is injected into the sample at a constant rate, displacing the
+other fluid. The core flooding protocol relies on a flux boundary condition
+to ensure that fluid is injected into the sample at a constant rate. The flux
+boundary condition implements a time-varying pressure boundary condition that
+adapts to ensure a constant volumetric flux. Details for the flux boundary
+condition are available
+(see: https://doi.org/10.1016/j.compfluid.2020.104670)
+
+.. code-block:: bash
+
+    protocol = "core flooding"
+    
+
+To match experimental conditions, it is usually important to match the capillary
+number, which is
+
+.. math::
+   \mbox{Ca} = \frac{\mu u_z}{\gamma}
+
+
+where :math:`\mu` is the dynamic viscosity, :math:`u_z` is the fluid
+(usually water) velocity and :math:`\gamma` is the interfacial tension.
+The volumetric flow rate is related to the fluid velocity based on
+
+.. math::
+   Q_z = \epsilon C_{xy} u_z
+
+where :math:`C_{xy}` is the cross-sectional area and :math:`\epsilon`
+is the porosity. Given a particular experimental system 
+self-similar conditions can be determined for the lattice Boltzmann
+simulation by matching the non-dimensional :math:`mbox{Ca}`. It is nearly
+awlays advantageous for the timestep to be as large as possible so
+that time-to-solution will be more favorable. This is accomplished by
+
+* use a high value for the numerical surface tension (e.g. ``alpha=1.0e-2``)
+* use a small value for the fluid viscosity (e.g. ``tau_w = 0.7`` and ``tau_n = 0.7`` )
+* determine the volumetric flow rate needed to match :math:`\mbox{Ca}`
+
+For the color LBM the interfacial tension is
+:math:`\gamma = 6 \alpha` and the dynamic viscosity is :math:`\mu =  \rho(\tau-1/2)/3`,
+where the units are relative to the lattice spacing, timestep and mass
+density. Agreemetn between the experimental and simulated values for
+:math:`\mbox{Ca}` is ensured by setting the volumetric flux
+
+.. math::
+   Q_z = \frac{\epsilon C_{xy} \gamma }{\mu} \mbox{Ca}
+
+where the LB units of the volumetric flux will be voxels per timestep.
+
+In some situations it may also be important to match other non-dimensional numbers,
+such as the viscosity ratio, density ratio, and/or Ohnesorge/Laplace number. This
+can be accomplished with an analogous procedure. Enforcing additional constraints
+will necessarily restrict the LB parameters that can be used, which are ultimately
+manifested as a constraint on the size of a timestep. 
+
+
+
+

docs/source/userGuide/models/color/protocols/fractionalFlow.rst

@@ -0,0 +1,17 @@
+==========================================
+Color model -- Fractional Flow Protocol
+==========================================
+
+The fractional flow protocol is designed to perform steady-state relative 
+permeability simulations by using an internal routine to change the fluid
+saturation by adding and subtracting mass to the fluid phases. The 
+mass density is updated for each fluid based on the locations where
+the local rate of flow is high. 
+
+.. code-block:: bash
+
+	protocol = "fractional flow"
+
+    
+    
+

docs/source/userGuide/models/color/protocols/imageSequence.rst

@@ -0,0 +1,27 @@
+======================================
+Color model -- Image Sequence Protocol
+======================================
+
+The image sequence protocol is designed to perform a set steady-state
+simulations based on a sequence of 3D (8-bit) images provided by the user.
+The images might be the output of a previous LBPM simulation, a sequence of
+(segmented) experimental data, or data generated from a custom routine.
+The image sequence protocol will apply the same set of flow conditions
+to all images in the sequence. This means
+
+* the image labels and any associated properties are the same
+* the external boundary conditions are the same
+* the physical simulation parameters are the same
+
+The image sequence protocol does not set boundary conditions by default.
+It is up to the user to determine the flow condition, with the understanding
+that the same set of will be applied to each image in the sequence.
+
+To enable the image sequence protocol, the following keys should be set
+within the ``Color`` section of the input database
+
+.. code-block:: bash
+
+   protocol = "image sequence"
+   image_sequence = "image1.raw", "image2.raw"
+    

docs/source/userGuide/models/color/protocols/shellAggregation.rst

@@ -0,0 +1,18 @@
+==========================================
+Color model -- Shell Aggregation Protocol
+==========================================
+
+The shell aggregation protocol is designed to perform steady-state relative 
+permeability simulations by using an internal routine to change the fluid
+saturation by moving the interface. The basic design for the shell aggregation
+protocol is described by Wang et al. ( https://doi.org/10.1016/j.jcp.2019.108966 ).
+
+
+.. code-block:: bash
+
+
+	protocol = "shell aggregation"
+
+    
+    
+

docs/source/userGuide/models/freeEnergy/freeEnergy.rst

@@ -0,0 +1,6 @@
+=============================================
+Free energy model
+=============================================
+
+The LBPM free energy model is constructed to solve the Allen-Cahn equations,
+which are typically used to describe liquid-gas systems. 

docs/source/userGuide/models/greyscale/greyscale.rst

@@ -0,0 +1,7 @@
+=============================================
+Greyscale model model
+=============================================
+
+The LBPM greyscale lattice Boltzmann model is constructed to approximate the
+solution of the Darcy-Brinkman equations in grey regions, coupled to a Navier-Stokes
+solution in open regions. 

docs/source/userGuide/models/index.rst

@@ -0,0 +1,23 @@
+###############################################################################
+LBPM model summary
+###############################################################################
+
+Currently supported lattice Boltzmann models
+
+.. toctree::
+   :glob:
+   :maxdepth: 2
+
+   color/*
+
+   mrt/*
+
+   nernstPlanck/*
+
+   PoissonBoltzmann/*
+   
+   greyscale/*
+
+   freeEnergy/*
+
+   

docs/source/userGuide/models/mrt/mrt.rst

@@ -0,0 +1,6 @@
+=============================================
+MRT model
+=============================================
+
+The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the
+solution of the Navier-Stokes equations.

docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst

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

docs/source/userGuide/visualization/visit.rst

@@ -0,0 +1,5 @@
+======================================
+Visualizing simulation data with visit
+======================================
+
+placeholder for visit

models/ColorModel.cpp

@@ -1997,188 +1997,197 @@ FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M){
 	Nz = M.Dm->Nz;
 	timestep=-1;
 	timestep_previous=-1;	
-	
+
 	phi.resize(Nx,Ny,Nz);         phi.fill(0);	    // phase indicator field
 	phi_t.resize(Nx,Ny,Nz);       phi_t.fill(0);	// time derivative for the phase indicator field
 }
 
 FlowAdaptor::~FlowAdaptor(){
-	
+
 }
 
-double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){	
+double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename){
+	int rank = M.rank;
+	int Nx = M.Nx;  int Ny = M.Ny;  int Nz = M.Nz;
+	if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
+	M.Mask->Decomp(Filename);
+	for (int i=0; i<Nx*Ny*Nz; i++) M.id[i] = M.Mask->id[i];  // save what was read
+	for (int i=0; i<Nx*Ny*Nz; i++) M.Dm->id[i] = M.Mask->id[i];  // save what was read
+
+	double *PhaseLabel;
+	PhaseLabel = new double[Nx*Ny*Nz];
+	M.AssignComponentLabels(PhaseLabel);
+
+	double Count = 0.0;
+	double PoreCount = 0.0;
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				if (M.id[Nx*Ny*k+Nx*j+i] == 2){
+					PoreCount++;
+					Count++;
+				}
+				else if (M.id[Nx*Ny*k+Nx*j+i] == 1){
+					PoreCount++;						
+				}
+			}
+		}
+	}
+
+	Count=M.Dm->Comm.sumReduce(  Count);
+	PoreCount=M.Dm->Comm.sumReduce(  PoreCount);
 	
-  	  double MASS_FRACTION_CHANGE = 0.0005;
-  	  if (M.db->keyExists( "FlowAdaptor" )){
-  		  auto flow_db = M.db->getDatabase( "FlowAdaptor" );
-  		  MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.0005);
-  	  }
-	  int Np = M.Np;
-	  double dA, dB, phi;
-	  double vx,vy,vz;
-	  double vax,vay,vaz;
-	  double vbx,vby,vbz;
-	  double vax_global,vay_global,vaz_global;
-	  double vbx_global,vby_global,vbz_global;
+	if (rank==0) printf("   new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
+	ScaLBL_CopyToDevice(M.Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
+	M.Dm->Comm.barrier();
+	
+	ScaLBL_D3Q19_Init(M.fq, M.Np);
+	ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np);
+	ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np);
+	M.Dm->Comm.barrier();
+	
+	ScaLBL_CopyToHost(M.Averages->Phi.data(),M.Phi,Nx*Ny*Nz*sizeof(double));
 
-	  double mass_a, mass_b, mass_a_global, mass_b_global;
-	  
-	  double *Aq_tmp, *Bq_tmp;
-	  double *Vel_x, *Vel_y, *Vel_z, *Phase;
-	  
-	  Aq_tmp = new double [7*Np];
-	  Bq_tmp = new double [7*Np];
-	  Phase = new double [Np];
-	  Vel_x = new double [Np];
-	  Vel_y = new double [Np];
-	  Vel_z = new double [Np];
+	double saturation = Count/PoreCount;
+	return saturation;
+}
 
-	  ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
-	  ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
-	  ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double));
-	  ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
-	  ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
-	  
-	  /* DEBUG STRUCTURES */
-	  int Nx = M.Nx;  int Ny = M.Ny;  int Nz = M.Nz;
-	  int N = Nx*Ny*Nz;
-	  double * DebugMassA, *DebugMassB;
-	  DebugMassA = new double[Np];
-	  DebugMassB = new double[Np];
-	  
-	  /* compute the total momentum */
-	  vax = vay = vaz = 0.0;
-	  vbx = vby = vbz = 0.0;
-	  mass_a = mass_b = 0.0;
-	  
-	  double maxSpeed = 0.0;
-	  double localMaxSpeed = 0.0;
-	  for (int k=1; k<Nz-1; k++){
-		  for (int j=1; j<Ny-1; j++){
-			  for (int i=1; i<Nx-1; i++){
-				  int n=M.Map(i,j,k);
-				  double distance = M.Averages->SDs(i,j,k);
-				  if (!(n<0) ){
-					  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
-					  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
-					  phi = (dA - dB) / (dA + dB);
-					  Phase[n] = phi;
-					  mass_a += dA;
-					  mass_b += dB;
-					  if (phi > 0.0){
-						  vax += Vel_x[n];
-						  vay += Vel_y[n];
-						  vaz += Vel_z[n];
-					  }
-					  else {
-						  vbx += Vel_x[n];
-						  vby += Vel_y[n];
-						  vbz += Vel_z[n];
-					  }
-					  double speed =  sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz);
-					  if (distance > 3.0 && speed > localMaxSpeed){
-						  localMaxSpeed = speed;
-					  }
-				  }
-			  }
-		  }
-	  }
-	  maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
-	  /*	  for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
-		  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
-		  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
-		  phi = (dA - dB) / (dA + dB);
-		  Phase[n] = phi;
-		  mass_a += dA;
-		  mass_b += dB;
-		  if (phi > 0.0){
-			  vax += Vel_x[n];
-			  vay += Vel_y[n];
-			  vaz += Vel_z[n];
-		  }
-		  else {
-			  vbx += Vel_x[n];
-			  vby += Vel_y[n];
-			  vbz += Vel_z[n];
-		  }
-	  }
-	  for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
-		    dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
-		    dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
-		    phi = (dA - dB) / (dA + dB);
-		    Phase[n] = phi;
-			  mass_a += dA;
-			  mass_b += dB;
-			  if (phi > 0.0){
-				  vax += Vel_x[n];
-				  vay += Vel_y[n];
-				  vaz += Vel_z[n];
-			  }
-			  else {
-				  vbx += Vel_x[n];
-				  vby += Vel_y[n];
-				  vbz += Vel_z[n];
-			  }
-	   }
-	   */
-	  mass_a_global = M.Dm->Comm.sumReduce(mass_a);
-	  mass_b_global = M.Dm->Comm.sumReduce(mass_b);
-	  vax_global = M.Dm->Comm.sumReduce(vax);
-	  vay_global = M.Dm->Comm.sumReduce(vay);
-	  vaz_global = M.Dm->Comm.sumReduce(vaz);
-	  vbx_global = M.Dm->Comm.sumReduce(vbx);
-	  vby_global = M.Dm->Comm.sumReduce(vby);
-	  vbz_global = M.Dm->Comm.sumReduce(vbz);	 
 
-	  double  total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
-	  double  total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
-	  /* compute the total mass change */
-	  double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
-	  if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
-		  TOTAL_MASS_CHANGE = 0.1*mass_a_global;
-	  if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
-		  TOTAL_MASS_CHANGE = 0.1*mass_b_global;
+double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){	
 
-	  double LOCAL_MASS_CHANGE = 0.0;
-	  for (int k=1; k<Nz-1; k++){
-		  for (int j=1; j<Ny-1; j++){
-			  for (int i=1; i<Nx-1; i++){
-				  int n=M.Map(i,j,k);
-				  if (!(n<0)){
-					  phi = Phase[n];
-					  vx = Vel_x[n];
-					  vy = Vel_y[n];
-					  vz = Vel_z[n];
-					  double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
-					  /* impose ceiling for spurious currents */
-					  if (local_momentum > maxSpeed) local_momentum =  maxSpeed;
-					  if (phi > 0.0){
-						  LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
-						  Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
-					  }
-					  else{
-						  LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
-						  Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-						  DebugMassB[n] = LOCAL_MASS_CHANGE;
-					  }
-				  }
-			  }
-		  }
-	  }
-/*	   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){	    
+	double MASS_FRACTION_CHANGE = 0.0005;
+	if (M.db->keyExists( "FlowAdaptor" )){
+		auto flow_db = M.db->getDatabase( "FlowAdaptor" );
+		MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.0005);
+	}
+	int Np = M.Np;
+	double dA, dB, phi;
+	double vx,vy,vz;
+	double vax,vay,vaz;
+	double vbx,vby,vbz;
+	double vax_global,vay_global,vaz_global;
+	double vbx_global,vby_global,vbz_global;
+
+	double mass_a, mass_b, mass_a_global, mass_b_global;
+
+	double *Aq_tmp, *Bq_tmp;
+	double *Vel_x, *Vel_y, *Vel_z, *Phase;
+
+	Aq_tmp = new double [7*Np];
+	Bq_tmp = new double [7*Np];
+	Phase = new double [Np];
+	Vel_x = new double [Np];
+	Vel_y = new double [Np];
+	Vel_z = new double [Np];
+
+	ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
+	ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
+	ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double));
+	ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
+	ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
+
+	/* DEBUG STRUCTURES */
+	int Nx = M.Nx;  int Ny = M.Ny;  int Nz = M.Nz;
+	//int N = Nx*Ny*Nz;
+	//double * DebugMassA, *DebugMassB;
+	//DebugMassA = new double[Np];
+	//DebugMassB = new double[Np];
+
+	/* compute the total momentum */
+	vax = vay = vaz = 0.0;
+	vbx = vby = vbz = 0.0;
+	mass_a = mass_b = 0.0;
+	double maxSpeed = 0.0;
+	double localMaxSpeed = 0.0;
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int n=M.Map(i,j,k);
+				double distance = M.Averages->SDs(i,j,k);
+				if (!(n<0) ){
+					dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+					dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+					phi = (dA - dB) / (dA + dB);
+					Phase[n] = phi;
+					mass_a += dA;
+					mass_b += dB;
+					if (phi > 0.0){
+						vax += Vel_x[n];
+						vay += Vel_y[n];
+						vaz += Vel_z[n];
+					}
+					else {
+						vbx += Vel_x[n];
+						vby += Vel_y[n];
+						vbz += Vel_z[n];
+					}
+					double speed =  sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz);
+					if (distance > 3.0 && speed > localMaxSpeed){
+						localMaxSpeed = speed;
+					}
+				}
+			}
+		}
+	}
+	maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
+	mass_a_global = M.Dm->Comm.sumReduce(mass_a);
+	mass_b_global = M.Dm->Comm.sumReduce(mass_b);
+	vax_global = M.Dm->Comm.sumReduce(vax);
+	vay_global = M.Dm->Comm.sumReduce(vay);
+	vaz_global = M.Dm->Comm.sumReduce(vaz);
+	vbx_global = M.Dm->Comm.sumReduce(vbx);
+	vby_global = M.Dm->Comm.sumReduce(vby);
+	vbz_global = M.Dm->Comm.sumReduce(vbz);	 
+
+	double  total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
+	double  total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
+	/* compute the total mass change */
+	double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
+	if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
+		TOTAL_MASS_CHANGE = 0.1*mass_a_global;
+	if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
+		TOTAL_MASS_CHANGE = 0.1*mass_b_global;
+
+	double LOCAL_MASS_CHANGE = 0.0;
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int n=M.Map(i,j,k);
+				if (!(n<0)){
+					phi = Phase[n];
+					vx = Vel_x[n];
+					vy = Vel_y[n];
+					vz = Vel_z[n];
+					double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
+					/* impose ceiling for spurious currents */
+					if (local_momentum > maxSpeed) local_momentum =  maxSpeed;
+					if (phi > 0.0){
+						LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
+						Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						//DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
+					}
+					else{
+						LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
+						Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						//DebugMassB[n] = LOCAL_MASS_CHANGE;
+					}
+				}
+			}
+		}
+	}
+	/*	   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){	    
 		   phi = Phase[n];
 		   vx = Vel_x[n];
 		   vy = Vel_y[n];
@@ -2237,10 +2246,10 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			   DebugMassB[n] = LOCAL_MASS_CHANGE;
 		   }
 	  }
- 	 */
-	  if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
+	 */
+	if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
 
-	  /* Print out debugging info with mass update  */
+	/* Print out debugging info with mass update  
 	  // initialize the array
 	  double value;
 	  char LocalRankFilename[40];
@@ -2257,7 +2266,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			  }
 		  }
 	  }
-	  
+
 	  FILE *AFILE;
 	  sprintf(LocalRankFilename,"dA.%05i.raw",M.rank);
 	  AFILE = fopen(LocalRankFilename,"wb");
@@ -2281,13 +2290,14 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 	  BFILE = fopen(LocalRankFilename,"wb");
 	  fwrite(regdata.data(),8,N,BFILE);
 	  fclose(BFILE);
+	 */
 
-	  // Need to initialize Aq, Bq, Den, Phi directly
-	  //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
-	  ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
-	  ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
+	// Need to initialize Aq, Bq, Den, Phi directly
+	//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
+	ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
+	ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
 
-	  return(TOTAL_MASS_CHANGE);
+	return(TOTAL_MASS_CHANGE);
 }
 
 void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){	

models/ColorModel.h

@@ -88,6 +88,8 @@ public:
 	double *ColorGrad;
 	double *Velocity;
 	double *Pressure;
+
+	void AssignComponentLabels(double *phase);
 		
 private:
 	Utilities::MPI comm;
@@ -101,7 +103,6 @@ private:
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);
-    void AssignComponentLabels(double *phase);
     double ImageInit(std::string filename);
     double MorphInit(const double beta, const double morph_delta);
     double SeedPhaseField(const double seed_water_in_oil);
@@ -113,6 +114,7 @@ public:
 	FlowAdaptor(ScaLBL_ColorModel &M);
 	~FlowAdaptor();
 	double MoveInterface(ScaLBL_ColorModel &M);
+	double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
 	double UpdateFractionalFlow(ScaLBL_ColorModel &M);
 	void Flatten(ScaLBL_ColorModel &M);
 	DoubleArray phi;

tests/lbpm_color_simulator.cpp

@@ -22,114 +22,138 @@
 int main( int argc, char **argv )
 {
 
-    // Initialize
-    Utilities::startup( argc, argv );
+	// Initialize
+	Utilities::startup( argc, argv );
 
-    { // Limit scope so variables that contain communicators will free before MPI_Finialize
+	{ // Limit scope so variables that contain communicators will free before MPI_Finialize
 
-        Utilities::MPI comm( MPI_COMM_WORLD );
-        int rank   = comm.getRank();
-        int nprocs = comm.getSize();
-        std::string SimulationMode = "production";
-        // Load the input database
-        auto db = std::make_shared<Database>( argv[1] );
-        if (argc > 2) {
-        	SimulationMode = "development";
-        }
-
-        if ( rank == 0 ) {
-            printf( "********************************************************\n" );
-            printf( "Running Color LBM	\n" );
-            printf( "********************************************************\n" );
-            if (SimulationMode == "development")
-            	printf("**** DEVELOPMENT MODE ENABLED *************\n");
-        }
-        // Initialize compute device
-        int device = ScaLBL_SetDevice( rank );
-        NULL_USE( device );
-        ScaLBL_DeviceBarrier();
-        comm.barrier();
-
-        PROFILE_ENABLE( 1 );
-        // PROFILE_ENABLE_TRACE();
-        // PROFILE_ENABLE_MEMORY();
-        PROFILE_SYNCHRONIZE();
-        PROFILE_START( "Main" );
-        Utilities::setErrorHandlers();
-
-        auto filename = argv[1];
-        ScaLBL_ColorModel ColorModel( rank, nprocs, comm );
-        ColorModel.ReadParams( filename );
-        ColorModel.SetDomain();
-        ColorModel.ReadInput();
-        ColorModel.Create();     // creating the model will create data structure to match the pore
-                                 // structure and allocate variables
-        ColorModel.Initialize(); // initializing the model will set initial conditions for variables
-       
-        if (SimulationMode == "development"){
-        	double MLUPS=0.0;
-        	int timestep = 0;
-		bool ContinueSimulation = true;
-        	int ANALYSIS_INTERVAL = ColorModel.timestepMax;
-        	/* flow adaptor keys to control */
-        	int SKIP_TIMESTEPS = 0;
-        	int MAX_STEADY_TIME = 1000000;
-		double ENDPOINT_THRESHOLD = 0.1;
-        	double FRACTIONAL_FLOW_INCREMENT = 0.05;
-        	if (ColorModel.db->keyExists( "FlowAdaptor" )){
-        		auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
-        		MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
-        		SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
-        		FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
-        	        ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
-        	}
-        	if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
-        		ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
-        	}
-        	/* Launch the simulation */
-        	FlowAdaptor Adapt(ColorModel);
-        	runAnalysis analysis(ColorModel);
-		while (ContinueSimulation){
-            	    /* this will run steady points */
-                    timestep += MAX_STEADY_TIME;
-		    MLUPS = ColorModel.Run(timestep);
-            	    if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
-            	    if (ColorModel.timestep > ColorModel.timestepMax){
-		      ContinueSimulation = false;
-		    }
-            	/* update the fractional flow by adding mass */
-            	int skip_time = 0;
-            	timestep = ColorModel.timestep;
-    			double SaturationChange = 0.0;
-    			double volB = ColorModel.Averages->gwb.V; 
-    			double volA = ColorModel.Averages->gnb.V; 
-    			double initialSaturation = volB/(volA + volB);
-
-			double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; 
-			double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; 
-			double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; 
-			double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; 
-			double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; 
-			double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; 			
-			double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
-			double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
-			if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
-		if (ContinueSimulation){
-            	while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
-            		timestep += ANALYSIS_INTERVAL;
-                	Adapt.UpdateFractionalFlow(ColorModel);
-                	MLUPS = ColorModel.Run(timestep);
-        			double volB = ColorModel.Averages->gwb.V; 
-        			double volA = ColorModel.Averages->gnb.V;
-        			SaturationChange = volB/(volA + volB) - initialSaturation;
-            		skip_time += ANALYSIS_INTERVAL;
+		Utilities::MPI comm( MPI_COMM_WORLD );
+		int rank   = comm.getRank();
+		int nprocs = comm.getSize();
+		std::string SimulationMode = "production";
+		// Load the input database
+		auto db = std::make_shared<Database>( argv[1] );
+		if (argc > 2) {
+			SimulationMode = "development";
 		}
-            	if (rank==0) printf("  *********************************************************************  \n");
-            	if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
-            	if (rank==0) printf("  *********************************************************************  \n");
+
+		if ( rank == 0 ) {
+			printf( "********************************************************\n" );
+			printf( "Running Color LBM	\n" );
+			printf( "********************************************************\n" );
+			if (SimulationMode == "development")
+				printf("**** DEVELOPMENT MODE ENABLED *************\n");
 		}
-            	/* apply timestep skipping algorithm to accelerate steady-state */
-            	/* skip_time = 0;
+		// Initialize compute device
+		int device = ScaLBL_SetDevice( rank );
+		NULL_USE( device );
+		ScaLBL_DeviceBarrier();
+		comm.barrier();
+
+		PROFILE_ENABLE( 1 );
+		// PROFILE_ENABLE_TRACE();
+		// PROFILE_ENABLE_MEMORY();
+		PROFILE_SYNCHRONIZE();
+		PROFILE_START( "Main" );
+		Utilities::setErrorHandlers();
+
+		auto filename = argv[1];
+		ScaLBL_ColorModel ColorModel( rank, nprocs, comm );
+		ColorModel.ReadParams( filename );
+		ColorModel.SetDomain();
+		ColorModel.ReadInput();
+		ColorModel.Create();     // creating the model will create data structure to match the pore
+		// structure and allocate variables
+		ColorModel.Initialize(); // initializing the model will set initial conditions for variables
+
+
+		if (SimulationMode == "development"){
+			double MLUPS=0.0;
+			int timestep = 0;
+			bool ContinueSimulation = true;
+			
+			/* Variables for simulation protocols */
+			auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "none" );
+			/* image sequence protocol */
+			int IMAGE_INDEX = 0;
+			int IMAGE_COUNT = 0;
+			std::vector<std::string> ImageList;
+			/* flow adaptor keys to control behavior */			
+			int SKIP_TIMESTEPS = 0;
+			int MAX_STEADY_TIME = 1000000;
+			double ENDPOINT_THRESHOLD = 0.1;
+			double FRACTIONAL_FLOW_INCREMENT = 0.05;
+			if (ColorModel.db->keyExists( "FlowAdaptor" )){
+				auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
+				MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
+				SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
+				FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
+				ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
+			}
+			/* analysis keys*/
+			int ANALYSIS_INTERVAL = ColorModel.timestepMax;
+			if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
+				ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
+			}
+			/* Launch the simulation */
+			FlowAdaptor Adapt(ColorModel);
+			runAnalysis analysis(ColorModel);
+			while (ContinueSimulation){
+				/* this will run steady points */
+				timestep += MAX_STEADY_TIME;
+				MLUPS = ColorModel.Run(timestep);
+				if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
+				if (ColorModel.timestep > ColorModel.timestepMax){
+					ContinueSimulation = false;
+				}
+				/* update the fluid configuration with the flow adapter */
+				int skip_time = 0;
+				timestep = ColorModel.timestep;
+				double SaturationChange = 0.0;
+				double volB = ColorModel.Averages->gwb.V; 
+				double volA = ColorModel.Averages->gnb.V; 
+				double initialSaturation = volB/(volA + volB);
+				double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; 
+				double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; 
+				double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; 
+				double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; 
+				double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; 
+				double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; 			
+				double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
+				double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
+				if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
+				if (ContinueSimulation){
+					while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
+						timestep += ANALYSIS_INTERVAL;
+						if (PROTOCOL == "fractional flow")	{							
+							Adapt.UpdateFractionalFlow(ColorModel);
+						}
+						else if (PROTOCOL == "image sequence"){
+							// Use image sequence
+							IMAGE_INDEX++;
+							if (IMAGE_INDEX < IMAGE_COUNT){
+								std::string next_image = ImageList[IMAGE_INDEX];
+								if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
+								ColorModel.color_db->putScalar<int>("image_index",IMAGE_INDEX);
+								Adapt.ImageInit(ColorModel, next_image);
+							}
+							else{
+								if (rank==0) printf("Finished simulating image sequence \n");
+								ColorModel.timestep =  ColorModel.timestepMax;
+							}
+						}
+						MLUPS = ColorModel.Run(timestep);
+						double volB = ColorModel.Averages->gwb.V; 
+						double volA = ColorModel.Averages->gnb.V;
+						SaturationChange = volB/(volA + volB) - initialSaturation;
+						skip_time += ANALYSIS_INTERVAL;
+					}
+					if (rank==0) printf("  *********************************************************************  \n");
+					if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
+					if (rank==0) printf("  *********************************************************************  \n");
+				}
+				/* apply timestep skipping algorithm to accelerate steady-state */
+				/* skip_time = 0;
             	timestep = ColorModel.timestep;
             	while (skip_time < SKIP_TIMESTEPS){
             		timestep += ANALYSIS_INTERVAL;
@@ -137,24 +161,24 @@ int main( int argc, char **argv )
             		Adapt.MoveInterface(ColorModel);
             		skip_time += ANALYSIS_INTERVAL;
             	}
-            	*/
-            }
-		  //ColorModel.WriteDebug();
-        }
+				 */
+			}
+			//ColorModel.WriteDebug();
+		}
 
-        else
-        	ColorModel.Run();        
+		else
+			ColorModel.Run();        
 
-        PROFILE_STOP( "Main" );
-        auto file  = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
-        auto level = db->getWithDefault<int>( "TimerLevel", 1 );
-	NULL_USE(level);
-        PROFILE_SAVE( file, level );
-        // ****************************************************
+		PROFILE_STOP( "Main" );
+		auto file  = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
+		auto level = db->getWithDefault<int>( "TimerLevel", 1 );
+		NULL_USE(level);
+		PROFILE_SAVE( file, level );
+		// ****************************************************
 
 
-    } // Limit scope so variables that contain communicators will free before MPI_Finialize
+	} // Limit scope so variables that contain communicators will free before MPI_Finialize
 
-    Utilities::shutdown();
-    return 0;
+	Utilities::shutdown();
+	return 0;
 }