NatLabRockies · marchdf · Jul 2, 2025 · Jun 27, 2025 · Jun 27, 2025 · Jun 27, 2025
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -14,17 +14,17 @@ concurrency:
 
 jobs:
   Formatting:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-24.04
     steps:
     - name: Clone
       uses: actions/checkout@v4
     - name: Check formatting
-      uses: DoozyX/clang-format-lint-action@v0.18.2
+      uses: DoozyX/clang-format-lint-action@v0.20
       with:
         source: './Source'
         exclude: '.'
         extensions: 'H,cpp'
-        clangFormatVersion: 18
+        clangFormatVersion: 20
   CPU-GNUmake:
     needs: Formatting
     runs-on: ubuntu-latest

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -64,10 +64,6 @@ endif()
 
 init_amrex()
 
-if(MARBLES_ENABLE_MPI)
-  find_package(MPI REQUIRED)
-endif()
-
 message(STATUS "Marbles Information:")
 message(STATUS "CMAKE_SYSTEM_NAME = ${CMAKE_SYSTEM_NAME}")
 message(STATUS "CMAKE_CXX_COMPILER_ID = ${CMAKE_CXX_COMPILER_ID}")

diff --git a/Source/BC.H b/Source/BC.H
@@ -485,6 +485,9 @@ struct BCFill
                             } else if (m_bc_type[bc_idx] == SLIPWALLZNORMAL) {
                                 bounce_back_zbc(iv, q, ev, data);
                             }
+                        } else {
+                            // to avoid spurious NaNs fill with a dummy value
+                            data(iv, q) = -1.0;
                         }
                     }
                 }

diff --git a/Source/FillPatchOps.H b/Source/FillPatchOps.H
@@ -22,6 +22,21 @@ public:
     virtual void
     fillpatch(const int lev, const amrex::Real time, amrex::MultiFab& mf) = 0;
 
+    virtual void fillpatch(
+        const int lev,
+        const amrex::Real time,
+        amrex::MultiFab& mf,
+        const amrex::IntVect& nghost) = 0;
+
+    virtual void
+    physbc(const int lev, const amrex::Real time, amrex::MultiFab& mf) = 0;
+
+    virtual void physbc(
+        const int lev,
+        const amrex::Real time,
+        amrex::MultiFab& mf,
+        const amrex::IntVect& nghost) = 0;
+
     virtual void fillpatch_from_coarse(
         const int lev, const amrex::Real time, amrex::MultiFab& mf) = 0;
 };
@@ -37,21 +52,31 @@ public:
         const amrex::Vector<amrex::IntVect>& ref_ratio,
         const amrex::Vector<amrex::BCRec>& bcs,
         BCOpCreator bc_op,
-        amrex::Vector<amrex::MultiFab>& f_,
-        bool is_energy_lattice = false)
+        amrex::Vector<amrex::MultiFab>& f_)
         : m_geom(geom)
         , m_ref_ratio(ref_ratio)
         , m_bcs(bcs)
         , m_op(std::move(bc_op))
         , m_f(f_)
-        , m_is_energy_lattice(is_energy_lattice)
     {}
 
     // compute a new multifab by coping in phi from valid region and filling
     // ghost cells works for single level and 2-level cases (fill fine grid
     // ghost by interpolating from coarse)
     void fillpatch(
         const int lev, const amrex::Real time, amrex::MultiFab& mf) override
+    {
+        fillpatch(lev, time, mf, mf.nGrowVect());
+    }
+
+    // compute a new multifab by coping in phi from valid region and filling
+    // ghost cells works for single level and 2-level cases (fill fine grid
+    // ghost by interpolating from coarse)
+    void fillpatch(
+        const int lev,
+        const amrex::Real time,
+        amrex::MultiFab& mf,
+        const amrex::IntVect& nghost) override
     {
         BL_PROFILE("FillPatchOps::fillpatch()");
         AMREX_ASSERT(mf.nComp() == constants::N_MICRO_STATES);
@@ -88,7 +113,7 @@ public:
             amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> physbc(
                 m_geom[lev], m_bcs, bc_functor());
             amrex::FillPatchSingleLevel(
-                mf, time, {&(m_f[lev])}, {time}, 0, 0, m_f[lev].nComp(),
+                mf, nghost, time, {&(m_f[lev])}, {time}, 0, 0, m_f[lev].nComp(),
                 m_geom[lev], physbc, 0);
         } else {
             amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> cphysbc(
@@ -97,14 +122,42 @@ public:
                 m_geom[lev], m_bcs, bc_functor());
 
             amrex::FillPatchTwoLevels(
-                mf, time, {&(m_f[lev - 1])}, {time}, {&(m_f[lev])}, {time}, 0,
-                0, m_f[lev].nComp(), m_geom[lev - 1], m_geom[lev], cphysbc, 0,
-                fphysbc, 0, m_ref_ratio[lev - 1], m_mapper, m_bcs, 0);
+                mf, nghost, time, {&(m_f[lev - 1])}, {time}, {&(m_f[lev])},
+                {time}, 0, 0, m_f[lev].nComp(), m_geom[lev - 1], m_geom[lev],
+                cphysbc, 0, fphysbc, 0, m_ref_ratio[lev - 1], m_mapper, m_bcs,
+                0);
         }
 
         amrex::Gpu::synchronize();
     }
 
+    // apply the physical boundary conditions
+    void
+    physbc(const int lev, const amrex::Real time, amrex::MultiFab& mf) override
+    {
+        physbc(lev, time, mf, mf.nGrowVect());
+    }
+
+    // apply the physical boundary conditions
+    void physbc(
+        const int lev,
+        const amrex::Real time,
+        amrex::MultiFab& mf,
+        const amrex::IntVect& nghost) override
+    {
+        BL_PROFILE("FillPatchOps::fillpatch()");
+        AMREX_ASSERT(mf.nComp() == constants::N_MICRO_STATES);
+        AMREX_ASSERT(mf.nComp() == m_f[0].nComp());
+        AMREX_ASSERT(m_bcs.size() == mf.nComp());
+
+        amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> physbc(
+            m_geom[lev], m_bcs, bc_functor());
+
+        physbc(mf, 0, m_f[lev].nComp(), nghost, time, 0);
+
+        amrex::Gpu::synchronize();
+    }
+
     // fill an entire multifab by interpolating from the coarser level
     // this comes into play when a new level of refinement appears
     void fillpatch_from_coarse(
@@ -135,7 +188,6 @@ protected:
     const BCOpCreator m_op;
     amrex::Vector<amrex::MultiFab>& m_f;
     amrex::Interpolater* m_mapper = &amrex::cell_cons_interp;
-    bool m_is_energy_lattice;
 };
 } // namespace lbm
 #endif
diff --git a/Source/IC.H b/Source/IC.H
@@ -419,7 +419,7 @@ struct SodTest
 
     explicit SodTest();
 
-    static std::string identifier() { return "ic_sod_test"; }
+    static std::string identifier() { return "ic_sod"; }
 
     DeviceType device_instance() const { return m_op; }
 

diff --git a/Source/LBM.H b/Source/LBM.H
@@ -9,6 +9,7 @@
 #include <AMReX_ParmParse.H>
 #include <AMReX_FluxRegister.H>
 #include <AMReX_PlotFileUtil.H>
+#include <AMReX_BCUtil.H>
 #include <AMReX_PhysBCFunct.H>
 #include <AMReX_ErrorList.H>
 #include "EB.H"
@@ -226,7 +227,7 @@ private:
     amrex::Vector<amrex::iMultiFab> m_mask;
 
     const int m_macrodata_nghost = 1;
-    const int m_f_nghost = 2;
+    const int m_f_nghost = 3;
     const int m_eq_nghost = 0;
     const int m_derived_nghost = 0;
 
@@ -237,7 +238,7 @@ private:
 
     // this is essentially a 2*DIM integer array storing the physical boundary
     // condition types at the lo/hi walls in each direction
-    amrex::Vector<amrex::BCRec> m_bcs; // 1-component
+    amrex::Vector<amrex::BCRec> m_bcs;
 
     std::unique_ptr<FillPatchOpsBase> m_fillpatch_op;
     std::unique_ptr<FillPatchOpsBase> m_fillpatch_g_op;

diff --git a/Source/LBM.cpp b/Source/LBM.cpp
@@ -353,24 +353,27 @@ void LBM::read_tagging_parameters()
             ppr.get("value_greater", value);
             std::string field;
             ppr.get("field_name", field);
-            m_err_tags.push_back(amrex::AMRErrorTag(
-                value, amrex::AMRErrorTag::GREATER, field, info));
+            m_err_tags.push_back(
+                amrex::AMRErrorTag(
+                    value, amrex::AMRErrorTag::GREATER, field, info));
             itexists = check_field_existence(field);
         } else if (ppr.countval("value_less") > 0) {
             amrex::Real value;
             ppr.get("value_less", value);
             std::string field;
             ppr.get("field_name", field);
-            m_err_tags.push_back(amrex::AMRErrorTag(
-                value, amrex::AMRErrorTag::LESS, field, info));
+            m_err_tags.push_back(
+                amrex::AMRErrorTag(
+                    value, amrex::AMRErrorTag::LESS, field, info));
             itexists = check_field_existence(field);
         } else if (ppr.countval("adjacent_difference_greater") > 0) {
             amrex::Real value;
             ppr.get("adjacent_difference_greater", value);
             std::string field;
             ppr.get("field_name", field);
-            m_err_tags.push_back(amrex::AMRErrorTag(
-                value, amrex::AMRErrorTag::GRAD, field, info));
+            m_err_tags.push_back(
+                amrex::AMRErrorTag(
+                    value, amrex::AMRErrorTag::GRAD, field, info));
             itexists = check_field_existence(field);
         } else if (realbox.ok()) {
             m_err_tags.push_back(amrex::AMRErrorTag(info));
@@ -498,6 +501,10 @@ void LBM::time_step(const int lev, const amrex::Real time, const int iteration)
         m_fillpatch_g_op->fillpatch(lev + 1, m_ts_new[lev + 1], m_g[lev + 1]);
 
         for (int i = 1; i <= m_nsubsteps[lev + 1]; ++i) {
+            m_fillpatch_op->physbc(lev + 1, m_ts_new[lev + 1], m_f[lev + 1]);
+
+            m_fillpatch_g_op->physbc(lev + 1, m_ts_new[lev + 1], m_g[lev + 1]);
+
             time_step(lev + 1, time + (i - 1) * m_dts[lev + 1], i);
         }
     }
@@ -1019,8 +1026,8 @@ void LBM::compute_eb_forces()
             m_f[lev], amrex::IntVect(0),
             [=] AMREX_GPU_DEVICE(
                 int nbx, int i, int j, int AMREX_D_PICK(, /*k*/, k)) noexcept
-            -> amrex::GpuTuple<AMREX_D_DECL(
-                amrex::Real, amrex::Real, amrex::Real)> {
+                -> amrex::GpuTuple<AMREX_D_DECL(
+                    amrex::Real, amrex::Real, amrex::Real)> {
                 const amrex::IntVect iv(AMREX_D_DECL(i, j, k));
                 amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> fs = {0.0};
                 if ((is_fluid_arrs[nbx](iv, 1) == 1) &&
@@ -1312,14 +1319,61 @@ void LBM::RemakeLevel(
     const amrex::DistributionMapping& dm)
 {
     BL_PROFILE("LBM::RemakeLevel()");
-    amrex::Abort("RemakeLevel not implemented");
-    amrex::MultiFab new_state(
-        ba, dm, m_macrodata[lev].nComp(), m_macrodata[lev].nGrow());
 
-    std::swap(new_state, m_macrodata[lev]);
+    if (Verbose() > 0) {
+        amrex::Print() << "Remaking level " << lev << std::endl;
+    }
+
+    m_factory[lev] = amrex::makeEBFabFactory(
+        Geom(lev), ba, dm, {5, 5, 5}, amrex::EBSupport::basic);
+    amrex::MultiFab new_f(
+        ba, dm, constants::N_MICRO_STATES, m_f_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+    amrex::MultiFab new_g(
+        ba, dm, constants::N_MICRO_STATES, m_f_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+
+    m_fillpatch_op->fillpatch(lev, time, new_f);
+    m_fillpatch_g_op->fillpatch(lev, time, new_g);
+
+    std::swap(new_f, m_f[lev]);
+    std::swap(new_g, m_g[lev]);
+
+    m_macrodata[lev].define(
+        ba, dm, constants::N_MACRO_STATES, m_macrodata_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+    m_is_fluid[lev].define(ba, dm, constants::N_IS_FLUID, m_f[lev].nGrow());
+    m_eq[lev].define(
+        ba, dm, constants::N_MICRO_STATES, m_eq_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+    m_eq_g[lev].define(
+        ba, dm, constants::N_MICRO_STATES, m_eq_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+    m_derived[lev].define(
+        ba, dm, constants::N_DERIVED, m_derived_nghost, amrex::MFInfo(),
+        *(m_factory[lev]));
+    m_mask[lev].define(ba, dm, 1, 0);
+
+    initialize_is_fluid(lev);
+    initialize_mask(lev);
+    fill_f_inside_eb(lev);
+    m_f[lev].FillBoundary(Geom(lev).periodicity());
+    m_g[lev].FillBoundary(Geom(lev).periodicity());
+    m_macrodata[lev].setVal(0.0);
+    m_eq[lev].setVal(0.0);
+    m_eq_g[lev].setVal(0.0);
+    m_derived[lev].setVal(0.0);
+
+    f_to_macrodata(lev);
+
+    macrodata_to_equilibrium(lev);
+
+    compute_derived(lev);
+
+    compute_q_corrections(lev);
 
     m_ts_new[lev] = time;
-    m_ts_old[lev] = constants::LOW_NUM;
+    m_ts_old[lev] = time - constants::SMALL_NUM;
 }
 
 // Delete level data
@@ -1357,7 +1411,7 @@ void LBM::set_bcs()
             VelBCOp(
                 m_mesh_speed, m_bc_type, m_g[0].nGrowVect(),
                 is_an_energy_lattice),
-            m_g, is_an_energy_lattice);
+            m_g);
 
     } else if (m_velocity_bc_type == "constant") {
 
@@ -1372,7 +1426,7 @@ void LBM::set_bcs()
             VelBCOp(
                 m_mesh_speed, m_bc_type, m_g[0].nGrowVect(),
                 is_an_energy_lattice),
-            m_g, is_an_energy_lattice);
+            m_g);
 
     } else if (m_velocity_bc_type == "channel") {
 
@@ -1387,7 +1441,7 @@ void LBM::set_bcs()
             VelBCOp(
                 m_mesh_speed, m_bc_type, m_g[0].nGrowVect(),
                 is_an_energy_lattice),
-            m_g, is_an_energy_lattice);
+            m_g);
 
     } else if (m_velocity_bc_type == "parabolic") {
 
@@ -1402,7 +1456,7 @@ void LBM::set_bcs()
             VelBCOp(
                 m_mesh_speed, m_bc_type, m_g[0].nGrowVect(),
                 is_an_energy_lattice),
-            m_g, is_an_energy_lattice);
+            m_g);
 
     } else {
         amrex::Abort("LBM::set_bcs(): Unknown velocity BC");
@@ -1425,7 +1479,7 @@ void LBM::set_ics()
         m_ic_op = std::make_unique<ic::Initializer<ic::ThermalDiffusivityTest>>(
             m_mesh_speed,
             ic::ThermalDiffusivityTest(ic::ThermalDiffusivityTest()), m_f, m_g);
-    } else if (m_ic_type == "sod_test") {
+    } else if (m_ic_type == "sod") {
         m_ic_op = std::make_unique<ic::Initializer<ic::SodTest>>(
             m_mesh_speed, ic::SodTest(ic::SodTest()), m_f, m_g);
     } else {
@@ -1504,6 +1558,15 @@ LBM::get_field(const std::string& name, const int lev, const int ngrow)
         amrex::Gpu::synchronize();
     }
 
+    amrex::Vector<amrex::BCRec> bcs(nc);
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        for (auto& bc : bcs) {
+            bc.setLo(idim, amrex::BCType::foextrap);
+            bc.setHi(idim, amrex::BCType::foextrap);
+        }
+    }
+    amrex::FillDomainBoundary(*mf, Geom(lev), bcs);
+
     return mf;
 }