diff --git a/Docs/sphinx/geometry/EB.rst b/Docs/sphinx/geometry/EB.rst index 60702e1c3..e0b25247f 100644 --- a/Docs/sphinx/geometry/EB.rst +++ b/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 ``ebfill`` function. 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. diff --git a/Exec/RegTests/CMakeLists.txt b/Exec/RegTests/CMakeLists.txt index ab78f9abb..17b5b2085 100644 --- a/Exec/RegTests/CMakeLists.txt +++ b/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) diff --git a/Exec/RegTests/EB-InflowBC/CMakeLists.txt b/Exec/RegTests/EB-InflowBC/CMakeLists.txt new file mode 100644 index 000000000..d7259efeb --- /dev/null +++ b/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) diff --git a/Exec/RegTests/EB-InflowBC/GNUmakefile b/Exec/RegTests/EB-InflowBC/GNUmakefile new file mode 100644 index 000000000..45c4fc26d --- /dev/null +++ b/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 diff --git a/Exec/RegTests/EB-InflowBC/Make.package b/Exec/RegTests/EB-InflowBC/Make.package new file mode 100644 index 000000000..bdc32b025 --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/Make.package @@ -0,0 +1,3 @@ +CEXE_headers += prob.H +CEXE_headers += prob_parm.H +CEXE_sources += prob.cpp diff --git a/Exec/RegTests/EB-InflowBC/eb-inflowbc.inp b/Exec/RegTests/EB-InflowBC/eb-inflowbc.inp new file mode 100644 index 000000000..e0969c400 --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/eb-inflowbc.inp @@ -0,0 +1,81 @@ +# ------------------ 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 diff --git a/Exec/RegTests/EB-InflowBC/example.inp b/Exec/RegTests/EB-InflowBC/example.inp new file mode 100644 index 000000000..e0969c400 --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/example.inp @@ -0,0 +1,81 @@ +# ------------------ 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 diff --git a/Exec/RegTests/EB-InflowBC/prob.H b/Exec/RegTests/EB-InflowBC/prob.H new file mode 100644 index 000000000..94e27e63d --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/prob.H @@ -0,0 +1,137 @@ +#ifndef PROB_H +#define PROB_H + +#include +#include +#include +#include +#include +#include + +#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 const& state, + amrex::GeometryData const& geomdata, + ProbParmDevice const& prob_parm) +{ + // Geometry + const amrex::Real* prob_lo = geomdata.ProbLo(); + const amrex::Real* prob_hi = geomdata.ProbHi(); + const amrex::Real* dx = geomdata.CellSize(); + + 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& /*turb_fluc*/) +{ +} + +void pc_prob_close(); + +struct MyProblemSpecificFunctions : public DefaultProblemSpecificFunctions +{ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + static bool ebfill( + 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(iv[0] + 0.5) * dx[0], + prob_lo[1] + static_cast(iv[1] + 0.5) * dx[1], + prob_lo[2] + static_cast(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 diff --git a/Exec/RegTests/EB-InflowBC/prob.cpp b/Exec/RegTests/EB-InflowBC/prob.cpp new file mode 100644 index 000000000..4fb96e8c4 --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/prob.cpp @@ -0,0 +1,70 @@ +#include "prob.H" + +void +pc_prob_close() +{ +} + +extern "C" { +void +amrex_probinit( + const int* init, + const int* name, + const int* namelen, + const amrex_real* problo, + const amrex_real* probhi) +{ + // Parse params + amrex::ParmParse pp("prob"); + pp.query("p_l", PeleC::h_prob_parm_device->p_l); + pp.query("T_l", PeleC::h_prob_parm_device->T_l); + pp.query("p_r", PeleC::h_prob_parm_device->p_r); + pp.query("T_r", PeleC::h_prob_parm_device->T_r); + pp.query("U_r", PeleC::h_prob_parm_device->U_r); + pp.get("left_gas", PeleC::prob_parm_host->gasL); + pp.get("right_gas", PeleC::prob_parm_host->gasR); + + if (PeleC::prob_parm_host->gasL == "N2") { + PeleC::h_prob_parm_device->left_gas_id = N2_ID; + } else { + PeleC::h_prob_parm_device->left_gas_id = O2_ID; + } + if (PeleC::prob_parm_host->gasR == "N2") { + PeleC::h_prob_parm_device->right_gas_id = N2_ID; + } else { + PeleC::h_prob_parm_device->right_gas_id = O2_ID; + } + amrex::Real e_l, e_r, cs, cp; + amrex::Real massfrac_l[NUM_SPECIES] = {0.0}; + amrex::Real massfrac_r[NUM_SPECIES] = {0.0}; + massfrac_l[PeleC::h_prob_parm_device->left_gas_id] = 1.0; + massfrac_r[PeleC::h_prob_parm_device->right_gas_id] = 1.0; + std::cout << "left " << PeleC::h_prob_parm_device->left_gas_id << std::endl; + std::cout << "right " << PeleC::h_prob_parm_device->right_gas_id << std::endl; + + auto eos = pele::physics::PhysicsType::eos(); + eos.PYT2RE( + PeleC::h_prob_parm_device->p_l, massfrac_l, PeleC::h_prob_parm_device->T_l, + PeleC::h_prob_parm_device->rho_l, e_l); + eos.PYT2RE( + PeleC::h_prob_parm_device->p_r, massfrac_r, PeleC::h_prob_parm_device->T_r, + PeleC::h_prob_parm_device->rho_r, e_r); + PeleC::h_prob_parm_device->rhoe_l = PeleC::h_prob_parm_device->rho_l * e_l; + PeleC::h_prob_parm_device->rhoe_r = PeleC::h_prob_parm_device->rho_r * e_r; +} +} + +void +PeleC::problem_post_timestep() +{ +} + +void +PeleC::problem_post_init() +{ +} + +void +PeleC::problem_post_restart() +{ +} diff --git a/Exec/RegTests/EB-InflowBC/prob_parm.H b/Exec/RegTests/EB-InflowBC/prob_parm.H new file mode 100644 index 000000000..268f61c4e --- /dev/null +++ b/Exec/RegTests/EB-InflowBC/prob_parm.H @@ -0,0 +1,28 @@ +#ifndef _PROB_PARM_H_ +#define _PROB_PARM_H_ + +#include +#include + +struct ProbParmDevice +{ + amrex::Real p_l = 1.0; // left pressure (erg/cc) + amrex::Real rho_l = 1.0; // left density (g/cc) + amrex::Real rhoe_l = 0.0; + amrex::Real T_l = 1.0; + amrex::Real p_r = 0.1; // right pressure (erg/cc) + amrex::Real rho_r = 0.125; // right density (g/cc) + amrex::Real rhoe_r = 0.0; + amrex::Real T_r = 1.0; + amrex::Real U_r = 1.0; + int left_gas_id = N2_ID; + int right_gas_id = O2_ID; +}; + +struct ProbParmHost +{ + std::string gasL = "N2"; + std::string gasR = "O2"; + ProbParmHost() {} +}; +#endif diff --git a/Source/Diffusion.cpp b/Source/Diffusion.cpp index bd6392818..1708d0fd0 100644 --- a/Source/Diffusion.cpp +++ b/Source/Diffusion.cpp @@ -437,8 +437,8 @@ PeleC::getMOLSrcTerm( amrex::Real* d_eb_flux_thdlocal = (nFlux > 0 ? eb_flux_thdlocal.dataPtr() : nullptr); pc_compute_hyp_mol_flux_eb( - cbox, qar, dx, d_sv_eb_bndry_geom, Ncut, d_eb_flux_thdlocal, - nFlux); + geom, cbox, qar, qauxar, dx, use_laxf_flux, vfrac.array(mfi), + d_sv_eb_bndry_geom, Ncut, d_eb_flux_thdlocal, nFlux); } } } diff --git a/Source/MOL.H b/Source/MOL.H index 875a755bc..96fe64ecd 100644 --- a/Source/MOL.H +++ b/Source/MOL.H @@ -130,9 +130,13 @@ void pc_compute_hyp_mol_flux( const amrex::Array4& flags); void pc_compute_hyp_mol_flux_eb( + amrex::Geometry const& geom, const amrex::Box& cbox, const amrex::Array4& q, + const amrex::Array4& qaux, const amrex::GpuArray& dx, + const bool use_laxf_flux, + const amrex::Array4& vfrac, const EBBndryGeom* ebg, const int Nebg, amrex::Real* ebflux, diff --git a/Source/MOL.cpp b/Source/MOL.cpp index eecbad085..9def07746 100644 --- a/Source/MOL.cpp +++ b/Source/MOL.cpp @@ -1,5 +1,6 @@ #include "MOL.H" #include "Godunov.H" +#include "prob.H" void pc_compute_hyp_mol_flux( @@ -205,9 +206,13 @@ pc_compute_hyp_mol_flux( void pc_compute_hyp_mol_flux_eb( + amrex::Geometry const& geom, const amrex::Box& cbox, const amrex::Array4& q, + const amrex::Array4& qaux, const amrex::GpuArray& dx, + const bool use_laxf_flux, + const amrex::Array4& vfrac, const EBBndryGeom* ebg, const int /*Nebg*/, amrex::Real* ebflux, @@ -216,8 +221,19 @@ pc_compute_hyp_mol_flux_eb( const int nextra = 0; + 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 int R_ADV = 5; + const int R_Y = R_ADV + NUM_ADV; + const int bc_test_val = 1; + const amrex::Real full_area = AMREX_D_PICK(1.0, dx[0], dx[0] * dx[1]); const amrex::Box bxg = amrex::grow(cbox, nextra - 1); + const auto geomdata = geom.data(); + ProbParmDevice const* prob_parm = PeleC::d_prob_parm_device; amrex::ParallelFor(nebflux, [=] AMREX_GPU_DEVICE(int L) { const amrex::IntVect& iv = ebg[L].iv; @@ -231,9 +247,137 @@ pc_compute_hyp_mol_flux_eb( } amrex::Real flux_tmp[NVAR] = {0.0}; - AMREX_D_TERM(flux_tmp[UMX] = -q(iv, QPRES) * ebnorm[0]; - , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm[1]; - , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm[2];) + + auto eos = pele::physics::PhysicsType::eos(); + amrex::Real qtempl[5 + NUM_SPECIES] = {0.0}; + amrex::Real spl[NUM_SPECIES] = {0.0}; + qtempl[R_UN] = -(AMREX_D_TERM( + q(iv, QU) * ebnorm[0], +q(iv, QV) * ebnorm[1], +q(iv, QW) * ebnorm[2])); + qtempl[R_UT1] = 0.0; + qtempl[R_UT2] = 0.0; + qtempl[R_P] = q(iv, QPRES); + qtempl[R_RHO] = q(iv, QRHO); + for (int n = 0; n < NUM_SPECIES; n++) { + qtempl[R_Y + n] = q(iv, QFS + n); + } + amrex::Real cavg = qaux(iv, QC); + + // Flip the velocity about the normal for the right state - will use + // left state for remainder of right state + amrex::Real qtempr[5 + NUM_SPECIES] = {0.0}; + qtempr[R_UN] = -1.0 * qtempl[R_UN]; + + amrex::Real eos_state_rho = qtempl[R_RHO]; + amrex::Real eos_state_p = qtempl[R_P]; + for (int n = 0; n < NUM_SPECIES; n++) { + spl[n] = qtempl[R_Y + n]; + } + amrex::Real eos_state_T; + eos.RYP2T(eos_state_rho, spl, eos_state_p, eos_state_T); + amrex::Real eos_state_e; + eos.RTY2E(eos_state_rho, eos_state_T, spl, eos_state_e); + amrex::Real rhoe_l = eos_state_rho * eos_state_e; + amrex::Real gamc_l; + eos.RTY2G(eos_state_rho, eos_state_T, spl, gamc_l); + + // Copy left state to right (default), except normal velocity which has + // already been flipped + qtempr[R_RHO] = qtempl[R_RHO]; + qtempr[R_UT1] = qtempl[R_UT1]; + qtempr[R_UT2] = qtempl[R_UT2]; + qtempr[R_P] = qtempl[R_P]; + amrex::Real spr[NUM_SPECIES] = {0.0}; + for (int n = 0; n < NUM_SPECIES; n++) { + spr[n] = spl[n]; + } + amrex::Real rhoe_r = rhoe_l; + amrex::Real gamc_r = gamc_l; + + const bool do_ebfill = ProblemSpecificFunctions::ebfill( + geomdata, vfrac(iv), iv, AMREX_D_DECL(ebnorm[0], ebnorm[1], ebnorm[2]), + qtempl, spl, rhoe_l, gamc_l, qtempr, spr, rhoe_r, gamc_r, prob_parm); + + if (do_ebfill) { + + amrex::Real ustar = 0.0; + if (!use_laxf_flux) { + amrex::Real qint_iu = 0.0, tmp1 = 0.0, tmp2 = 0.0, tmp3 = 0.0, + tmp4 = 0.0; + riemann( + qtempl[R_RHO], qtempl[R_UN], qtempl[R_UT1], qtempl[R_UT2], + qtempl[R_P], spl, qtempr[R_RHO], qtempr[R_UN], qtempr[R_UT1], + qtempr[R_UT2], qtempr[R_P], spr, bc_test_val, cavg, ustar, + flux_tmp[URHO], &flux_tmp[UFS], flux_tmp[UMX], flux_tmp[UMY], + flux_tmp[UMZ], flux_tmp[UEDEN], flux_tmp[UEINT], qint_iu, tmp1, + tmp2, tmp3, tmp4); +#if NUM_ADV > 0 + for (int n = 0; n < NUM_ADV; n++) { + pc_cmpflx_passive( + ustar, flux_tmp[URHO], qtempl[R_ADV + n], qtempr[R_ADV + n], + flux_tmp[UFA + n]); + } +#endif +#if NUM_AUX > 0 + const int R_AUX = R_Y + NUM_SPECIES; + for (int n = 0; n < NUM_AUX; n++) { + pc_cmpflx_passive( + ustar, flux_tmp[URHO], qtempl[R_AUX + n], qtempr[R_AUX + n], + flux_tmp[UFX + n]); + } +#endif +#if NUM_LIN > 0 + const int R_LIN = R_Y + NUM_SPECIES + NUM_AUX; + for (int n = 0; n < NUM_LIN; n++) { + pc_cmpflx_passive( + ustar, qint_iu, qtempl[R_LIN + n], qtempr[R_LIN + n], + flux_tmp[ULIN + n]); + } +#endif + } else { + amrex::Real maxeigval = 0.0; + laxfriedrich_flux( + qtempl[R_RHO], qtempl[R_UN], qtempl[R_UT1], qtempl[R_UT2], + qtempl[R_P], spl, qtempr[R_RHO], qtempr[R_UN], qtempr[R_UT1], + qtempr[R_UT2], qtempr[R_P], spr, bc_test_val, cavg, ustar, + maxeigval, flux_tmp[URHO], &flux_tmp[UFS], flux_tmp[UMX], + flux_tmp[UMY], flux_tmp[UMZ], flux_tmp[UEDEN], flux_tmp[UEINT]); +#if NUM_ADV > 0 + for (int n = 0; n < NUM_ADV; n++) { + pc_lax_cmpflx_passive( + qtempl[R_UN], qtempr[R_UN], qtempl[R_RHO], qtempr[R_RHO], + qtempl[R_ADV + n], qtempr[R_ADV + n], maxeigval, + flux_tmp[UFA + n]); + } +#endif +#if NUM_AUX > 0 + const int R_AUX = R_Y + NUM_SPECIES; + for (int n = 0; n < NUM_AUX; n++) { + pc_lax_cmpflx_passive( + qtempl[R_UN], qtempr[R_UN], qtempl[R_RHO], qtempr[R_RHO], + qtempl[R_AUX + n], qtempr[R_AUX + n], maxeigval, + flux_tmp[UFX + n]); + } +#endif +#if NUM_LIN > 0 + const int R_LIN = R_Y + NUM_SPECIES + NUM_AUX; + for (int n = 0; n < NUM_LIN; n++) { + pc_lax_cmpflx_passive( + qtempl[R_UN], qtempr[R_UN], 1., 1., qtempl[R_LIN + n], + qtempr[R_LIN + n], maxeigval, flux_tmp[ULIN + n]); + } +#endif + } + + const amrex::Real tmp_flx_umx = flux_tmp[UMX]; + AMREX_D_TERM(flux_tmp[UMX] = -tmp_flx_umx * ebnorm[0]; + , flux_tmp[UMY] = -tmp_flx_umx * ebnorm[1]; + , flux_tmp[UMZ] = -tmp_flx_umx * ebnorm[2];) + + } else { + AMREX_D_TERM(flux_tmp[UMX] = -q(iv, QPRES) * ebnorm[0]; + , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm[1]; + , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm[2];) + } // Copy result into ebflux vector. Being a bit chicken here and only // copy values where ebg % iv is within box diff --git a/Source/ProblemSpecificFunctions.H b/Source/ProblemSpecificFunctions.H index b7582ee15..e6bfb8c2b 100644 --- a/Source/ProblemSpecificFunctions.H +++ b/Source/ProblemSpecificFunctions.H @@ -98,6 +98,29 @@ struct DefaultProblemSpecificFunctions Twall(i, j, k) = Twall_in; } + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + static bool ebfill( + 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*/) + { + return false; + } + static void problem_modify_ext_sources( amrex::Real /*time*/, amrex::Real /*dt*/, diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 70cf22c5c..18959c4ff 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -249,6 +249,7 @@ add_test_rv(eb-c11 EB-C11) add_test_rv(eb-c12 EB-C12) # add_test_r(eb-c14 EB-C14) # disable due to FPE in ghost cells add_test_r(eb-converging-nozzle EB-ConvergingNozzle) +add_test_r(eb-inflowbc EB-InflowBC) if(PELE_DIM GREATER 2) add_test_r(shock-cylinder Sod) # can run in 2D but needs input file change endif()