From f7f70b58da535c964712b5abc9be63d19b07984f Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 14 Jun 2021 11:47:18 -0400 Subject: [PATCH 1/4] added sphinx rtd --- docs/Makefile | 20 ++++ docs/README.md | 12 +++ docs/source/conf.py | 70 ++++++++++++++ .../buildingModels/overview.rst | 20 ++++ .../designOverview/dataStructures.rst | 6 ++ docs/source/developerGuide/index.rst | 18 ++++ .../testingModels/unitTests.rst | 9 ++ docs/source/examples/running.rst | 9 ++ docs/source/index.rst | 25 +++++ docs/source/install.rst | 10 ++ docs/source/publications/publications.rst | 13 +++ docs/source/userGuide/IO/fileformat.rst | 13 +++ docs/source/userGuide/analysis/analysis.rst | 5 + docs/source/userGuide/index.rst | 17 ++++ .../PoissonBoltzmann/PoissonBoltzmann.rst | 6 ++ docs/source/userGuide/models/color/index.rst | 91 +++++++++++++++++++ .../models/color/protocols/centrifuge.rst | 12 +++ .../models/color/protocols/coreFlooding.rst | 63 +++++++++++++ .../models/color/protocols/fractionalFlow.rst | 17 ++++ .../models/color/protocols/imageSequence.rst | 27 ++++++ .../color/protocols/shellAggregation.rst | 18 ++++ .../models/freeEnergy/freeEnergy.rst | 6 ++ .../userGuide/models/greyscale/greyscale.rst | 7 ++ docs/source/userGuide/models/index.rst | 23 +++++ docs/source/userGuide/models/mrt/mrt.rst | 6 ++ .../models/nernstPlanck/nerstPlanck.rst | 6 ++ docs/source/userGuide/visualization/visit.rst | 5 + 27 files changed, 534 insertions(+) create mode 100644 docs/Makefile create mode 100644 docs/README.md create mode 100644 docs/source/conf.py create mode 100644 docs/source/developerGuide/buildingModels/overview.rst create mode 100644 docs/source/developerGuide/designOverview/dataStructures.rst create mode 100644 docs/source/developerGuide/index.rst create mode 100644 docs/source/developerGuide/testingModels/unitTests.rst create mode 100644 docs/source/examples/running.rst create mode 100644 docs/source/index.rst create mode 100644 docs/source/install.rst create mode 100644 docs/source/publications/publications.rst create mode 100644 docs/source/userGuide/IO/fileformat.rst create mode 100644 docs/source/userGuide/analysis/analysis.rst create mode 100644 docs/source/userGuide/index.rst create mode 100644 docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst create mode 100644 docs/source/userGuide/models/color/index.rst create mode 100644 docs/source/userGuide/models/color/protocols/centrifuge.rst create mode 100644 docs/source/userGuide/models/color/protocols/coreFlooding.rst create mode 100644 docs/source/userGuide/models/color/protocols/fractionalFlow.rst create mode 100644 docs/source/userGuide/models/color/protocols/imageSequence.rst create mode 100644 docs/source/userGuide/models/color/protocols/shellAggregation.rst create mode 100644 docs/source/userGuide/models/freeEnergy/freeEnergy.rst create mode 100644 docs/source/userGuide/models/greyscale/greyscale.rst create mode 100644 docs/source/userGuide/models/index.rst create mode 100644 docs/source/userGuide/models/mrt/mrt.rst create mode 100644 docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst create mode 100644 docs/source/userGuide/visualization/visit.rst diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 00000000..7fad3c1f --- /dev/null +++ b/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) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 00000000..f3ae3c75 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,12 @@ +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 diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 00000000..fa5df96e --- /dev/null +++ b/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") diff --git a/docs/source/developerGuide/buildingModels/overview.rst b/docs/source/developerGuide/buildingModels/overview.rst new file mode 100644 index 00000000..223b7c99 --- /dev/null +++ b/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 + + + diff --git a/docs/source/developerGuide/designOverview/dataStructures.rst b/docs/source/developerGuide/designOverview/dataStructures.rst new file mode 100644 index 00000000..bea4f91b --- /dev/null +++ b/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. diff --git a/docs/source/developerGuide/index.rst b/docs/source/developerGuide/index.rst new file mode 100644 index 00000000..d8448fa1 --- /dev/null +++ b/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/* + + diff --git a/docs/source/developerGuide/testingModels/unitTests.rst b/docs/source/developerGuide/testingModels/unitTests.rst new file mode 100644 index 00000000..4a4e84bf --- /dev/null +++ b/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.) diff --git a/docs/source/examples/running.rst b/docs/source/examples/running.rst new file mode 100644 index 00000000..83993f59 --- /dev/null +++ b/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. diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 00000000..fc5225d1 --- /dev/null +++ b/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` diff --git a/docs/source/install.rst b/docs/source/install.rst new file mode 100644 index 00000000..232da8d0 --- /dev/null +++ b/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 + diff --git a/docs/source/publications/publications.rst b/docs/source/publications/publications.rst new file mode 100644 index 00000000..aea8e5b6 --- /dev/null +++ b/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 + + diff --git a/docs/source/userGuide/IO/fileformat.rst b/docs/source/userGuide/IO/fileformat.rst new file mode 100644 index 00000000..faec04a6 --- /dev/null +++ b/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 -- + diff --git a/docs/source/userGuide/analysis/analysis.rst b/docs/source/userGuide/analysis/analysis.rst new file mode 100644 index 00000000..10d5282a --- /dev/null +++ b/docs/source/userGuide/analysis/analysis.rst @@ -0,0 +1,5 @@ +=========================== +Internal Analysis Framework +=========================== + +placeholder for analysis diff --git a/docs/source/userGuide/index.rst b/docs/source/userGuide/index.rst new file mode 100644 index 00000000..49025202 --- /dev/null +++ b/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/* diff --git a/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst b/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst new file mode 100644 index 00000000..93dbaff0 --- /dev/null +++ b/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. diff --git a/docs/source/userGuide/models/color/index.rst b/docs/source/userGuide/models/color/index.rst new file mode 100644 index 00000000..23bb5949 --- /dev/null +++ b/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 + diff --git a/docs/source/userGuide/models/color/protocols/centrifuge.rst b/docs/source/userGuide/models/color/protocols/centrifuge.rst new file mode 100644 index 00000000..21cf4eee --- /dev/null +++ b/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" + + diff --git a/docs/source/userGuide/models/color/protocols/coreFlooding.rst b/docs/source/userGuide/models/color/protocols/coreFlooding.rst new file mode 100644 index 00000000..12d1e799 --- /dev/null +++ b/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. + + + + diff --git a/docs/source/userGuide/models/color/protocols/fractionalFlow.rst b/docs/source/userGuide/models/color/protocols/fractionalFlow.rst new file mode 100644 index 00000000..7f5a26f0 --- /dev/null +++ b/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" + + + + diff --git a/docs/source/userGuide/models/color/protocols/imageSequence.rst b/docs/source/userGuide/models/color/protocols/imageSequence.rst new file mode 100644 index 00000000..25c91e6a --- /dev/null +++ b/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" + diff --git a/docs/source/userGuide/models/color/protocols/shellAggregation.rst b/docs/source/userGuide/models/color/protocols/shellAggregation.rst new file mode 100644 index 00000000..1ab54131 --- /dev/null +++ b/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" + + + + diff --git a/docs/source/userGuide/models/freeEnergy/freeEnergy.rst b/docs/source/userGuide/models/freeEnergy/freeEnergy.rst new file mode 100644 index 00000000..f53cfaad --- /dev/null +++ b/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. diff --git a/docs/source/userGuide/models/greyscale/greyscale.rst b/docs/source/userGuide/models/greyscale/greyscale.rst new file mode 100644 index 00000000..47c0621f --- /dev/null +++ b/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. diff --git a/docs/source/userGuide/models/index.rst b/docs/source/userGuide/models/index.rst new file mode 100644 index 00000000..6a158aa6 --- /dev/null +++ b/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/* + + diff --git a/docs/source/userGuide/models/mrt/mrt.rst b/docs/source/userGuide/models/mrt/mrt.rst new file mode 100644 index 00000000..1e232ea6 --- /dev/null +++ b/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. diff --git a/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst b/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst new file mode 100644 index 00000000..3842e2e6 --- /dev/null +++ b/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. diff --git a/docs/source/userGuide/visualization/visit.rst b/docs/source/userGuide/visualization/visit.rst new file mode 100644 index 00000000..f9ff2c8a --- /dev/null +++ b/docs/source/userGuide/visualization/visit.rst @@ -0,0 +1,5 @@ +====================================== +Visualizing simulation data with visit +====================================== + +placeholder for visit From 7f4fb3afe1dc6c4acf38029655c93ec3ce1116e2 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 14 Jun 2021 11:50:10 -0400 Subject: [PATCH 2/4] update docs README --- docs/README.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/README.md b/docs/README.md index f3ae3c75..bd31ba05 100644 --- a/docs/README.md +++ b/docs/README.md @@ -10,3 +10,10 @@ pip install sphinx-rtd-theme 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 +# \ No newline at end of file From 792310074de0e33bb419ef40fb4f073f67afbe8e Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Mon, 14 Jun 2021 13:08:28 -0400 Subject: [PATCH 3/4] remove debugging code from flow adapter --- models/ColorModel.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8e98f48a..98e6e002 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -2025,10 +2025,10 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ /* 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]; + //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; @@ -2145,7 +2145,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ 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; + //DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE; } else{ LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B; @@ -2156,7 +2156,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ 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; + //DebugMassB[n] = LOCAL_MASS_CHANGE; } } } @@ -2224,7 +2224,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ */ 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]; @@ -2247,7 +2247,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ AFILE = fopen(LocalRankFilename,"wb"); fwrite(regdata.data(),8,N,AFILE); fclose(AFILE); - + regdata.fill(0.f); for (int k=0; k Date: Mon, 14 Jun 2021 16:33:24 -0400 Subject: [PATCH 4/4] add image sequence to FlowAdaptor --- models/ColorModel.cpp | 371 +++++++++++++++++---------------- models/ColorModel.h | 4 +- tests/lbpm_color_simulator.cpp | 262 ++++++++++++----------- 3 files changed, 336 insertions(+), 301 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 98e6e002..3d726296 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -1981,188 +1981,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; iid[i]; // save what was read + for (int i=0; iid[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; kComm.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( "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; kSDs(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 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( "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; kSDs(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 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]; @@ -2221,10 +2230,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]; @@ -2241,13 +2250,13 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ } } } - + FILE *AFILE; sprintf(LocalRankFilename,"dA.%05i.raw",M.rank); AFILE = fopen(LocalRankFilename,"wb"); fwrite(regdata.data(),8,N,AFILE); fclose(AFILE); - + regdata.fill(0.f); for (int k=0; k 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); @@ -97,6 +98,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; diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 56b6d3dd..3591bbcd 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/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( 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( "max_steady_timesteps", 1000000 ); - SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 100000 ); - FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault( "fractional_flow_increment", 0.05); - ENDPOINT_THRESHOLD = flow_db->getWithDefault( "endpoint_threshold", 0.1); - } - if (ColorModel.analysis_db->keyExists( "analysis_interval" )){ - ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar( "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( 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( "protocol", "none" ); + /* image sequence protocol */ + int IMAGE_INDEX = 0; + int IMAGE_COUNT = 0; + std::vector 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( "max_steady_timesteps", 1000000 ); + SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 100000 ); + FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault( "fractional_flow_increment", 0.05); + ENDPOINT_THRESHOLD = flow_db->getWithDefault( "endpoint_threshold", 0.1); + } + /* analysis keys*/ + int ANALYSIS_INTERVAL = ColorModel.timestepMax; + if (ColorModel.analysis_db->keyExists( "analysis_interval" )){ + ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar( "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("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( "TimerFile", "lbpm_color_simulator" ); - auto level = db->getWithDefault( "TimerLevel", 1 ); - NULL_USE(level); - PROFILE_SAVE( file, level ); - // **************************************************** + PROFILE_STOP( "Main" ); + auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); + auto level = db->getWithDefault( "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; }