diff --git a/.travis.yml b/.travis.yml
index a981ea2544538cc96de76b958daf9fbeb25d7c82..b9ec906e73f428ac9ab8bb514f0d300356e07960 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,14 +1,12 @@
+dist: xenial
+
language: generic
python:
- 3.6
-install:
- - pip install --user sphinx # sphinx_rtd_theme
- # pypi's sphinx_rtd_theme is outdated. Internet recommends using current
- # master. Note: at some point, it will be OK to add `sphinx_rtd_theme` to
- # the `pip install` above => remove line below:
- - python -m pip install --user https://github.com/rtfd/sphinx_rtd_theme/archive/master.zip -U
+install:
+ - python -m pip install --user sphinx sphinx-rtd-theme
script: bash ./build_and_deploy.sh
@@ -24,7 +22,7 @@ addons:
deploy:
provider: pages
target-branch: gh-pages
- local-dir: docs/build
+ local-dir: docs/webroot
skip-cleanup: true
github-token: $GITHUB_TOKEN # Set in travis-ci.org dashboard, marked secure
keep-history: true
diff --git a/build_and_deploy.sh b/build_and_deploy.sh
index 4abeb2f926b549d7d3c3b1db78813c8a0d07ee09..a6a794b046b76bb22f8bf2db07dfc46ea9c2345b 100755
--- a/build_and_deploy.sh
+++ b/build_and_deploy.sh
@@ -1,4 +1,4 @@
-#!/bin/bash
+#!/bin/bash -ex
set -e # Exit with nonzero exit code if anything fails
# Then we build and deploy the sphinx / doxygen documentation
@@ -16,43 +16,5 @@ fi
#SSH_REPO=${REPO/https:\/\/github.com\//git@github.com:}
#SHA=`git rev-parse --verify HEAD`
-# Change working directory to /docs
-cd docs # ..................................... now we're in /docs
-
-
-############################# DOXYGEN ########################################
-
-
-cd doxygen # ..................................... now we're in /docs/doxygen
-
-# create temporary clone of MFiX-Exa SOURCE_BRANCH
-git clone https://gitlab-ci-token:$GITLAB_TOKEN@mfix.netl.doe.gov/gitlab/exa/mfix.git
-
-# build the Doxygen documentation
-doxygen doxygen.conf
-
-cd .. # ..................................... now we're in /docs
-
-# move it to the deploy directory: /docs/build
-mkdir build
-mv doxygen/html build/doxygen
-
-
-############################# SPHINX ########################################
-
-
-# now do sphinx
-make html
-cd build # ..................................... now we're in /docs/build
-mv html docs_html
-
-
-############################# WEB SERVER STUFF ################################
-
-
-# Sphinx is set up to treat build/html (and by mv, build/docs_html) as the
-# webroot. Hence the .nojekyll file is in the wrong place
-mv docs_html/.nojekyll .
-
-# PWD is currently /docs/build, we want to move /docs/webroot into /docs/build
-mv ../webroot/* .
+# build Doxygen docs and Sphinx docs
+make -C docs
diff --git a/docs/Makefile b/docs/Makefile
index 8422dc9b04f1f5207ee2ce7cb8277a2461aff0d1..a3c288c0ba74ba62a0e895895d177b4ea2ebf4fb 100644
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -1,14 +1,15 @@
-# Minimal makefile for Sphinx documentation
-#
-
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = python -msphinx
-SPHINXPROJ = amrex
SOURCEDIR = source
-BUILDDIR = build
+BUILDDIR = webroot
+
+all: html webroot/doxygen
+ mv webroot/html webroot/docs_html
+# Sphinx is set up to treat its build directory as the webroot, hence the
+# .nojekyll file is in the wrong place
+ mv webroot/docs_html/.nojekyll webroot
-# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
@@ -17,4 +18,9 @@ help:
# 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)
\ No newline at end of file
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+webroot/doxygen:
+ git clone "https://gitlab-ci-token:$(GITLAB_TOKEN)@mfix.netl.doe.gov/gitlab/exa/mfix.git" doxygen/mfix
+ cd doxygen && doxygen doxygen.conf
+ mv doxygen/html $@
diff --git a/docs/source/BuildingMacVelocities.rst b/docs/source/BuildingMacVelocities.rst
new file mode 100644
index 0000000000000000000000000000000000000000..5ed22ad3b065499105f91003e8bd4fe905d71c63
--- /dev/null
+++ b/docs/source/BuildingMacVelocities.rst
@@ -0,0 +1,62 @@
+Creating the MAC velocities
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To create the normal velocities on faces, we first extrapolate from the cell centers on each side using the
+slopes as computed earlier, and upwind the face value to define :math:`U^{pred}`
+
+. To compute the x-velocity on the x-faces of regular (ie not cut) cells, we call
+
+ .. code:: shell
+
+ AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,
+ {
+ // X-faces
+ Real upls = ccvel_fab(i ,j,k,0) - 0.5 * xslopes_fab(i ,j,k,0);
+ Real umns = ccvel_fab(i-1,j,k,0) + 0.5 * xslopes_fab(i-1,j,k,0);
+ if ( umns < 0.0 && upls > 0.0 ) {
+ umac_fab(i,j,k) = 0.0;
+ } else {
+ Real avg = 0.5 * ( upls + umns );
+ if ( std::abs(avg) < small_vel) { umac_fab(i,j,k) = 0.0;
+ } else if (avg >= 0) { umac_fab(i,j,k) = umns;
+ } else { umac_fab(i,j,k) = upls;
+ }
+ }
+ });
+
+For cut cells we test on whether the area fraction is non-zero:
+
+ .. code:: shell
+
+ AMREX_CUDA_HOST_DEVICE_FOR_3D(ubx, i, j, k,
+ {
+ // X-faces
+ if (ax_fab(i,j,k) > 0.0)
+ {
+ Real upls = ccvel_fab(i ,j,k,0) - 0.5 * xslopes_fab(i ,j,k,0);
+ Real umns = ccvel_fab(i-1,j,k,0) + 0.5 * xslopes_fab(i-1,j,k,0);
+ if ( umns < 0.0 && upls > 0.0 ) {
+ umac_fab(i,j,k) = 0.0;
+ } else {
+ Real avg = 0.5 * ( upls + umns );
+ if ( std::abs(avg) < small_vel) { umac_fab(i,j,k) = 0.0;
+ } else if (avg >= 0) { umac_fab(i,j,k) = umns;
+ } else { umac_fab(i,j,k) = upls;
+ }
+ }
+ } else {
+ umac_fab(i,j,k) = huge_vel;
+ }
+ });
+
+We then perform a MAC projection on the face-centered velocities to enforce that they satisfy
+
+.. math:: \nabla \cdot (\varepsilon_g U^{MAC}) = 0
+
+We do this by solving
+
+.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi^{MAC} = \nabla \cdot \left( \varepsilon_g U^{pred} \right)
+
+then defining
+
+.. math:: U^{MAC} = U^{pred} - \frac{1}{\rho_g} \nabla \phi^{MAC}
diff --git a/docs/source/CITests.rst b/docs/source/CITests.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6a1eb505d0ff39b6f03991ecc2865c6b0e61b1b4
--- /dev/null
+++ b/docs/source/CITests.rst
@@ -0,0 +1,107 @@
+.. _Chap:CITesting :
+
+Continuous Integration
+======================
+
+The following regression tests are run every time a commit is pushed to the main
+MFiX-Exa repository on the NETL gitlab.
+
+For each of the tests in the chart below, there are
+three directional variations; these are identified in the repository as,
+for example, FLD01-x, FLD01-y, and FLD01-z.
+
+For each direction, where appropriate, there are multiple versions, with the following notations:
+
+ * SGS: single grid serial
+
+ * MGS: multiple grid serial
+
+ * TGS: tiled grid serial
+
+ * MGP: multiple grid parallel
+
+Below Ng = number of grids, Npa = number of particles, Np = number of MPI ranks.
+
+All the FLD cases are fluid-only and steady state.
+All the DEM cases are particle-only except for DEM06 and DEM07 which are fluid and particles;
+these both use the "BVK2" drag type.
+In all cases the particle data were read in from "particle_input.dat"
+
+None of these tests have non-rectangular geometry.
+
+"NSW" means "No Slip Wall" and "Per" is "periodic."
+"MI/PO" refers to Mass Inflow at the low end of the domain and Pressure Outflow at the high end.
+"PI/PO" refers to Pressure Inflow at the low end of the domain and Pressure Outflow at the high end.
+
+Additional detail about these problems is given in tests/README.md
+
+Single-grid, single-process (SGS) particle-only tests:
+
++-------+----+----+----+------+------+-------+-----+--------------------+
+| Test | nx | ny | nz | bc_x | bc_y | bc_z | Npa | Description |
++=======+====+====+====+======+======+=======+=====+====================+
+| DEM01 | 2 | 5 | 5 | NSW | Per | Per | 1 | Freely falling |
+| | | | | | | | | particle with |
+| | | | | | | | | wall collision |
++-------+----+----+----+------+------+-------+-----+--------------------+
+| DEM02 | 2 | 5 | 5 | NSW | Per | Per | 1 | Multiple bounces |
+| | | | | | | | | with bounce height |
+| | | | | | | | | measured |
++-------+----+----+----+------+------+-------+-----+--------------------+
+| DEM03 | 2 | 5 | 5 | NSW | Per | Per | 2 | Two stacked |
+| | | | | | | | | compressed |
+| | | | | | | | | particles |
++-------+----+----+----+------+------+-------+-----+--------------------+
+| DEM04 | 4 | 4 | 4 | NSW | Per | Per | 1 | Single particle |
+| | | | | | | | | slipping on a |
+| | | | | | | | | rough surface |
++-------+----+----+----+------+------+-------+-----+--------------------+
+| DEM05 | 5 | 2 | 5 | Per | Per | Per | 93 | Oblique particle |
+| | | | | | | | | collisions |
+| | | | | | | | | |
++-------+----+----+----+------+------+-------+-----+--------------------+
+
+
+Steady-state fluid-only tests:
+
++-------+-----+----+----+----+-------+------+------+-----+----------------------+
+| Test | | nx | ny | nz | bc_x | bc_y | bc_z | Ng | Np |
++=======+=====+====+====+====+=======+======+======+=====+======================+
+| FLD01 | | 8 | 8 | 4 | Per | NSW | Per | Poiseuille flow |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+| | SGS | | | | | | | 1 | 1 | |
+| | MGS | | | | | | | 4 | 1 | |
+| | MGP | | | | | | | 4 | 8 | |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+| FLD02 | | 80 |16 | 16 | MI/PO | NSW | NSW | Couette flow |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+| | SGS | | | | | | | 1 | 1 | |
+| | MGS | | | | | | | 40 | 1 | |
+| | MGP | | | | | | | 40 | 8 | |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+| FLD03 | | 8 | 8 | 4 | PI/PO | NSW | Per | Poiseuille flow |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+| | SGS | | | | | | | 1 | 1 | |
+| | MGS | | | | | | | 4 | 1 | |
+| | MGP | | | | | | | 4 | 8 | |
++-------+-----+----+----+----+-------+------+------+-----+----+-----------------+
+
+Coupled particle/fluid tests:
+
++-------+-----+----+----+----+------+------+------+------+----+--------------------+
+| Test | | nx | ny | nz | bc_x | bc_y | bc_z | Npa | Ng | Np |
++=======+=====+====+====+====+======+======+======+======+====+====================+
+| DEM06 | | 50 | 5 | 5 | NSW | NSW | NSW | 1 | Single particle falling |
+| | | | | | | | | | under gravity |
++-------+-----+----+----+----+------+------+------+------+----+----+---------------+
+| | SGS | | | | | | | | 1 | 1 | |
+| | MGS | | | | | | | | 10 | 1 | |
+| | MGP | | | | | | | | 10 | 8 | |
++-------+-----+----+----+----+------+------+------+------+----+----+---------------+
+| DEM07 | | 20 | 20 | 20 | Per | Per | Per | 1222 | Homogeneous cooling |
+| | | | | | | | | | system |
++-------+-----+----+----+----+------+------+------+------+----+----+---------------+
+| | SGS | | | | | | | | 1 | 1 | |
+| | MGS | | | | | | | | 8 | 1 | |
+| | MGP | | | | | | | | 8 | 8 | |
++-------+-----+----+----+----+------+------+------+------+----+----+---------------+
diff --git a/docs/source/Debugging.rst b/docs/source/Debugging.rst
new file mode 100644
index 0000000000000000000000000000000000000000..ff7742fb08a42c29454ded79abba731f95b57d32
--- /dev/null
+++ b/docs/source/Debugging.rst
@@ -0,0 +1,65 @@
+.. role:: cpp(code)
+ :language: c++
+
+.. _Chap:Debugging:
+
+Debugging
+=========
+
+Debugging is an art. Everyone has their own favorite method. Here we
+offer a few tips we have found to be useful.
+
+Compiling in debug mode (e.g., :cpp:`make DEBUG=TRUE` for gmake users;
+:cpp:`cmake -DDEBUG` for cmake users) and running with
+:cpp:`amrex.fpe_trap_invalid=1` in the inputs file can be helpful.
+In debug mode, many compiler debugging flags are turned on and all
+:cpp:`MultiFab` s are initialized to signaling NaNs. The
+:cpp:`amrex.fpe_trap_invalid` parameter will result in backtrace files
+when a floating point exception occurs. One can then examine those
+files to track down the origin of the issue.
+
+Several other ways to look at the data include:
+
+1) Writing a :cpp:`MultiFab` to disk with
+
+.. highlight:: c++
+
+::
+
+ VisMF::Write(const FabArray<FArrayBox>& mf, const std::string& name);
+
+and examining it with ``Amrvis`` (section :ref:`sec:amrvis` in the AMReX documentation).
+
+2) You can also use the :cpp:`print_state` routine:
+
+.. highlight:: c++
+
+::
+
+ void print_state(const MultiFab& mf, const IntVect& cell, const int n=-1);
+
+which outputs the data for a single cell.
+
+3) If you want to compare old and new plotfiles,
+
+.. highlight:: c++
+
+::
+
+ fcompare --infile1 plt00000_run1 --infile2 plt00000_run2 --diffvar u_g
+
+will print out the maximum absolute and relative differences between the two plotfiles
+for each variable and will also create a new plotfile "diffs" that contains the difference
+in u_g (in this case) between the two plotfiles.
+
+The :cpp:`fcompare` executable can be built in AMReX (go to amrex/Tools/Plotfile and type "make").
+
+4) Valgrind is another useful debugging tool. Note that for runs using
+more than one MPI process, one can tell valgrind to output to different
+files for different processes. For example,
+
+.. highlight:: console
+
+::
+
+ mpiexec -n 4 valgrind --leak-check=yes --track-origins=yes --log-file=vallog.%p ./mfix.exe ...
diff --git a/docs/source/DualGrid.rst b/docs/source/DualGrid.rst
new file mode 100644
index 0000000000000000000000000000000000000000..de11d73448aba6d45dff19bed4657d9145f58cc8
--- /dev/null
+++ b/docs/source/DualGrid.rst
@@ -0,0 +1,25 @@
+.. role:: cpp(code)
+ :language: c++
+
+.. role:: fortran(code)
+ :language: fortran
+
+.. _ss:dual_grid:
+
+Dual Grid Approach
+------------------
+
+In MFiX-Exa the mesh work and particle work have very different requirements for load balancing.
+
+Rather than using a combined work estimate to create the same grids for mesh and particle
+data, we have the option to pursue a "dual grid" approach.
+
+With this approach the mesh (:cpp:`MultiFab`) and particle (:cpp:`ParticleContainer`) data
+are allocated on different :cpp:`BoxArrays` with different :cpp:`DistributionMappings`.
+
+This enables separate load balancing strategies to be used for the mesh and particle work.
+
+The cost of this strategy, of course, is the need to copy mesh data onto temporary
+:cpp:`MultiFabs` defined on the particle :cpp:`BoxArrays` when mesh-particle communication
+is required.
+
diff --git a/docs/source/EB.rst b/docs/source/EB.rst
new file mode 100644
index 0000000000000000000000000000000000000000..4786a36a69f11de6c0bce4edbb73ed14f4b50ed7
--- /dev/null
+++ b/docs/source/EB.rst
@@ -0,0 +1,16 @@
+.. _Chap:EB:
+
+Embedded Boundaries
+===================
+
+MFiX-Exa uses AMReX's Embedded Boundaries to represent container walls. This
+allows MFiX-Exa to simulate a wide range of shapes and boundary conditions using
+constructive solid geometry. MFiX-Exa also has the capability of locally
+refining the computational meshes near walls.
+
+
+.. toctree::
+ :maxdepth: 1
+ :caption: Contents:
+
+ EBWalls
diff --git a/docs/source/EBWalls.rst b/docs/source/EBWalls.rst
new file mode 100644
index 0000000000000000000000000000000000000000..2f5373e4c01c4c7eb4ab104177574865e16b8f8e
--- /dev/null
+++ b/docs/source/EBWalls.rst
@@ -0,0 +1,308 @@
+.. role:: cpp(code)
+ :language: c++
+
+.. role:: fortran(code)
+ :language: fortran
+
+
+.. _sec:EB-basics:
+
+
+Constructing Embedded Boundaries in MFiX-Exa
+============================================
+
+MFiX uses AMReX's constructive solid geometry framework defined in the namespace
+:cpp:`amrex::EB2`. See the `AMReX EB documentation`_ for more details. These are
+defined in ``src/eb/mfix_eb.cpp``. A the function :cpp:`mfix::make_eb_geometry`
+(also defined in ``src/eb/mfix_eb.cpp``) selects :cpp:one of the following
+geometries depending on the value of the :cpp:``mfix.geometry`` setting in the
+``inputs`` file.
+
++------------------------------+----------------------+-------------------------+
+| Description | ``mfix.geometry`` | Implementation Satus |
++==============================+======================+=========================+
+| Planar walls from mfix.dat | Don't specify | Fully implemented |
+| (on by default) | | |
++------------------------------+----------------------+-------------------------+
+| Box (up to six walls) | ``box`` | Fully implemented |
++------------------------------+----------------------+-------------------------+
+| Cylinder | ``cylinder`` | Fully implemented |
++------------------------------+----------------------+-------------------------+
+| Hopper | ``hopper`` | Fully implemented |
++------------------------------+----------------------+-------------------------+
+| Cyclone | ``cyclone`` | Fully implemented |
++------------------------------+----------------------+-------------------------+
+| General | ``general`` | Partially implemented |
+| | cf. note 1 below | |
++------------------------------+----------------------+-------------------------+
+| Hourglass | ``hourglass`` | Deprecated |
+| | cf. note 1 | cf. note 2 below |
++------------------------------+----------------------+-------------------------+
+| CLR (chemical looping | ``clr`` | Deprecated |
+| reactor) | cf. note 1 | cf. note 2 |
++------------------------------+----------------------+-------------------------+
+| CLR Riser | ``clr_riser`` | Deprecated |
+| | cf. note 1 | cf. note 2 |
++------------------------------+----------------------+-------------------------+
+
+1. Older (legacy) alternative settings are:
+
+ +-----------------------------+-------------------------------+
+ | Value of ``mfix.geometry`` | Alternative |
+ +=============================+===============================+
+ | ``general`` | ``mfix.use_poly2 = true`` |
+ | | or ``mfix.use_walls = true`` |
+ +-----------------------------+-------------------------------+
+ | ``hourglass`` | ``mfix.hourglass = true`` |
+ +-----------------------------+-------------------------------+
+ | ``clr`` | ``mfix.clr = true`` |
+ +-----------------------------+-------------------------------+
+ | ``clr_riser`` | ``mfix.clr_riser = true`` |
+ +-----------------------------+-------------------------------+
+
+2. These geometries where not ported from AMReX's old :cpp:``EB`` system to the
+ new :cpp:``EB2``.
+
+Also note that planar boundary conditions can be specified in the ``mfix.dat``
+file. Even if the user does not specify an ``mfix.geometry`` in the ``inputs``,
+any no-slip or free-slip boundary conditions are expressed as EB walls.
+
+There are also two parameters (specified in the ``inputs`` file) that influence
+the level-set creation:
+
++-------------------------------+------------------------------------------------+
+| Parameter | Description |
++===============================+================================================+
+| ``amr.max_level`` | If greater than 0, MFiX operates in multi-level|
+| | mode. The level-set grids follow all other |
+| | grids. If equal to 1, the level-set has two |
+| | levels (one additional level with higher |
+| | refinement). |
++-------------------------------+------------------------------------------------+
+| ``mfix.levelset__refinement`` | If ``amr.max_level > 1`` this parameter is |
+| | ignored. Otherwise it sets the maximum |
+| | refinement of the level-set |
++-------------------------------+------------------------------------------------+
+
+
+How MFiX-Exa Constructs the EB Geometry
+---------------------------------------
+
+Once a geometry is selected by :cpp:`mfix::make_eb_geometry`, the procedure is
+the same for (almost) all geometries. Also see the `AMReX geometry
+documentation`_ for information on how to construct new geometries:
+
+1. Construct an implicit function representing the geometry (using the language
+ of constructive solid geometry). For example
+
+.. highlight:: c++
+
+::
+
+ EB2::CylinderIF my_cyl(radius, height, direction, center, inside);
+ auto gshop_cyl = EB2::makeShop(my_cyl);
+
+2. (Optional) Construct the implicit function representing the EB seen by the
+ particles. This might deviate from the "standard" EB depending on the
+ specific application. In all standard cases used by mfix, this step is
+ omitted.
+
+3. Call :cpp:`mfix::build_eb_levels(gshop)` this function builds the EB levels
+ and fills the implicit function :cpp:`MultiFab` (the later being used to
+ construct the level-set function). Note that this also makes the particle EB
+ levels point to the fluid eb levels.
+
+4. (Optional, do this if you did 2.) Call
+ :cpp:`mfix::build_particle_eb_levels(gshop_part)` **last**. This will update
+ the particle EB levels.
+
+
+MFiX's EB Data Structures
+-------------------------
+
+The :cpp:`mfix` class stores the following EB data:
+
+.. highlight:: c++
+
+::
+
+ //! EB levels representing fluid boundary conditions
+ Vector<const EB2::Level *> eb_levels;
+ //! EB levels representing particle boundary conditions (same as
+ //! `mfix::eb_levels` but might include additional walls).
+ Vector<const EB2::Level *> particle_eb_levels;
+
+ //! EB factory that lives on the fluid grids
+ Vector< std::unique_ptr<amrex::EBFArrayBoxFactory> > ebfactory;
+ //! EB factory that lives on the particle grids
+ Vector< std::unique_ptr<amrex::EBFArrayBoxFactory> > particle_ebfactory;
+
+As discussed in the previous sub-section, the difference between
+:cpp:`mfix::eb_levels` and :cpp:`mfix::particle_eb_levels` enables the user to
+specify a modfied EB geometry for particles only. Whereas the fluid sees the EB
+geometry in :cpp:`mfix::eb_levels`. If no addition particle EB geometry is
+specified (point 4 in the previous section), then
+:cpp:`mfix::particle_eb_levels` points to :cpp:`mfix::eb_levels`.
+
+In the same spirit, the :cpp:`mfix::ebfactory` is constructed over the fluid
+grid and using the fluid EB levels, whereas :cpp:`mfix::particle_ebfactory` is
+constructed over the particle grid using the particle EB levels.
+
+
+A note about constructing EB Levels
+-----------------------------------
+
+MFiX-Exa builds EB levels in :cpp:`mfix::build_eb_levels` (via
+:cpp:`LSCore<F>::BuildEBLevel`)
+
+.. highlight:: c++
+
+::
+
+ EB2::Build(gshop, geom[lev], required_crse_lev, max_crse_level);
+ const EB2::IndexSpace & ebis = EB2::IndexSpace::top();
+
+
+When building an EB level, the maximum coarsening level (:cpp:`int
+max_crse_level`) and the required coarsening level (:cpp:`int
+required_crse_lev`) need to be specified. The reason for this is that we need to
+specify to which level of coarseness the EB is still defined. It might not be
+immediately obvious, but the Poisson solver (used in the fluid solve) also
+depends indirectly on these parameters. Thus changing these during EB level
+creation might restrict how many levels the MLMG solver can use, and therefore
+give slightly different answers in the fluid solve.
+
+
+
+Local Mesh Refinement at Walls
+==============================
+
+MFiX-Exa has the capability of locally refining the computational grid near EBs.
+This is done by tagging (in :cpp:`mfix::ErrorEst`) any cells with volume
+fraction between 0 and 1. To enable local mesh refinement, set ``amr.max_level``
+to a value greater than 1. Note that the parameter ``mfix.levelset__refinement``
+is ignored on all cases except when ``amr.max_level = 1``.
+
+
+MFiX-Exa Initialization Process
+-------------------------------
+
+Since MFiX requires the volume fraction when building grids (because this is
+needed by :cpp:`mfix::ErrorEst`), the EB geometries need to be built before
+calling :cpp:`mfix::Init`. The recommended procedure therefore is
+
+.. highlight:: c++
+
+::
+
+ // Default constructor (geom[lev] is defined here)
+ mfix my_mfix;
+
+ // Initialize internals from ParamParse database
+ my_mfix.InitParams(solve_fluid, solve_dem, call_udf);
+
+ // Initialize memory for data-array internals
+ my_mfix.ResizeArrays();
+
+ // Construct EB (must be done _before_ mfix::Init)
+ my_mfix.make_eb_geometry();
+
+ // Initialize derived internals. Grids are create here.
+ my_mfix.Init(dt, time);
+
+ // Create EB factories on new grids
+ my_mfix.make_eb_factories();
+
+ if (solve_dem)
+ {
+ // Fill level-sets on each level (must be done _after_ mfix::Init)
+ my_mfix.fill_eb_levelsets();
+ }
+
+ // Finish constructing levels
+ my_mfix.InitLevelData(dt,time);
+
+ // Regrid (ensure all MultiFabs are on their correct grids)
+ my_mfix.Regrid();
+
+
+Also note that mfix defines boundary conditions in Fortran also (via the
+mfix.dat). Since these are potentially needed to build EB walls,
+:cpp:`mfix::make_eb_geometry` also calls :cpp:`mfix_set_bc_type`.
+
+The grids for each level are build in the :cpp:`mfix::Init` by invoking the
+initialization functions inherited from :cpp:`amrex::AmrCore`.
+
+.. highlight:: c++
+
+::
+
+ // This tells the AmrMesh class not to iterate when creating the initial
+ // grid hierarchy
+ SetIterateToFalse();
+
+ // This tells the Cluster routine to use the new chopping routine which
+ // rejects cuts if they don't improve the efficiency
+ SetUseNewChop();
+
+ // This Builds the new Grids
+ InitFromScratch(0.);
+
+
+
+The Level-Set Function
+======================
+
+MFiX-Exa uses a level-set function to resolve particle-wall collisions. See the
+`AMReX Level-Set documentation`_ for more details. The level-set function is
+stored on the nodal :cpp:`Vector<std::unique_ptr<MultiFab>> mfix::level_sets`.
+The level-set data is always stored on the particle grids. Depending on the
+input ``amr.max_level`` The level-set can be in one of two modes:
+
+1. MFiX-Exa is running in single-level mode (:cpp:`nlev == 1`). Then
+ :cpp:`mfix::level_sets[0]` will be at the same resolution as the fluid
+ (except that it is stored on the particle grid). Even though :cpp:`nlev == 1`,
+ there is a second level, :cpp:`level_sets[1]`. This level is the same as
+ :cpp:`level_sets[0]` but refined by :cpp:`mfix::levelset__refinement`. This
+ way the level-set always has the appropriate resolution to resolve structures
+ in the EB, even if the fluid is defined on a fairly coarse grid.
+
+2. MFiX-Exa is running in multi-level mode (:cpp:`nlev > 1`). The the parameter
+ :cpp:`mfix::levelset__refinement` is ignored. :cpp:`mfix::level_sets` then
+ follows the rest of MFiX, i.e. it is defined on the particle grids on all
+ levels.
+
+The level-set is used in two places:
+
+1. The function :cpp:`MFIXParticleContainer::EvolveParticles` interpolates the
+ level-set onto each particle's position in order to resolve collisions with
+ the EBs. If :cpp:`nlev == 1`, :cpp:`level_sets[1]` is used to evolve the
+ particle positions. Otherwise :cpp:`level_sets[lev]` is used for each level.
+
+2. The fluid-particle coupling can sometimes rely on neighbor stencils where one
+ or more cell is covered by an EB. In order to avoid values that do not
+ conform with the boundary conditions, the fluid velocity is reconstructed in
+ those cells. The algorithm relies on the level-set, and uses
+ :cpp:`level_sets[lev]` on each level.
+
+
+Special Cases Involving Level-Sets
+----------------------------------
+
+The level-set function is filled by the `mfix::fill_eb_levelsets()` function.
+There are two special cases involving level-sets:
+
+1. Mass-Inflow boundary conditions are not given EB walls. However, we don't
+ want particles to fall out of a MI either, so at the very end of the
+ `mfix::fill_eb_levelsets()` function we call `mfix::intersect_ls_walls()`.
+ This performs an intersection operation with the level-set representing a
+ wall at each MI.
+
+2. Box geometries and regular geometries are comprised entirely out of planar
+ surfaces. Therefore the levelset is not construction out of an EB factory (as
+ would be the case for all other geometries). But out of an intersection with
+ all planar surfaces. This has the advantage of correctly describing corners.
+
+.. _AMReX EB documentation: https://amrex-codes.github.io/amrex/docs_html/EB_Chapter.html
+.. _AMReX Level-Set documentation: https://amrex-codes.github.io/amrex/docs_html/EB.html#level-sets
+.. _AMReX geometry documentation: https://amrex-codes.github.io/amrex/docs_html/EB.html#initializing-the-geometric-database
diff --git a/docs/source/FluidEquations.rst b/docs/source/FluidEquations.rst
index 125b3ffacd2ae087c8b9746991cacdebeea18d85..3b814b4b605d8d61d39063444ed9017ac0ef1795 100644
--- a/docs/source/FluidEquations.rst
+++ b/docs/source/FluidEquations.rst
@@ -22,13 +22,62 @@ Conservation of fluid mass:
.. math:: \frac{\partial (\varepsilon_g \rho_g)}{\partial t} + \nabla \cdot (\varepsilon_g \rho_g U_g) = 0
-Conservation of fluid momentum:
+Unlike two-fluid modeling, :math:`\varepsilon_g = 1 - \varepsilon_s`, is not a solution
+variable. Rather, :math:`\varepsilon_s` is evaluated from the particle field through
+volume filtering,
-.. math:: \frac{ \partial (\varepsilon_g \rho_g U)}{\partial t} + \nabla \cdot (\varepsilon_g \rho_g U_g U_g) + \varepsilon_g \nabla p_g = \nabla \cdot \tau
- + \sum_p \beta_p (V_p - U_g) + \rho_g g
+.. math::
+ (1 - \varepsilon_g) A (\boldsymbol{x},t) \approx
+ \sum_{i=1}^{N_p} A_i(\boldsymbol{X}_i,t) {\mathcal G}
+ (\left| \boldsymbol{x} - \boldsymbol{X}_i \right|) {\mathcal V}_i
-where :math:`\sum_p \beta_p (V_p - U_g)` is the drag term in which :math:`V_p` represents the particle velocity and :math:`\beta_p` is the drag coefficient associated with that particle
+where :math:`A_i` is a general particle property, :math:`{\mathcal V}_i` is the particle
+volume and :math:`\mathcal G` is a transfer kernel with compact support--here linear hat.
+Setting :math:`A_i = 1` gives the local void fraction.
-Conservation of fluid volume:
+
+Assuming the fluid phase is incompressible, :math:`\frac{D \rho_g}{Dt} = 0`, the conservation of fluid mass is equivalent to the conservation of fluid volume:
.. math:: \frac{\partial \varepsilon_g}{\partial t} + \nabla \cdot (\varepsilon_g U_g) = 0
+
+The conservation of fluid momentum is:
+
+.. math:: \frac{ \partial (\varepsilon_g \rho_g U_g)}{\partial t}
+ + \nabla \cdot (\varepsilon_g \rho_g U_g U_g) + \varepsilon_g \nabla p_g
+ = \nabla \cdot \tau + M_{sg} + \varepsilon_g \rho_g g
+
+where :math:`M_{sg} = - M_{gs}` is the generalized interfacial momentum transfer from
+the solid particles to the fluid-phase. Like :math:`\varepsilon_s`,
+:math:`M_{gs}` is determined from the L-E transfer kernel by setting :math:`A_i = F_{gi}`,
+where :math:`F_{gi}` is the force due to the fluid-phase on the ith particle. Following
+MFiX classic (and many other CFD-DEM codes designed for high density ratio gas-solids
+flows), only buoyancy (pressure gradient) and steady drag are considered:
+
+.. math::
+ F_{gi} = - \mathcal{V}_i \nabla p_g
+ - \frac{1}{2} C_D \rho_g \boldsymbol{V}_{ig} \left|\boldsymbol{V}_{ig}\right| A_i^{(proj)}
+
+where :math:`\boldsymbol{V}_{ig} = \boldsymbol{V}_i - \boldsymbol{U}_g ( \boldsymbol{X}_i )`
+is the velocity of ith particle relative to the fluid-phase (at the particle position
+:math:`\boldsymbol{X}_i`). :math:`F_{gi}` is closed by the specification of a drag
+coefficient, :math:`C_D`. Currently, MFiX-Exa includes Wen-Yu, Gidaspow and BVK2 drag laws.
+
+.. wdf todo
+ gidaspow form: C_D^{Gidaspow} = \chi C_D^{(Wen-Yu)} + \left(1 - \chi \right) C_D^{(Ergun)}
+ \chi = \frac{\arctan 150 \left( \varepsilon_g - 0.8 \right)}{\pi} + \frac{1}{2}
+ wen-yu: C_D^{(Wen-Yu)} = \max \left[\frac{24}{Re_i}\left(1 + 0.15 Re_i^{0.687}\right),\ 0.44 \right] \left( 1 - \varepsilon_g \right)^{-1.65}
+ ergun: C_D^{(Ergun)} = \frac{200 \left(1 - \varepsilon_g\right)}{Re_i} + \frac{7}{3}
+ BVK2: C_D^{BVK2} = ...
+ particle reynolds number: Re_i = \frac{\rho_g \left( 1 - \varepsilon_g \right) d_i \left| \boldsymbol{V}_{ig} \right| }{\mu_g}
+
+
+In chemical engineering literature, it is common to lump all drag-related terms of
+:math:`M_{gs}` into :math:`\beta`. With this simplification and some re-arrangement,
+the fluid momentum takes the more convenient form:
+
+.. math:: \frac{ \partial (\varepsilon_g \rho_g U_g)}{\partial t}
+ + \nabla \cdot (\varepsilon_g \rho_g U_g U_g) + \varepsilon_g \nabla p_g
+ = \nabla \cdot \tau + \sum_p \beta_p (V_p - U_g) + \varepsilon_g \rho_g g
+
+
+
diff --git a/docs/source/FluidTimeDiscretization.rst b/docs/source/FluidTimeDiscretization.rst
index 74f31a7f08f6a989c191eeadf93ab3615978c147..2654ad8a6d92b116147367621306ec87a04a869b 100644
--- a/docs/source/FluidTimeDiscretization.rst
+++ b/docs/source/FluidTimeDiscretization.rst
@@ -2,51 +2,10 @@
Time Discretization
===================
-In the absence of reactions, we assume that the fluid density is unchanged.
-We compute the fluid volume fraction directly from the particle locations.
+.. toctree::
+ :maxdepth: 1
-Thus here we focus on the discretization of the momentum equation
-
-In the predictor
-
-- Define :math:`U^{MAC,n}`, the face-centered (staggered) MAC velocity which is used for advection, using :math:`U^n`
-
-- Define an approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\ast}` by setting
-
-.. math:: (\varepsilon_g \rho_g U)^{\ast} &= (\varepsilon_g \rho_g U)^n -
- \Delta t \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g) + \varepsilon_g \nabla {p_g}^{n-1/2} \right) \\ &+
- \Delta t \left( \nabla \cdot \tau^n + \sum_p \beta_p (V_p - {U_g}^{\ast}) + \rho_g \varepsilon_g g \right)
-
-- Project :math:`U^{\ast}` by solving
-
-.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g U)^{\ast}+ {\varepsilon_g}{\rho_g} \nabla {p_g}^{n-1/2} \right)
-
-then defining
-
-.. math:: U^{\ast \ast} = U^{\ast} - \frac{1}{\rho_g} \nabla \phi
-
-and
-
-.. math:: {p_g}^{n+1/2, \ast} = \phi
-
-
-In the corrector
-
-- Define :math:`U^{MAC,\ast \ast}` at the "new" time using :math:`U^{\ast \ast}`
-
-- Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U)^{\ast \ast \ast}` by setting
-
-.. math:: (\varepsilon_g \rho_g U)^{\ast \ast \ast} &= (\varepsilon_g \rho_g U)^n - \frac{\Delta t}{2} \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^n + \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^{\ast \ast}\right) + \\ &+ \frac{\Delta t}{2} \left( \nabla \cdot \tau^n + \nabla \cdot \tau^{\ast \ast} \right) + \Delta t \left( - \varepsilon_g \nabla {p_g}^{n+1/2,\ast} + \sum_p \beta_p (V_p - {U_g}^{\ast \ast \ast}) + \varepsilon_g \rho_g g \right)
-
-- Project :math:`U^{\ast \ast \ast}` by solving
-
-.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g U)^{\ast \ast \ast} + \frac{\varepsilon_g}{\rho_g} \nabla {p_g}^{n+1/2,\ast} \right)
-
-then defining
-
-.. math:: U^{n+1} = U^{\ast \ast \ast} - \frac{1}{\rho_g} \nabla \phi
-
-and
-
-.. math:: {p_g}^{n+1/2} = \phi
+ FluidTimeStep
+ Slopes
+ BuildingMacVelocities
diff --git a/docs/source/FluidTimeStep.rst b/docs/source/FluidTimeStep.rst
new file mode 100644
index 0000000000000000000000000000000000000000..5239d5c799d254a0efb6d22f5ed06df65f5e0c7d
--- /dev/null
+++ b/docs/source/FluidTimeStep.rst
@@ -0,0 +1,52 @@
+
+Fluid Time Step
+~~~~~~~~~~~~~~~
+
+In the absence of reactions, we assume that the fluid density is unchanged.
+
+We compute the fluid volume fraction directly from the particle locations.
+
+Thus here we focus on the discretization of the momentum equation
+
+In the predictor
+
+- Define :math:`U^{MAC,n}`, the face-centered (staggered) MAC velocity which is used for advection, using :math:`U_g^n`
+
+- Define an approximation to the new-time state, :math:`(\varepsilon_g \rho_g U_g)^{\ast}` by setting
+
+.. math:: (\varepsilon_g \rho_g U_g)^{\ast} &= (\varepsilon_g \rho_g U_g)^n -
+ \Delta t \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g) + \varepsilon_g \nabla {p_g}^{n-1/2} \right) \\ &+
+ \Delta t \left( \nabla \cdot \tau^n + \sum_p \beta_p (V_p - {U_g}^{\ast}) + \rho_g \varepsilon_g g \right)
+
+- Project :math:`U_g^{\ast}` by solving
+
+.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g U_g)^{\ast}+ {\varepsilon_g}{\rho_g} \nabla {p_g}^{n-1/2} \right)
+
+then defining
+
+.. math:: U_g^{\ast \ast} = U_g^{\ast} - \frac{1}{\rho_g} \nabla \phi
+
+and
+
+.. math:: {p_g}^{n+1/2, \ast} = \phi
+
+
+In the corrector
+
+- Define :math:`U^{MAC,\ast \ast}` at the "new" time using :math:`U_g^{\ast \ast}`
+
+- Define a new approximation to the new-time state, :math:`(\varepsilon_g \rho_g U_g)^{\ast \ast \ast}` by setting
+
+.. math:: (\varepsilon_g \rho_g U_g)^{\ast \ast \ast} &= (\varepsilon_g \rho_g U_g)^n - \frac{\Delta t}{2} \left( \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^n + \nabla \cdot (\varepsilon_g \rho_g U^{MAC} U_g)^{\ast \ast}\right) + \\ &+ \frac{\Delta t}{2} \left( \nabla \cdot \tau^n + \nabla \cdot \tau^{\ast \ast \ast} \right) + \Delta t \left( - \varepsilon_g \nabla {p_g}^{n+1/2,\ast} + \sum_p \beta_p (V_p - {U_g}^{\ast \ast \ast}) + \varepsilon_g \rho_g g \right)
+
+- Project :math:`U_g^{\ast \ast \ast}` by solving
+
+.. math:: \nabla \cdot \frac{\varepsilon_g}{\rho_g} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} (\varepsilon_g U_g)^{\ast \ast \ast} + \frac{\varepsilon_g}{\rho_g} \nabla {p_g}^{n+1/2,\ast} \right)
+
+then defining
+
+.. math:: U_g^{n+1} = U_g^{\ast \ast \ast} - \frac{1}{\rho_g} \nabla \phi
+
+and
+
+.. math:: {p_g}^{n+1/2} = \phi
diff --git a/docs/source/Fluids.rst b/docs/source/Fluids_Chapter.rst
similarity index 93%
rename from docs/source/Fluids.rst
rename to docs/source/Fluids_Chapter.rst
index 5a8a100fc5c6143f537684d08c8c904a9b29cf65..f65ccf2ad1ed485b83a1e411ed0bf2f592258c88 100644
--- a/docs/source/Fluids.rst
+++ b/docs/source/Fluids_Chapter.rst
@@ -7,7 +7,6 @@ Solving the Fluid Equations
===========================
.. toctree::
- :maxdepth: 1
FluidEquations
FluidTimeDiscretization
diff --git a/docs/source/GridCreation.rst b/docs/source/GridCreation.rst
new file mode 100644
index 0000000000000000000000000000000000000000..97ed0dbda21848e5aa2d658e747d6fcddbcc7fd0
--- /dev/null
+++ b/docs/source/GridCreation.rst
@@ -0,0 +1,64 @@
+.. role:: cpp(code)
+ :language: c++
+
+.. role:: fortran(code)
+ :language: fortran
+
+.. _sec:grid_creation:
+
+Grid Creation
+-------------
+
+To run MFiX-Exa you must specify :cpp:`n_cell` in the inputs file --
+this is the number of cells spanning the domain in each coordinate direction at level 0.
+
+Users often specify :cpp:`max_grid_size` as well. The default load balancing algorithm then divides the
+domain in every direction so that each grid is no longer than :cpp:`max_grid_size` in that direction.
+If not specified by the user, :cpp:`max_grid_size` defaults to 128 in 2D and 32 in 3D (in each coordinate direction).
+
+Another popular input is :cpp:`blocking_factor`. The value of :cpp:`blocking_factor`
+constrains grid creation in that in that each grid must be divisible by :cpp:`blocking_factor`.
+Note that both the domain (at each level) and :cpp:`max_grid_size` must be divisible by :cpp:`blocking_factor`
+If not specified by the user, :cpp:`blocking_factor` defaults to 8 in each coordinate direction.
+The typical purpose of :cpp:`blocking_factor` is to ensure that the grids will be
+sufficiently coarsenable for good multigrid performance.
+
+There is one more default behavior to be aware of. There is a boolean :cpp:`refine_grid_layout`
+that defaults to true but can be over-ridden at run-time.
+If :cpp:`refine_grid_layout` is true and the number of grids created is less than the number of processors
+(Ngrids < Nprocs), then grids will be further subdivided until Ngrids >= Nprocs.
+
+Caveat: if subdividing the grids to achieve Ngrids >= Nprocs would violate the
+:cpp:`blocking_factor` criterion then additional grids are not created and the
+number of grids will remain less than the number of processors
+
+Note that :cpp:`n_cell` must be given as three separate integers, one for each coordinate direction.
+
+However, :cpp:`max_grid_size` and :cpp:`blocking_factor` can be specified as a single value
+applying to all coordinate directions, or as separate values for each direction.
+
+ - if :cpp:`max_grid_size` (or :cpp:`blocking_factor`) is specified as multiple integers then the first
+ integer applies to level 0, the second to level 1, etc. If you don't specify as many
+ integers as there are levels, the final value will be used for the remaining levels.
+
+ - if different values of :cpp:`max_grid_size` (or :cpp:`blocking_factor`) are wanted for each coordinate direction,
+ then :cpp:`max_grid_size_x`, :cpp:`max_grid_size_y` and :cpp:`max_grid_size_z`
+ (or :cpp:`blocking_factor_x`, :cpp:`blocking_factor_y` and :cpp:`blocking_factor_z`) must be used.
+ If you don't specify as many integers as there are levels, the final value will be used for the remaining levels.
+
+ - note that :cpp:`max_grid_size` is just an upper bound; with :cpp:`n_cell = 48`
+ and :cpp:`max_grid_size = 32`, we will typically have one grid of length 32 and one of length 16.
+
+The grid creation process at level 0 proceeds as follows (if not using the KD-tree approach):
+
+#. The domain is initially defined by a single grid of size :cpp:`n_cell`.
+
+#. If :cpp:`n_cell` is greater than :cpp:`max_grid_size` then the grids are subdivided until
+ each grid is no longer than :cpp:`max_grid_size` cells on each side. The :cpp:`blocking_factor` criterion
+ (ie that the length of each side of each grid is divisible by :cpp:`blocking_factor` in that direction)
+ is satisfied during this process.
+
+#. Next, if :cpp:`refine_grid_layout = true` and there are more processors than grids
+ at this level, then the grids at this level are further divided until Ngrids >= Nprocs
+ (unless doing so would violate the :cpp:`blocking_factor` criterion).
+
diff --git a/docs/source/Inputs.rst b/docs/source/Inputs.rst
deleted file mode 100644
index 3f01cca91eb2b78cde69d4a4b35e943c828add03..0000000000000000000000000000000000000000
--- a/docs/source/Inputs.rst
+++ /dev/null
@@ -1,68 +0,0 @@
-.. _Chap:Inputs:
-
-Run-time Inputs
-===============
-
-The following inputs must be preceded by "amr."
-
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| File | Description | Type | Default |
-+=================+=======================================================================+=============+===========+
-| max_step | Maximum number of time steps to take | Int | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| stop_time | Maximum time to reach | Real | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| max_level | Maximum level of refinement allowed (0 when single-level) | Int | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| n_cell | Number of cells at level 0 in each coordinate direction | Int | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| plot_int | Frequency of plotfile output | Int | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| plot_file | Prefix to use for plotfile output | String | plt |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| check_int | Frequency of checkpoint output | Int | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| chk_file | Prefix to use for checkpoint output | String | chk |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| restart | Name of file from which to restart the run | String | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| max_grid_size_x | Maximum number of cells at level 0 in each grid in x-direction | Int | 32 |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| max_grid_size_y | Maximum number of cells at level 0 in each grid in y-direction | Int | 32 |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| max_grid_size_z | Maximum number of cells at level 0 in each grid in z-direction | Int | 32 |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-
-
-The following inputs must be preceded by "geometry."
-
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| File | Description | Type | Default |
-+=================+=======================================================================+=============+===========+
-| coord_sys | 0 for Cartesian | Int | 0 |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| is_periodic | 1 for true, 0 for false (one value for each coordinate direction) | Ints | 0 0 0 |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| prob_lo | Low corner of physical domain (physical not index space) | Reals | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-| prob_hi | High corner of physical domain (physical not index space) | Reals | None |
-+-----------------+-----------------------------------------------------------------------+-------------+-----------+
-
-
-The following inputs must be preceded by "mfix."
-
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| File | Description | Type | Default |
-+====================+=======================================================================+=============+===========+
-| fixed_dt | Should we use a fixed timestep? | Int | 0 |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| dt_max | Maximum value of dt if calculating with cfl | Real | 1.e14 |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| cfl | CFL constraint (dt < cfl * dx / u) if fixed_dt not 1 | Real | 0.5 |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| geometry | Which type of geometry are we using? | String | |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| explicit_diffusion | Should we use explicit or implicit diffusion? | Int | 1 |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
-| particle_init_type | How do we initialize the particles? "Auto" vs AsciiFile | String | AsciiFile |
-+--------------------+-----------------------------------------------------------------------+-------------+-----------+
diff --git a/docs/source/InputsCheckpoint.rst b/docs/source/InputsCheckpoint.rst
new file mode 100644
index 0000000000000000000000000000000000000000..13bae0e467fa43a1ed0247f203748e315aad72c5
--- /dev/null
+++ b/docs/source/InputsCheckpoint.rst
@@ -0,0 +1,18 @@
+.. _Chap:InputsCheckpoint:
+
+Checkpoint/Restart
+==================
+
+The following inputs must be preceded by "amr" and control checkpoint/restart.
+
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++==================+=======================================================================+=============+===========+
+| restart | If present, then the name of file to restart from | String | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| check_int | Frequency of checkpoint output; | Int | -1 |
+| | if -1 then no checkpoints will be written | | |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| check_file | Prefix to use for checkpoint output | String | chk |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+
diff --git a/docs/source/InputsDrag.rst b/docs/source/InputsDrag.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6dbb7cf4d048db6349703f2ca742499a27d6d47b
--- /dev/null
+++ b/docs/source/InputsDrag.rst
@@ -0,0 +1,128 @@
+Drag Types
+==========
+
+The following inputs must be preceded by "mfix."
+
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++===================+=======================================================================+=============+===========+
+| drag_type | Which drag model to use | String | None |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The options currently supported in mfix are :cpp:`WenYu`, :cpp:`Gidaspow`, :cpp:`BVK2`, or :cpp:`UserDrag`.
+
+If one of these is not specified, the code will abort with
+
+.. highlight:: c++
+
+::
+
+ amrex::Abort::0::"Don't know this drag type!!!
+
+The drag models are defined in src/src_des/des_drag_K.H
+
+If the user wishes to use their own drag model, they must
+
+ * specify :cpp:`mfix.drag_type = UserDrag` in the inputs file
+
+ * provide the code in the ComputeDragUser routine in a local usr_drag.cpp file.
+ An example can be found in tests/DEM06-x.
+
+With the variables defined as follows:
+
+ .. code:: shell
+
+ * \brief Returns: the calculated drag coefficient.
+ *
+ * Inputs:
+ * EPg - gas volume fraction
+ * Mug - gas laminar viscosity
+ * ROpg - gas density * EP_g
+ * vrel - magnitude of gas-solids relative velocity
+ * DPM - particle diamater of solids phase M
+ * DPA - average particle diameter
+ * PHIS - solids volume fraction of solids phases
+ * fvelx - x component of the fluid velocity at the particle position
+ * fvely - y component of the fluid velocity at the particle position
+ * fvelz - z component of the fluid velocity at the particle position
+ * i, j, k - particle cell indices
+ * pid - particle id number
+ */
+
+The WenYu model is defined as
+
+ .. code:: shell
+
+ RE = (Mug > 0.0) ? DPM*vrel*ROPg/Mug : DEMParams::large_number;
+
+ if (RE <= 1000.0)
+ {
+ C_d = (24.0/(RE+DEMParams::small_number)) * (1.0 + 0.15*std::pow(RE, 0.687));
+ }
+ else
+ {
+ C_d = 0.44;
+ }
+
+ if (RE < DEMParams::eps) return 0.0;
+ return 0.75 * C_d * vrel * ROPg * std::pow(EPg, -2.65) / DPM;
+
+The Gidaspow model is defined as
+
+ .. code:: shell
+
+ ROg = ROPg / EPg;
+
+ RE = (Mug > 0.0) ? DPM*vrel*ROPg/Mug : DEMParams::large_number;
+
+ // Dense phase - EPg <= 0.8
+ Ergun = 150.0*(1.0 - EPg)*Mug / (EPg*DPM*DPM) + 1.75*ROg*vrel/DPM;
+
+ // Dilute phase - EPg > 0.8
+ if (RE <= 1000.0)
+ {
+ C_d = (24.0/(RE+DEMParams::small_number)) * (1.0 + 0.15*std::pow(RE, 0.687));
+ }
+ else
+ {
+ C_d = 0.44;
+ }
+
+ WenYu = 0.75*C_d*vrel*ROPg*std::pow(EPg, -2.65) / DPM;
+
+ // switch function
+ PHI_gs = atan(150.0*1.75*(EPg - 0.8))/M_PI / DPM;
+
+ // blend the models
+ if (RE < DEMParams::eps) return 0.0;
+ return (1.0 - PHI_gs)*Ergun + PHI_gs*WenYu;
+
+The Gidaspow model is defined as
+
+ .. code:: shell
+
+ amrex::Real RE = (Mug > 0.0) ? DPA*vrel*ROPg/Mug : DEMParams::large_number;
+
+ if (RE > DEMParams::eps)
+ {
+ oEPgfour = 1.0 / EPg / EPg / EPg / EPg;
+
+ // eq(9) BVK J. fluid. Mech. 528, 2005
+ // (this F_Stokes is /= of Koch_Hill by a factor of ep_g)
+ F_Stokes = 18.0*Mug*EPg/DPM/DPM;
+
+ F = 10.0*PHIS/EPg/EPg + EPg*EPg*(1.0 + 1.5*sqrt(PHIS));
+
+ F += RE*(0.11*PHIS*(1.0+PHIS) - 4.56e-3*oEPgfour +
+ std::pow(RE, -0.343)*(0.169*EPg + 6.44e-2*oEPgfour));
+
+ // F += 0.413*RE/(24.0*EPg*EPg) *
+ // (1.0/EPg + 3.0*EPg*PHIS + 8.4/std::pow(RE, 0.343)) /
+ // (1.0 + std::pow(10.0, 3.0*PHIS)/std::pow(RE, 0.5 + 2.0*PHIS));
+
+ return F*F_Stokes;
+ }
+ else
+ {
+ return 0.0;
+ }
diff --git a/docs/source/InputsInitialization.rst b/docs/source/InputsInitialization.rst
new file mode 100644
index 0000000000000000000000000000000000000000..296b3860c13bf498e534391e9e0ae55e4bfabbe9
--- /dev/null
+++ b/docs/source/InputsInitialization.rst
@@ -0,0 +1,42 @@
+.. _Chap:InputsInitialization:
+
+Initialization
+==============
+
+The following inputs must be preceded by "amr" and determine how we initialize a calculation:
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
+| input_deck | Read physical data from this file | String | mfix.dat |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| restart | If set, then restart from this file rather than from scratch | String | None |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| repl_x | Replicate initial data by this factor in the x-direction | Int | 1 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| repl_y | Replicate initial data by this factor in the y-direction | Int | 1 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| repl_z | Replicate initial data by this factor in the z-direction | Int | 1 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+
+
+The following inputs must be preceded by "mfix" and determine how we initialize a calculation:
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| particle_init_type | How do we initialize the particles? "Auto" vs AsciiFile | String | AsciiFile |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| do_initial_proj | Should we do the initial projection? | Bool | True |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| initial_iterations | How many pressure iterations before starting the first timestep | Int | 3 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+
+The following inputs must be preceded by "particles."
+
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++====================+===========================================================================+=============+===========+
+| removeOutOfRange | Should we remove particles at initialization that are touching the wall | Int | 1 |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
diff --git a/docs/source/InputsLoadBalancing.rst b/docs/source/InputsLoadBalancing.rst
new file mode 100644
index 0000000000000000000000000000000000000000..25a8d322a5c698f24ba1024f1ec355b28ab01a1d
--- /dev/null
+++ b/docs/source/InputsLoadBalancing.rst
@@ -0,0 +1,59 @@
+.. _Chap:InputsLoadBalancing:
+
+Gridding and Load Balancing
+===========================
+
+The following inputs must be preceded by "amr" and determine how we create the grids and how often we regrid.
+
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+===========+
+| regrid_int | How often to regrid (in number of steps at level 0) | Int | -1 |
+| | if regrid_int = -1 then no regridding will occur | | |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_grid_size_x | Maximum number of cells at level 0 in each grid in x-direction | Int | 32 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_grid_size_y | Maximum number of cells at level 0 in each grid in y-direction | Int | 32 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_grid_size_z | Maximum number of cells at level 0 in each grid in z-direction | Int | 32 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| blocking_factor_x | Each grid must be divisible by blocking_factor_x in x-direction | Int | 8 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| blocking_factor_y | Each grid must be divisible by blocking_factor_y in y-direction | Int | 8 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+| blocking_factor_z | Each grid must be divisible by blocking_factor_z in z-direction | Int | 8 |
++----------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "particles"
+
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++===================+=======================================================================+=============+===========+
+| max_grid_size_x | Maximum number of cells at level 0 in each grid in x-direction | Int | 32 |
+| | for grids in the ParticleBoxArray if dual_grid is true | | |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_grid_size_y | Maximum number of cells at level 0 in each grid in y-direction | Int | 32 |
+| | for grids in the ParticleBoxArray if dual_grid is true | | |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_grid_size_z | Maximum number of cells at level 0 in each grid in z-direction | Int | 32 |
+| | for grids in the ParticleBoxArray if dual_grid is true. | | |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "mfix" and determine how we load balance:
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
+| dual_grid | If true then use the "dual_grid" approach for load balancing | Bool | False |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| load_balance_fluid | Only relevant if (dual_grid); if so do we also regrid mesh data | Int | 1 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| load_balance_type | What strategy to use for load balancing | String | FixedSize |
+| | Options are "FixedSize", "KDTree" or "Knapsack" | | |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| knapsack_weight_type | What weighting function to use if using Knapsack load balancing | String | RunTimeCosts |
+| | Options are "RunTimeCosts" or "NumParticles"" | | |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| knapsack_nmax | Maximum number of grids per MPI process if using knapsack algorithm | Int | 128 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+
diff --git a/docs/source/InputsMonitors.rst b/docs/source/InputsMonitors.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3d54a4a1d20c83d2984c9b314323d715603c44d2
--- /dev/null
+++ b/docs/source/InputsMonitors.rst
@@ -0,0 +1,41 @@
+.. _Chap:InputsMonitors:
+
+Monitors
+========
+
+The following inputs must be preceded by "amr" and control whether to compute
+monitors, i.e., spatial averages, and how often to output the results.
+n is the number of monitors implicitly defined by the size of avg_region_x_w.
+
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++==================+=======================================================================+=============+===========+
+| avg_int | Interval, in number of CFD dt's, to write output | Int | -1 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_file | Base file name which is appended with the data type (vel_p_, p_g_, | String | avg_region|
+| | ep_g_ or vel_g_), the number of this type of averaging, and the .csv | | |
+| | file extension | | |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_vel_p | Average and save particle velocity (if 1) | n*Int | 0 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_p_g | Average and save fluid-phase pressure (if 1) | n*Int | 0 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_ep_g | Average and save fluid-phase volume fraction (if 1) | n*Int | 0 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_vel_g | Average and save fluid-phase velocity (if 1) | n*Int | 0 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_x_w | Lower bound of averaging region in x-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_x_e | Upper bound of averaging region in x-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_y_s | Lower bound of averaging region in y-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_y_n | Upper bound of averaging region in y-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_z_b | Lower bound of averaging region in z-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| avg_region_z_t | Upper bound of averaging region in z-direction | n*Real | None |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+
+
diff --git a/docs/source/InputsMultigrid.rst b/docs/source/InputsMultigrid.rst
new file mode 100644
index 0000000000000000000000000000000000000000..7ed2a9fb4dbe85276443f783b0327d4c6cceb9bc
--- /dev/null
+++ b/docs/source/InputsMultigrid.rst
@@ -0,0 +1,76 @@
+.. _Chap:InputsMultigrid:
+
+Multigrid Inputs
+================
+
+The following inputs can be set directly in the AMReX solver classes but we
+set them via the MFiX-Exa routines because we may want different inputs for the
+different solvers called by MFiX-Exa.
+
+These control the nodal projection and must be preceded by "mfix":
+
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_verbose | Verbosity of multigrid solver in nodal projection | Int | 0 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_cg_verbose | Verbosity of BiCGStab solver in nodal projection | Int | 0 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_rtol | Relative tolerance in nodal projection | Real | 1.e-11 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_atol | Absolute tolerance in nodal projection | Real | 1.e-14 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_maxiter | Maximum number of iterations in the nodal projection | Int | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_cg_maxiter | Maximum number of iterations in the nodal projection | Int | |
+| | bottom solver if using bicg, cg, bicgcg or cgbicg | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_max_coarsening_level | Maximum number of coarser levels to allowin the nodal projection | Int | |
+| | If set to 0, the bottom solver will be called at the current level | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| bottom_solver_type | Which bottom solver to use in the nodal projection | String | bicgcg |
+| | Options are bicgstab, cg, cgbicg, smoother or hypre | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+
+These control the MAC projection and must be preceded by "mac":
+
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++=========================+=======================================================================+=============+==============+
+| mg_verbose | Verbosity of multigrid solver in MAC projection | Int | 0 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_cg_verbose | Verbosity of BiCGStab solver in MAC projection | Int | 0 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_rtol | Relative tolerance in MAC projection | Real | 1.e-11 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_atol | Absolute tolerance in MAC projection | Real | 1.e-14 |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_maxiter | Maximum number of iterations in the MAC projection | Int | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_cg_maxiter | Maximum number of iterations in the MAC projection | Int | |
+| | bottom solver if using bicg, cg, bicgcg or cgbicg | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_max_coarsening_level | Maximum number of coarser levels to allow in the nodal projection | Int | |
+| | If set to 0, the bottom solver will be called at the current level | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+| bottom_solver_type | Which bottom solver to use in the MAC projection | String | bicgcg |
+| | Options are bicgstab, cg, cgbicg, smoother or hypre | | |
++-------------------------+-----------------------------------------------------------------------+-------------+--------------+
+
+These control the diffusion solver and must be preceded by "diff":
+The following inputs must be preceded by "diff"
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
+| mg_verbose | Verbosity of linear solver for diffusion solve | Int | 0 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_cg_verbose | Verbosity of BiCGStab solver in diffusion solve | Int | 0 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_rtol | Relative tolerance in diffusion solve | Real | 1.e-11 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| mg_atol | Absolute tolerance in diffusion solve | Real | 1.e-14 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| bottom_solver_type | Which bottom solver to use in the diffusion solve | String | bicgcg |
+| | Options are bicgstab, cg, cgbicg, smoother or hypre | | |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
diff --git a/docs/source/InputsPlotFiles.rst b/docs/source/InputsPlotFiles.rst
new file mode 100644
index 0000000000000000000000000000000000000000..7bc928c3c789012704ffda859f1f57be35e39133
--- /dev/null
+++ b/docs/source/InputsPlotFiles.rst
@@ -0,0 +1,79 @@
+.. _Chap:InputsPlotfiles:
+
+Plotfiles and Other Output
+==========================
+
+The following inputs must be preceded by "amr" and control frequency and naming of plotfile generation as well
+as whether the EB geometry or level set should be written out, and if the particles should be written out in Ascii
+format (for debugging).
+
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++=====================+=======================================================================+=============+===========+
+| plot_int | Frequency of plotfile output; | Int | -1 |
+| | if -1 then no plotfiles will be written | | |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plotfile_on_restart | Should we write a plotfile when we restart (only used if plot_int>0) | Bool | False |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plot_file | Prefix to use for plotfile output | String | plt |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| write_ls | Should we write a plotfile holding the level set and volfrac? | Bool | False |
+| | If true, it will only be written once,after initialization or restart | | |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| write_eb_surface | Should we write out the EB geometry in vtp format | Bool | False |
+| | If true, it will only be written once,after initialization or restart | | |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| par_ascii_file | Prefix to use for ascii particle output | String | par |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| par_ascii_int | Frequency of ascii particle output; | Int | -1 |
+| | if -1 then no plotfiles will be written | | |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "amr" and control what variables will be written in plotfiles.
+
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++=====================+=======================================================================+=============+===========+
+| plt_regtest | Save all variables to plot file (overrides all other IO flags) | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_vel_g | Save fluid velocity data to plot file | Int | 1 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_ep_g | Save fluid volume fraction to plot file | Int | 1 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_p_g | Save fluid pressure to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_ro_g | Save fluid density to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_mu_g | Save fluid viscosity to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_diveu | Save div(ep_g . u) to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_volfrac | Save Eulerian grid volume fraction (from cut cells) to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_gradp_g | Save gradient of pressure filed to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_vort | Save vorticity to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_vel_p | Save particle velocity to plot file | Int | 1 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_radius | Save particle radius to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_volume | Save particle volume to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_volume | Save particle volume to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_mass | Save particle mass to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_ro_p | Save particle density to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_omoi | Save (one divided by the) particle momentum of inertia to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_mass | Save particle mass to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_omega_p | Save particle angular velocity to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_drag_p | Save particle drag force to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
+| plt_phase | Save particle type to plot file | Int | 0 |
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
++---------------------+-----------------------------------------------------------------------+-------------+-----------+
diff --git a/docs/source/InputsProblemDefinition.rst b/docs/source/InputsProblemDefinition.rst
new file mode 100644
index 0000000000000000000000000000000000000000..bd6c14c7f0dca0719e8e66913405d0b933efce85
--- /dev/null
+++ b/docs/source/InputsProblemDefinition.rst
@@ -0,0 +1,68 @@
+Problem Definition
+==================
+
+The following inputs must be preceded by "amr."
+
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++===================+=======================================================================+=============+===========+
+| n_cell | Number of cells at level 0 in each coordinate direction | Int Int Int | None |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+| max_level | Maximum level of refinement allowed (0 when single-level) | Int | None |
++-------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "geometry."
+
++-----------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++=================+=======================================================================+=============+===========+
+| coord_sys | 0 for Cartesian | Int | 0 |
++-----------------+-----------------------------------------------------------------------+-------------+-----------+
+| is_periodic | 1 for true, 0 for false (one value for each coordinate direction) | Ints | 0 0 0 |
++-----------------+-----------------------------------------------------------------------+-------------+-----------+
+| prob_lo | Low corner of physical domain (physical not index space) | Reals | None |
++-----------------+-----------------------------------------------------------------------+-------------+-----------+
+| prob_hi | High corner of physical domain (physical not index space) | Reals | None |
++-----------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "mfix."
+
++----------------------+-------------------------------------------------------------------------+-------------+--------+
+| | Description | Type | Default |
++======================+=========================================================================+==========+===========+
+| geometry | Which type of EB geometry are we using? | String | |
++----------------------+-------------------------------------------------------------------------+----------+-----------+
+| levelset__refinement | Refinement factor of levelset resolution relative to level 0 resolution | Int | 1 !
++----------------------+-------------------------------------------------------------------------+----------+-----------+
+
+Setting basic EB walls can be specified by inputs preceded by "xlo", "xhi", "ylo", "yhi", "zlo", and "zhi"
+
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++====================+===========================================================================+=============+===========+
+| type | Used to define boundary type. Available options include: | String | None |
+| | | | |
+| | * 'pi' or 'pressure_inflow' | | |
+| | * 'po' or 'pressure_outflow' | | |
+| | * 'mi' or 'mass_inflow' | | |
+| | * 'nsw' or 'no_slip_wall' | | |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| pressure | Sets boundary pressure for pressure inflows, outflows and mass inflows | Real | None |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| velocity | Sets boundary velocity for mass inflows | Real | None |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+| location | Specifies an offset from the domain boundary for no-slip walls | Real | None |
++--------------------+---------------------------------------------------------------------------+-------------+-----------+
+
+To specify multiple mass inflows (e.g., define a jet and uniform background flow), provide multiple velocities for the region and define the physical extents of the sub-region. The first velocity is applied to the entire flow plane. Subsequent velocities are successively applied to the specified sub-regions. If multiple sub-regions overlap, the velocity of last specified region is used. An example of a uniform mass inflow with a square-jet centered at (0.5x0.5) is given below.
+
+
+.. code-block:: none
+
+ xlo.type = "mi"
+ xlo.velocity = 0.01 0.10
+
+ xlo.ylo = 0.25
+ xlo.yhi = 0.75
+ xlo.zlo = 0.25
+ xlo.zhi = 0.75
diff --git a/docs/source/InputsTimeStepping.rst b/docs/source/InputsTimeStepping.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b201e05ff3cb55a99e0d12200c5b124d6c98bfda
--- /dev/null
+++ b/docs/source/InputsTimeStepping.rst
@@ -0,0 +1,98 @@
+.. sec:InputsTimeStepping:
+
+Time Stepping
+=============
+
+The following inputs must be preceded by "amr." Note that if both are specified, both criteria
+are used and the simulation still stop when the first criterion is hit. In the case of unsteady flow,
+the simulation will stop when either the number of steps reaches max_step or time reaches stop_time.
+In the case of unsteady flow, the simulation will stop when either the tolerance (difference between
+subsequent steps) is reached or the number of iterations reaches the maximum number specified.
+
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| | Description | Type | Default |
++==================+=======================================================================+=============+===========+
+| max_step | Maximum number of time steps to take | Int | -1 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+| stop_time | Maximum time to reach | Real | -1.0 |
++------------------+-----------------------------------------------------------------------+-------------+-----------+
+
+The following inputs must be preceded by "mfix."
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
+| fixed_dt | Should we use a fixed timestep? | Int | 0 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| dt_min | Abort if dt gets smaller than this value | Real | 1.e-6 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| dt_max | Maximum value of dt if calculating with cfl | Real | 1.e14 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| cfl | CFL constraint (dt < cfl * dx / u) if fixed_dt not 1 | Real | 0.5 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+
+The following inputs must be preceded by "mfix" and are only relevant if running a problem to steady state.
+Currently, the criterion for setting "steady_state" to true is if "dt" is undefined in mfix.dat
+
++-----------------------+-----------------------------------------------------------------------+-------------+------------+
+| | Description | Type | Default |
++=======================+=======================================================================+=============+============+
+| steady_state | Are we running a steady-state calculation? | Int | 0 |
++-----------------------+-----------------------------------------------------------------------+-------------+------------+
+| steady_state_tol | Tolerance for checking if we have reached steady state | Real | None |
+| | | | |
+| | (Must be set if steady_state_tol = 1) | | |
++-----------------------+-----------------------------------------------------------------------+-------------+------------+
+| steady_state_max_iter | Maximum number of allowed iterations to converge to steady state | Int | 100000000 |
++-----------------------+-----------------------------------------------------------------------+-------------+------------+
+
+Setting the Time Step
+---------------------
+
+There are several ways that the inputs are used to determine what time step
+is used in the evolution of the fluid-particle system in MFiX-Exa.
+
+1) In a pure particle case, the :cpp:`mfix.fixed_dt`, if specified, is only used to determine the frequency
+of outputs, it has no effect on the "dtsolid" used in the particle evaluation. If you do not specify a positive
+value of :cpp:`mfix.fixed_dt` then the code will abort.
+
+.. highlight:: c++
+
+::
+
+ amrex::Abort::0::If running particle-only must specify fixed_dt in the inputs file !!!
+
+The particle time step "dtsolid" is determined by computing the collision time "tcoll" from particle properties,
+then setting "dtsolid" to be "tcoll / 50".
+
+2) In a pure fluid case, there are two options:
+
+ * If you want to fix the dt, simply set :cpp:`mfix.fixed_dt = XXX` and the fluid time
+ step will always be that number.
+
+ * If you want to let the code determine the appropriate time step using the advective CFL
+ condition, then set :cpp:`mfix.cfl = 0.7` for example, and the fluid time step will
+ be computed to be dt = 0.5 * dx / max(vel).
+
+ * If dt as computed in the compute_dt routine is smaller than the user-specified
+ ```mfix.dt_min``` then the code will abort:
+
+ .. highlight:: c++
+
+ ::
+
+ amrex::Abort::0::"Current dt is smaller than dt_min !!!
+
+ * If dt as computed in the compute_dt routine is larger than the user-specified
+ ```mfix.dt_max``` then dt will be set to the minimum of its computed value and dt_max
+
+ * Note that the cfl defaults to 0.5 so it does not have to be set in the inputs file. If neither
+ :cpp:`mfix.cfl` nor :cpp:`fixed_dt` is set, then default value of cfl will be used.
+ If :cpp:`mfix.fixed_dt` is set, then it will override the cfl option whether
+ :cpp:`mfix.cfl` is set or not.
+
+These options apply to steady state calculations as well as unsteady runs.
+
+3) In a coupled particle-fluid case, dt is determined as in the pure-fluid case. In this case
+ the particle time step "subdt" is first computed as in the particle-only case ("dtsolid"),
+ then is adjusted so that an integral number of particle steps fit into a single fluid time step.
diff --git a/docs/source/InputsVerbosity.rst b/docs/source/InputsVerbosity.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a5f693a39ba3d059aebaedf607be7c7156550bcc
--- /dev/null
+++ b/docs/source/InputsVerbosity.rst
@@ -0,0 +1,15 @@
+.. _Chap:InputsVerbosity:
+
+Verbosity
+=========
+
+The following inputs must be preceded by "mfix."
+
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| | Description | Type | Default |
++======================+=======================================================================+=============+==============+
+| verbose | Verbosity in MFiX-Exa routines | Int | 0 |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+| ooo_debug | If true then print the name of the routine we are in | Bool | False |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
+
diff --git a/docs/source/Inputs_Chapter.rst b/docs/source/Inputs_Chapter.rst
new file mode 100644
index 0000000000000000000000000000000000000000..20de996627a8f99e2ad3708bb090c77c7496036e
--- /dev/null
+++ b/docs/source/Inputs_Chapter.rst
@@ -0,0 +1,18 @@
+.. _Chap:Inputs:
+
+Run-time Inputs
+===============
+
+.. toctree::
+ :maxdepth: 1
+
+ InputsProblemDefinition
+ InputsDrag
+ InputsTimeStepping
+ InputsInitialization
+ InputsLoadBalancing
+ InputsMultigrid
+ InputsPlotFiles
+ InputsCheckpoint
+ InputsMonitors
+ InputsVerbosity
diff --git a/docs/source/Introduction.rst b/docs/source/Introduction.rst
index fd7c5a6abf611f4d4534e3b516e9d05bdda555ed..cd89521dfae345fd6c87dc9874021d3956349cbc 100644
--- a/docs/source/Introduction.rst
+++ b/docs/source/Introduction.rst
@@ -30,5 +30,9 @@ time discretizations differ. Specifically,
- Plotfile format supported by AmrVis, VisIt, ParaView, and yt.
-MFiX-Exa is being developed at NETL, LBNL, and CU as part of DOE's Exascale Computing Project.
+MFiX-Exa is being developed at NETL and LBNL as part of the U.S. Department of Energy's
+Exascale Computing Project (ECP).
+
+MFiX-Exa heavily leverages AMReX (see https://amrex-codes.github.io/) which is also supported by
+ECP as part of the AMReX Co-Design Center.
diff --git a/docs/source/LoadBalancing.rst b/docs/source/LoadBalancing.rst
new file mode 100644
index 0000000000000000000000000000000000000000..7496e3934c3b7ea16081b9e1f1e97057350dc60d
--- /dev/null
+++ b/docs/source/LoadBalancing.rst
@@ -0,0 +1,31 @@
+.. role:: cpp(code)
+ :language: c++
+
+.. role:: fortran(code)
+ :language: fortran
+
+.. _sec:load_balancing:
+
+Load Balancing
+--------------
+
+The process of load balancing is typically independent of the process of grid creation;
+the inputs to load balancing are a given set of grids with a set of weights
+assigned to each grid. (The exception to this is the KD-tree approach in which the
+grid creation process is governed by trying to balance the work in each grid.)
+
+Single-level load balancing algorithms are sequentially applied to each AMR level independently,
+and the resulting distributions are mapped onto the ranks taking into account the weights
+already assigned to them (assign heaviest set of grids to the least loaded rank)
+
+Options supported by AMReX include:
+
+- Knapsack: the default weight of a grid in the knapsack algorithm is the number of grid cells,
+ but AMReX supports the option to pass an array of weights – one per grid – or alternatively
+ to pass in a MultiFab of weights per cell which is used to compute the weight per grid
+
+- SFC: enumerate grids with a space-filling Z-morton curve, then partition the
+ resulting ordering across ranks in a way that balances the load
+
+- Round-robin: sort grids and assign them to ranks in round-robin fashion -- specifically
+ FAB *i* is owned by CPU *i*%N where N is the total number of MPI ranks.
diff --git a/docs/source/ManagingGridHierarchy_Chapter.rst b/docs/source/ManagingGridHierarchy_Chapter.rst
new file mode 100644
index 0000000000000000000000000000000000000000..deb241fe679a48711cb3cfffb4f9af6f739cd863
--- /dev/null
+++ b/docs/source/ManagingGridHierarchy_Chapter.rst
@@ -0,0 +1,37 @@
+.. role:: cpp(code)
+ :language: c++
+
+Gridding and Load Balancing
+===========================
+
+MFiX-Exa has a great deal of flexibility when it comes to how to decompose the
+computational domain into individual rectangular grids, and how to distribute
+those grids to MPI ranks. There can be grids of different sizes,
+more than one grid per MPI rank, and different strategies for distributing the grids to MPI ranks.
+
+We use the phrase "load balancing" here to refer to the combined process
+of grid creation (and re-creation when regridding) and distribution of grids to MPI ranks.
+
+See :ref:`sec:grid_creation` for grids are created, i.e. how the :cpp:`BoxArray` on which
+:cpp:`MultiFabs` will be built is defined at each level.
+
+See :ref:`sec:load_balancing` for the strategies AMReX supports for distributing
+grids to MPI ranks, i.e. defining the :cpp:`DistributionMapping` with which
+:cpp:`MultiFabs` at that level will be built.
+
+MFiX-Exa also allows for the "dual grid approach", in which mesh and particle data are allocated
+on different box layouts with different mappings to MPI ranks. This option is enabled
+by setting :cpp:`amr.dual_grid = 1` in the inputs file.
+See :ref:`sec:dual_grid` for more about this approach.
+
+When running on multicore machines with OpenMP, we can also control the distribution of
+work by setting the size of grid tiles (by defining :cpp:`fabarray_mfiter.tile_size`), and if relevant, of
+particle tiles (by defining :cpp:`particle.tile_size`). We can also specify the strategy for assigning
+tiles to OpenMP threads. See :ref:`sec:basics:mfiter:tiling:` for more about tiling.
+
+.. toctree::
+ :maxdepth: 1
+
+ GridCreation
+ DualGrid
+ LoadBalancing
diff --git a/docs/source/NightlyTests.rst b/docs/source/NightlyTests.rst
new file mode 100644
index 0000000000000000000000000000000000000000..aab9423c5efbed349589a80a662682bcf2c01802
--- /dev/null
+++ b/docs/source/NightlyTests.rst
@@ -0,0 +1,114 @@
+.. _Chap:NightlyTesting :
+
+Nightly Tests
+=============
+
+The following regression tests are run nightly with MFiX-Exa. The plotfiles generated in each night's test
+are compared with the benchmark plotfiles using the AMReX :cpp:`fcompare` utility to compare the mesh data
+and :cpp:`particle_compare` to compare the particle data.
+
+The results of these tests can be found at https://ccse.lbl.gov/pub/RegressionTesting/MFIX-Exa/
+
+Below Ng = number of grids, Npa = number of particles, Np = number of MPI ranks.
+
+"Auto" means the particles were generated automatically with the random number
+generator; if "Auto" is not specified the particle data were read in from "particle_input.dat"
+
+These first tests have both fluid and particles and are run in rectangular geometries;
+all tests except DEM06 use drag type "BVK2".
+
+"NSW" means "No Slip Wall" and "Per" is "periodic."
+"MI/PO" refers to Mass Inflow at the low end of the domain and Pressure Outflow at the high end.
+
++-------------------+----+--------+------+-------+----+----+----------------------+
+| Test | nx | bc_x | EB | Npa | Ng | Np | What does this test? |
+| | ny | bc_y | | | | | |
+| | nz | bc_z | | | | | |
++===================+====+========+======+=======+====+====+======================+
+| BENCH01 | 32 | Per | None | 5005 | 1 | 1 | Triply periodic |
+| Size0001 | 32 | Per | | | | | |
+| | 32 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH01 | 64 | Per | None | 40040 | 8 | 4 | Replicate |
+| Size0001 | 64 | Per | | | | | |
+| replicate | 64 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH01 | 32 | Per | None | 5005 | 8 | 4 | Restart |
+| Size0001 | 32 | Per | | | | | |
+| restart | 32 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH02 | 10 | Per | None | 1611 | 1 | 1 | Mixed NSW / Per |
+| Size0001 | 10 | NSW | | | | | |
+| | 10 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH02 | 10 | NSW | None | 1611 | 1 | 1 | NSW on all faces |
+| Size0001 | 10 | NSW | | | | | |
+| walls | 10 | NSW | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH03 | 4 | Per | None | 2500 | 1 | 1 | Mixed MI/PO + Per |
+| Size0001 | 50 | MI/PO | | | | | |
+| | 4 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| BENCH04 | 4 | Per | None | 224 | 1 | 1 | Triply periodic |
+| Size0001 | 50 | Per | | | | | |
+| | 4 | Per | | | | | |
++-------------------+----+--------+------+-------+----+----+----------------------+
+| DEM06 | 5 | Per | None | 1 | 10 | 4 | Single particle |
+| z multiple | 5 | Per | | | | | falling in fluid |
+| | 50 | MI/PO | | | | | (user_drag) |
++-------------------+----+--------+------+-------+----+----+----------------------+
+
+This second set of tests have both fluid and particles and are run in cylindrial geometries
+interior to the domain boundaries; they also use drag type "BVK2". Here "IGN" means
+those domain boundaries should be ignored because they are outside the EB boundary.
+
++-------------------+----+-------+------+--------+----+----+----------------------+
+| Test | nx | bc_x | EB | Npa | Ng | Np | What does this test? |
+| | ny | bc_y | | | | | |
+| | nz | bc_z | | | | | |
++===================+====+=======+======+========+====+====+======================+
+| BENCH05 | 40 | MI/PO | Cyl | 7949 | 4 | 4 | EB in parallel |
+| Size0008 | 10 | IGN | | Auto | | | |
+| | 10 | IGN | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| BENCH05 | 40 | MI/PO | Cyl | 7968 | 4 | 1 | EB in serial |
+| Size0008 | 10 | IGN | | Auto | | | |
+| serial | 10 | IGN | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| BENCH05 | 40 | MI/PO | Cyl | 36672 | 16 | 4 | Regrid & dual grid |
+| Size0008 | 20 | IGN | | Auto | | | |
+| medium | 20 | IGN | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| BENCH05 | 40 | MI/PO | Cyl | 157106 | 16 | 4 | Regrid & dual grid |
+| Size0008 | 40 | IGN | | Auto | | | |
+| wide | 40 | IGN | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| BENCH06 | 40 | Per | Cyl | 627 | 4 | 1 | EB |
+| Size0008 | 10 | IGN | | Auto | | | with periodic |
+| serial | 10 | IGN | | | | | serial |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| BENCH06 | 40 | Per | Cyl | 624 | 4 | 4 | EB |
+| Size0008 | 10 | IGN | | Auto | | | with periodic |
+| | 10 | IGN | | | | | parallel |
++-------------------+----+-------+------+--------+----+----+----------------------+
+
+
+This third set of tests is particles-only in rectangular geometries.
+
++-------------------+----+-------+------+--------+----+----+----------------------+
+| Test | nx | bc_x | EB | Npa | Ng | Np | What does this test? |
+| | ny | bc_y | | | | | |
+| | nz | bc_z | | | | | |
++===================+====+=======+======+========+====+====+======================+
+| DEM01 | 4 | NSW | None | 1 | 1 | 1 | Particle only |
+| x single | 4 | NSW | | | | | |
+| | 4 | NSW | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| DEM03 | 5 | Per | None | 2 | 1 | 1 | Particles only |
+| z single | 5 | Per | | | | | |
+| | 2 | NSW | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
+| DEM04 | 4 | NSW | None | 1 | 1 | 1 | Particles only |
+| z single | 4 | Per | | | | | |
+| | 4 | Per | | | | | |
++-------------------+----+-------+------+--------+----+----+----------------------+
diff --git a/docs/source/ParticlesOnGpus.rst b/docs/source/ParticlesOnGpus.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a20bd0561004c02ceeda83ae51aa13924c0efdff
--- /dev/null
+++ b/docs/source/ParticlesOnGpus.rst
@@ -0,0 +1,66 @@
+Particles on GPUs
+==========================
+
+The particle components of MFIX-Exa are a natural candidate for offloading to the GPU.
+The particle kernels are compute-intensive and can in principle be processed asynchronously with parts of the fluid advance.
+The core components of the particle method in MFIX-Exa are:
+
+- Neighbor List Construction
+- Particle-Particle Collisions
+- Particle-Wall Collisions
+
+Of these operations, the neighbor list construction requires the most care.
+A neighbor list is a pre-computed list of all the neighbors a given particle can interact with over the next *n* timesteps.
+Neighbor lists are usually constructed by binning the particles by an interaction distance,
+and then performing the N\ :sup:`2` distance check only on the particles in neighboring bins. In detail, the CPU version of the neighbor list algorithm is as follows:
+
+- For each tile on each level, loop over the particles, identifying the bin it belongs to.
+- Add the particle to a linked-list for the cell that `owns` it.
+- For each cell, loop over all the particles, and then loop over all potential collisions partners in the neighboring cells.
+- If a collision partner is close enough, add it to that particle's neighbor list.
+
+To port this algorithm to the GPU, we use the parallel algorithms library Thrust, distributed as part of the CUDA Toolkit. Thrust provides parallel sorting, searching, and prefix summing algorithms that are particularly useful in porting particle algorithms. To construct the neighbor list on the GPU, we follow the basic approach used by Canaba, a product of the Particle Co-Design Center:
+
+- Sort the particles on each grid by bin, using a parallel counting sort. We use Thrust's `exclusive\_scan` function to implement the prefix sum phase of the sort, and hand-coded kernels for the rest. This step does not actually involving rearranging the particle data - rather, we compute a permutation that would put the particles in order without actually reordering them.
+- Once the particles are sorted by bin, we can loop over the particles in neighboring bins. We make two passes over the particles. First, we launch a kernel to count the number of collision partners for each particle.
+- Then, we sum these numbers and allocate space for our neighbor list.
+- Finally, we make a another pass over the particles, putting them into to list at the appropriate place.
+
+Note that we build a \emph{full} neighbor list, meaning that if particle $i$ appears in particle $j$'s list, then particle $j$ also appears in particle $i$'s list. This simplifies the force-computation step when using these lists, since the forces and torques for a given particle can be updated without atomics.
+
+The final on-grid neighbor list data structure consists of two arrays. First, we have the neighbor list itself, stored as a big, 1D array of particle indices. Then, we have an `offsets` array that stores, for each particle, where in the neighbor list array to look. The details of this data structure have been hidden inside an iterator, so that user code can look like:
+
+.. code-block:: c
+
+ // now we loop over the neighbor list and compute the forces
+ AMREX_FOR_1D ( np, i,
+ {
+ ParticleType& p1 = pstruct[i];
+ p1.rdata(PIdx::ax) = 0.0;
+ p1.rdata(PIdx::ay) = 0.0;
+ p1.rdata(PIdx::az) = 0.0;
+
+ for (const auto& p2 : nbor_data.getNeighbors(i))
+ {
+ Real dx = p1.pos(0) - p2.pos(0);
+ Real dy = p1.pos(1) - p2.pos(1);
+ Real dz = p1.pos(2) - p2.pos(2);
+
+ ...
+ }
+
+Note that, because of our use of managed memory to store the particle data and the neighbor list, the above code will work when compiled for either CPU or GPU.
+
+The above algorithm deals with constructing a neighbor list for the particles on a single grid. When domain decomposition is used, one must also make copies of particles on adjacent grids, potentially performing the necessary MPI communication for grids associated with other processes. The routines `fillNeighbors`, which computes which particles needed to be ghosted to which grid, and `updateNeighbors`, which copies up-to-date data for particles that have already been ghosted, have also been offloaded to the GPU, using techniques similar to AMReX's `Redistribute` routine. The important thing for users is that calling these functions does not trigger copying data off the GPU.
+
+Once the neighbor list has been constructed, collisions with both particles and walls can easily be processed.
+
+We have created a GPU branch of MFIX that is capable of running with GPU support. As of this writing, the following operations in MFIX have been offloaded:
+
+- Neighbor particles / neighbor list construction
+- Particle-particle collisions
+- Particle-wall collisions
+- PIC Deposition (used in putting the drag force and solids volume fraction on the grid)
+
+
+
diff --git a/docs/source/Particles.rst b/docs/source/Particles_Chapter.rst
similarity index 98%
rename from docs/source/Particles.rst
rename to docs/source/Particles_Chapter.rst
index 0fa2e7a13bbdf5b81120b5df9da2cb148bcc54c4..2465434ac9c48086aed61eee7be5aa8bfe391ef2 100644
--- a/docs/source/Particles.rst
+++ b/docs/source/Particles_Chapter.rst
@@ -26,6 +26,4 @@ particle-particle, particle-fluid, and particle-wall interactions.
ParticleBasics
ParticleFluid
ParticleWalls
-
-
-
+ ParticlesOnGpus
diff --git a/docs/source/Slopes.rst b/docs/source/Slopes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..ea366a0534a6eca1a66bcd783f1fa70ed04ccbc2
--- /dev/null
+++ b/docs/source/Slopes.rst
@@ -0,0 +1,45 @@
+Computing slopes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Slopes are computed one direction at a time for each scalar and each
+velocity component.
+
+We use the second order Monotonized Central (MC)
+limiter (van Leer, 1977). The scheme is described below for the u-velocity.
+
+The limiter computes the slope at cell "i" by combining the left, central
+and right u-variation "du":
+
+.. code:: shell
+
+ du_l = u(i) - u(i-1) = left variation
+ du_c = 0.5 * ( u(i+1) - u(i-1) ) = central (umlimited) variation
+ du_r = u(i+1) - u(i) = right variation
+
+Finally, the u-variation at cell "i" is given by :
+
+.. code:: shell
+
+ du(i) = sign(du_c) * min(2|du_l|, |du_c|, 2|du_r|)) if du_l*du_r > 0
+ du(i) = 0 otherwise
+
+The above procedure is applied direction by direction.
+
+BOUNDARY CONDITIONS
+When periodic or Neumann's BCs are imposed, the scheme can be applied
+without any change since the ghost cells at the boundary are filled
+by either periodicity or by extrapolation.
+For Dirichlet's BCs in the transversal direction, the scheme can again
+be applied as is since the velocity is known at the first ghost cell
+out of the domain.
+However, for Dirichlet's BCs in the longitudinal direction, the velocity
+is not known outside the domain since the BC is applied directly at the first
+valid node which lies on the boundary itself. Therefore, the scheme must be
+arranged as follows to use ONLY values from inside the domain.
+For a left boundary (i=0), the u-variations are:
+
+.. code:: shell
+
+ du_l = 0 Dont use values on the left
+ du_c = -1.5*u(0) + 2*u(1) -0.5*u(2) 2nd order right-biased
+ du_r = u(1) - u(0) Right variation
diff --git a/docs/source/Structure.rst b/docs/source/Structure.rst
index c106641007e2d2976f492fc60162a42bc6d4cb7e..98195a68bf0d8f12804a606f48f73f9e8c0905e1 100644
--- a/docs/source/Structure.rst
+++ b/docs/source/Structure.rst
@@ -12,17 +12,7 @@ Directory overview
+---------------+--------------------------------------------------+
| exec | Directory for building with gmake (optional) |
+---------------+--------------------------------------------------+
-| exec_cc | Directory for building with gmake (optional) |
-+---------------+--------------------------------------------------+
-| src_des | Source files for particle-only operations |
-+---------------+--------------------------------------------------+
-| src_ebs | Source files for EB-only operations |
-+---------------+--------------------------------------------------+
-| src_staggered | Source files for SIMPLE and projection algorithm |
-| | with face-centered velocity components |
-+---------------+--------------------------------------------------+
-| src_cc | Source files for projection algorithms with |
-| | cell-centered velocity components |
+| src | Source files |
+---------------+--------------------------------------------------+
| tests | Regression tests (see tests/README.md) |
+---------------+--------------------------------------------------+
diff --git a/docs/source/_static/theme_overrides.css b/docs/source/_static/theme_overrides.css
index db92babc73e090e6fd1ee73203f6b74493522c83..7c1a5202231086af91e2bd4b388d8de4b064f27f 100644
--- a/docs/source/_static/theme_overrides.css
+++ b/docs/source/_static/theme_overrides.css
@@ -8,10 +8,3 @@
max-width: 100%;
overflow: visible;
}
-
-/* rtd_theme currently colours code-blocks in the pygments style (ie. green).
- * This overrides this design choice with white. Note: future versions of the
- * sphinx_rtd_theme might not need this HACK. */
-.highlight {
- background: #ffffff;
-}
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 69abf9d8865ac26cf0d68ce0395c2b39135f6015..f1c57b7933c948eaaaab69d6803945a8caf48834 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -22,9 +22,14 @@ the master branch at the beginning of each month.
Introduction
GettingStarted
- Inputs
- Fluids
- Particles
+ Inputs_Chapter
+ ManagingGridHierarchy_Chapter
+ Fluids_Chapter
+ Particles_Chapter
+ EB
+ CITests
+ NightlyTests
+ Debugging
Notice
------