Potential fix to reflections at AMR levels

Source/LBM.cpp

@@ -1019,13 +1019,41 @@
                 amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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 (q >= 0 && q < constants::N_MICRO_STATES &&
+                            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]]);
+
+                            // Bounds checking. Debug
+                            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 +1163,11 @@
 
     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 +1219,11 @@
 
     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 +1381,12 @@
 
     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 +1935,11 @@
 
         f_to_macrodata(lev);
 
+        compute_q_corrections(lev);
+
         macrodata_to_equilibrium(lev);
 
         compute_derived(lev);
-
-        compute_q_corrections(lev);
     }
 }
 

Tests/CMakeLists.txt

@@ -96,6 +96,8 @@ 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)
+  add_test_r(cylinder_turek_2d2_amr)
 endif()
 
 add_test_re(sod)

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