From 304eb2a1467330014b435ca9c415226a4236985d Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 18 May 2026 09:38:00 -0700 Subject: [PATCH 1/4] more precision changes --- src/embedded_boundaries/eb_annulus.cpp | 2 +- src/embedded_boundaries/eb_box.cpp | 18 +++++++++++------- src/embedded_boundaries/eb_cylinder.cpp | 2 +- src/embedded_boundaries/eb_sphere.cpp | 2 +- src/setup/init.cpp | 2 +- src/utilities/io.cpp | 4 ++-- 6 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/embedded_boundaries/eb_annulus.cpp b/src/embedded_boundaries/eb_annulus.cpp index f8c96a83a..b7eed7152 100644 --- a/src/embedded_boundaries/eb_annulus.cpp +++ b/src/embedded_boundaries/eb_annulus.cpp @@ -39,7 +39,7 @@ void incflo::make_eb_annulus() // Compute distance between cylinder centres Real offset = 0.0; for(int i = 0; i < AMREX_SPACEDIM; i++) - offset += std::pow(outer_center[i] - inner_center[i], 2); + offset += amrex::Math::powi<2>(outer_center[i] - inner_center[i]); offset = std::sqrt(offset); // Check that the inner cylinder is fully contained in the outer one diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 3c361e1d8..9a372ae09 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -32,7 +32,11 @@ void incflo::make_eb_box() ************************************************************************/ Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); - Real offset = 1.0e-15; +#ifdef AMREX_USE_FLOAT + Real offset = Real(1e-8); +#else + Real offset = Real(1e-15); +#endif bool inside = true; for(int i = 0; i < AMREX_SPACEDIM; i++) @@ -54,8 +58,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(0)) { - xlo = 2.0 * geom[0].ProbLo(0) - geom[0].ProbHi(0); - xhi = 2.0 * geom[0].ProbHi(0) - geom[0].ProbLo(0); + xlo = Real(2) * geom[0].ProbLo(0) - geom[0].ProbHi(0); + xhi = Real(2) * geom[0].ProbHi(0) - geom[0].ProbLo(0); } Real ylo = boxLo[1] + offset; @@ -65,8 +69,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(1)) { - ylo = 2.0 * geom[0].ProbLo(1) - geom[0].ProbHi(1); - yhi = 2.0 * geom[0].ProbHi(1) - geom[0].ProbLo(1); + ylo = Real(2) * geom[0].ProbLo(1) - geom[0].ProbHi(1); + yhi = Real(2) * geom[0].ProbHi(1) - geom[0].ProbLo(1); } #if (AMREX_SPACEDIM > 2) @@ -77,8 +81,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(2)) { - zlo = 2.0 * geom[0].ProbLo(2) - geom[0].ProbHi(2); - zhi = 2.0 * geom[0].ProbHi(2) - geom[0].ProbLo(2); + zlo = Real(2) * geom[0].ProbLo(2) - geom[0].ProbHi(2); + zhi = Real(2) * geom[0].ProbHi(2) - geom[0].ProbLo(2); } #endif diff --git a/src/embedded_boundaries/eb_cylinder.cpp b/src/embedded_boundaries/eb_cylinder.cpp index 02afa9815..cc5ba9c03 100644 --- a/src/embedded_boundaries/eb_cylinder.cpp +++ b/src/embedded_boundaries/eb_cylinder.cpp @@ -33,7 +33,7 @@ void incflo::make_eb_cylinder() pp.getarr("center", centervec, 0, 3); Array center = {AMREX_D_DECL(centervec[0], centervec[1], centervec[2])}; - rotation = (rotation/180.)*M_PI; + rotation = (rotation/180.)*Real(M_PI); // Print info about cylinder amrex::Print() << " " << "\n"; diff --git a/src/embedded_boundaries/eb_sphere.cpp b/src/embedded_boundaries/eb_sphere.cpp index cdf735e06..09ef3563b 100644 --- a/src/embedded_boundaries/eb_sphere.cpp +++ b/src/embedded_boundaries/eb_sphere.cpp @@ -16,7 +16,7 @@ void incflo::make_eb_sphere() { // Initialise sphere parameters bool inside = true; - Real radius = 0.0002; + Real radius = Real(0.0002); Vector centervec(3); // Get sphere information from inputs file. * diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 2a3db0709..db7259667 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -222,7 +222,7 @@ void incflo::ReadParameters () amrex::Real tol_deg(0.); pp_eb_flow.query("normal_tol", tol_deg); - m_eb_flow.normal_tol = tol_deg*M_PI/amrex::Real(180.); + m_eb_flow.normal_tol = tol_deg*amrex::Real(M_PI)/amrex::Real(180); } if (m_advect_tracer && m_eb_flow.enabled && m_eb_flow.tracer.empty()) { diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 801011e10..4f3096bb9 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -175,7 +175,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_lo[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_lo[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } @@ -186,7 +186,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_hi[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_hi[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } From 71cf392c18e7ff13aa7e2e97b9c647cd8626f497 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 21 Sep 2026 15:35:54 -0700 Subject: [PATCH 2/4] Fix issues #250-264 Correctness and physics: - #250: delete the dead set_inflow_velocity / prob_set_inflow_velocity, which have had no caller since PR #145 and wrote velocity component 2 with no AMREX_SPACEDIM guard for probtypes 33/333. - #251: refuse incflo.use_mac_phi_in_godunov, which only drops grad p from the MAC face prediction and is never read back as a pressure, and zero mac_phi in the LevelData constructor so the "macphi"/"error_mac_p" plot variables are defined on regridded and restarted levels. - #252: use_tensor_correction now forces godunov_include_diff_in_forcing off, because divtau_o then holds the (tensor - scalar) difference rather than the full explicit viscous term the edge-state forcing expects. - #253: InitialPressureProjection passes set_inflow_bc = false; the field being projected is the body force (rho-rho0)/rho*g, not a velocity. - #254: ApplyCCProjection uses its own face scratch instead of the caller's u_mac/v_mac/w_mac, which the caller still needs for particle advection. - #256: refuse probtypes 1100/1102 in a 2D build and guard the z-direction code in init_jump, which otherwise wrote past the end of a 2-component FAB. Defensive gaps and diagnostics: - #255: reject a mixed BC combined with use_tensor_solve at read time (neither MLTensorOp nor MLEBTensorOp implements Robin BCs), add the BC::mixed case to get_diffuse_tensor_bc, and fix get_diffuse_velocity_bc's abort message. - #257: every EB builder now requires max_level coarsenings, as the STL builder already did, so coarse AMR levels get EB data in their domain ghost cells. - #258: an unrecognised incflo.geometry aborts instead of silently building a regular geometry; only an empty name or all_regular selects make_eb_regular. - #259: refuse a restart whose checkpoint prob_lo/prob_hi differ from the inputs the EB was built from. - #264: refuse the unimplemented "divu" plot variable when the inputs are read rather than at the first plotfile, and refuse a non-zero incflo.ic_p, which only sets m_p000 and would otherwise be silently dropped. Build coverage, docs and shipped decks: - #260: build particle support in one GCC and one CUDA CI job, run the 2D smoke test with tracer particles, and add USE_PARTICLES to the four GNUmakefiles. - #261: drop the stale quarter-cell seeding note removed by PR #213 and the tracer_particles_mass_density plot variable that nothing can produce. - #262: document the input keys the code actually queries (ro_0, incflo.write_eb_surface, amr.refine_particles, mg_rtol/mg_atol, the scalar_diffusion/tensor_diffusion prefixes and the mg_ iteration keys). - #263: rename incflo.v, amr.plot_p and incflo.use_godunov in the shipped decks, drop the unread cylinder.height, and mark cylinder.speed (and the two decks it was the only flow driver for) as not implemented. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/cuda.yml | 1 + .github/workflows/gcc.yml | 3 +- .../source/InputsAlgorithm.rst | 2 +- .../source/InputsMultigrid.rst | 15 +-- .../source/InputsPlotFiles.rst | 3 +- .../sphinx_documentation/source/Particles.rst | 10 +- .../benchmark.taylor_vortex_decaying | 3 +- convergence_2d/todo | 11 ++- .../benchmark.taylor_vortex_decaying | 3 +- convergence_3d/todo | 9 +- src/boundary_conditions/CMakeLists.txt | 1 - src/boundary_conditions/Make.package | 1 - .../boundary_conditions.cpp | 6 ++ .../incflo_set_velocity_bcs.cpp | 39 -------- src/diffusion/incflo_diffusion.cpp | 10 +- src/embedded_boundaries/eb_annulus.cpp | 5 +- src/embedded_boundaries/eb_box.cpp | 5 +- src/embedded_boundaries/eb_chkptfile.cpp | 5 +- src/embedded_boundaries/eb_csg.cpp | 5 +- src/embedded_boundaries/eb_cyl_tuscan.cpp | 5 +- src/embedded_boundaries/eb_cylinder.cpp | 5 +- src/embedded_boundaries/eb_regular.cpp | 5 +- src/embedded_boundaries/eb_sphere.cpp | 5 +- src/embedded_boundaries/eb_spherecube.cpp | 5 +- src/embedded_boundaries/eb_tuscan.cpp | 5 +- src/embedded_boundaries/eb_twocylinders.cpp | 5 +- .../embedded_boundaries.cpp | 22 ++++- src/incflo.H | 5 - src/prob/prob_bc.cpp | 98 ------------------- src/prob/prob_init_fluid.cpp | 17 +++- src/projection/incflo_apply_cc_projection.cpp | 22 ++++- src/setup/incflo_arrays.cpp | 5 + src/setup/init.cpp | 51 +++++++++- src/utilities/io.cpp | 25 ++++- test_2d/GNUmakefile | 4 + test_2d/benchmark.inviscid_rotated_planes_x | 1 - test_2d/benchmark.rotated_planes_x | 1 - test_2d/benchmark.rotated_planes_y | 1 - test_2d/benchmark.tracer_adv_diff | 1 - test_2d/benchmark.tracer_adv_diff_cn | 1 - test_2d/benchmark.tracer_advection | 1 - test_2d/benchmark.vortex_in_sphere | 2 +- test_3d/GNUmakefile | 4 + test_3d/benchmark.rotated_cylinder_x | 1 - test_3d/benchmark.rotated_cylinder_y | 1 - test_3d/benchmark.rotated_cylinder_z | 1 - test_3d/benchmark.slanted_cylinder | 1 - test_3d/benchmark.tracer_adv_diff | 1 - test_3d/benchmark.tracer_adv_diff_cn | 1 - test_3d/benchmark.tracer_advection | 1 - test_3d/inputs.channel_sphere | 2 +- test_3d/inputs.couette_cylinder | 4 +- test_3d/inputs.rotating_cylinder | 9 +- test_3d/inputs.taylor_couette | 9 +- .../benchmark.inviscid_rotated_cylinder_x | 1 - .../benchmark.inviscid_rotated_cylinder_y | 1 - .../benchmark.inviscid_rotated_cylinder_z | 1 - .../benchmark.rotated_cylinder_x | 1 - .../benchmark.rotated_cylinder_y | 1 - .../benchmark.rotated_cylinder_z | 1 - test_no_eb_2d/GNUmakefile | 4 + test_no_eb_3d/GNUmakefile | 4 + 62 files changed, 261 insertions(+), 217 deletions(-) delete mode 100644 src/boundary_conditions/incflo_set_velocity_bcs.cpp diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 337f3de6a..063191ebb 100644 --- a/.github/workflows/cuda.yml +++ b/.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 \ diff --git a/.github/workflows/gcc.yml b/.github/workflows/gcc.yml index 57f48f9b9..846905e21 100644 --- a/.github/workflows/gcc.yml +++ b/.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 diff --git a/Docs/sphinx_documentation/source/InputsAlgorithm.rst b/Docs/sphinx_documentation/source/InputsAlgorithm.rst index ce7814f8a..0bf17516f 100644 --- a/Docs/sphinx_documentation/source/InputsAlgorithm.rst +++ b/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) | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ diff --git a/Docs/sphinx_documentation/source/InputsMultigrid.rst b/Docs/sphinx_documentation/source/InputsMultigrid.rst index 07bfcc6ca..eb1f97724 100644 --- a/Docs/sphinx_documentation/source/InputsMultigrid.rst +++ b/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 | diff --git a/Docs/sphinx_documentation/source/InputsPlotFiles.rst b/Docs/sphinx_documentation/source/InputsPlotFiles.rst index 42e0d1779..c2b0c8682 100644 --- a/Docs/sphinx_documentation/source/InputsPlotFiles.rst +++ b/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 | +---------------------+-----------------------------------------------------------------------+-------------+-----------+ diff --git a/Docs/sphinx_documentation/source/Particles.rst b/Docs/sphinx_documentation/source/Particles.rst index 3b96b9ed8..efed0f514 100644 --- a/Docs/sphinx_documentation/source/Particles.rst +++ b/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. diff --git a/convergence_2d/benchmark.taylor_vortex_decaying b/convergence_2d/benchmark.taylor_vortex_decaying index 589a97869..335dd8bd2 100644 --- a/convergence_2d/benchmark.taylor_vortex_decaying +++ b/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 diff --git a/convergence_2d/todo b/convergence_2d/todo index 8e93c64c2..d1d85bdcf 100644 --- a/convergence_2d/todo +++ b/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 diff --git a/convergence_3d/benchmark.taylor_vortex_decaying b/convergence_3d/benchmark.taylor_vortex_decaying index c7c6ccbd3..ab3e60a6e 100644 --- a/convergence_3d/benchmark.taylor_vortex_decaying +++ b/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 diff --git a/convergence_3d/todo b/convergence_3d/todo index 9539e9df3..a296a96bd 100644 --- a/convergence_3d/todo +++ b/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 diff --git a/src/boundary_conditions/CMakeLists.txt b/src/boundary_conditions/CMakeLists.txt index 7bf0c1725..ce9bd6251 100644 --- a/src/boundary_conditions/CMakeLists.txt +++ b/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 ) diff --git a/src/boundary_conditions/Make.package b/src/boundary_conditions/Make.package index 8e382a0ff..9de2ae895 100644 --- a/src/boundary_conditions/Make.package +++ b/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 diff --git a/src/boundary_conditions/boundary_conditions.cpp b/src/boundary_conditions/boundary_conditions.cpp index f4afacb76..2939884b7 100644 --- a/src/boundary_conditions/boundary_conditions.cpp +++ b/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"; diff --git a/src/boundary_conditions/incflo_set_velocity_bcs.cpp b/src/boundary_conditions/incflo_set_velocity_bcs.cpp deleted file mode 100644 index 31cb74bd5..000000000 --- a/src/boundary_conditions/incflo_set_velocity_bcs.cpp +++ /dev/null @@ -1,39 +0,0 @@ -#ifdef AMREX_USE_EB -#include -#endif - -#include - -using namespace amrex; - -void -incflo::set_inflow_velocity (int lev, amrex::Real time, MultiFab& vel, int nghost) -{ - Geometry const& gm = Geom(lev); - Box const& domain = gm.growPeriodicDomain(nghost); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - Orientation olo(dir,Orientation::low); - Orientation ohi(dir,Orientation::high); - if (m_bc_type[olo] == BC::mass_inflow || m_bc_type[ohi] == BC::mass_inflow) { - Box dlo = (m_bc_type[olo] == BC::mass_inflow) ? amrex::adjCellLo(domain,dir,nghost) : Box(); - Box dhi = (m_bc_type[ohi] == BC::mass_inflow) ? amrex::adjCellHi(domain,dir,nghost) : Box(); -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(vel); mfi.isValid(); ++mfi) { - Box const& gbx = amrex::grow(mfi.validbox(),nghost); - Box blo = gbx & dlo; - Box bhi = gbx & dhi; - Array4 const& v = vel[mfi].array(); - int gid = mfi.index(); - if (blo.ok()) { - prob_set_inflow_velocity(gid, olo, blo, v, lev, time); - } - if (bhi.ok()) { - prob_set_inflow_velocity(gid, ohi, bhi, v, lev, time); - } - } - } - } - vel.EnforcePeriodicity(gm.periodicity()); -} diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp index d4d5cb5e4..3e049a1d9 100644 --- a/src/diffusion/incflo_diffusion.cpp +++ b/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"); }; } } diff --git a/src/embedded_boundaries/eb_annulus.cpp b/src/embedded_boundaries/eb_annulus.cpp index b7eed7152..7f0530537 100644 --- a/src/embedded_boundaries/eb_annulus.cpp +++ b/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); } diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 9a372ae09..f9db545a4 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/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); } diff --git a/src/embedded_boundaries/eb_chkptfile.cpp b/src/embedded_boundaries/eb_chkptfile.cpp index e7608930d..1176a71ad 100644 --- a/src/embedded_boundaries/eb_chkptfile.cpp +++ b/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); } diff --git a/src/embedded_boundaries/eb_csg.cpp b/src/embedded_boundaries/eb_csg.cpp index 77c7d46c7..b3404f62f 100644 --- a/src/embedded_boundaries/eb_csg.cpp +++ b/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); } diff --git a/src/embedded_boundaries/eb_cyl_tuscan.cpp b/src/embedded_boundaries/eb_cyl_tuscan.cpp index ff9103b18..685fd5c3e 100644 --- a/src/embedded_boundaries/eb_cyl_tuscan.cpp +++ b/src/embedded_boundaries/eb_cyl_tuscan.cpp @@ -76,7 +76,10 @@ void incflo::make_eb_cyl_tuscan() auto gshop = EB2::makeShop(twocylinders); // 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); } diff --git a/src/embedded_boundaries/eb_cylinder.cpp b/src/embedded_boundaries/eb_cylinder.cpp index 0c1eff0ef..43613b4c5 100644 --- a/src/embedded_boundaries/eb_cylinder.cpp +++ b/src/embedded_boundaries/eb_cylinder.cpp @@ -58,7 +58,10 @@ void incflo::make_eb_cylinder() auto gshop = EB2::makeShop(my_cyl_rot); // 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); } diff --git a/src/embedded_boundaries/eb_regular.cpp b/src/embedded_boundaries/eb_regular.cpp index a2efc82d1..3c3d63a0f 100644 --- a/src/embedded_boundaries/eb_regular.cpp +++ b/src/embedded_boundaries/eb_regular.cpp @@ -10,5 +10,8 @@ void incflo::make_eb_regular() { EB2::AllRegularIF my_regular; auto gshop = EB2::makeShop(my_regular); - EB2::Build(gshop, geom.back(), 0, 100); + // 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. + EB2::Build(gshop, geom.back(), max_level, max_level + 100); } diff --git a/src/embedded_boundaries/eb_sphere.cpp b/src/embedded_boundaries/eb_sphere.cpp index 09ef3563b..2db1cddc8 100644 --- a/src/embedded_boundaries/eb_sphere.cpp +++ b/src/embedded_boundaries/eb_sphere.cpp @@ -44,7 +44,10 @@ void incflo::make_eb_sphere() auto gshop = EB2::makeShop(my_sphere); // 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); } diff --git a/src/embedded_boundaries/eb_spherecube.cpp b/src/embedded_boundaries/eb_spherecube.cpp index 502d0dc41..ca5924113 100644 --- a/src/embedded_boundaries/eb_spherecube.cpp +++ b/src/embedded_boundaries/eb_spherecube.cpp @@ -26,7 +26,10 @@ void incflo::make_eb_spherecube() auto gshop = EB2::makeShop(cubesphere); // 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); } diff --git a/src/embedded_boundaries/eb_tuscan.cpp b/src/embedded_boundaries/eb_tuscan.cpp index 5b7cc0d25..a5693c66d 100644 --- a/src/embedded_boundaries/eb_tuscan.cpp +++ b/src/embedded_boundaries/eb_tuscan.cpp @@ -118,7 +118,10 @@ void incflo::make_eb_tuscan() auto gshop = EB2::makeShop(boxes_inside); // 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); } diff --git a/src/embedded_boundaries/eb_twocylinders.cpp b/src/embedded_boundaries/eb_twocylinders.cpp index 0665d4c91..b6299598e 100644 --- a/src/embedded_boundaries/eb_twocylinders.cpp +++ b/src/embedded_boundaries/eb_twocylinders.cpp @@ -62,7 +62,10 @@ void incflo::make_eb_twocylinders() EB2::CylinderIF cyl1(radius1, direction1, center1, false); EB2::CylinderIF cyl2(radius2, direction2, center2, false); // 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; // NOTE: this must not be written as a ternary -- a conditional expression has diff --git a/src/embedded_boundaries/embedded_boundaries.cpp b/src/embedded_boundaries/embedded_boundaries.cpp index 7731abc56..05ff22681 100644 --- a/src/embedded_boundaries/embedded_boundaries.cpp +++ b/src/embedded_boundaries/embedded_boundaries.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -86,12 +87,31 @@ void incflo::MakeEBGeometry() make_eb_csg(csg_file); } #endif - else + else if (geom_type.empty() || geom_type == "all_regular") { amrex::Print() << "\n No EB geometry declared in inputs => " << " Will build all regular geometry." << "\n"; make_eb_regular(); } + else + { + // An unrecognised name (a typo, a 3D-only geometry in a 2D build, or "csg" + // in a build without CSG support) must not silently become a regular + // geometry: the user asked for an embedded boundary and would get a plain + // box with no diagnostic. + std::string msg = "incflo.geometry = " + geom_type + " is not a known EB geometry"; +#ifndef CSG_EB + if (geom_type == "csg") { + msg += " in this build (csg requires a build with CSG support)"; + } +#endif +#if (AMREX_SPACEDIM == 2) + if (geom_type == "twocylinders" || geom_type == "spherecube" || geom_type == "tuscan") { + msg += " in this build (it is only available in 3D)"; + } +#endif + amrex::Abort(msg); + } amrex::Print() << "Done making the EB geometry index space.\n" << "\n"; if (m_write_geom_chk) { diff --git a/src/incflo.H b/src/incflo.H index 72e84d46b..9484469f4 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -125,8 +125,6 @@ public: amrex::iMultiFab make_nodalBC_mask (int lev); amrex::Vector make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr); - void set_inflow_velocity (int lev, amrex::Real time, amrex::MultiFab& vel, int nghost); - #ifdef AMREX_USE_EB void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost); void set_eb_density (int lev, amrex::Real time, amrex::MultiFab& eb_density, int nghost); @@ -289,9 +287,6 @@ public: amrex::Array4 const& bcval, int lev); - void prob_set_inflow_velocity (int grid_id, amrex::Orientation ori, amrex::Box const& bx, - amrex::Array4 const& v, int lev, amrex::Real time); - #include "incflo_prob_I.H" #include "incflo_prob_usr_I.H" diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp index a82184022..e28d37c11 100644 --- a/src/prob/prob_bc.cpp +++ b/src/prob/prob_bc.cpp @@ -311,101 +311,3 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx, +std::to_string(m_probtype)); } } - -void incflo::prob_set_inflow_velocity (int /*grid_id*/, Orientation ori, Box const& bx, - Array4 const& vel, int lev, Real /*time*/) -{ - if (6 == m_probtype) - { - AMREX_D_TERM(Real u = m_ic_u;, - Real v = m_ic_v;, - Real w = m_ic_w;); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - AMREX_D_TERM(vel(i,j,k,0) = u;, - vel(i,j,k,1) = v;, - vel(i,j,k,2) = w;); - }); - } - else if (31 == m_probtype) - { - Real dyinv = Real(1) / Geom(lev).Domain().length(1); - Real u = m_ic_u; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real y = (j+Real(0.5))*dyinv; - vel(i,j,k,0) = Real(6) * u * y * (Real(1)-y); - }); - } -#if (AMREX_SPACEDIM == 3) - else if (311 == m_probtype) - { - Real dzinv = Real(1) / Geom(lev).Domain().length(2); - Real u = m_ic_u; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = (k+Real(0.5))*dzinv; - vel(i,j,k,0) = Real(6) * u * z * (Real(1)-z); - }); - } - else if (41 == m_probtype) - { - Real dzinv = Real(1) / Geom(lev).Domain().length(2); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = (k+Real(0.5))*dzinv; - vel(i,j,k,0) = Real(0.5)*z; - }); - } - else if (32 == m_probtype) - { - Real dzinv = Real(1) / Geom(lev).Domain().length(2); - Real v = m_ic_v; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real z = (k+Real(0.5))*dzinv; - vel(i,j,k,1) = Real(6) * v * z * (Real(1)-z); - }); - } -#endif - else if (322 == m_probtype) - { - Real dxinv = Real(1) / Geom(lev).Domain().length(0); - Real v = m_ic_v; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = (i+Real(0.5))*dxinv; - vel(i,j,k,1) = Real(6) * v * x * (Real(1)-x); - }); - } - else if (33 == m_probtype) - { - Real dxinv = Real(1) / Geom(lev).Domain().length(0); - Real w = m_ic_w; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = (i+Real(0.5))*dxinv; - vel(i,j,k,2) = Real(6) * w * x * (Real(1)-x); - }); - } - else if (333 == m_probtype) - { - Real dyinv = Real(1) / Geom(lev).Domain().length(1); - Real w = m_ic_w; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real y = (j+Real(0.5))*dyinv; - vel(i,j,k,2) = Real(6) * w * y * (Real(1)-y); - }); - } - else - { - const int dir = ori.coordDir(); - const Real bcv = m_bc_velocity[ori][dir]; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - vel(i,j,k,dir) = bcv; - }); - }; -} diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 5d1e79ca4..73f631ae1 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -97,6 +97,12 @@ void incflo::prob_init_fluid (int lev) } else if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype) { +#if (AMREX_SPACEDIM == 2) + if (1100 == m_probtype || 1102 == m_probtype) { + amrex::Abort("prob_init_fluid: probtypes 1100 and 1102 involve the z direction " + "and need a 3D build"); + } +#endif init_jump(vbx, gbx, ld.velocity.array(mfi), ld.density.array(mfi), @@ -717,32 +723,41 @@ void incflo::init_jump (Box const& vbx, Box const& /*gbx*/, GpuArray const& /*problo*/, GpuArray const& /*probhi*/) const { + // Only 1101 (jump along y, flipping u) is meaningful in 2D; 1100 and 1102 flip w + // / index the z direction and are refused for a 2D build in prob_init_fluid. int direction = 0; if (1101 == m_probtype) { direction = 1; } +#if (AMREX_SPACEDIM == 3) else if (1102 == m_probtype) { direction = 2; } +#endif int half_num_cells = domain.length(direction) / 2; ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { +#if (AMREX_SPACEDIM == 3) if (direction == 0) { if (i <= half_num_cells) { vel(i,j,k,2) = -vel(i,j,k,2); } } - else if (direction == 1) { + else +#endif + if (direction == 1) { if (j <= half_num_cells) { vel(i,j,k,0) = -vel(i,j,k,0); } } +#if (AMREX_SPACEDIM == 3) else if (direction == 2) { if (k <= half_num_cells) { vel(i,j,k,1) = -vel(i,j,k,1); } } +#endif }); } diff --git a/src/projection/incflo_apply_cc_projection.cpp b/src/projection/incflo_apply_cc_projection.cpp index 5db553a39..46d09a8b7 100644 --- a/src/projection/incflo_apply_cc_projection.cpp +++ b/src/projection/incflo_apply_cc_projection.cpp @@ -351,12 +351,26 @@ void incflo::ApplyCCProjection (Vector density, cc_gphi[lev].setVal(0.); } + // Scratch face velocities for the projection. The u_mac/v_mac/w_mac passed in + // are the caller's MAC-projected half-time velocities, which the caller still + // needs after this routine (tracer particle advection, small-cell correction), + // so they must not be overwritten here. + Vector > umac_proj(finest_level+1); Vector > mac_vec(finest_level+1); for (int lev=0; lev <= finest_level; ++lev) { - AMREX_D_TERM(mac_vec[lev][0] = u_mac[lev];, - mac_vec[lev][1] = v_mac[lev];, - mac_vec[lev][2] = w_mac[lev];); + AMREX_D_TERM(umac_proj[lev][0].define(u_mac[lev]->boxArray(), dmap[lev], 1, + u_mac[lev]->nGrow(), MFInfo(), Factory(lev));, + umac_proj[lev][1].define(v_mac[lev]->boxArray(), dmap[lev], 1, + v_mac[lev]->nGrow(), MFInfo(), Factory(lev));, + umac_proj[lev][2].define(w_mac[lev]->boxArray(), dmap[lev], 1, + w_mac[lev]->nGrow(), MFInfo(), Factory(lev));); + AMREX_D_TERM(umac_proj[lev][0].setVal(0.);, + umac_proj[lev][1].setVal(0.);, + umac_proj[lev][2].setVal(0.);); + AMREX_D_TERM(mac_vec[lev][0] = &umac_proj[lev][0];, + mac_vec[lev][1] = &umac_proj[lev][1];, + mac_vec[lev][2] = &umac_proj[lev][2];); } // Compute velocity on faces @@ -367,7 +381,7 @@ void incflo::ApplyCCProjection (Vector density, vel[lev]->FillBoundary(geom[lev].periodicity()); #if 1 MOL::ExtrapVelToFaces(*vel[lev], - AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]), + AMREX_D_DECL(*mac_vec[lev][0], *mac_vec[lev][1], *mac_vec[lev][2]), geom[lev], get_velocity_bcrec(), get_velocity_bcrec_device_ptr()); diff --git a/src/setup/incflo_arrays.cpp b/src/setup/incflo_arrays.cpp index 77af5aaeb..5ecf7f53b 100644 --- a/src/setup/incflo_arrays.cpp +++ b/src/setup/incflo_arrays.cpp @@ -23,6 +23,11 @@ incflo::LevelData::LevelData (amrex::BoxArray const& ba, conv_density_o (ba, dm, 1 , 0, MFInfo(), fact), conv_tracer_o (ba, dm, my_incflo->m_ntrac, 0, MFInfo(), fact) { + // mac_phi is only written by the MAC projection when use_mac_phi_in_godunov is + // on, but it can be plotted ("macphi", "error_mac_p") on any level, including + // levels created by regrid or restart, so it must not hold uninitialized memory. + mac_phi.setVal(0.0); + if (my_incflo->m_use_cc_proj) { p_cc.define(ba , dm, 1, 1, MFInfo(), fact); } else { diff --git a/src/setup/init.cpp b/src/setup/init.cpp index b379805ab..60453907b 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -4,6 +4,8 @@ #include #endif +#include + using namespace amrex; void incflo::ReadParameters () @@ -83,6 +85,14 @@ void incflo::ReadParameters () pp.query("godunov_use_forces_in_trans" , m_godunov_use_forces_in_trans); pp.query("godunov_include_diff_in_forcing" , m_godunov_include_diff_in_forcing); pp.query("use_mac_phi_in_godunov" , m_use_mac_phi_in_godunov); + if (m_use_mac_phi_in_godunov) { + // Only half of this option exists: it drops grad p from the forcing used + // to predict the MAC velocities, and mac_phi is solved for and doubled + // (see compute_MAC_projected_velocities) but never read back as the + // pressure anywhere in the Godunov forcing. + amrex::Abort("incflo.use_mac_phi_in_godunov is not implemented: mac_phi is " + "never used as the pressure in the Godunov forcing"); + } pp.query("use_cc_proj" , m_use_cc_proj); // What type of redistribution algorithm; @@ -128,6 +138,19 @@ void incflo::ReadParameters () amrex::Abort("We cannot have use_tensor_correction be true and diffusion type not Implicit"); } + // With use_tensor_correction, compute_divtau redefines divtau as the + // (tensor - scalar) difference, because the scalar part is handled by the + // implicit solve in update_velocity. ld.divtau_o is then *not* the full + // explicit viscous term that the Godunov edge-state forcing expects, so + // including it there would give the predictor essentially no viscous + // forcing instead of div(eta grad u)/rho. + if (use_tensor_correction && m_godunov_include_diff_in_forcing) { + m_godunov_include_diff_in_forcing = false; + amrex::Print() << "WARNING: incflo.use_tensor_correction = 1 sets " + "godunov_include_diff_in_forcing = 0, because divtau_o then holds " + "only the tensor-minus-scalar correction, not the full viscous term\n"; + } + if (m_advection_type == "MOL" && m_cfl > 0.5) { amrex::Abort("We currently require cfl <= 0.5 when using the MOL advection scheme"); } @@ -149,6 +172,12 @@ void incflo::ReadParameters () pp.query("ic_v", m_ic_v); pp.query("ic_w", m_ic_w); pp.query("ic_p", m_ic_p); + if (std::abs(m_ic_p) > Real(0.0)) { + // set_background_pressure copies ic_p into m_p000, and nothing reads + // m_p000, so a non-zero initial pressure would be silently dropped. + amrex::Abort("incflo.ic_p is not implemented: it only sets m_p000, which no code " + "reads, so a non-zero initial pressure would be silently ignored"); + } if ( !pp.queryarr("ic_t", m_ic_t, 0, m_ntrac) ) { m_ic_t.resize(m_ntrac, 0.); } @@ -446,6 +475,21 @@ void incflo::ReadIOParameters() m_smallplotVars.clear(); pp.queryarr("smallplotVariables", m_smallplotVars); } + + // "divu" is accepted by the parsing above (amr.plt_divu, amr.plotVariables or + // amr.smallplotVariables), but WritePlotVariables has no implementation for it + // and only aborts. Refuse it here rather than after the whole initialization, + // at the first plotfile. + for (auto const& v : m_plotVars) { + if (v == "divu") { + amrex::Abort("plotfile variable 'divu' (amr.plt_divu) is not implemented"); + } + } + for (auto const& v : m_smallplotVars) { + if (v == "divu") { + amrex::Abort("smallplotfile variable 'divu' is not implemented"); + } + } } // @@ -624,10 +668,15 @@ void incflo::InitialPressureProjection() // Always zero this here Vector Source(finest_level+1, nullptr); + // The field being projected is a body force, (rho-rho0)/rho * g, not a velocity, + // so its ghost cells must not be filled with the inflow velocity (and there is + // nothing meaningful for enforceInOutSolvability to rescale). Leave the ghost + // cells at zero, as the incremental/small-dt projections do, so that the + // boundary nodes see grad(phi).n = F.n, i.e. hydrostatic balance. // FIXME FIXME FIXME - THIS ONLY WORKS RIGHT FOR NODAL PROJ ApplyProjection(get_density_new_const(), GetVecOfPtrs(vel), Source, m_cur_time, dummy_dt, false /*incremental*/, - true /*set_inflow_bc*/); + false /*set_inflow_bc*/); } #ifdef AMREX_USE_EB diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index e536cb1ee..e31d69e9b 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -2,6 +2,11 @@ #include #include #include +#ifdef AMREX_USE_EB +#include +#endif + +#include using namespace amrex; @@ -198,6 +203,21 @@ void incflo::ReadCheckpointFile() // Set up problem domain RealBox rb(prob_lo, prob_hi); +#ifdef AMREX_USE_EB + // The EB index space was built in the incflo constructor from the inputs' + // geometry.prob_lo/prob_hi; the physical domain is about to be reset from the + // checkpoint header, so refuse a mismatch instead of silently leaving the EB at + // the wrong physical location. + if (!EB2::IndexSpace::top().getLevel(Geom(0)).isAllRegular()) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (std::abs(prob_lo[idim] - Geom(0).ProbLo(idim)) > Real(1.e-12) || + std::abs(prob_hi[idim] - Geom(0).ProbHi(idim)) > Real(1.e-12)) { + amrex::Abort("ReadCheckpointFile: geometry.prob_lo/prob_hi differ from the " + "checkpoint, but the EB geometry was already built from the inputs"); + } + } + } +#endif Geometry::ResetDefaultProbDomain(rb); for (int lev = 0; lev <= max_level; ++lev) { SetGeometry(lev, Geometry(Geom(lev).Domain(), rb, Geom(lev).CoordInt(), @@ -679,9 +699,8 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo ++icomp; } else if (vars[n] == "divu") { - amrex::Abort("plt_divu: xxxxx TODO"); - pltscaVarsName.push_back("divu"); - ++icomp; + // Refused in ReadIOParameters; kept here so the list stays exhaustive. + amrex::Abort("plotfile variable 'divu' (amr.plt_divu) is not implemented"); } else if (vars[n] == "particle_count") { #ifdef INCFLO_USE_PARTICLES diff --git a/test_2d/GNUmakefile b/test_2d/GNUmakefile index 67768a03a..58e06a97d 100644 --- a/test_2d/GNUmakefile +++ b/test_2d/GNUmakefile @@ -27,4 +27,8 @@ USE_EB = TRUE USE_CSG = FALSE +# Set to TRUE to build the tracer-particle module (src/particles, and every +# INCFLO_USE_PARTICLES block); needed for incflo.use_tracer_particles = 1. +USE_PARTICLES = FALSE + include ../src/Make.incflo diff --git a/test_2d/benchmark.inviscid_rotated_planes_x b/test_2d/benchmark.inviscid_rotated_planes_x index 35802d299..4e0adaf1b 100644 --- a/test_2d/benchmark.inviscid_rotated_planes_x +++ b/test_2d/benchmark.inviscid_rotated_planes_x @@ -70,7 +70,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 45 cylinder.rotation_axe = 2 diff --git a/test_2d/benchmark.rotated_planes_x b/test_2d/benchmark.rotated_planes_x index ab54672c7..7fecab0e0 100644 --- a/test_2d/benchmark.rotated_planes_x +++ b/test_2d/benchmark.rotated_planes_x @@ -74,7 +74,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_2d/benchmark.rotated_planes_y b/test_2d/benchmark.rotated_planes_y index bf1f9fde1..c4e2cbee6 100644 --- a/test_2d/benchmark.rotated_planes_y +++ b/test_2d/benchmark.rotated_planes_y @@ -72,7 +72,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_2d/benchmark.tracer_adv_diff b/test_2d/benchmark.tracer_adv_diff index cd04b0a45..1d0cef7f7 100644 --- a/test_2d/benchmark.tracer_adv_diff +++ b/test_2d/benchmark.tracer_adv_diff @@ -63,7 +63,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_2d/benchmark.tracer_adv_diff_cn b/test_2d/benchmark.tracer_adv_diff_cn index d48df307a..de9a1e1b6 100644 --- a/test_2d/benchmark.tracer_adv_diff_cn +++ b/test_2d/benchmark.tracer_adv_diff_cn @@ -63,7 +63,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_2d/benchmark.tracer_advection b/test_2d/benchmark.tracer_advection index 0acd27d1e..362866ed3 100644 --- a/test_2d/benchmark.tracer_advection +++ b/test_2d/benchmark.tracer_advection @@ -58,7 +58,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_2d/benchmark.vortex_in_sphere b/test_2d/benchmark.vortex_in_sphere index d0451d4b4..b067a80ae 100644 --- a/test_2d/benchmark.vortex_in_sphere +++ b/test_2d/benchmark.vortex_in_sphere @@ -32,7 +32,7 @@ incflo.mu = 0.0 # Dynamic viscosity coefficient amr.n_cell = 64 64 # Grid cells at coarsest AMRlevel amr.max_level = 0 # Max AMR level in hierarchy amr.v = 1 -incflo.v = 1 +incflo.verbose = 1 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # GEOMETRY # diff --git a/test_3d/GNUmakefile b/test_3d/GNUmakefile index df7e424ea..231489204 100644 --- a/test_3d/GNUmakefile +++ b/test_3d/GNUmakefile @@ -26,4 +26,8 @@ USE_EB = TRUE USE_CSG = FALSE +# Set to TRUE to build the tracer-particle module (src/particles, and every +# INCFLO_USE_PARTICLES block); needed for incflo.use_tracer_particles = 1. +USE_PARTICLES = FALSE + include ../src/Make.incflo diff --git a/test_3d/benchmark.rotated_cylinder_x b/test_3d/benchmark.rotated_cylinder_x index a291c218f..c96c072b2 100644 --- a/test_3d/benchmark.rotated_cylinder_x +++ b/test_3d/benchmark.rotated_cylinder_x @@ -86,7 +86,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_3d/benchmark.rotated_cylinder_y b/test_3d/benchmark.rotated_cylinder_y index 832ec7660..baf55a064 100644 --- a/test_3d/benchmark.rotated_cylinder_y +++ b/test_3d/benchmark.rotated_cylinder_y @@ -85,7 +85,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 0 diff --git a/test_3d/benchmark.rotated_cylinder_z b/test_3d/benchmark.rotated_cylinder_z index fad813ac2..2e4bba828 100644 --- a/test_3d/benchmark.rotated_cylinder_z +++ b/test_3d/benchmark.rotated_cylinder_z @@ -85,7 +85,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 1 diff --git a/test_3d/benchmark.slanted_cylinder b/test_3d/benchmark.slanted_cylinder index 3073fb480..998121db4 100644 --- a/test_3d/benchmark.slanted_cylinder +++ b/test_3d/benchmark.slanted_cylinder @@ -71,7 +71,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_3d/benchmark.tracer_adv_diff b/test_3d/benchmark.tracer_adv_diff index aba1f31e5..c1529cbe9 100644 --- a/test_3d/benchmark.tracer_adv_diff +++ b/test_3d/benchmark.tracer_adv_diff @@ -63,7 +63,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_3d/benchmark.tracer_adv_diff_cn b/test_3d/benchmark.tracer_adv_diff_cn index f9d311bb2..be4e1833a 100644 --- a/test_3d/benchmark.tracer_adv_diff_cn +++ b/test_3d/benchmark.tracer_adv_diff_cn @@ -68,7 +68,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_3d/benchmark.tracer_advection b/test_3d/benchmark.tracer_advection index ddd6e855b..8124a5632 100644 --- a/test_3d/benchmark.tracer_advection +++ b/test_3d/benchmark.tracer_advection @@ -62,7 +62,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.49 -cylinder.height = -1.0 cylinder.direction = 0 cylinder.center = 0.5 1.5 1.5 diff --git a/test_3d/inputs.channel_sphere b/test_3d/inputs.channel_sphere index 03d4d124a..4191acccb 100644 --- a/test_3d/inputs.channel_sphere +++ b/test_3d/inputs.channel_sphere @@ -17,7 +17,7 @@ incflo.cfl = 0.45 # CFL factor amr.plot_int = 100 # Steps between plot files amr.check_int = 1000 # Steps between checkpoint files amr.restart = "" # Checkpoint to restart from -amr.plot_p = 1 +amr.plt_p = 1 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # PHYSICS # diff --git a/test_3d/inputs.couette_cylinder b/test_3d/inputs.couette_cylinder index 4fea7da7a..325b03dcb 100644 --- a/test_3d/inputs.couette_cylinder +++ b/test_3d/inputs.couette_cylinder @@ -50,7 +50,9 @@ cylinder.internal_flow = false cylinder.radius = 0.5 cylinder.direction = 2 cylinder.center = 3. 3. .0 -cylinder.speed = 0.15 +# NOTE: a rotating EB wall is not implemented -- nothing in incflo reads +# cylinder.speed, so this key is kept only as a record of the intent. +# cylinder.speed = 0.15 # Boundary conditions ylo.type = "nsw" diff --git a/test_3d/inputs.rotating_cylinder b/test_3d/inputs.rotating_cylinder index 082b8267d..11d950493 100644 --- a/test_3d/inputs.rotating_cylinder +++ b/test_3d/inputs.rotating_cylinder @@ -1,3 +1,8 @@ +# WARNING: as shipped this deck has no flow driver (fully periodic, probtype 0, +# zero initial velocity, zero gravity and no delp), because the rotating +# cylinder that cylinder.speed below was meant to prescribe is not +# implemented. It therefore simulates a fluid at rest. + #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # SIMULATION STOP # #.......................................# @@ -48,7 +53,9 @@ cylinder.internal_flow = false cylinder.radius = 0.5 cylinder.direction = 2 cylinder.center = 3. 3. .0 -cylinder.speed = 0.1 +# NOTE: a rotating EB wall is not implemented -- nothing in incflo reads +# cylinder.speed, so this key is kept only as a record of the intent. +# cylinder.speed = 0.1 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # INITIAL CONDITIONS # diff --git a/test_3d/inputs.taylor_couette b/test_3d/inputs.taylor_couette index 0bd4988f3..2784f2bd0 100644 --- a/test_3d/inputs.taylor_couette +++ b/test_3d/inputs.taylor_couette @@ -1,3 +1,8 @@ +# WARNING: as shipped this deck has no flow driver (fully periodic, zero initial +# velocity, zero gravity and zero delp), because the rotating cylinders +# that cylinder.speed below was meant to prescribe are not implemented. +# It therefore simulates a fluid at rest. + #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # SIMULATION STOP # #.......................................# @@ -47,7 +52,9 @@ annulus.outer_radius = 4.0 # Outer cylinder radius annulus.inner_radius = 1.0 # Inner cylinder radius annulus.outer_center = 5.0 5.0 0. # Outer cylinder center annulus.inner_center = 5.0 5.0 0. # Inner cylinder center -cylinder.speed = 0.1 +# NOTE: a rotating EB wall is not implemented -- nothing in incflo reads +# cylinder.speed, so this key is kept only as a record of the intent. +# cylinder.speed = 0.1 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # NUMERICAL PARAMETERS # diff --git a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_x b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_x index 78ef20091..73f0c65fe 100644 --- a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_x +++ b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_x @@ -90,7 +90,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_y b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_y index 4ee87ae43..9c0464552 100644 --- a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_y +++ b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_y @@ -91,7 +91,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 0 diff --git a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_z b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_z index e1bc4e030..a7d4f5060 100644 --- a/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_z +++ b/test_3d/test_redistribution/benchmark.inviscid_rotated_cylinder_z @@ -89,7 +89,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 1 diff --git a/test_3d/test_redistribution/benchmark.rotated_cylinder_x b/test_3d/test_redistribution/benchmark.rotated_cylinder_x index 694b3b974..a84944c92 100644 --- a/test_3d/test_redistribution/benchmark.rotated_cylinder_x +++ b/test_3d/test_redistribution/benchmark.rotated_cylinder_x @@ -90,7 +90,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 2 diff --git a/test_3d/test_redistribution/benchmark.rotated_cylinder_y b/test_3d/test_redistribution/benchmark.rotated_cylinder_y index f7e4c0df2..5cc0aecff 100644 --- a/test_3d/test_redistribution/benchmark.rotated_cylinder_y +++ b/test_3d/test_redistribution/benchmark.rotated_cylinder_y @@ -89,7 +89,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 0 diff --git a/test_3d/test_redistribution/benchmark.rotated_cylinder_z b/test_3d/test_redistribution/benchmark.rotated_cylinder_z index 6697e90e7..b3d9d93ca 100644 --- a/test_3d/test_redistribution/benchmark.rotated_cylinder_z +++ b/test_3d/test_redistribution/benchmark.rotated_cylinder_z @@ -89,7 +89,6 @@ incflo.geometry = "cylinder" cylinder.internal_flow = true cylinder.radius = 0.172 -cylinder.height = -1.0 cylinder.rotation = 30 cylinder.rotation_axe = 1 diff --git a/test_no_eb_2d/GNUmakefile b/test_no_eb_2d/GNUmakefile index ef9332d27..2b5781afd 100644 --- a/test_no_eb_2d/GNUmakefile +++ b/test_no_eb_2d/GNUmakefile @@ -26,4 +26,8 @@ USE_EB = FALSE USE_CSG = FALSE +# Set to TRUE to build the tracer-particle module (src/particles, and every +# INCFLO_USE_PARTICLES block); needed for incflo.use_tracer_particles = 1. +USE_PARTICLES = FALSE + include ../src/Make.incflo diff --git a/test_no_eb_3d/GNUmakefile b/test_no_eb_3d/GNUmakefile index 0299d3eba..6dae0cb2b 100644 --- a/test_no_eb_3d/GNUmakefile +++ b/test_no_eb_3d/GNUmakefile @@ -26,4 +26,8 @@ USE_EB = FALSE USE_CSG = FALSE +# Set to TRUE to build the tracer-particle module (src/particles, and every +# INCFLO_USE_PARTICLES block); needed for incflo.use_tracer_particles = 1. +USE_PARTICLES = FALSE + include ../src/Make.incflo From e354cf7ef540ae8d5078b15caa7c447d281f4b67 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 21 Sep 2026 16:13:56 -0700 Subject: [PATCH 3/4] Fix -Werror=null-dereference in the particle_count / particle-refinement paths incflo::ErrorEst and incflo::WritePlotVariables both dereferenced the result of ParticleData::operator[] without checking it, and that operator returns nullptr when no species of the given name exists. With INCFLO_PARTICLES=ON GCC inlines ParticleContainer::Increment into both call sites and reports the potential null dereference, which fails the GCC NO EB 2D CI job now that it builds particle support. Look the container up once outside the level loop, abort if it is missing, and use Increment rather than IncrementWithTotal in ErrorEst: the total was discarded, so it was a global reduction per level for nothing. Co-Authored-By: Claude Opus 5 (1M context) --- src/incflo_tagging.cpp | 7 ++++++- src/utilities/io.cpp | 12 ++++++++++-- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/incflo_tagging.cpp b/src/incflo_tagging.cpp index 23aed27c3..4e7ebcf0c 100644 --- a/src/incflo_tagging.cpp +++ b/src/incflo_tagging.cpp @@ -156,13 +156,18 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/) // then come back at the next regridding // const auto& particles_namelist( particleData.getNames() ); + auto* pc = particleData[particles_namelist[0]]; + if (pc == nullptr) { + amrex::Abort("incflo::ErrorEst: no particle container named " + +particles_namelist[0]); + } std::unique_ptr mf = std::make_unique(grids[levc], dmap[levc], 1, 0); mf->setVal(0.0); IntVect rr = IntVect::TheUnitVector(); for (int lev = levc; lev <= finest_level; lev++) { MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0); - particleData[particles_namelist[0]]->IncrementWithTotal(temp_dat, lev); + pc->Increment(temp_dat, lev); MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0); diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index e31d69e9b..1dc7da43f 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -705,11 +705,19 @@ void incflo::WritePlotVariables(Vector vars, const std::string& plo else if (vars[n] == "particle_count") { #ifdef INCFLO_USE_PARTICLES const auto& particles_namelist( particleData.getNames() ); + incflo_PC* pc = nullptr; + if (!particles_namelist.empty()) { + pc = particleData[particles_namelist[0]]; + if (pc == nullptr) { + amrex::Abort("incflo::WritePlotVariables: no particle container named " + +particles_namelist[0]); + } + } for (int lev = 0; lev <= finest_level; ++lev) { MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0); temp_dat.setVal(0); - if (!particles_namelist.empty()) { - particleData[particles_namelist[0]]->Increment(temp_dat, lev); + if (pc != nullptr) { + pc->Increment(temp_dat, lev); } MultiFab::Copy(mf[lev], temp_dat, 0, icomp, 1, 0); } From 129c6582063c8c6622a074b45e1c51d9f26e261e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 21 Sep 2026 17:00:44 -0700 Subject: [PATCH 4/4] Fix the clang-tidy errors in the newly-linted particle sources The GCC NO EB 2D job runs clang-tidy with --warnings-as-errors=*, and now that it builds with particles the src/particles files are analysed for the first time. They are clean at the compiler level but tripped eight checks: - clang-analyzer-optin.cplusplus.VirtualCall: both constructors called the virtual readInputs() during construction; qualify the calls so they are static, which is what actually happens anyway. - readability-inconsistent-declaration-parameter-name: the MAC-velocity arguments of EvolveParticles/AdvectWithFlow and the box argument of initializeParticlesUniformDistributionInBox are named differently in the header and in the definition; the header now matches the definitions. - bugprone-narrowing-conversions: make the Long -> int conversion of numPts() explicit. - readability-qualified-auto, modernize-use-auto, modernize-use-bool-literals, readability-redundant-string-cstr, readability-redundant-control-flow and performance-unnecessary-copy-initialization: mechanical cleanups. Also guard the two remaining unchecked dereferences of ParticleData's operator[] in incflo_Tracers.cpp, matching the two fixed in the previous commit, and make all three null branches end in a return so that neither GCC nor the clang static analyzer treats the code after the Abort as reachable with a null pointer. Co-Authored-By: Claude Opus 5 (1M context) --- src/incflo_tagging.cpp | 1 + src/particles/incflo_PC.H | 22 +++++++++++-------- src/particles/incflo_PCEvolve.cpp | 1 - src/particles/incflo_PCInit.cpp | 7 ++----- src/particles/incflo_PCUtils.cpp | 2 -- src/particles/incflo_Tracers.cpp | 35 +++++++++++++++++++------------ 6 files changed, 38 insertions(+), 30 deletions(-) diff --git a/src/incflo_tagging.cpp b/src/incflo_tagging.cpp index 4e7ebcf0c..1e82a0a42 100644 --- a/src/incflo_tagging.cpp +++ b/src/incflo_tagging.cpp @@ -160,6 +160,7 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/) if (pc == nullptr) { amrex::Abort("incflo::ErrorEst: no particle container named " +particles_namelist[0]); + return; // not reached: Abort() does not return } std::unique_ptr mf = std::make_unique(grids[levc], dmap[levc], 1, 0); mf->setVal(0.0); diff --git a/src/particles/incflo_PC.H b/src/particles/incflo_PC.H index 2ab623414..8b31c8725 100644 --- a/src/particles/incflo_PC.H +++ b/src/particles/incflo_PC.H @@ -72,7 +72,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: { BL_PROFILE("incflo_PC::incflo_PC()"); m_name = a_name; - readInputs(); + // Qualified so that this is a static call: a derived class is not yet + // constructed here, so virtual dispatch would not reach its override. + incflo_PC::readInputs(); } /*! Constructor */ @@ -88,7 +90,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: { BL_PROFILE("incflo_PC::incflo_PC()"); m_name = a_name; - readInputs(); + // Qualified so that this is a static call: a derived class is not yet + // constructed here, so virtual dispatch would not reach its override. + incflo_PC::readInputs(); } /*! Initialize particles in domain */ @@ -102,9 +106,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: * responsible for calling Redistribute() once all levels are done */ virtual void EvolveParticles (int, amrex::Real, - AMREX_D_DECL(const amrex::MultiFab* a_u_mac, - const amrex::MultiFab* a_v_mac, - const amrex::MultiFab* a_w_mac)); + AMREX_D_DECL(const amrex::MultiFab* a_umac, + const amrex::MultiFab* a_vmac, + const amrex::MultiFab* a_wmac)); /*! Get real-type particle attribute names */ virtual amrex::Vector varNames () const @@ -123,9 +127,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: /*! Uses midpoint method to advance particles using flow velocity. */ virtual void AdvectWithFlow ( int, amrex::Real, - AMREX_D_DECL(const amrex::MultiFab* a_u_mac, - const amrex::MultiFab* a_v_mac, - const amrex::MultiFab* a_w_mac)); + AMREX_D_DECL(const amrex::MultiFab* a_umac, + const amrex::MultiFab* a_vmac, + const amrex::MultiFab* a_wmac)); /*! Compute mass density */ virtual void massDensity ( amrex::MultiFab&, const int&, const int& a_comp = 0) const; @@ -153,7 +157,7 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: // public due to CUDA extended lambda capture rules /*! Default particle initialization */ - void initializeParticlesUniformDistributionInBox ( const amrex::RealBox& particle_box + void initializeParticlesUniformDistributionInBox ( const amrex::RealBox& particle_init_domain #ifdef AMREX_USE_EB ,amrex::EBFArrayBoxFactory const& ebfact #endif diff --git a/src/particles/incflo_PCEvolve.cpp b/src/particles/incflo_PCEvolve.cpp index 8283feab1..c06153af8 100644 --- a/src/particles/incflo_PCEvolve.cpp +++ b/src/particles/incflo_PCEvolve.cpp @@ -24,7 +24,6 @@ void incflo_PC::EvolveParticles ( int a_l // hand a particle that crossed from level lev into a level lev+1 grid // to level lev+1, which then advects it a second time in the same step. // - return; } // diff --git a/src/particles/incflo_PCInit.cpp b/src/particles/incflo_PCInit.cpp index 4ef53e3c7..f8f3bef78 100644 --- a/src/particles/incflo_PCInit.cpp +++ b/src/particles/incflo_PCInit.cpp @@ -45,8 +45,6 @@ void incflo_PC::readInputs () m_advect_w_flow = (m_name == incfloParticleNames::tracers ? true : false); pp.query("advect_with_flow", m_advect_w_flow); - - return; } /*! Initialize particles in domain */ @@ -71,7 +69,6 @@ void incflo_PC::InitializeParticles ( << m_name << " particle species.\n"; Error("See error message!"); } - return; } /*! Uniform distribution: the number of particles per grid cell is specified @@ -132,7 +129,7 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par int np = 0; { - int ncell = num_particles[mfi].numPts(); + int ncell = static_cast(num_particles[mfi].numPts()); const int* in = num_particles[mfi].dataPtr(); int* out = offsets[mfi].dataPtr(); np = Scan::PrefixSum( ncell, @@ -145,7 +142,7 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par auto& particle_tile = DefineAndReturnParticleTile(lev, mfi); particle_tile.resize(np); - auto aos = &particle_tile.GetArrayOfStructs()[0]; + auto* aos = &particle_tile.GetArrayOfStructs()[0]; auto& soa = particle_tile.GetStructOfArrays(); AMREX_D_TERM(auto* vx_ptr = soa.GetRealData(incflo_ParticlesRealIdxSoA::vx).data();, auto* vy_ptr = soa.GetRealData(incflo_ParticlesRealIdxSoA::vy).data();, diff --git a/src/particles/incflo_PCUtils.cpp b/src/particles/incflo_PCUtils.cpp index 0dea62e19..6123cf702 100644 --- a/src/particles/incflo_PCUtils.cpp +++ b/src/particles/incflo_PCUtils.cpp @@ -34,8 +34,6 @@ void incflo_PC::massDensity ( MultiFab& a_mf, return mass*inv_cell_volume; }); }); - - return; } #endif diff --git a/src/particles/incflo_Tracers.cpp b/src/particles/incflo_Tracers.cpp index 19af633d7..78b5e311b 100644 --- a/src/particles/incflo_Tracers.cpp +++ b/src/particles/incflo_Tracers.cpp @@ -11,14 +11,13 @@ void incflo::readTracerParticlesParams () { ParmParse pp("incflo"); - m_use_tracer_particles = 0; + m_use_tracer_particles = false; - pp.query(std::string("use_"+incfloParticleNames::tracers).c_str(), m_use_tracer_particles); + pp.query("use_"+incfloParticleNames::tracers, m_use_tracer_particles); if (m_use_tracer_particles) { particleData.addName(incfloParticleNames::tracers); } - return; } /*! Initialize tracer particles */ @@ -32,12 +31,12 @@ void incflo::initializeTracerParticles ( ParGDBBase* a_gdb for (auto it = namelist_unalloc.begin(); it != namelist_unalloc.end(); ++it) { - std::string species_name( *it ); + const std::string& species_name( *it ); if (species_name == incfloParticleNames::tracers) { AMREX_ASSERT(m_use_tracer_particles); - incflo_PC* pc = new incflo_PC(a_gdb, incfloParticleNames::tracers); + auto* pc = new incflo_PC(a_gdb, incfloParticleNames::tracers); #ifdef AMREX_USE_EB pc->InitializeParticles(ebfact); #else @@ -55,17 +54,21 @@ void incflo::initializeTracerParticles ( ParGDBBase* a_gdb // just created; the restart path gets this from ParticleContainer::Restart. const auto& particles_namelist( particleData.getNames() ); for (auto it = particles_namelist.begin(); it != particles_namelist.end(); ++it) { - std::string species_name( *it ); + const std::string& species_name( *it ); if (species_name == incfloParticleNames::tracers) { - if (!particleData[incfloParticleNames::tracers]->OK()) { - particleData[incfloParticleNames::tracers]->resizeData(); + auto* pc = particleData[incfloParticleNames::tracers]; + if (pc == nullptr) { + amrex::Abort("incflo::initializeTracerParticles: no particle container named " + +incfloParticleNames::tracers); + return; // not reached: Abort() does not return } - particleData[incfloParticleNames::tracers]->Redistribute(); + if (!pc->OK()) { + pc->resizeData(); + } + pc->Redistribute(); } } - - return; } /*! Evolve tracers particles for one time step*/ @@ -74,10 +77,16 @@ void incflo::evolveTracerParticles (AMREX_D_DECL(Vector const& Vector const& w_mac)) { if (m_use_tracer_particles) { + auto* pc = particleData[incfloParticleNames::tracers]; + if (pc == nullptr) { + amrex::Abort("incflo::evolveTracerParticles: no particle container named " + +incfloParticleNames::tracers); + return; // not reached: Abort() does not return + } for (int lev = 0; lev <= finest_level; ++lev) { - particleData[incfloParticleNames::tracers]->EvolveParticles(lev, m_dt, - AMREX_D_DECL(u_mac[lev],v_mac[lev],w_mac[lev])); + pc->EvolveParticles(lev, m_dt, + AMREX_D_DECL(u_mac[lev],v_mac[lev],w_mac[lev])); } //