Fix the incompressible option. (#138)

* Cannot access I_R if the run is purely incompressible.

* Do not avgDown I_R if pure incompressible.

* Add a 3D flow over bump input.

Exec/RegTests/EB_FlowPastCylinder/input.3d-regt_WallBump

@@ -0,0 +1,70 @@
+#----------------------DOMAIN DEFINITION------------------------
+geometry.is_periodic = 0 0 1               # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     = -0.02 0.0 0.0        # x_lo y_lo (z_lo)
+geometry.prob_hi     =  0.06 0.04 0.02  # x_hi y_hi (z_hi)
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Inflow  NoSlipWallAdiab Interior
+peleLM.hi_bc = Outflow Symmetry Interior
+
+#-------------------------AMR CONTROL----------------------------
+amr.n_cell          = 64 32 16         # Level 0 number of cells in each direction   
+amr.v               = 1                # AMR verbose
+amr.max_level       = 2                # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2          # refinement ratio
+amr.regrid_int      = 5                # how often to regrid
+amr.n_error_buf     = 2 2 2 2          # number of buffer cells in error est
+amr.grid_eff        = 0.7              # what constitutes an efficient grid
+amr.blocking_factor = 16               # block factor in grid generation (min box size)
+amr.max_grid_size   = 64               # max box size
+
+#--------------------------- Problem -------------------------------
+prob.T_mean = 298.0
+prob.P_mean = 101325.0
+prob.type   = ConvectedVortex
+prob.meanFlowMag = 5.0
+prob.meanFlowDir = 1
+
+#------------  INPUTS TO CONSTANT TRANSPORT -----------------
+transport.const_viscosity        = 0.0001 
+transport.const_bulk_viscosity   = 0.0 
+transport.const_conductivity     = 0.0 
+transport.const_diffusivity      = 0.0 
+
+#-------------------------PeleLM CONTROL----------------------------
+peleLM.v = 3
+peleLM.incompressible = 1
+peleLM.rho = 1.17
+peleLM.mu = 0.0001
+
+amr.plot_int = 20
+amr.plot_file = "plt"
+amr.max_step = 20
+amr.dt_shrink = 1.0
+amr.stop_time = 0.02
+amr.cfl = 0.7
+amr.derive_plot_vars = avg_pressure mag_vort
+
+#------------------------- EB SETUP -----------------------------
+eb2.geom_type = cylinder
+eb2.cylinder_radius = 0.006
+eb2.cylinder_direction = 2
+eb2.cylinder_center = 0.0 0.0 0.0
+eb2.cylinder_has_fluid_inside = 0
+#eb2.small_volfrac = 1.0e-4
+
+#--------------------REFINEMENT CONTROL------------------------
+#amr.refinement_indicators = vort_low vort_high
+#amr.vort_low.max_level     = 2
+#amr.vort_low.value_less    = -2000
+#amr.vort_low.field_name    = mag_vort
+#amr.vort_high.max_level     = 3
+#amr.vort_high.value_greater = 2000
+#amr.vort_high.field_name    = mag_vort
+
+#amrex.fpe_trap_invalid = 1
+#amrex.fpe_trap_zero = 1
+#amrex.fpe_trap_overflow = 1

Source/PeleLMEvolve.cpp

@@ -24,7 +24,7 @@ void PeleLM::Evolve() {
       if ( (m_regrid_int > 0) && (m_nstep > 0) && (m_nstep%m_regrid_int == 0) ) {
          if (m_verbose > 0) amrex::Print() << " Regridding...\n";
          // Average down I_R to have proper values in newly uncovered areas
-         averageDownReaction();
+         if (!m_incompressible) averageDownReaction();
          regrid(0, m_cur_time);
          resetMacProjector();
          resetCoveredMask();

Source/PeleLMUtils.cpp

@@ -582,7 +582,7 @@ PeleLM::derive(const std::string &a_name,
           const Box& bx = mfi.growntilebox(nGrow);
           FArrayBox& derfab = (*mf)[mfi];
           FArrayBox const& statefab = (*statemf)[mfi];
-          FArrayBox const& reactfab = ldataR_p->I_R[mfi];
+          FArrayBox const& reactfab = (m_incompressible) ? ldata_p->press[mfi] : ldataR_p->I_R[mfi];
           FArrayBox const& pressfab = ldata_p->press[mfi];
           rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, reactfab, pressfab, geom[lev], a_time, stateBCs, lev);
       }
@@ -634,7 +634,7 @@ PeleLM::deriveComp(const std::string &a_name,
           const Box& bx = mfi.growntilebox(nGrow);
           FArrayBox& derfab = derTemp[mfi];
           FArrayBox const& statefab = (*statemf)[mfi];
-          FArrayBox const& reactfab = ldataR_p->I_R[mfi];
+          FArrayBox const& reactfab = (m_incompressible) ? ldata_p->press[mfi] : ldataR_p->I_R[mfi];
           FArrayBox const& pressfab = ldata_p->press[mfi];
           rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, reactfab, pressfab, geom[lev], a_time, stateBCs, lev);
       }