diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ff039a5..6f2692e 100644 --- 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 index c1ab220..ffcaa5a 100644 --- 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 index f333e3f..09a7b3a 100644 --- 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 index 74b0362..7eeeff7 100644 --- 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,14 +52,12 @@ public: const amrex::Vector& ref_ratio, const amrex::Vector& bcs, BCOpCreator bc_op, - amrex::Vector& f_, - bool is_energy_lattice = false) + amrex::Vector& 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 @@ -52,6 +65,18 @@ public: // 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> 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> 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> 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& 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 index 80db6db..67b22e5 100644 --- 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 index 0b15322..6f05dbb 100644 --- a/Source/LBM.H +++ b/Source/LBM.H @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include "EB.H" @@ -226,7 +227,7 @@ private: amrex::Vector 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 m_bcs; // 1-component + amrex::Vector m_bcs; std::unique_ptr m_fillpatch_op; std::unique_ptr m_fillpatch_g_op; diff --git a/Source/LBM.cpp b/Source/LBM.cpp index 1d3cb31..04285ae 100644 --- 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::GpuTuple { const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); amrex::GpuArray 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>( 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>( 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 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; } diff --git a/Source/Utilities.H b/Source/Utilities.H index 6402690..8c56ac5 100644 --- a/Source/Utilities.H +++ b/Source/Utilities.H @@ -142,15 +142,15 @@ AMREX_GPU_DEVICE AMREX_INLINE void set_extended_grad_expansion_generic( ux * ((jy)-rho * uy) - uy * ((jx)-rho * ux)) * one_bystheta0 * one_bystheta0); - AMREX_3D_ONLY(const amrex::Real a2xz = - ((pxz)-0 - rho * ux * uz - ux * ((jz)-rho * uz) - - uz * ((jx)-rho * ux)) * - one_bystheta0 * one_bystheta0;); - - AMREX_3D_ONLY(const amrex::Real a2yz = - ((pyz)-0 - rho * uy * uz - uy * ((jz)-rho * uz) - - uz * ((jy)-rho * uy)) * - one_bystheta0 * one_bystheta0;); + AMREX_3D_ONLY( + const amrex::Real a2xz = ((pxz)-0 - rho * ux * uz - + ux * ((jz)-rho * uz) - uz * ((jx)-rho * ux)) * + one_bystheta0 * one_bystheta0;); + + AMREX_3D_ONLY( + const amrex::Real a2yz = ((pyz)-0 - rho * uy * uz - + uy * ((jz)-rho * uz) - uz * ((jy)-rho * uy)) * + one_bystheta0 * one_bystheta0;); f = rho + AMREX_D_TERM(a1x * ev[0], +a1y * ev[1], +a1z * ev[2]); diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 5bf555a..4883ca0 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -88,15 +88,18 @@ add_test_r(channel_regular) add_test_r(tg) add_test_r(single_rotated_box) -if(PELEC_DIM EQUAL 2) +if(MARBLES_DIM EQUAL 2) add_test_r(single_cylinder_2d) endif() -if(PELEC_DIM EQUAL 3) +if(MARBLES_DIM EQUAL 3) add_test_r(single_cylinder) add_test_r(channel_cylinder) endif() +add_test_re(sod) +add_test_re(sod_amr) + #============================================================================= # Regression tests excluded from CI #============================================================================= diff --git a/Tests/test_files/sod_test/sod.inp b/Tests/test_files/sod/sod.inp similarity index 53% rename from Tests/test_files/sod_test/sod.inp rename to Tests/test_files/sod/sod.inp index 7643cd0..10f7951 100644 --- a/Tests/test_files/sod_test/sod.inp +++ b/Tests/test_files/sod/sod.inp @@ -1,8 +1,8 @@ # How to use this sod test tube case? -# set ic_sod_test.initial_temperature=lbm.initial_temperature to a value <= 0.3333 +# set ic_sod.initial_temperature=lbm.initial_temperature to a value <= 0.3333 # set lbm.nu to a value <= 0.0100 -# Further instructions in sod_tester.ipynb -# Note. Always: geometry.prob_hi[0]=amr.n_cell[0]=amr.max_grid_size=2.0*ic_sod_test.xDiscontinuity +# Further instructions in soder.ipynb +# Note. Always: geometry.prob_hi[0]=amr.n_cell[0]=2.0*ic_sod.xDiscontinuity max_step = 1500 @@ -14,8 +14,7 @@ geometry.is_periodic = 0 1 1 # timestepping amr.n_cell = 3000 2 2 amr.max_level = 0 -amr.max_grid_size = 3000 -amr.blocking_factor = 1 +amr.blocking_factor = 2 amr.plot_int = 100 amr.chk_int = 100 amr.file_name_digits = 5 @@ -29,20 +28,20 @@ lbm.alpha = 0.010 lbm.save_streaming = 0 lbm.compute_forces = 0 -lbm.ic_type = "sod_test" -ic_sod_test.density = 0.50 -ic_sod_test.velocity = 0.0 0.0 0.0 -ic_sod_test.mach_components = 0.0 0.0 0.0 +lbm.ic_type = "sod" +ic_sod.density = 0.50 +ic_sod.velocity = 0.0 0.0 0.0 +ic_sod.mach_components = 0.0 0.0 0.0 -# Set ic_sod_test.xDiscontinuity equal to half length to place +# Set ic_sod.xDiscontinuity equal to half length to place # discontinuity in the middle -ic_sod_test.x_discontinuity = 1500.0 +ic_sod.x_discontinuity = 1500.0 -ic_sod_test.initial_temperature = 0.20 -ic_sod_test.adiabatic_exponent = 2.0 -ic_sod_test.mean_molecular_mass = 28.96 -ic_sod_test.density_ratio = 4.00 -ic_sod_test.temperature_ratio = 0.1250 +ic_sod.initial_temperature = 0.20 +ic_sod.adiabatic_exponent = 2.0 +ic_sod.mean_molecular_mass = 28.96 +ic_sod.density_ratio = 4.00 +ic_sod.temperature_ratio = 0.1250 lbm.initial_temperature = 0.20 lbm.adiabatic_exponent = 2.0 lbm.mean_molecular_mass = 28.96 diff --git a/Tests/test_files/sod_amr/sod_amr.inp b/Tests/test_files/sod_amr/sod_amr.inp new file mode 100644 index 0000000..96ee73a --- /dev/null +++ b/Tests/test_files/sod_amr/sod_amr.inp @@ -0,0 +1,64 @@ +# How to use this sod test tube case? +# set ic_sod.initial_temperature=lbm.initial_temperature to a value <= 0.3333 +# set lbm.nu to a value <= 0.0100 + +max_step = 500 + +# geometry parameters +geometry.prob_lo = 0.0 -2.0 -2.0 +geometry.prob_hi = 1000.0 2.0 2.0 +geometry.is_periodic = 0 1 1 + +# timestepping +amr.n_cell = 1000 4 4 +amr.regrid_int = 1 +amr.max_level = 1 +amr.blocking_factor = 4 +amr.plot_int = 100 +amr.chk_int = 100 +amr.file_name_digits = 5 +amr.n_error_buf = 8 + +lbm.bc_lo = 5 0 0 +lbm.bc_hi = 5 0 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 +lbm.nu = 0.010 +lbm.alpha = 0.010 +lbm.save_streaming = 1 +lbm.compute_forces = 0 + +lbm.ic_type = "sod" +ic_sod.density = 0.50 +ic_sod.velocity = 0.0 0.0 0.0 +ic_sod.mach_components = 0.0 0.0 0.0 + +# Set ic_sod.xDiscontinuity equal to half length to place +# discontinuity in the middle +ic_sod.x_discontinuity = 500.0 + +ic_sod.initial_temperature = 0.20 +ic_sod.adiabatic_exponent = 2.0 +ic_sod.mean_molecular_mass = 28.96 +ic_sod.density_ratio = 4.00 +ic_sod.temperature_ratio = 0.1250 +lbm.initial_temperature = 0.20 +lbm.adiabatic_exponent = 2.0 +lbm.mean_molecular_mass = 28.96 + +eb2.geom_type = "all_regular" + +tagging.refinement_indicators = rho bc +tagging.rho.max_level = 2 +tagging.rho.adjacent_difference_greater = 0.1 +tagging.rho.field_name = rho + +tagging.bc.in_box_lo = -10 -10 -10 +tagging.bc.in_box_hi = 20 10 10 + + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 1