From acc0f7e3ae79fa4b12a93058b92aed63b3e9d9df Mon Sep 17 00:00:00 2001 From: Bruce Perry <53018946+baperry2@users.noreply.github.com> Date: Mon, 23 Oct 2023 20:18:14 -0600 Subject: [PATCH] Some improvements and testing for the EB covered state (#704) --- Docs/sphinx/geometry/EB.rst | 16 ++++++ Exec/RegTests/EB-C11/eb-c11.inp | 7 +-- Exec/RegTests/EB-C11/example.inp | 7 +-- Exec/RegTests/EB-C3/eb-c3.inp | 17 +++--- Exec/RegTests/EB-C4-5/eb-c4.inp | 8 +-- Source/EB.H | 6 ++ Source/EB.cpp | 56 +++++++++++++++++-- Source/InitEB.cpp | 39 +++++-------- Source/Params/_cpp_parameters | 3 + Source/Params/param_includes/pelec_defaults.H | 1 + Source/Params/param_includes/pelec_params.H | 1 + Source/Params/param_includes/pelec_queries.H | 1 + Source/PeleC.H | 2 - Source/PeleC.cpp | 9 +++ Tests/CMakeLists.txt | 10 +++- 15 files changed, 130 insertions(+), 53 deletions(-) diff --git a/Docs/sphinx/geometry/EB.rst b/Docs/sphinx/geometry/EB.rst index 5cd26c2bf..60702e1c3 100644 --- a/Docs/sphinx/geometry/EB.rst +++ b/Docs/sphinx/geometry/EB.rst @@ -162,6 +162,12 @@ Applying boundary and face stencils ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When processing geometry cells, the cached datastructures can be applied efficiently, for example, to interpolate fluxes from face centers to face centroids in cut cells. +Two stencil types are available for computing the gradients at the EB for diffusive fluxes, and they are selected as follows: + +* ``ebd.boundary_grad_stencil_type = 0``: Quadratic stencil (default). On poorly resolved geometries, this stencil may reach into covered cells, in which case the simulation will + fail with a warning. See `Johansen and Collela `_ for further details. + +* ``ebd.boundary_grad_stencil_type = 1``: Least-squares stencil. See `Anderson and Bonhaus `_ for further details. .. include:: /geometry/geometry_init.rst @@ -182,3 +188,13 @@ These are the input options to read from a saved EB geometry: eb2.geom_type="chkfile" eb2.chkfile="chk_geom" # optional, defaults to "chk_geom" eb2.max_grid_size=32 # optional, defaults to 64, must match the max_grid_size used to generate the EB in the first place + +Setting the Covered State +------------------------- + +By default, the state for all cells completely covered by the EB is set to the value of the initial condition of the first valid fluid cell on the base grid. +The values in the covered region do not affect values in the fluid region, but should still have valid values because some basic operations are still +carried out in the covered region, and invalid operations or NaNs may be detected during these operations if the specified values are not valid. +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. diff --git a/Exec/RegTests/EB-C11/eb-c11.inp b/Exec/RegTests/EB-C11/eb-c11.inp index 083b3d810..d5804cbc5 100644 --- a/Exec/RegTests/EB-C11/eb-c11.inp +++ b/Exec/RegTests/EB-C11/eb-c11.inp @@ -74,8 +74,7 @@ prob.right_gas = HE pelec.eb_boundary_T = 300. pelec.eb_isothermal = 1 -eb2.geom_type = box -eb2.box_lo = 1.0 0.0 0.0 -eb2.box_hi = 1.4 0.1 0.1 -eb2.box_has_fluid_inside = 0 +eb2.geom_type = plane +eb2.plane_point = 1.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-C11/example.inp b/Exec/RegTests/EB-C11/example.inp index 083b3d810..d5804cbc5 100644 --- a/Exec/RegTests/EB-C11/example.inp +++ b/Exec/RegTests/EB-C11/example.inp @@ -74,8 +74,7 @@ prob.right_gas = HE pelec.eb_boundary_T = 300. pelec.eb_isothermal = 1 -eb2.geom_type = box -eb2.box_lo = 1.0 0.0 0.0 -eb2.box_hi = 1.4 0.1 0.1 -eb2.box_has_fluid_inside = 0 +eb2.geom_type = plane +eb2.plane_point = 1.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-C3/eb-c3.inp b/Exec/RegTests/EB-C3/eb-c3.inp index 78c6a089f..7c45baee9 100644 --- a/Exec/RegTests/EB-C3/eb-c3.inp +++ b/Exec/RegTests/EB-C3/eb-c3.inp @@ -5,7 +5,7 @@ stop_time = 3e-5 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 0 0 0 geometry.coord_sys = 0 # 0 => cart -geometry.prob_lo = 0 0 0 +geometry.prob_lo = 0 0 0 geometry.prob_lo = 0.0 -0.05 -0.05 geometry.prob_hi = 0.1 0.05 0.05 @@ -28,10 +28,10 @@ pelec.do_mol = 1 pelec.do_react = 1 pelec.chem_integrator = "ReactorRK64" pelec.allow_negative_energy = 0 -pelec.diffuse_temp = 0 -pelec.diffuse_vel = 0 -pelec.diffuse_spec = 0 -pelec.diffuse_enth = 0 +pelec.diffuse_temp = 1 +pelec.diffuse_vel = 1 +pelec.diffuse_spec = 1 +pelec.diffuse_enth = 1 # TIME STEP CONTROL pelec.dt_cutoff = 5.e-20 # level 0 timestep below which we halt @@ -63,13 +63,14 @@ amr.plot_file = plt amr.plot_int = 100 amr.derive_plot_vars = density xmom ymom zmom rho_E rho_e Temp rho_H2 rho_O2 rho_H2O rho_H rho_O rho_OH rho_HO2 rho_H2O2 rho_N2 rho_omega_H2 rho_omega_O2 rho_omega_H2O rho_omega_H rho_omega_O rho_omega_OH rho_omega_HO2 rho_omega_H2O2 rho_omega_N2 pressure kineng enstrophy soundspeed MachNumber entropy magvort divu eint_E eint_e logden x_velocity y_velocity z_velocity magvel radvel magmom vfrac cp cv viscosity bulk_viscosity conductivity D(H2) D(O2) D(H2O) D(H) D(O) D(OH) D(HO2) D(H2O2) D(N2) -#extruded triangles lets the user create a maximum of 5 triangles +#extruded triangles lets the user create a maximum of 5 triangles #in 2D that will be extruded in the z direction -#make sure the coordinates are in anti-clockwise direction +#make sure the coordinates are in anti-clockwise direction #eb2.geom_type = "all_regular" eb2.geom_type = "ICE_PistonBowl" -ebd.boundary_grad_stencil_type = 0 +ebd.boundary_grad_stencil_type = 1 +pelec.eb_zero_body_state = 1 #fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM PARAMETERS diff --git a/Exec/RegTests/EB-C4-5/eb-c4.inp b/Exec/RegTests/EB-C4-5/eb-c4.inp index 671b36312..2273bf916 100644 --- a/Exec/RegTests/EB-C4-5/eb-c4.inp +++ b/Exec/RegTests/EB-C4-5/eb-c4.inp @@ -6,7 +6,7 @@ stop_time = 3e-3 # PROBLEM SIZE & GEOMETRY geometry.is_periodic = 0 0 0 geometry.coord_sys = 0 # 0 => cart -geometry.prob_lo = 0 0 0 +geometry.prob_lo = 0 0 0 geometry.prob_lo = 0.0 -0.05 -0.05 geometry.prob_hi = 0.1 0.05 0.05 @@ -68,13 +68,13 @@ amr.plot_int = 100 #amr.plot_int = 1 amr.derive_plot_vars=ALL -#extruded triangles lets the user create a maximum of 5 triangles +#extruded triangles lets the user create a maximum of 5 triangles #in 2D that will be extruded in the z direction -#make sure the coordinates are in anti-clockwise direction +#make sure the coordinates are in anti-clockwise direction #eb2.geom_type = "all_regular" eb2.geom_type = "ICE_PistonBowl" -ebd.boundary_grad_stencil_type = 0 +ebd.boundary_grad_stencil_type = 1 #fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM PARAMETERS diff --git a/Source/EB.H b/Source/EB.H index 08e082de3..c7a9768bc 100644 --- a/Source/EB.H +++ b/Source/EB.H @@ -222,6 +222,12 @@ void pc_fill_bndry_grad_stencil_ls( const amrex::Array4& flags, EBBndrySten* /*grad_stencil*/); +void pc_check_bndry_grad_stencil( + const amrex::Box& /*bx*/, + const int /*Nsten*/, + const amrex::Array4& flags, + const EBBndrySten* /*grad_stencil*/); + void pc_fill_bndry_grad_stencil_quadratic( const amrex::Box& /*bx*/, const amrex::Real /*dx*/, diff --git a/Source/EB.cpp b/Source/EB.cpp index b689a4fb9..3c78ae3af 100644 --- a/Source/EB.cpp +++ b/Source/EB.cpp @@ -389,6 +389,44 @@ pc_fill_bndry_grad_stencil_ls( }); } +/* Because the quadratic EB boundary flux stencil can unintentionally access + covered cells, here we check whether any of the accessed cells are covered + and raise an error if they are */ +void +pc_check_bndry_grad_stencil( + const amrex::Box& bx, + const int Nsten, + const amrex::Array4& flags, + const EBBndrySten* sten) +{ + amrex::ParallelFor(Nsten, [=] AMREX_GPU_DEVICE(int L) { + const amrex::IntVect iv = sten[L].iv; + if (bx.contains(iv)) { + for (int ii = 0; ii < 3; ii++) { + for (int jj = 0; jj < 3; jj++) { +#if AMREX_SPACEDIM > 2 + for (int kk = 0; kk < 3; kk++) { +#endif + const amrex::IntVect ivp = amrex::IntVect(AMREX_D_DECL( + sten[L].iv_base[0] + ii, sten[L].iv_base[1] + jj, + sten[L].iv_base[2] + kk)); + if ( + flags(ivp).isCovered() && + std::abs(sten[L].val AMREX_D_TERM([ii], [jj], [kk])) > 1e-14) { + amrex::Abort( + "EB Boundary Stencil for noslip/isothermal EB accesses a " + "covered cell. Try again with alternate option for" + " ebd.boundary_grad_stencil_type or refine your base grid."); + } +#if AMREX_SPACEDIM > 2 + } +#endif + } + } + } + }); +} + void pc_fill_flux_interp_stencil( const amrex::Box& bx, @@ -455,8 +493,9 @@ pc_apply_face_stencil( #endif const amrex::IntVect ivd = amrex::IntVect{ AMREX_D_DECL(iv[0], iv[1] - 1 + t0, iv[2] - 1 + t1)}; + amrex::Real stenval = sten[L].val AMREX_D_TERM(, [t0], [t1]); d_newval[L] += - sten[L].val AMREX_D_TERM(, [t0], [t1]) * vout(ivd, n); + (std::abs(stenval) > 1e-14) ? stenval * vout(ivd, n) : 0.0; #if AMREX_SPACEDIM > 2 } #endif @@ -468,8 +507,9 @@ pc_apply_face_stencil( #endif const amrex::IntVect ivd = amrex::IntVect{ AMREX_D_DECL(iv[0] - 1 + t0, iv[1], iv[2] - 1 + t1)}; + amrex::Real stenval = sten[L].val AMREX_D_TERM(, [t0], [t1]); d_newval[L] += - sten[L].val AMREX_D_TERM(, [t0], [t1]) * vout(ivd, n); + (std::abs(stenval) > 1e-14) ? stenval * vout(ivd, n) : 0.0; #if AMREX_SPACEDIM > 2 } #endif @@ -480,8 +520,9 @@ pc_apply_face_stencil( for (int t1 = 0; t1 < 3; t1++) { const amrex::IntVect ivd = amrex::IntVect{ AMREX_D_DECL(iv[0] - 1 + t0, iv[1] - 1 + t1, iv[2])}; + amrex::Real stenval = sten[L].val AMREX_D_TERM(, [t0], [t1]); d_newval[L] += - sten[L].val AMREX_D_TERM(, [t0], [t1]) * vout(ivd, n); + (std::abs(stenval) > 1e-14) ? stenval * vout(ivd, n) : 0.0; } } #endif @@ -652,8 +693,10 @@ pc_apply_eb_boundry_visc_flux_stencil( for (int kk = 0; kk < 3; kk++) { #endif for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - sum[idir] += sten[L].val AMREX_D_TERM([ii], [jj], [kk]) * - Ut AMREX_D_TERM([ii], [jj], [kk])[idir]; + amrex::Real stenval = sten[L].val AMREX_D_TERM([ii], [jj], [kk]); + sum[idir] += (std::abs(stenval) > 1e-14) + ? stenval * Ut AMREX_D_TERM([ii], [jj], [kk])[idir] + : 0.0; } #if AMREX_SPACEDIM > 2 } @@ -708,8 +751,9 @@ pc_apply_eb_boundry_flux_stencil( const amrex::IntVect ivp = amrex::IntVect(AMREX_D_DECL( sten[L].iv_base[0] + ii, sten[L].iv_base[1] + jj, sten[L].iv_base[2] + kk)); + amrex::Real stenval = sten[L].val AMREX_D_TERM([ii], [jj], [kk]); sum += - sten[L].val AMREX_D_TERM([ii], [jj], [kk]) * s(ivp, scomp + n); + (std::abs(stenval) > 1e-14) ? stenval * s(ivp, scomp + n) : 0.0; #if AMREX_SPACEDIM > 2 } #endif diff --git a/Source/InitEB.cpp b/Source/InitEB.cpp index eb252a854..91919a2b8 100644 --- a/Source/InitEB.cpp +++ b/Source/InitEB.cpp @@ -147,6 +147,13 @@ PeleC::initialize_eb2_structs() amrex::Abort(); } + if (eb_noslip or eb_isothermal) { + amrex::Box sbox = amrex::grow(tbox, -3); + pc_check_bndry_grad_stencil( + sbox, ncutcells, flags.array(mfi), + sv_eb_bndry_grad_stencil[iLocal].data()); + } + sv_eb_flux[iLocal].define(sv_eb_bndry_grad_stencil[iLocal], NVAR); sv_eb_bcval[iLocal].define(sv_eb_bndry_grad_stencil[iLocal], QVAR); @@ -237,6 +244,14 @@ PeleC::define_body_state() return; } + if (eb_zero_body_state) { + for (int n = 0; n < NVAR; ++n) { + body_state[n] = 0.0; + } + body_state_set = true; + return; + } + // Scan over data and find a point in the fluid to use to // set computable values in cells outside the domain // We look for the lexicographically first valid point @@ -303,30 +318,6 @@ PeleC::set_body_state(amrex::MultiFab& S) amrex::Gpu::synchronize(); } -void -PeleC::zero_in_body(amrex::MultiFab& S) const -{ - BL_PROFILE("PeleC::zero_in_body()"); - - if (!eb_in_domain) { - return; - } - - amrex::GpuArray zeros = {0.0}; - auto const& fact = dynamic_cast(Factory()); - auto const& flags = fact.getMultiEBCellFlagFab(); - - auto const& sarrs = S.arrays(); - auto const& flagarrs = flags.const_arrays(); - const amrex::IntVect ngs(0); - amrex::ParallelFor( - S, ngs, NVAR, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { - pc_set_body_state(i, j, k, n, flagarrs[nbx], zeros, sarrs[nbx]); - }); - amrex::Gpu::synchronize(); -} - // Sets up implicit function using EB2 infrastructure void initialize_EB2( diff --git a/Source/Params/_cpp_parameters b/Source/Params/_cpp_parameters index 1c9af9023..28083ccd5 100644 --- a/Source/Params/_cpp_parameters +++ b/Source/Params/_cpp_parameters @@ -222,6 +222,9 @@ eb_srd_max_order int 0 # Weight types for ML redistribution (0 = 1.0, 1 = energy, 2 = density, 3 = vfrac) eb_weights_type int 2 +# Set body state in covered cells to 0 +eb_zero_body_state bool false + #----------------------------------------------------------------------------- # category: method of manufactured solution #----------------------------------------------------------------------------- diff --git a/Source/Params/param_includes/pelec_defaults.H b/Source/Params/param_includes/pelec_defaults.H index ce27e7fc0..c0718366e 100644 --- a/Source/Params/param_includes/pelec_defaults.H +++ b/Source/Params/param_includes/pelec_defaults.H @@ -66,6 +66,7 @@ bool PeleC::eb_clean_massfrac = true; amrex::Real PeleC::eb_clean_massfrac_threshold = 0.0; int PeleC::eb_srd_max_order = 0; int PeleC::eb_weights_type = 2; +bool PeleC::eb_zero_body_state = false; bool PeleC::do_mms = false; std::string PeleC::masa_solution_name = "ad_cns_3d_les"; amrex::Real PeleC::fixed_dt = -1.0; diff --git a/Source/Params/param_includes/pelec_params.H b/Source/Params/param_includes/pelec_params.H index 0a1701737..52174fe3b 100644 --- a/Source/Params/param_includes/pelec_params.H +++ b/Source/Params/param_includes/pelec_params.H @@ -64,6 +64,7 @@ static bool eb_clean_massfrac; static amrex::Real eb_clean_massfrac_threshold; static int eb_srd_max_order; static int eb_weights_type; +static bool eb_zero_body_state; static bool do_mms; static std::string masa_solution_name; static amrex::Real fixed_dt; diff --git a/Source/Params/param_includes/pelec_queries.H b/Source/Params/param_includes/pelec_queries.H index dca2d50a8..865d5f27a 100644 --- a/Source/Params/param_includes/pelec_queries.H +++ b/Source/Params/param_includes/pelec_queries.H @@ -82,6 +82,7 @@ pp.query("eb_clean_massfrac", eb_clean_massfrac); pp.query("eb_clean_massfrac_threshold", eb_clean_massfrac_threshold); pp.query("eb_srd_max_order", eb_srd_max_order); pp.query("eb_weights_type", eb_weights_type); +pp.query("eb_zero_body_state", eb_zero_body_state); pp.query("do_mms", do_mms); pp.query("masa_solution_name", masa_solution_name); pp.query("fixed_dt", fixed_dt); diff --git a/Source/PeleC.H b/Source/PeleC.H index 5ab21fdfd..989aecfee 100644 --- a/Source/PeleC.H +++ b/Source/PeleC.H @@ -232,8 +232,6 @@ public: void set_body_state(amrex::MultiFab& S); - void zero_in_body(amrex::MultiFab& S) const; - void initialize_signed_distance(); void eb_distance(const int lev, amrex::MultiFab& signDistLev); diff --git a/Source/PeleC.cpp b/Source/PeleC.cpp index fdfbe7e53..9a31d92c9 100644 --- a/Source/PeleC.cpp +++ b/Source/PeleC.cpp @@ -1123,6 +1123,7 @@ PeleC::post_timestep(int iteration) amrex::MultiFab& S_new = get_new_data(State_Type); int ng_pts = 0; computeTemp(S_new, ng_pts); + set_body_state(S_new); #ifdef PELEC_USE_SOOT clipSootMoments(S_new, ng_pts); @@ -1407,6 +1408,14 @@ PeleC::okToContinue() << std::endl; } + if ( + get_new_data(State_Type).contains_nan() || + get_new_data(State_Type).contains_inf()) { + test = 0; + amrex::Print() << " Signalling a stop because NaNs detected in the Solution" + << std::endl; + } + return test; } diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 6214bdacc..7fb6a059e 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -113,6 +113,14 @@ function(add_test_rf TEST_NAME TEST_EXE_DIR) set_tests_properties(${TEST_NAME} PROPERTIES WILL_FAIL TRUE) endfunction(add_test_rf) +# Regression tests no FPE trapping ever +function(add_test_rn TEST_NAME TEST_EXE_DIR) + setup_test() + set(RUNTIME_OPTIONS "max_step=10 ${RUNTIME_OPTIONS} amrex.signal_handling=0") + add_test(${TEST_NAME} sh -c "${MPI_COMMANDS} ${CURRENT_TEST_EXE} ${MPIEXEC_POSTFLAGS} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.inp ${RUNTIME_OPTIONS} > ${TEST_NAME}.log ${SAVE_GOLDS_COMMAND} ${FCOMPARE_COMMAND}") + set_tests_properties(${TEST_NAME} PROPERTIES TIMEOUT 18000 PROCESSORS ${PELEC_NP} WORKING_DIRECTORY "${CURRENT_TEST_BINARY_DIR}/" LABELS "regression" ATTACHED_FILES_ON_FAIL "${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.log") +endfunction(add_test_rn) + # Verification test with 1 resolution function(add_test_v1 TEST_NAME TEST_EXE_DIR) setup_test() @@ -229,6 +237,7 @@ add_test_r(sod-2 Sod) add_test_rv(sod-3 Sod) add_test_rv(sod-4 Sod) add_test_r(channel-1 ChannelFlow) +add_test_rn(eb-c3 EB-C3) add_test_r(eb-c4 EB-C4-5) add_test_r(eb-c5 EB-C4-5) add_test_rv(eb-c9 EB-C9) @@ -268,7 +277,6 @@ add_test_re(pmf-lidryer-cvode PMF) add_test_re(sedov-1 Sedov) add_test_re(shu-osher-1 Shu-Osher) add_test_re(zerod-1 zeroD) -add_test_re(eb-c3 EB-C3) add_test_re(soot-flame Soot-Flame) add_test_re(eb-c8 EB-C8) add_test_re(eb-c8-rere EB-C8)