diff --git a/Source/LBM.cpp b/Source/LBM.cpp index f6d7d11..3dc17be 100644 --- a/Source/LBM.cpp +++ b/Source/LBM.cpp @@ -1019,13 +1019,39 @@ void LBM::compute_eb_forces() amrex::GpuArray fs = {0.0}; if ((is_fluid_arrs[nbx](iv, 1) == 1) && (mask_arrs[nbx](iv) == 0)) { + + const auto f_arr = f_arrs[nbx]; + const auto is_fluid_arr = is_fluid_arrs[nbx]; + + const auto& lb = amrex::lbound(f_arr); + const auto& ub = amrex::ubound(f_arr); + const amrex::Box fbox( + amrex::IntVect(AMREX_D_DECL(lb.x, lb.y, lb.z)), + amrex::IntVect(AMREX_D_DECL(ub.x, ub.y, ub.z))); + + const auto& is_lb = amrex::lbound(is_fluid_arr); + const auto& is_ub = amrex::ubound(is_fluid_arr); + const amrex::Box is_box( + amrex::IntVect(AMREX_D_DECL(is_lb.x, is_lb.y, is_lb.z)), + amrex::IntVect( + AMREX_D_DECL(is_ub.x, is_ub.y, is_ub.z))); + for (int q = 0; q < constants::N_MICRO_STATES; q++) { - const auto& ev = evs[q]; - const amrex::IntVect ivr(iv + evs[bounce_dirs[q]]); - for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - fs[idir] += 2.0 * ev[idir] * f_arrs[nbx](ivr, q) * - is_fluid_arrs[nbx](ivr, 0); + if (bounce_dirs[q] >= 0 && + bounce_dirs[q] < constants::N_MICRO_STATES) { + + const auto& ev = evs[q]; + const amrex::IntVect ivr(iv + evs[bounce_dirs[q]]); + + if (fbox.contains(ivr) && is_box.contains(ivr)) { + for (int idir = 0; idir < AMREX_SPACEDIM; + idir++) { + fs[idir] += 2.0 * ev[idir] * + f_arrs[nbx](ivr, q) * + is_fluid_arrs[nbx](ivr, 0); + } + } } } } @@ -1135,11 +1161,11 @@ void LBM::MakeNewLevelFromCoarse( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } // Make a new level from scratch using provided BoxArray and @@ -1191,11 +1217,11 @@ void LBM::MakeNewLevelFromScratch( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } void LBM::initialize_f(const int lev) @@ -1353,12 +1379,12 @@ void LBM::RemakeLevel( f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - compute_q_corrections(lev); - m_ts_new[lev] = time; m_ts_old[lev] = time - constants::SMALL_NUM; } @@ -1907,11 +1933,11 @@ void LBM::read_checkpoint_file() f_to_macrodata(lev); + compute_q_corrections(lev); + macrodata_to_equilibrium(lev); compute_derived(lev); - - compute_q_corrections(lev); } } diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index f0ddd60..f22bfa9 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -95,7 +95,7 @@ endif() if(MARBLES_DIM EQUAL 3) add_test_r(single_cylinder) add_test_r(channel_cylinder) - # add_test_r(channel_cylinder_amr) + add_test_r(cylinder_turek_2d2) endif() add_test_re(sod) diff --git a/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp b/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp new file mode 100644 index 0000000..3a444e9 --- /dev/null +++ b/Tests/test_files/cylinder_turek_2d2_amr/cylinder_turek_2d2_amr.inp @@ -0,0 +1,65 @@ +max_step = 100000 + +# geometry parameters +geometry.prob_lo = 0.0 0.0 -2.0 +geometry.prob_hi = 660.0 123.0 2.0 +geometry.is_periodic = 0 0 1 + +# timestepping +amr.n_cell = 660 123 4 +amr.max_level = 1 +amr.max_grid_size = 32 +amr.blocking_factor_x = 2 +amr.blocking_factor_y = 1 +amr.blocking_factor_z = 1 +amr.plot_int = 100 +amr.chk_int = 1000 +amr.file_name_digits = 5 + +lbm.bc_lo = 2 1 0 +lbm.bc_hi = 5 1 0 +lbm.dx_outer = 1.0 +lbm.dt_outer = 1.0 + +#nu = u D / Re +#nu = M a D / Re +#nu = M sqrt(gamma R T) D / Re +#nu = 0.1 sqrt(1.4 1 0.03333) 30 / 100 + +lbm.nu = 0.00648 +lbm.alpha = 0.00648 +lbm.initial_temperature = 0.03333 +lbm.adiabatic_exponent = 1.4 +lbm.save_streaming = 0 +lbm.compute_forces = 1 + +lbm.velocity_bc_type = "parabolic" +velocity_bc_parabolic.Mach_ref = 0.1 +velocity_bc_parabolic.normal_dir = 1 +velocity_bc_parabolic.tangential_dir = 0 +velocity_bc_parabolic.initial_density = 1.0 +velocity_bc_parabolic.initial_temperature = 0.03333 +velocity_bc_parabolic.adiabatic_exponent = 1.4 + +lbm.ic_type = "constant" +ic_constant.mach_components = 0.1 0.0 0.0 +ic_constant.density = 1.0 +ic_constant.initial_temperature = 0.03333 +ic_constant.adiabatic_exponent = 1.4 + +eb2.geom_type = "cylinder" +eb2.cylinder_radius = 15.0 +eb2.cylinder_center = 60.0 60.0 0.0 +eb2.cylinder_has_fluid_inside = 0 +eb2.cylinder_height = 1000.0 +eb2.cylinder_direction = 2 + +tagging.refinement_indicators = box +tagging.box.in_box_lo = 30.0 30.0 -2.0 +tagging.box.in_box_hi = 90.0 90.0 2.0 + +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 \ No newline at end of file