Implement the state prescription at EB for the Riemann solver (#774)

Docs/sphinx/geometry/EB.rst

@@ -198,3 +198,10 @@ carried out in the covered region, and invalid operations or NaNs may be detecte
 For debugging purposes, one may specify ``pelec.eb_zero_body_state = true``, in which case all state variables in the covered region will be set to zero.
 This will lead to NaNs when primitive variables are computed (dividing by density), but these should not propagate into fluid cells. The ``EB-C3`` RegTest
 is run with this option to ensure that covered cells do not affect fluid cells.
+
+Problem specific inflow conditions on an EB
+-------------------------------------------
+
+It is possible for the user to define problem specific conditions on an EB surface. This is done by defining an ``problem_eb_state`` function and then including `pelec.eb_problem_state = 1` in the input file. An example of this is found in the ``EB-InflowBC`` case.
+
+.. warning:: This is a beta feature. Currently this will only affect the calculation of the hydrodynamic fluxes so this works best for advection dominated EB conditions.

Exec/RegTests/CMakeLists.txt

@@ -29,6 +29,7 @@ add_subdirectory(EB-C10)
 add_subdirectory(EB-C11)
 add_subdirectory(EB-C12)
 add_subdirectory(EB-ConvergingNozzle)
+add_subdirectory(EB-InflowBC)
 add_subdirectory(Soot-ZeroD)
 if(NOT PELE_EXCLUDE_BUILD_IN_CI)
   add_subdirectory(Soot-Flame)

Exec/RegTests/EB-InflowBC/CMakeLists.txt

@@ -0,0 +1,7 @@
+set(PELE_PHYSICS_EOS_MODEL Fuego)
+set(PELE_PHYSICS_CHEMISTRY_MODEL air)
+set(PELE_PHYSICS_TRANSPORT_MODEL Simple)
+set(PELE_PHYSICS_ENABLE_SOOT OFF)
+set(PELE_PHYSICS_ENABLE_SPRAY OFF)
+set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
+include(BuildExeAndLib)

Exec/RegTests/EB-InflowBC/GNUmakefile

@@ -0,0 +1,38 @@
+# AMReX
+DIM = 3
+COMP = gnu
+PRECISION = DOUBLE
+
+# Profiling
+PROFILE = FALSE
+TINY_PROFILE = FALSE
+COMM_PROFILE = FALSE
+TRACE_PROFILE = FALSE
+MEM_PROFILE = FALSE
+USE_GPROF = FALSE
+
+# Performance
+USE_MPI = TRUE
+USE_OMP = FALSE
+USE_CUDA = FALSE
+USE_HIP = FALSE
+USE_SYCL = FALSE
+
+# Debugging
+DEBUG = FALSE
+FSANITIZER = FALSE
+THREAD_SANITIZER = FALSE
+
+# PeleC
+PELE_CVODE_FORCE_YCORDER = FALSE
+PELE_USE_MAGMA = FALSE
+PELE_COMPILE_AJACOBIAN = FALSE
+Eos_Model := Fuego
+Transport_Model := Simple
+Chemistry_Model := air
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+PELE_HOME := ../../..
+include $(PELE_HOME)/Exec/Make.PeleC

Exec/RegTests/EB-InflowBC/Make.package

@@ -0,0 +1,3 @@
+CEXE_headers += prob.H
+CEXE_headers += prob_parm.H
+CEXE_sources += prob.cpp

Exec/RegTests/EB-InflowBC/eb-inflowbc.inp

@@ -0,0 +1,82 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000
+#stop_time =  1.959e-6 #final time is 0.2*L*sqrt(rhoL/pL)
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0  0  1
+geometry.coord_sys   = 0  # 0 => cart, 1 => RZ  2=>spherical
+
+geometry.prob_lo     = -1.0  -1.0  -1.0
+geometry.prob_hi     =  1.0   1.0   1.0
+amr.n_cell           =  32    32    32
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+pelec.lo_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+pelec.hi_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 1
+pelec.diffuse_vel = 1
+pelec.diffuse_temp = 1
+pelec.diffuse_spec = 1
+pelec.do_react = 0
+pelec.diffuse_enth = 0
+pelec.chem_integrator = "ReactorRK64"
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.2     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.8     # scale back initial timestep
+pelec.change_max     = 1.05    # scale back initial timestep
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in Castro.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+amr.data_log         = datlog
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 1       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 16
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000        # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt        # root name of plotfile
+amr.plot_int        = 100       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+prob.p_l = 1e6
+prob.T_l = 300
+prob.p_r = 1e6
+prob.T_r = 300
+prob.U_r = 1000
+prob.left_gas = O2
+prob.right_gas = N2
+
+# Problem setup
+pelec.eb_boundary_T = 300.
+pelec.eb_isothermal = 0
+
+# TAGGING
+tagging.denerr = 1e20
+tagging.dengrad = 4e-5
+tagging.max_denerr_lev = 3
+tagging.max_dengrad_lev = 3
+
+eb2.geom_type = plane
+eb2.plane_point = 0.0 0.0 0.0
+eb2.plane_normal = 1.0 0.0 0.0
+ebd.boundary_grad_stencil_type = 0
+pelec.eb_problem_state = 1

Exec/RegTests/EB-InflowBC/example.inp

@@ -0,0 +1,82 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000
+#stop_time =  1.959e-6 #final time is 0.2*L*sqrt(rhoL/pL)
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0  0  1
+geometry.coord_sys   = 0  # 0 => cart, 1 => RZ  2=>spherical
+
+geometry.prob_lo     = -1.0  -1.0  -1.0
+geometry.prob_hi     =  1.0   1.0   1.0
+amr.n_cell           =  32    32    32
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+pelec.lo_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+pelec.hi_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 1
+pelec.diffuse_vel = 1
+pelec.diffuse_temp = 1
+pelec.diffuse_spec = 1
+pelec.do_react = 0
+pelec.diffuse_enth = 0
+pelec.chem_integrator = "ReactorRK64"
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.2     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.8     # scale back initial timestep
+pelec.change_max     = 1.05    # scale back initial timestep
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in Castro.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+amr.data_log         = datlog
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 1       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 16
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000        # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt        # root name of plotfile
+amr.plot_int        = 100       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+prob.p_l = 1e6
+prob.T_l = 300
+prob.p_r = 1e6
+prob.T_r = 300
+prob.U_r = 1000
+prob.left_gas = O2
+prob.right_gas = N2
+
+# Problem setup
+pelec.eb_boundary_T = 300.
+pelec.eb_isothermal = 0
+
+# TAGGING
+tagging.denerr = 1e20
+tagging.dengrad = 4e-5
+tagging.max_denerr_lev = 3
+tagging.max_dengrad_lev = 3
+
+eb2.geom_type = plane
+eb2.plane_point = 0.0 0.0 0.0
+eb2.plane_normal = 1.0 0.0 0.0
+ebd.boundary_grad_stencil_type = 0
+pelec.eb_problem_state = 1

Exec/RegTests/EB-InflowBC/prob.H

@@ -0,0 +1,132 @@
+#ifndef PROB_H
+#define PROB_H
+
+#include <AMReX_Print.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_REAL.H>
+#include <AMReX_GpuMemory.H>
+
+#include "mechanism.H"
+
+#include "PeleC.H"
+#include "IndexDefines.H"
+#include "PelePhysics.H"
+#include "Tagging.H"
+#include "ProblemSpecificFunctions.H"
+#include "prob_parm.H"
+#include "Utilities.H"
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+pc_initdata(
+  int i,
+  int j,
+  int k,
+  amrex::Array4<amrex::Real> const& state,
+  amrex::GeometryData const& /*geomdata*/,
+  ProbParmDevice const& prob_parm)
+{
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, UFS + n) = 0.0;
+  }
+
+  // Set the states
+  state(i, j, k, URHO) = prob_parm.rho_l;
+  state(i, j, k, UMX) = 0.0;
+  state(i, j, k, UMY) = 0.0;
+  state(i, j, k, UMZ) = 0.0;
+  state(i, j, k, UEDEN) = prob_parm.rhoe_l;
+  state(i, j, k, UEINT) = prob_parm.rhoe_l;
+  state(i, j, k, UTEMP) = prob_parm.T_l;
+  state(i, j, k, UFS + prob_parm.left_gas_id) = state(i, j, k, URHO);
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+bcnormal(
+  const amrex::Real* /*x[AMREX_SPACEDIM]*/,
+  const amrex::Real* /*s_int[NVAR]*/,
+  amrex::Real* /*s_ext[NVAR]*/,
+  const int /*idir*/,
+  const int /*sgn*/,
+  const amrex::Real /*time*/,
+  amrex::GeometryData const& /*geomdata*/,
+  ProbParmDevice const& /*prob_parm*/,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& /*turb_fluc*/)
+{
+}
+
+void pc_prob_close();
+
+struct MyProblemSpecificFunctions : public DefaultProblemSpecificFunctions
+{
+  AMREX_GPU_DEVICE
+  AMREX_FORCE_INLINE
+  static bool problem_eb_state(
+    amrex::GeometryData const& geomdata,
+    const amrex::Real /*vf*/,
+    const amrex::IntVect iv,
+    AMREX_D_DECL(
+      const amrex::Real /*nx*/,
+      const amrex::Real /*ny*/,
+      const amrex::Real /*nz*/),
+    const amrex::Real ql[5 + NUM_SPECIES],
+    const amrex::Real* /*spl[NUM_SPECIES]*/,
+    const amrex::Real /*rhoe_l*/,
+    const amrex::Real /*gamc_l*/,
+    amrex::Real qr[5 + NUM_SPECIES],
+    amrex::Real spr[NUM_SPECIES],
+    amrex::Real& rhoe_r,
+    amrex::Real& gamc_r,
+    ProbParmDevice const* prob_parm)
+  {
+    bool do_ebfill = false;
+    const amrex::Real* prob_lo = geomdata.ProbLo();
+    const amrex::Real* dx = geomdata.CellSize();
+    const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL(
+      prob_lo[0] + static_cast<amrex::Real>(iv[0] + 0.5) * dx[0],
+      prob_lo[1] + static_cast<amrex::Real>(iv[1] + 0.5) * dx[1],
+      prob_lo[2] + static_cast<amrex::Real>(iv[2] + 0.5) * dx[2])};
+
+    const int R_RHO = 0;
+    const int R_UN = 1;
+    const int R_UT1 = 2;
+    const int R_UT2 = 3;
+    const int R_P = 4;
+
+    const amrex::Real max_r = 0.1;
+    const amrex::Real radius =
+      std::sqrt(AMREX_D_TERM(, x[1] * x[1], +x[2] * x[2]));
+    if (radius < max_r) {
+      qr[R_RHO] = prob_parm->rho_r;
+      qr[R_UN] = -prob_parm->U_r;
+      qr[R_UT1] = ql[R_UT1];
+      qr[R_UT2] = ql[R_UT2];
+      for (int n = 0; n < NUM_SPECIES; n++) {
+        spr[n] = 0.0;
+      }
+      spr[prob_parm->right_gas_id] = 1.0;
+
+      amrex::Real eos_state_rho = qr[R_RHO];
+      amrex::Real eos_state_p = qr[R_P];
+      amrex::Real eos_state_e;
+      auto eos = pele::physics::PhysicsType::eos();
+      eos.RYP2E(eos_state_rho, spr, eos_state_p, eos_state_e);
+      rhoe_r = eos_state_rho * eos_state_e;
+      amrex::Real eos_state_T;
+      eos.RYP2T(eos_state_rho, spr, eos_state_p, eos_state_T);
+      eos.TY2G(eos_state_T, spr, gamc_r);
+      do_ebfill = true;
+    }
+
+    return do_ebfill;
+  }
+};
+
+using ProblemSpecificFunctions = MyProblemSpecificFunctions;
+
+#endif