Fix issues #250-264

.github/workflows/cuda.yml

@@ -129,6 +129,7 @@ jobs:
            -DINCFLO_MPI=ON                                           \
            -DINCFLO_OMP=OFF                                          \
            -DINCFLO_EB=OFF                                           \
+           -DINCFLO_PARTICLES=ON                                     \
            -DINCFLO_CUDA=ON                                          \
            -DAMReX_GPU_BACKEND=CUDA                                  \
            -DAMReX_CUDA_ARCH=70                                      \

.github/workflows/gcc.yml

@@ -302,6 +302,7 @@ jobs:
            -DINCFLO_MPI=OFF                                          \
            -DINCFLO_OMP=ON                                           \
            -DINCFLO_EB=OFF                                           \
+           -DINCFLO_PARTICLES=ON                                     \
            -DAMREX_HOME=${{ github.workspace }}/amrex                \
            -DAMREX_HYDRO_HOME=${{ github.workspace }}/AMReX-Hydro    \
            -DCMAKE_VERBOSE_MAKEFILE=ON                               \
@@ -321,4 +322,4 @@ jobs:
       run: |
         export OMP_NUM_THREADS=2
         cd ${{ github.workspace }}/incflo
-        build/incflo.ex test_no_eb_2d/benchmark.bouss_bubble_god max_step=10 incflo.verbose=1 mac_proj.verbose=1 nodal_proj.verbose=1
+        build/incflo.ex test_no_eb_2d/benchmark.bouss_bubble_god max_step=10 incflo.verbose=1 mac_proj.verbose=1 nodal_proj.verbose=1 incflo.use_tracer_particles=1 amr.refine_particles=0

Docs/sphinx_documentation/source/InputsAlgorithm.rst

@@ -22,7 +22,7 @@ The following inputs must be preceded by "incflo."
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  constant_density    |  Only evolve the continuity equation if false                         |  bool       |  true        |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
-|  rho_0               |  density (if constant)                                                |  Real       |  1.0         |
+|  ro_0                |  density (if constant)                                                |  Real       |  1.0         |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  diffusion_type      |  Diffusion type (0 = Explicit, 1 = Crank-Nicholson, 2 = Implicit)     |       int   | 2 (Implicit) |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+

Docs/sphinx_documentation/source/InputsMultigrid.rst

@@ -5,7 +5,10 @@ Multigrid Inputs
 
 Below is a list of the most commonly used multigrid settings options.
 To control the nodal projection precede with "nodal_proj", for the MAC projection use "mac_proj", and
-for the diffusion solver use "diffusion"
+for the diffusion solvers use "scalar_diffusion" (tracers, temperature and the component-wise
+velocity solve) or "tensor_diffusion" (the tensor velocity solve).  The diffusion solvers read
+the verbosity and iteration keys with an "mg_" prefix: mg_verbose, mg_bottom_verbose,
+mg_max_iter, mg_bottom_maxiter, mg_maxorder.
 
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
 |                         |  Description                                                          |   Type      | Default        |
@@ -14,18 +17,18 @@ for the diffusion solver use "diffusion"
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
 | bottom_verbose          |  Verbosity of BiCGStab solver                                         |    Int      |   0            |
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
-| rtol                    |  Relative tolerance                                                   |  | Real     |  | 1.e-11      |
+| mg_rtol                 |  Relative tolerance                                                   |  | Real     |  | 1.e-11      |
 |                         |                                                                       |  | float    |  | 1.e-4       |
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
-| atol                    |  Absolute tolerance                                                   |  | Real     |  | 1.e-14      |
+| mg_atol                 |  Absolute tolerance                                                   |  | Real     |  | 1.e-14      |
 |                         |                                                                       |  | float    |  | 1.e-7       |
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
-| maxiter                 |  Maximum number of iterations                                         |    Int      | nodal 100      |
+| maxiter                 |  Maximum number of iterations (diffusion: mg_max_iter)                |    Int      | nodal 100      |
 |                         |                                                                       |             | MAC   200      |
 |                         |                                                                       |             | diffusion 100  |
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
-| bottom_maxiter          |  Maximum number of iterations in the                                  |    Int      | nodal 100      |
-|                         |  bottom solver if using bicg, cg, bicgcg or cgbicg                    |             | MAC   200      |
+| bottom_maxiter          |  Maximum number of iterations in the bottom solver                    |    Int      | nodal 100      |
+|                         |  if using bicg, cg, bicgcg or cgbicg (diffusion: mg_bottom_maxiter)   |             | MAC   200      |
 |                         |                                                                       |             | diffusion 100  |
 +-------------------------+-----------------------------------------------------------------------+-------------+----------------+
 | mg_max_coarsening_level |  Maximum number of coarser levels to allow.                           |    Int      |   100          |

Docs/sphinx_documentation/source/InputsPlotFiles.rst

@@ -23,6 +23,7 @@ as whether the EB geometry should be written out.
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+
 | 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 |             |           |
+|                     | NOTE: this key takes the "incflo" prefix (incflo.write_eb_surface)    |             |           |
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+
 
 The following inputs must be preceded by "amr" and control what variables will be written in plotfiles.
@@ -66,7 +67,7 @@ and if using EB, volume fraction.
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+
 | plt_strainrate      | Save strain rate to plot file                                         |    Int      | 0         |
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+
-| plt_divu            | Save velocity divergence to plot file                                 |    Int      | 0         |
+| plt_divu            | Not implemented: requesting divu is refused when the inputs are read  |    Int      | 0         |
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+
 | plt_vfrac           | Save EB volume fraction to plot file                                  |    Int      | 1         |
 +---------------------+-----------------------------------------------------------------------+-------------+-----------+

Docs/sphinx_documentation/source/Particles.rst

@@ -106,10 +106,6 @@ distribution in a box (``initial_distribution_type = box``), which is implemente
 - If the cylinder inputs are present, particles lying outside the cylinder are removed immediately after
   initialization.
 
-Note that, as currently implemented, particles are placed only in cells whose indices satisfy
-:math:`i \bmod 2 = 0` and :math:`j \bmod 2 = 1`, so only one quarter of the cells in the box receive
-particles.
-
 Advection
 ---------
 
@@ -141,8 +137,8 @@ following mesh variables can be derived from the particles:
 
 - A particle count per cell (computed with :cpp:`incflo_PC::Increment()`), which is included in plotfiles
   as ``particle_count`` and is used in :cpp:`incflo::ErrorEst()` as a refinement criterion: when
-  ``incflo.refine_particles`` is true (the default), cells containing at least one particle are tagged
+  ``amr.refine_particles`` is true (the default), cells containing at least one particle are tagged
   for refinement.
 - A mass density computed by depositing the particle mass onto the mesh with linear interpolation in
-  :cpp:`incflo_PC::massDensity()`, which appears in plotfiles as the variable
-  ``tracer_particles_mass_density``.
+  :cpp:`incflo_PC::massDensity()`. This is not yet connected to the plotfile writer, so there is no
+  ``tracer_particles_mass_density`` plot variable at present.

convergence_2d/benchmark.taylor_vortex_decaying

@@ -5,13 +5,12 @@ stop_time               =   0.2         # Max (simulated) time to evolve
 max_step                =  -1           # Max number of time steps
 steady_state            =   0           # Steady-state solver? 
 
-incflo.use_godunov = true               #
+incflo.advection_type = "Godunov"  #
 incflo.use_ppm     = true               #
 incflo.use_mac_phi_in_godunov = false   #
 
 incflo.initial_iterations = 10
 
-incflo.v = 4
 incflo.verbose = 4
 
 amr.max_grid_size = 1024

convergence_2d/todo

@@ -4,9 +4,12 @@ rm -rf out* plt* chk*
 ../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64   incflo.use_mac_phi_in_godunov = false incflo.fixed_dt = 0.004 > out.64.false
 ../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 128 128  incflo.use_mac_phi_in_godunov = false incflo.fixed_dt = 0.002 > out.128.false
 
-../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 16  16   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.016 > out.16.true
-../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 32  32   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.008 > out.32.true
-../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.004 > out.64.true
-../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 128 128  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.002 > out.128.true
+# incflo.use_mac_phi_in_godunov is not implemented (it drops grad p from the MAC
+# face prediction but nothing ever reads mac_phi back), so the runs below are
+# disabled; the code aborts if the option is set to true.
+#../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 16  16   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.016 > out.16.true
+#../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 32  32   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.008 > out.32.true
+#../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64   incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.004 > out.64.true
+#../test_no_eb_2d/incflo2d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 128 128  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.002 > out.128.true
 
 source todo_print

convergence_3d/benchmark.taylor_vortex_decaying

@@ -5,13 +5,12 @@ stop_time               =   0.2         # Max (simulated) time to evolve
 max_step                =  -1           # Max number of time steps
 steady_state            =   0           # Steady-state solver? 
 
-incflo.use_godunov = true               #
+incflo.advection_type = "Godunov"  #
 incflo.use_ppm     = true               #
 incflo.use_mac_phi_in_godunov = false   #
 
 incflo.initial_iterations = 10
 
-incflo.v = 4
 incflo.verbose = 4
 
 amr.max_grid_size = 1024

convergence_3d/todo

@@ -4,9 +4,12 @@ rm -rf out* plt* chk*
 ../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64  64  incflo.use_mac_phi_in_godunov = false incflo.fixed_dt = 0.004 > out.64.false
 #../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 128 128 128  incflo.use_mac_phi_in_godunov = false incflo.fixed_dt = 0.002 > out.128.false
 
-../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 16  16 16  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.016 > out.16.true
-../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 32  32 32  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.008 > out.32.true
-../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64 64  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.004 > out.64.true
+# incflo.use_mac_phi_in_godunov is not implemented (it drops grad p from the MAC
+# face prediction but nothing ever reads mac_phi back), so the runs below are
+# disabled; the code aborts if the option is set to true.
+#../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 16  16 16  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.016 > out.16.true
+#../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 32  32 32  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.008 > out.32.true
+#../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 64  64 64  incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.004 > out.64.true
 #../test_no_eb_3d/incflo3d.gnu.MPI.ex benchmark.taylor_vortex_decaying amr.n_cell = 128 128 128 incflo.use_mac_phi_in_godunov = true incflo.fixed_dt = 0.002 > out.128.true
 
 source todo_print

src/boundary_conditions/CMakeLists.txt

@@ -4,5 +4,4 @@ target_sources(incflo
    incflo_fillpatch.cpp
    incflo_fillphysbc.cpp
    incflo_set_bcs.cpp
-   incflo_set_velocity_bcs.cpp
    )

src/boundary_conditions/Make.package

@@ -2,4 +2,3 @@
 CEXE_sources += boundary_conditions.cpp
 CEXE_sources += incflo_fillpatch.cpp incflo_fillphysbc.cpp
 CEXE_sources += incflo_set_bcs.cpp
-CEXE_sources += incflo_set_velocity_bcs.cpp

src/boundary_conditions/boundary_conditions.cpp

@@ -203,6 +203,12 @@ void incflo::init_bcs ()
 #ifdef AMREX_USE_EB
             // ReadParameters() already called
             if (m_advection_type != "Godunov") { amrex::Abort("mixed BCs require Godunov"); }
+            // MLTensorOp/MLEBTensorOp have no Robin BC, so the mixed BC is only
+            // supported by the component-wise velocity solve.  Catch this here rather
+            // than at the first diffusion solve of the first time step.
+            if (use_tensor_solve) {
+                amrex::Abort("mixed BCs require incflo.use_tensor_solve = false");
+            }
 
             ParmParse ipp("incflo");
             std::string eb_geom = "null";

src/diffusion/incflo_diffusion.cpp

@@ -178,6 +178,14 @@ incflo::get_diffuse_tensor_bc (Orientation::Side side) const noexcept
                 r[dir][dir] = LinOpBCType::Dirichlet;
                 break;
             }
+            case BC::mixed:
+            {
+                // Neither MLTensorOp nor MLEBTensorOp implements Robin BCs, so the
+                // mixed BC can only be done by the component-wise velocity solve.
+                amrex::Abort("get_diffuse_tensor_bc: mixed BCs are not supported by the "
+                             "tensor solve; set incflo.use_tensor_solve = false");
+                break;
+            }
             default:
                 amrex::Abort("get_diffuse_tensor_bc: undefined BC type");
             };
@@ -236,7 +244,7 @@ incflo::get_diffuse_velocity_bc (Orientation::Side side, int comp) const noexcep
                 break;
             }
             default:
-                amrex::Abort("get_diffuse_tensor_bc: undefined BC type");
+                amrex::Abort("get_diffuse_velocity_bc: undefined BC type");
             };
         }
     }

src/embedded_boundaries/eb_annulus.cpp

@@ -78,7 +78,10 @@ void incflo::make_eb_annulus()
     auto gshop = EB2::makeShop(annulus);
 
     // Build index space
-    int max_level_here = 0;
+    // geom.back() is the finest AMR level; requiring max_level coarsenings
+    // makes AMReX build EB data (including domain ghost cells) for every
+    // AMR level, as documented for EB2::Build.
+    int max_level_here = max_level;
     int max_coarsening_level = 100;
     EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level);
 }

src/embedded_boundaries/eb_box.cpp

@@ -95,7 +95,10 @@ void incflo::make_eb_box()
         auto gshop = EB2::makeShop(my_box);
 
         // Build index space
-        int max_level_here = 0;
+        // geom.back() is the finest AMR level; requiring max_level coarsenings
+        // makes AMReX build EB data (including domain ghost cells) for every
+        // AMR level, as documented for EB2::Build.
+        int max_level_here = max_level;
         int max_coarsening_level = 100;
         EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level);
     }

src/embedded_boundaries/eb_chkptfile.cpp

@@ -7,7 +7,10 @@ using namespace amrex;
 void incflo::make_eb_chkptfile()
 {
    // Build index space
-   int max_level_here = 0;
+   // geom.back() is the finest AMR level; requiring max_level coarsenings
+   // makes AMReX build EB data (including domain ghost cells) for every
+   // AMR level, as documented for EB2::Build.
+   int max_level_here = max_level;
    int max_coarsening_level = 100;
    EB2::BuildFromChkptFile("geom_chk", geom.back(), max_level_here, max_level_here + max_coarsening_level);
 }

src/embedded_boundaries/eb_csg.cpp

@@ -45,7 +45,10 @@ void incflo::make_eb_csg(const std::string& geom_file)
     auto gshop = EB2::makeShop(final_csg_if);
 
     // Build index space
-    int max_level_here = 0;
+    // geom.back() is the finest AMR level; requiring max_level coarsenings
+    // makes AMReX build EB data (including domain ghost cells) for every
+    // AMR level, as documented for EB2::Build.
+    int max_level_here = max_level;
     int max_coarsening_level = 100;
     EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level);
 }