From 63548adf1a49c34b766875f57b18589548615fe8 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 17 Jan 2025 20:25:34 -0800 Subject: [PATCH 1/3] Update eb (#2073) * update EB to use new data structures * update cmake * add CI fixes * address CI issues * fix for gpu * more fixes for new EB stuff * fix for gpu --- CMake/BuildERFExe.cmake | 9 +- CMakeLists.txt | 2 +- Exec/DevTests/EB_Test/ERF_Prob.cpp | 17 +- Exec/DevTests/EB_Test/inputs | 3 +- Source/EB/CMakeLists.txt | 10 + Source/EB/ERF_EB.H | 73 ++ Source/EB/ERF_EB.cpp | 101 +++ Source/EB/ERF_EBAux.H | 48 ++ Source/EB/ERF_EBAux.cpp | 304 ++++++++ Source/EB/ERF_EBBox.cpp | 248 +++--- Source/EB/ERF_EBCutCell.H | 723 ++++++++++++++++++ Source/EB/ERF_EBCutCell.cpp | 43 ++ Source/EB/ERF_EBIF.H | 213 +----- .../EB/{ERF_TerrainIF.H => ERF_EBIFTerrain.H} | 3 +- Source/EB/ERF_EBPolygon.H | 196 +++++ ...edistribute.cpp => ERF_EBRedistribute.cpp} | 27 +- Source/EB/ERF_EBRegular.cpp | 14 - Source/EB/ERF_EBUtils.H | 60 ++ ...teEBSurface.cpp => ERF_EBWriteSurface.cpp} | 0 Source/EB/ERF_InitEB.cpp | 63 -- Source/EB/Make.package | 23 +- Source/ERF.H | 26 +- Source/ERF.cpp | 11 +- .../TimeIntegration/ERF_TI_no_substep_fun.H | 65 +- Source/main.cpp | 1 - 25 files changed, 1825 insertions(+), 458 deletions(-) create mode 100644 Source/EB/CMakeLists.txt create mode 100644 Source/EB/ERF_EB.H create mode 100644 Source/EB/ERF_EB.cpp create mode 100644 Source/EB/ERF_EBAux.H create mode 100644 Source/EB/ERF_EBAux.cpp create mode 100644 Source/EB/ERF_EBCutCell.H create mode 100644 Source/EB/ERF_EBCutCell.cpp rename Source/EB/{ERF_TerrainIF.H => ERF_EBIFTerrain.H} (93%) create mode 100644 Source/EB/ERF_EBPolygon.H rename Source/EB/{ERF_Redistribute.cpp => ERF_EBRedistribute.cpp} (78%) delete mode 100644 Source/EB/ERF_EBRegular.cpp create mode 100644 Source/EB/ERF_EBUtils.H rename Source/EB/{ERF_WriteEBSurface.cpp => ERF_EBWriteSurface.cpp} (100%) delete mode 100644 Source/EB/ERF_InitEB.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 856d87a47..ef225b3c6 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -125,11 +125,12 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Diffusion/ERF_ComputeStrain_N.cpp ${SRC_DIR}/Diffusion/ERF_ComputeStrain_T.cpp ${SRC_DIR}/Diffusion/ERF_ComputeTurbulentViscosity.cpp - ${SRC_DIR}/EB/ERF_InitEB.cpp + ${SRC_DIR}/EB/ERF_EBAux.cpp ${SRC_DIR}/EB/ERF_EBBox.cpp - ${SRC_DIR}/EB/ERF_EBRegular.cpp - ${SRC_DIR}/EB/ERF_InitEB.cpp - ${SRC_DIR}/EB/ERF_WriteEBSurface.cpp + ${SRC_DIR}/EB/ERF_EB.cpp + ${SRC_DIR}/EB/ERF_EBCutCell.cpp + ${SRC_DIR}/EB/ERF_EBRedistribute.cpp + ${SRC_DIR}/EB/ERF_EBWriteSurface.cpp ${SRC_DIR}/Initialization/ERF_InitBCs.cpp ${SRC_DIR}/Initialization/ERF_InitCustom.cpp ${SRC_DIR}/Initialization/ERF_InitFromHSE.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 99c439199..fc22f53ce 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -245,7 +245,7 @@ target_link_libraries(erf_api PUBLIC erf_srclib) add_library(${PROJECT_NAME}::erf_api ALIAS erf_srclib) # Collect all headers and make them installable with the target -set(ERF_INCLUDES "Source/ERF.H;Source/ERF_Constants.H;Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H;Source/WindFarmParametrization/EWP/ERF_EWP.H;Source/WindFarmParametrization/Null/ERF_NullWindFarm.H;Source/WindFarmParametrization/ERF_WindFarm.H;Source/WindFarmParametrization/Fitch/ERF_Fitch.H;Source/BoundaryConditions/ERF_PhysBCFunct.H;Source/BoundaryConditions/ERF_MOSTAverage.H;Source/BoundaryConditions/ERF_MOSTRoughness.H;Source/BoundaryConditions/ERF_ABLMost.H;Source/BoundaryConditions/ERF_FillPatcher.H;Source/BoundaryConditions/ERF_MOSTStress.H;Source/BoundaryConditions/ERF_TimeInterpolatedData.H;Source/Utils/ERF_InitTerrain.H;Source/Utils/ERF_Interpolation.H;Source/Utils/ERF_TileNoZ.H;Source/Utils/ERF_PlaneAverage.H;Source/Utils/ERF_Interpolation_WENO.H;Source/Utils/ERF_DirectionSelector.H;Source/Utils/ERF_ParFunctions.H;Source/Utils/ERF_Wstar.H;Source/Utils/ERF_MicrophysicsUtils.H;Source/Utils/ERF_SatMethods.H;Source/Utils/ERF_Interpolation_1D.H;Source/Utils/ERF_Interpolation_UPW.H;Source/Utils/ERF_TerrainMetrics.H;Source/Utils/ERF_Interpolation_WENO_Z.H;Source/Utils/ERF_Thetav.H;Source/Utils/ERF_WaterVaporSaturation.H;Source/Utils/ERF_Utils.H;Source/Utils/ERF_Orbit.H;Source/Utils/ERF_EOS.H;Source/Utils/ERF_HSEUtils.H;Source/EB/ERF_TerrainIF.H;Source/EB/ERF_EBIF.H;Source/Particles/ERFPC.H;Source/Particles/ERF_ParticleData.H;Source/Prob/ERF_InitDensityHSEDry.H;Source/Prob/ERF_InitRayleighDamping.H;Source/Prob/ERF_InitConstantDensityHSE.H;Source/ERF_ProbCommon.H;Source/ERF_Derive.H;Source/Radiation/ERF_Mam4Constituents.H;Source/Radiation/ERF_Mam4Aero.H;Source/Radiation/ERF_Optics.H;Source/Radiation/ERF_ModalAeroWaterUptake.H;Source/Radiation/ERF_CloudRadProps.H;Source/Radiation/ERF_PhysProp.H;Source/Radiation/ERF_Radiation.H;Source/Radiation/ERF_Albedo.H;Source/Radiation/ERF_Parameterizations.H;Source/Radiation/ERF_RadConstants.H;Source/Radiation/ERF_AeroRadProps.H;Source/Radiation/ERF_M2005EffRadius.H;Source/Radiation/ERF_LinearInterpolate.H;Source/Radiation/ERF_Slingo.H;Source/Radiation/ERF_RRTMGP.H;Source/Radiation/ERF_EbertCurry.H;Source/SourceTerms/ERF_NumericalDiffusion.H;Source/SourceTerms/ERF_SrcHeaders.H;Source/SourceTerms/ERF_ForestDrag.H;Source/IO/ERF_SampleData.H;Source/IO/ERF_NCInterface.H;Source/IO/ERF_NCWpsFile.H;Source/IO/ERF_NCPlotFile.H;Source/IO/ERF_ReadBndryPlanes.H;Source/IO/ERF_WriteBndryPlanes.H;Source/PBL/ERF_MYNNStruct.H;Source/PBL/ERF_PBLModels.H;Source/PBL/ERF_PBLHeight.H;Source/TimeIntegration/ERF_TI_substep_fun.H;Source/TimeIntegration/ERF_TI_slow_headers.H;Source/TimeIntegration/ERF_TI_slow_rhs_fun.H;Source/TimeIntegration/ERF_TI_fast_headers.H;Source/TimeIntegration/ERF_TI_utils.H;Source/TimeIntegration/ERF_MRI.H;Source/TimeIntegration/ERF_TI_no_substep_fun.H;Source/LandSurfaceModel/Null/ERF_NullSurf.H;Source/LandSurfaceModel/ERF_LandSurface.H;Source/LandSurfaceModel/MM5/ERF_MM5.H;Source/LandSurfaceModel/SLM/ERF_SLM.H;Source/ERF_IndexDefines.H;Source/Advection/ERF_AdvectionSrcForMom_N.H;Source/Advection/ERF_AdvectionSrcForScalars.H;Source/Advection/ERF_AdvectionSrcForMom_T.H;Source/Advection/ERF_Advection.H;Source/MultiBlock/ERF_MultiBlockContainer.H;Source/Initialization/ERF_MetgridUtils.H;Source/Diffusion/ERF_EddyViscosity.H;Source/Diffusion/ERF_Diffusion.H;Source/Microphysics/Null/ERF_NullMoistLagrangian.H;Source/Microphysics/Null/ERF_NullMoist.H;Source/Microphysics/ERF_Microphysics.H;Source/Microphysics/ERF_LagrangianMicrophysics.H;Source/Microphysics/ERF_EulerianMicrophysics.H;Source/Microphysics/Kessler/ERF_Kessler.H;Source/Microphysics/SAM/ERF_SAM.H;Source/Microphysics/SatAdj/ERF_SatAdj.H;Source/DataStructs/ERF_InputSpongeData.H;Source/DataStructs/ERF_TurbPertStruct.H;Source/DataStructs/ERF_SpongeStruct.H;Source/DataStructs/ERF_AdvStruct.H;Source/DataStructs/ERF_DataStruct.H;Source/DataStructs/ERF_InputSoundingData.H;Source/DataStructs/ERF_DiffStruct.H;Source/DataStructs/ERF_TurbStruct.H;Source/LinearSolvers/ERF_TerrainPoisson.H;Source/LinearSolvers/ERF_FFTUtils.H") +set(ERF_INCLUDES "Source/ERF.H;Source/ERF_Constants.H;Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H;Source/WindFarmParametrization/EWP/ERF_EWP.H;Source/WindFarmParametrization/Null/ERF_NullWindFarm.H;Source/WindFarmParametrization/ERF_WindFarm.H;Source/WindFarmParametrization/Fitch/ERF_Fitch.H;Source/BoundaryConditions/ERF_PhysBCFunct.H;Source/BoundaryConditions/ERF_MOSTAverage.H;Source/BoundaryConditions/ERF_MOSTRoughness.H;Source/BoundaryConditions/ERF_ABLMost.H;Source/BoundaryConditions/ERF_FillPatcher.H;Source/BoundaryConditions/ERF_MOSTStress.H;Source/BoundaryConditions/ERF_TimeInterpolatedData.H;Source/Utils/ERF_InitTerrain.H;Source/Utils/ERF_Interpolation.H;Source/Utils/ERF_TileNoZ.H;Source/Utils/ERF_PlaneAverage.H;Source/Utils/ERF_Interpolation_WENO.H;Source/Utils/ERF_DirectionSelector.H;Source/Utils/ERF_ParFunctions.H;Source/Utils/ERF_Wstar.H;Source/Utils/ERF_MicrophysicsUtils.H;Source/Utils/ERF_SatMethods.H;Source/Utils/ERF_Interpolation_1D.H;Source/Utils/ERF_Interpolation_UPW.H;Source/Utils/ERF_TerrainMetrics.H;Source/Utils/ERF_Interpolation_WENO_Z.H;Source/Utils/ERF_Thetav.H;Source/Utils/ERF_WaterVaporSaturation.H;Source/Utils/ERF_Utils.H;Source/Utils/ERF_Orbit.H;Source/Utils/ERF_EOS.H;Source/Utils/ERF_HSEUtils.H;Source/EB/ERF_EBIFTerrain.H;Source/EB/ERF_EBIF.H;Source/EB/ERF_EBAux.H;Source/EB/ERF_EBCutCell.H;Source/EB/ERF_EB.H/Source/EB/ERF_EBPolygon.H;Source/EB/ERF_EBUtils.H;Source/Particles/ERFPC.H;Source/Particles/ERF_ParticleData.H;Source/Prob/ERF_InitDensityHSEDry.H;Source/Prob/ERF_InitRayleighDamping.H;Source/Prob/ERF_InitConstantDensityHSE.H;Source/ERF_ProbCommon.H;Source/ERF_Derive.H;Source/Radiation/ERF_Mam4Constituents.H;Source/Radiation/ERF_Mam4Aero.H;Source/Radiation/ERF_Optics.H;Source/Radiation/ERF_ModalAeroWaterUptake.H;Source/Radiation/ERF_CloudRadProps.H;Source/Radiation/ERF_PhysProp.H;Source/Radiation/ERF_Radiation.H;Source/Radiation/ERF_Albedo.H;Source/Radiation/ERF_Parameterizations.H;Source/Radiation/ERF_RadConstants.H;Source/Radiation/ERF_AeroRadProps.H;Source/Radiation/ERF_M2005EffRadius.H;Source/Radiation/ERF_LinearInterpolate.H;Source/Radiation/ERF_Slingo.H;Source/Radiation/ERF_RRTMGP.H;Source/Radiation/ERF_EbertCurry.H;Source/SourceTerms/ERF_NumericalDiffusion.H;Source/SourceTerms/ERF_SrcHeaders.H;Source/SourceTerms/ERF_ForestDrag.H;Source/IO/ERF_SampleData.H;Source/IO/ERF_NCInterface.H;Source/IO/ERF_NCWpsFile.H;Source/IO/ERF_NCPlotFile.H;Source/IO/ERF_ReadBndryPlanes.H;Source/IO/ERF_WriteBndryPlanes.H;Source/PBL/ERF_MYNNStruct.H;Source/PBL/ERF_PBLModels.H;Source/PBL/ERF_PBLHeight.H;Source/TimeIntegration/ERF_TI_substep_fun.H;Source/TimeIntegration/ERF_TI_slow_headers.H;Source/TimeIntegration/ERF_TI_slow_rhs_fun.H;Source/TimeIntegration/ERF_TI_fast_headers.H;Source/TimeIntegration/ERF_TI_utils.H;Source/TimeIntegration/ERF_MRI.H;Source/TimeIntegration/ERF_TI_no_substep_fun.H;Source/LandSurfaceModel/Null/ERF_NullSurf.H;Source/LandSurfaceModel/ERF_LandSurface.H;Source/LandSurfaceModel/MM5/ERF_MM5.H;Source/LandSurfaceModel/SLM/ERF_SLM.H;Source/ERF_IndexDefines.H;Source/Advection/ERF_AdvectionSrcForMom_N.H;Source/Advection/ERF_AdvectionSrcForScalars.H;Source/Advection/ERF_AdvectionSrcForMom_T.H;Source/Advection/ERF_Advection.H;Source/MultiBlock/ERF_MultiBlockContainer.H;Source/Initialization/ERF_MetgridUtils.H;Source/Diffusion/ERF_EddyViscosity.H;Source/Diffusion/ERF_Diffusion.H;Source/Microphysics/Null/ERF_NullMoistLagrangian.H;Source/Microphysics/Null/ERF_NullMoist.H;Source/Microphysics/ERF_Microphysics.H;Source/Microphysics/ERF_LagrangianMicrophysics.H;Source/Microphysics/ERF_EulerianMicrophysics.H;Source/Microphysics/Kessler/ERF_Kessler.H;Source/Microphysics/SAM/ERF_SAM.H;Source/Microphysics/SatAdj/ERF_SatAdj.H;Source/DataStructs/ERF_InputSpongeData.H;Source/DataStructs/ERF_TurbPertStruct.H;Source/DataStructs/ERF_SpongeStruct.H;Source/DataStructs/ERF_AdvStruct.H;Source/DataStructs/ERF_DataStruct.H;Source/DataStructs/ERF_InputSoundingData.H;Source/DataStructs/ERF_DiffStruct.H;Source/DataStructs/ERF_TurbStruct.H;Source/LinearSolvers/ERF_TerrainPoisson.H;Source/LinearSolvers/ERF_FFTUtils.H") set_target_properties( erf_srclib PROPERTIES PUBLIC_HEADER "${ERF_INCLUDES}") diff --git a/Exec/DevTests/EB_Test/ERF_Prob.cpp b/Exec/DevTests/EB_Test/ERF_Prob.cpp index 773ca4eab..ff85f96e3 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.cpp +++ b/Exec/DevTests/EB_Test/ERF_Prob.cpp @@ -82,6 +82,8 @@ Problem::init_custom_pert( const bool use_moisture = (sc.moisture_type != MoistureType::None); + AMREX_ALWAYS_ASSERT(sc.terrain_type == TerrainType::EB); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); // ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -169,17 +171,10 @@ Problem::init_custom_pert( dxInv[2] = 1. / dx[2]; // Set the z-velocity from impenetrable condition - if (sc.terrain_type == TerrainType::EB) { - ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel_pert(i, j, k) = 0.0; - }); - } else { - ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv); - }); - } + ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + z_vel_pert(i, j, k) = 0.0; + }); amrex::Gpu::streamSynchronize(); } diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index 99ed8fa5e..c91d19a00 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -57,7 +57,6 @@ erf.plot_vars_1 = density rhotheta theta rhoadv_0 x_velocity y_velocity z_ve # SOLVER CHOICE erf.use_gravity = false erf.use_coriolis = false -erf.les_type = "None" erf.les_type = "None" erf.molec_diff_type = "None" @@ -67,6 +66,8 @@ erf.alpha_C = 0.0 erf.init_type = "uniform" +erf.terrain_type = EB + # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 diff --git a/Source/EB/CMakeLists.txt b/Source/EB/CMakeLists.txt new file mode 100644 index 000000000..5f0c0ae2a --- /dev/null +++ b/Source/EB/CMakeLists.txt @@ -0,0 +1,10 @@ +target_include_directories( simplexacore PUBLIC + $) + +target_sources( simplexacore PRIVATE + eb.cpp + eb_box.cpp + eb_terrain.cpp + eb_cut_cell.cpp + eb_aux.cpp + ) diff --git a/Source/EB/ERF_EB.H b/Source/EB/ERF_EB.H new file mode 100644 index 000000000..107c60c10 --- /dev/null +++ b/Source/EB/ERF_EB.H @@ -0,0 +1,73 @@ +#ifndef ERF_EB_H_ +#define ERF_EB_H_ + +#include +#include +#include + +#include +#include + +#include + +class eb_ { + + public: + + ~eb_ (); + + eb_ (amrex::Geometry const& a_geom, + amrex::FArrayBox const& terrain_fab, + amrex::Gpu::DeviceVector& a_dz_stretched, + bool is_anelastic); + + void make_factory ( amrex::Geometry const& a_geom, + amrex::DistributionMapping const& a_dmap, + amrex::BoxArray const& a_grids); + + int nghost_basic () const { return 2; } + int nghost_volume () const { return 2; } + int nghost_full () const { return 2; } + + amrex::EBFArrayBoxFactory const* get_const_factory () noexcept + { return m_factory.get(); } + + amrex::EB2::Level const* get_level () const noexcept + { return m_eb_level; } + + private: + + int m_has_eb; + std::string m_type; + + amrex::EBSupport m_support_level; + + int m_write_eb_surface; + + //! EB level constructed from building GeometryShop + amrex::EB2::Level const* m_eb_level; + + std::unique_ptr m_factory; + + eb_aux_ m_u_factory; + eb_aux_ m_v_factory; + eb_aux_ m_w_factory; + + void make_box (amrex::Geometry const& a_geom); + void make_terrain (amrex::Geometry const& a_geom); + + //! Construct EB levels from Geometry shop. + template + void build_level (amrex::Geometry const& a_geom, + amrex::EB2::GeometryShop a_gshop) + { + int const req_lev(0); + int const max_lev(2); + + amrex::EB2::Build(a_gshop, a_geom, req_lev, max_lev); + const amrex::EB2::IndexSpace& ebis = amrex::EB2::IndexSpace::top(); + m_eb_level = &(ebis.getLevel(a_geom)); + } + +}; +#endif diff --git a/Source/EB/ERF_EB.cpp b/Source/EB/ERF_EB.cpp new file mode 100644 index 000000000..46f8b5139 --- /dev/null +++ b/Source/EB/ERF_EB.cpp @@ -0,0 +1,101 @@ +#include +#include +#include + +#include +#include + +#include + +using namespace amrex; + +eb_::~eb_() +{ + if (m_factory) { m_factory.reset(nullptr); } +} + +eb_:: eb_ ( Geometry const& a_geom, amrex::FArrayBox const& terrain_fab, + amrex::Gpu::DeviceVector& a_dz_stretched, + bool is_anelastic) + : m_has_eb(0), + m_support_level(EBSupport::full), + m_write_eb_surface(0) +{ + m_type = "terrain"; + + int max_coarsening_level; + if (is_anelastic) { + max_coarsening_level = 100; + } else { + max_coarsening_level = 0; + } + + int max_level_here = 0; + + if (m_type == "terrain") + { + Print() << "\nBuilding EB geometry based on terrain.\n"; + + TerrainIF ebterrain(terrain_fab, a_geom, a_dz_stretched); + + auto gshop = EB2::makeShop(ebterrain); + + EB2::Build(gshop, a_geom, max_level_here, max_level_here+max_coarsening_level); + + m_has_eb = 1; + + } else if (m_type == "box") { + + Print() << "\nBuilding box geometry.\n"; + make_box(a_geom); + m_has_eb = 1; + + } else { + + EB2::AllRegularIF regular; + auto gshop = EB2::makeShop(regular); + build_level(a_geom, gshop); + } +} + +void +eb_:: +make_factory ( Geometry const& a_geom, + DistributionMapping const& a_dmap, + BoxArray const& a_grids) +{ + Print() << "making EB factory\n"; + m_factory = std::make_unique(*m_eb_level, a_geom, a_grids, a_dmap, + Vector{nghost_basic(), nghost_volume(), nghost_full()}, m_support_level); + + if (m_write_eb_surface) { + WriteEBSurface(a_grids, a_dmap, a_geom, m_factory.get()); + } + + { int const idim(0); + + Print() << "making EB staggered u-factory\n"; + //m_u_factory.set_verbose(); + m_u_factory.define(idim, a_geom, a_grids, a_dmap, + Vector{nghost_basic(), nghost_volume(), nghost_full()}, + m_factory.get()); + } + + { int const idim(1); + Print() << "making EB staggered v-factory\n"; + //m_v_factory.set_verbose(); + m_v_factory.define(idim, a_geom, a_grids, a_dmap, + Vector{nghost_basic(), nghost_volume(), nghost_full()}, + m_factory.get()); + } + + { int const idim(2); + Print() << "making EB staggered w-factory\n"; + //m_w_factory.set_verbose(); + m_w_factory.define(idim, a_geom, a_grids, a_dmap, + Vector{nghost_basic(), nghost_volume(), nghost_full()}, + m_factory.get()); + } + + Print() << "\nDone making EB factory.\n\n"; +} diff --git a/Source/EB/ERF_EBAux.H b/Source/EB/ERF_EBAux.H new file mode 100644 index 000000000..ead438395 --- /dev/null +++ b/Source/EB/ERF_EBAux.H @@ -0,0 +1,48 @@ +#ifndef ERF_EB_AUX_H_ +#define ERF_EB_AUX_H_ + +#include +#include +#include + +#include +#include + +class eb_aux_ +{ + public: + + eb_aux_ (); + ~eb_aux_ (); + + void define ( int const& a_idim, + amrex::Geometry const& a_geom, + amrex::BoxArray const& a_grids, + amrex::DistributionMapping const& a_dmap, + amrex::Vector const& a_ngrow, + amrex::EBFArrayBoxFactory const* a_factory); + + void set_verbose ( ) { m_verbose = 1; } + + private: + + int m_verbose; + + // int m_defined; + + amrex::FabArray* m_cellflags = nullptr; + + amrex::MultiFab* m_volfrac = nullptr; + + // amrex::MultiCutFab* m_centroid = nullptr; + // amrex::MultiCutFab* m_bndrycent = nullptr; + // amrex::MultiCutFab* m_bndryarea = nullptr; + // amrex::MultiCutFab* m_bndrynorm = nullptr; + + // amrex::Array m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; + // amrex::Array m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; + // amrex::Array m_edgecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; + +}; + +#endif diff --git a/Source/EB/ERF_EBAux.cpp b/Source/EB/ERF_EBAux.cpp new file mode 100644 index 000000000..1cccdd73d --- /dev/null +++ b/Source/EB/ERF_EBAux.cpp @@ -0,0 +1,304 @@ +#include +#include + +using namespace amrex; + +eb_aux_:: +~eb_aux_ () +{ +} + +eb_aux_:: +eb_aux_ () + : m_verbose(0) +// ,m_defined(0) +{} + +void +eb_aux_:: +define( int const& a_idim, + Geometry const& /*a_geom*/, + BoxArray const& a_grids, + DistributionMapping const& a_dmap, + Vector const& a_ngrow, + EBFArrayBoxFactory const* a_factory) +{ + // Box dbox(a_geom.Domain()); + + const IntVect vdim(IntVect::TheDimensionVector(a_idim)); + + const BoxArray& grids = amrex::convert(a_grids, vdim); + + m_cellflags = new FabArray(grids, a_dmap, 1, a_ngrow[0], MFInfo(), + DefaultFabFactory()); + + m_cellflags->setVal(EBCellFlag::TheDefaultCell()); + + m_volfrac = new MultiFab(grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory()); + + +#if 0 + m_centroid = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_ngrow[1], *m_cellflags); + m_bndrycent = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_grow[2], *m_cellflags); + + m_bndryarea = new MultiCutFab(a_ba, a_dm, 1, m_grow[2], *m_cellflags); + m_bndrynorm = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, m_grow[2], *m_cellflags); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const BoxArray& faceba = amrex::convert(a_ba, IntVect::TheDimensionVector(idim)); + m_areafrac[idim] = new MultiCutFab(faceba, a_dm, 1, m_grow[2], *m_cellflags); + m_facecent[idim] = new MultiCutFab(faceba, a_dm, AMREX_SPACEDIM-1, m_grow[2], *m_cellflags); + } +#endif + + const auto& FlagFab = a_factory->getMultiEBCellFlagFab(); + + for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.validbox(); + + if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) { + + Array4 const& flag = FlagFab.const_array(mfi); + + Array4 const& vfrac = (a_factory->getVolFrac()).const_array(mfi); + // Array4 const& ccent = (a_factory->getCentroid()).const_array(mfi); + + Array4 const& afrac = (a_factory->getAreaFrac()[a_idim])->const_array(mfi); + + // EB normal and face centroid + Array4 const& bnorm = a_factory->getBndryNormal()[mfi].const_array(); + Array4 const& bcent = a_factory->getBndryCent()[mfi].const_array(); + + Array4 const& aux_flag = m_cellflags->array(mfi); + Array4 const& aux_vfrac = m_volfrac->array(mfi); + + ParallelFor(bx, [ +#ifndef AMREX_USE_GPU + verbose=m_verbose, +#endif + vfrac, afrac, bnorm, bcent, flag, + aux_flag, aux_vfrac, vdim, idim=a_idim ] + AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + aux_flag(i,j,k).setCovered(); + aux_flag(i,j,k).setDisconnected(); + + aux_vfrac(i,j,k) = 0.0; + + IntVect const iv_hi(i,j,k); + IntVect const iv_lo(iv_hi - vdim); + + if ( flag(iv_lo).isRegular() && flag(iv_hi).isRegular()) { + + aux_flag(i,j,k).setRegular(); + aux_flag(i,j,k).setConnected(vdim); + + aux_vfrac(i,j,k) = 1.0; + + } else if ( flag(iv_lo).isCovered() && flag(iv_hi).isCovered()) { + + // defaults to covered and disconnected. + + } else { + + +#ifndef AMREX_USE_GPU + if (verbose) { Print() << "\ncell: " << amrex::IntVect(i,j,k) << "\n"; } +#endif + Array lo_arr = {-0.5,-0.5,-0.5}; + Array hi_arr = { 0.5, 0.5, 0.5}; + + // plane point and normal + RealVect lo_point(bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2)); + RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2)); + + // High side of low cell + lo_arr[idim] = 0.0; + hi_arr[idim] = 0.5; + RealBox lo_rbx(lo_arr.data(), hi_arr.data()); + + eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal); + + // cell iv_lo covered (regular) imples lo_eb_cc is covered (regular) + // The inverse is not always true. + AMREX_ASSERT( !flag(iv_lo).isCovered() || lo_eb_cc.isCovered() ); + AMREX_ASSERT( !flag(iv_lo).isRegular() || lo_eb_cc.isRegular() ); + + // plane point and normal + RealVect hi_point(bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2)); + RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2)); + + // Low side of high cell + lo_arr[idim] = -0.5; + hi_arr[idim] = 0.0; + RealBox hi_rbx(lo_arr.data(), hi_arr.data()); + + eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal); + + // cell iv_hi covered (regular) imples hi_eb_cc is covered (regular) + // The inverse is not always true. + AMREX_ASSERT( !flag(iv_hi).isCovered() || hi_eb_cc.isCovered() ); + AMREX_ASSERT( !flag(iv_hi).isRegular() || hi_eb_cc.isRegular() ); + +#if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1 + + { /***************************** SANITY CHECK ***********************\ + * Perform some basic sanity checks to verify that what we computed * + * for cell (i,j,k) compares to what we know to be true. * + \******************************************************************/ + + // Compute the cut-cell for the high side of the high cell. This is + // only needed for sanity checks. + + eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal); + + // cell iv_hi covered (regular) imples hi_hi_eb_cc is covered (regular) + // The inverse is not always true. +#ifndef AMREX_USE_GPU + if ( !(!flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular()) || + !(!flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered()) ) { + Print() << "flag(iv_hi) and hi_hi_eb_cc flags do not agree\n" + << "\n isRegular() " << flag(iv_hi).isRegular() << " " << hi_hi_eb_cc.isRegular() + << "\n isCovered() " << flag(iv_hi).isCovered() << " " << hi_hi_eb_cc.isCovered() + << "\n"; + } +#endif + // If cell iv_hi is regular or covered, then hi_hi_eb_cc must also + // be regular or covered. The inverse is not true. + AMREX_ALWAYS_ASSERT( !flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular() ); + AMREX_ALWAYS_ASSERT( !flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered() ); + + // The area and volume fractions that are computed for the scalar grid + // are slightly different than those we compute from the geometric + // reconstruction using the EB point and normal. However, we expect + // that the area fractions computed here will give back the same + // normal we used to compute them. + if ( flag(iv_hi).isSingleValued() ) { + + Real const adx = (idim == 0) + ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0)) + : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0)) + - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0)); + + Real const ady = (idim == 1) + ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1)) + : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1)) + - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1)); + + Real const adz = (idim == 2) + ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2)) + : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2)) + - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2)); + + Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz); + + // EB normal + Real const apnorminv = 1. / apnorm; + RealVect const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv); + Real const dot_normals = normal.dotProduct(hi_normal); + +#ifndef AMREX_USE_GPU + if ( !amrex::almostEqual(dot_normals, 1.0) ) { + Print() << "\nFail: check-1 dot_normals " << dot_normals + << '\n'; + + hi_eb_cc.debug(); + hi_hi_eb_cc.debug(); + + } else if (verbose) { + Print() << "Pass: dot_normals = 1.0\n"; + + } +#endif + AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) ); + } + + // The idim area of hi_eb_cc.areaHi() should equal hi_hi_eb_cc.areaLo() + { +#ifndef AMREX_USE_GPU + Real const abs_err = std::abs( hi_eb_cc.areaHi(idim) - hi_hi_eb_cc.areaLo(idim) ); + Real machine_tol = 10.0*std::numeric_limits::epsilon(); + if ( abs_err >= machine_tol ) { + Print() << "\nFail: check-4 area abs_err: " << abs_err + << "\n hi_eb_cc.areaHi " << hi_eb_cc.areaHi(idim) + << "\n hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(idim) + << '\n'; + } else if (verbose) { + Print() << "Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo" + << " abs_err: " << abs_err << "\n"; + } + AMREX_ALWAYS_ASSERT( abs_err < machine_tol ); +#endif + } + + // The low-side area of hi_eb_cc should equal idim afrac. + { Real const abs_err = amrex::min(std::abs(lo_eb_cc.areaHi(idim) - afrac(iv_hi)), + std::abs(hi_eb_cc.areaLo(idim) - afrac(iv_hi))); + Real compare_tol = 5.0e-6; +#ifndef AMREX_USE_GPU + if ( abs_err >= compare_tol ) { + //hi_eb_cc.debug(); + Print() << "\nFail: check-2 area abs_err " << abs_err + << "\n hi_eb_cc.areaLo(" << idim << ") = " << hi_eb_cc.areaLo(idim) + << "\n lo_eb_cc.areaHi(" << idim << ") = " << lo_eb_cc.areaHi(idim) + << "\n afrac" << iv_hi << " = " << afrac(iv_hi) + << '\n'; + } else if (verbose) { + Print() << "Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi) + << " abs_err: " << abs_err << "\n"; + } +#endif + AMREX_ALWAYS_ASSERT( abs_err < compare_tol ); + } + + // The combined volumes of hi_eb_cc.areaHi() and hi_hi_eb_cc should + // equal vfrac(iv_hi). + { Real const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume(); + Real const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol); + Real compare_tol = 5.0e-6; +#ifndef AMREX_USE_GPU + if ( abs_err >= compare_tol ) { + hi_eb_cc.debug(); + hi_hi_eb_cc.debug(); + amrex::Print() << "\nFail: check-4 volume abs_err: " << abs_err + << "\n point: " << hi_point + << "\n normal: " << hi_normal + << "\n hi_eb_cc.volume() " << hi_eb_cc.volume() + << "\n hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume() + << "\n vfrac: " << vfrac(iv_hi) + << '\n'; + } else if (verbose) { + Print() << "Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi) + << " abs_err: " << abs_err << "\n"; + } +#endif + AMREX_ALWAYS_ASSERT( abs_err < compare_tol ); + } + } +#endif + + if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) { + + // defaults to covered and disconnected. + + } else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) { + + aux_flag(i,j,k).setRegular(); + aux_flag(i,j,k).setConnected(vdim); + + aux_vfrac(i,j,k) = 1.0; + + } else { + + if (lo_eb_cc.isCovered()) { } + if (hi_eb_cc.isCovered()) { } + + } + } + }); + + } + + } +} diff --git a/Source/EB/ERF_EBBox.cpp b/Source/EB/ERF_EBBox.cpp index e48ef0b76..d2598a565 100644 --- a/Source/EB/ERF_EBBox.cpp +++ b/Source/EB/ERF_EBBox.cpp @@ -1,137 +1,125 @@ #include -#include -#include +#include +#include +#include +#include -#include -#include +#include using namespace amrex; -/**************************************************************************** - * Function to create a simple rectangular box with EB walls. * - * * - ****************************************************************************/ -void ERF::make_eb_box() +/************************************************************************ + * * + * Function to create a simple box geometry: * + * -> box.{Lo,Hi} vector storing box lo/hi * + * -> box.offset vector storing box offset * + * * + * NOTE: walls are placed _outside_ domain for periodic directions. * + * * + ***********************************************************************/ +void +eb_:: +make_box ( Geometry const& a_geom) { - // Get box information from inputs file - ParmParse pp("box"); - - if(geom[0].isAllPeriodic()) - { - make_eb_regular(); - } - else - { - /************************************************************************ - * * - * Define Box geometry: * - * -> box.{Lo,Hi} vector storing box lo/hi * - * -> box.offset vector storing box offset * - * NOTE: walls are placed _outside_ domain for periodic directions. * - * * - ************************************************************************/ - - Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); - Real offset = 1.0e-15; - - for(int i = 0; i < AMREX_SPACEDIM; i++) - { - boxLo[i] = geom[0].ProbLo(i); - boxHi[i] = geom[0].ProbHi(i); - } - - pp.queryarr("Lo", boxLo, 0, AMREX_SPACEDIM); - pp.queryarr("Hi", boxHi, 0, AMREX_SPACEDIM); - - pp.query("offset", offset); - - Real xlo = boxLo[0] + offset; - Real xhi = boxHi[0] - offset; - - // This ensures that the walls won't even touch the ghost cells. By - // putting them one domain width away - if(geom[0].isPeriodic(0)) - { - xlo = 2.0 * geom[0].ProbLo(0) - geom[0].ProbHi(0); - xhi = 2.0 * geom[0].ProbHi(0) - geom[0].ProbLo(0); - } - - Real ylo = boxLo[1] + offset; - Real yhi = boxHi[1] - offset; - - // This ensures that the walls won't even touch the ghost cells. By - // putting them one domain width away - if(geom[0].isPeriodic(1)) - { - ylo = 2.0 * geom[0].ProbLo(1) - geom[0].ProbHi(1); - yhi = 2.0 * geom[0].ProbHi(1) - geom[0].ProbLo(1); - } - -#if (AMREX_SPACEDIM == 2) - Array point_lox{xlo, 0.0}; - Array normal_lox{-1.0, 0.0}; - Array point_hix{xhi, 0.0}; - Array normal_hix{1.0, 0.0}; - - Array point_loy{0.0, ylo}; - Array normal_loy{0.0, -1.0}; - Array point_hiy{0.0, yhi}; - Array normal_hiy{0.0, 1.0}; - - EB2::PlaneIF plane_lox(point_lox, normal_lox); - EB2::PlaneIF plane_hix(point_hix, normal_hix); - - EB2::PlaneIF plane_loy(point_loy, normal_loy); - EB2::PlaneIF plane_hiy(point_hiy, normal_hiy); - - // Generate GeometryShop - auto gshop = EB2::makeShop(EB2::makeUnion(plane_lox, plane_hix, - plane_loy, plane_hiy)); -#else - Real zlo = boxLo[2] + offset; - Real zhi = boxHi[2] - offset; - - // This ensures that the walls won't even touch the ghost cells. By - // putting them one domain width away - if(geom[0].isPeriodic(2)) - { - zlo = 2.0 * geom[0].ProbLo(2) - geom[0].ProbHi(2); - zhi = 2.0 * geom[0].ProbHi(2) - geom[0].ProbLo(2); - } - - Array point_lox{xlo, 0.0, 0.0}; - Array normal_lox{-1.0, 0.0, 0.0}; - Array point_hix{xhi, 0.0, 0.0}; - Array normal_hix{1.0, 0.0, 0.0}; - - Array point_loy{0.0, ylo, 0.0}; - Array normal_loy{0.0, -1.0, 0.0}; - Array point_hiy{0.0, yhi, 0.0}; - Array normal_hiy{0.0, 1.0, 0.0}; - - Array point_loz{0.0, 0.0, zlo}; - Array normal_loz{0.0, 0.0, -1.0}; - Array point_hiz{0.0, 0.0, zhi}; - Array normal_hiz{0.0, 0.0, 1.0}; - - EB2::PlaneIF plane_lox(point_lox, normal_lox); - EB2::PlaneIF plane_hix(point_hix, normal_hix); - - EB2::PlaneIF plane_loy(point_loy, normal_loy); - EB2::PlaneIF plane_hiy(point_hiy, normal_hiy); - - EB2::PlaneIF plane_loz(point_loz, normal_loz); - EB2::PlaneIF plane_hiz(point_hiz, normal_hiz); - - // Generate GeometryShop - auto gshop = EB2::makeShop(EB2::makeUnion(plane_lox, plane_hix, - plane_loy, plane_hiy, - plane_loz, plane_hiz)); -#endif - - // Build index space - int max_level_here = 0; - int max_coarsening_level = 100; - EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level); - } + ParmParse pp_box("eb.box"); + + bool inside = true; + + Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); + + Real const* dx = a_geom.CellSize(); + + if ( !almostEqual(dx[0],dx[1]) || !almostEqual(dx[1],dx[2])) { + amrex::Error(" EB Error: Mesh spacing must be uniform!.\n"); + } + + Real offset = 0.01*dx[0]; + + for (int i = 0; i < AMREX_SPACEDIM; i++) { + boxLo[i] = a_geom.ProbLo(i); + boxHi[i] = a_geom.ProbHi(i); + } + + pp_box.queryarr("lo", boxLo, 0, AMREX_SPACEDIM); + pp_box.queryarr("hi", boxHi, 0, AMREX_SPACEDIM); + + pp_box.query("internal_flow", inside); + pp_box.query("offset", offset); + + Real xlo = boxLo[0] + offset; + Real xhi = boxHi[0] - offset; + + // This ensures that the walls won't even touch the ghost cells. By + // putting them one domain width away + if (a_geom.isPeriodic(0)) + { + xlo = 2.0*a_geom.ProbLo(0) - a_geom.ProbHi(0); + xhi = 2.0*a_geom.ProbHi(0) - a_geom.ProbLo(0); + } + + + Real ylo = boxLo[1] + offset; + Real yhi = boxHi[1] - offset; + + // This ensures that the walls won't even touch the ghost cells. By + // putting them one domain width away + if (a_geom.isPeriodic(1)) + { + ylo = 2.0*a_geom.ProbLo(1) - a_geom.ProbHi(1); + yhi = 2.0*a_geom.ProbHi(1) - a_geom.ProbLo(1); + } + + Real zlo = boxLo[2] + offset; + Real zhi = boxHi[2] - offset; + + // This ensures that the walls won't even touch the ghost cells. By + // putting them one domain width away + if (a_geom.isPeriodic(2)) + { + zlo = 2.0*a_geom.ProbLo(2) - a_geom.ProbHi(2); + zhi = 2.0*a_geom.ProbHi(2) - a_geom.ProbLo(2); + } + + Array point_lox{ xlo, 0.0, 0.0}; + Array normal_lox{-1.0, 0.0, 0.0}; + Array point_hix{ xhi, 0.0, 0.0}; + Array normal_hix{ 1.0, 0.0, 0.0}; + + Array point_loy{0.0, ylo, 0.0}; + Array normal_loy{0.0,-1.0, 0.0}; + Array point_hiy{0.0, yhi, 0.0}; + Array normal_hiy{0.0, 1.0, 0.0}; + + Array point_loz{0.0, 0.0, zlo}; + Array normal_loz{0.0, 0.0,-1.0}; + Array point_hiz{0.0, 0.0, zhi}; + Array normal_hiz{0.0, 0.0, 1.0}; + + EB2::PlaneIF plane_lox(point_lox,normal_lox); + EB2::PlaneIF plane_hix(point_hix,normal_hix); + + EB2::PlaneIF plane_loy(point_loy,normal_loy); + EB2::PlaneIF plane_hiy(point_hiy,normal_hiy); + + EB2::PlaneIF plane_loz(point_loz,normal_loz); + EB2::PlaneIF plane_hiz(point_hiz,normal_hiz); + + auto box = EB2::makeUnion(plane_lox, plane_hix, plane_loy, + plane_hiy, plane_loz, plane_hiz ); + + + if( inside ) { + + auto gshop = EB2::makeShop(box); + build_level(a_geom, gshop); + + } else { + + auto xob = EB2::makeComplement(box); + auto gshop = EB2::makeShop(xob); + + build_level(a_geom, gshop); + } + } + diff --git a/Source/EB/ERF_EBCutCell.H b/Source/EB/ERF_EBCutCell.H new file mode 100644 index 000000000..f465a7da4 --- /dev/null +++ b/Source/EB/ERF_EBCutCell.H @@ -0,0 +1,723 @@ +#ifndef ERF_EB_CUT_CELL_H_ +#define ERF_EB_CUT_CELL_H_ + +#include +#include +#include +#include + +#include "ERF_EBPolygon.H" +#include "ERF_EBUtils.H" + +class eb_cut_cell_ { + + public: + + AMREX_GPU_HOST_DEVICE + eb_cut_cell_ ( amrex::EBCellFlag const& a_flag, + amrex::RealBox const& a_rbox, + amrex::RealVect const& a_point, + amrex::RealVect const& a_normal ); + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool isCovered () const noexcept { return m_flag.isCovered(); } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool isRegular () const noexcept { return m_flag.isRegular(); } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool isSingleValued () const noexcept { return m_flag.isSingleValued(); } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void set_covered () { m_flag.setCovered(); } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void set_regular () { + + m_flag.setRegular(); + + m_F1.set_area( m_rbox.length(0)*m_rbox.length(1) ); + m_F3.set_area( m_rbox.length(0)*m_rbox.length(1) ); + + m_F2.set_area( m_rbox.length(1)*m_rbox.length(2) ); + m_F4.set_area( m_rbox.length(1)*m_rbox.length(2) ); + + m_F5.set_area( m_rbox.length(0)*m_rbox.length(2) ); + m_F6.set_area( m_rbox.length(0)*m_rbox.length(2) ); + + m_F7.set_area(0.); + + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real volume () { + + if (m_flag.isCovered() ) { return 0.; } + if (m_flag.isRegular() ) { return m_rbox.volume(); } + + amrex::Real volume(0.); + + amrex::Real const* lo = m_rbox.lo(); + amrex::RealVect v0(lo[0], lo[1], lo[2]); + + if( m_F2.ok() ) { volume += m_F2.area() * m_F2.distance(v0); } + if (m_F3.ok() ) { volume += m_F3.area() * m_F3.distance(v0); } + if (m_F6.ok() ) { volume += m_F6.area() * m_F6.distance(v0); } + if (m_F7.ok() ) { volume += m_F7.area() * m_F7.distance(v0); } + + volume /= 3.; + + return m_invert*volume + (1.-m_invert)*(m_rbox.volume()-volume); + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real areaLo ( int const a_idim ) const noexcept { + AMREX_ASSERT( a_idim >0 && a_idim < AMREX_SPACEDIM ); + if (m_flag.isCovered() ) { return 0.; } + if (m_flag.isRegular() ) { return 1.; } + amrex::Real const area(m_lo_faces[a_idim]->area()); + return m_invert*area + (1.-m_invert)*(m_rbox_area[a_idim] - area); + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real areaHi ( int const a_idim ) const noexcept { + AMREX_ASSERT( a_idim > 0 && a_idim < AMREX_SPACEDIM ); + if (m_flag.isCovered() ) { return 0.; } + if (m_flag.isRegular() ) { return 1.; } + amrex::Real const area(m_hi_faces[a_idim]->area()); + return m_invert*area + (1.-m_invert)*(m_rbox_area[a_idim] - area); + } + + void debug ( int const a_face = -1 ); + + private: + + amrex::RealBox const m_rbox; + amrex::RealVect const m_eb_point; + amrex::RealVect const m_eb_normal; + + amrex::Real m_invert; + + amrex::RealVect m_rbox_area; + + amrex::EBCellFlag m_flag; + + // Cell faces + + polygon_ m_F1; + polygon_ m_F2; + polygon_ m_F3; + polygon_ m_F4; + polygon_ m_F5; + polygon_ m_F6; + + // cut face + polygon_ m_F7; + + amrex::Array m_lo_faces {&m_F4, &m_F5, &m_F1}; + amrex::Array m_hi_faces {&m_F2, &m_F6, &m_F3}; + + AMREX_GPU_HOST_DEVICE + void calc_edge_intersections ( int const a_dry_run = 0 ); + +}; + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +eb_cut_cell_:: +eb_cut_cell_ ( amrex::EBCellFlag const& a_flag, + amrex::RealBox const& a_rbox, + amrex::RealVect const& a_point, + amrex::RealVect const& a_normal ) + : m_rbox(a_rbox) + , m_eb_point(a_point) + , m_eb_normal(a_normal) + , m_invert(0.0) + , m_F1(a_point, a_normal) + , m_F2(a_point, a_normal) + , m_F3(a_point, a_normal) + , m_F4(a_point, a_normal) + , m_F5(a_point, a_normal) + , m_F6(a_point, a_normal) +{ + using namespace amrex; + + m_rbox_area[0] = m_rbox.length(1)*m_rbox.length(2); + m_rbox_area[1] = m_rbox.length(0)*m_rbox.length(2); + m_rbox_area[2] = m_rbox.length(0)*m_rbox.length(1); + + RealVect v0(a_rbox.lo(0), a_rbox.lo(1), a_rbox.lo(2)); + RealVect v7(a_rbox.hi(0), a_rbox.hi(1), a_rbox.hi(2)); + + if (a_flag.isCovered() ) { + + set_covered(); + + } else if (a_flag.isRegular() ) { + + set_regular(); + + } else { // Check that the box and plane intersect. + + RealVect c = 0.5*(v0 + v7); + RealVect e = v7 - c; + + Real r = e[0]*amrex::Math::abs(a_normal[0]) + + e[1]*amrex::Math::abs(a_normal[1]) + + e[2]*amrex::Math::abs(a_normal[2]); + + Real s = amrex::Math::abs(c.dotProduct(a_normal) + - a_point.dotProduct(a_normal)); + + if (s > r) { + if (a_normal.dotProduct(v0 - a_point) > 0.) + { set_covered(); } else { set_regular(); } + } else { m_flag.setSingleValued(); } + } + + if ( m_flag.isSingleValued() ) { + + m_invert = ((m_eb_normal.dotProduct(v0 - m_eb_point)) > 0.) ? 0.0 : 1.0; + + calc_edge_intersections(); + + } // end singleValued + + m_F1.define(); + m_F2.define(); + m_F3.define(); + m_F4.define(); + m_F5.define(); + m_F6.define(); + m_F7.define(); + +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +eb_cut_cell_:: +calc_edge_intersections ( int const a_dry_run ) +{ + using namespace amrex; + + RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2)); + RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2)); + RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2)); + RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2)); + RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2)); + RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2)); + RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2)); + RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2)); + + if (!a_dry_run) { + m_F1.add_vertex(v0); + m_F4.add_vertex(v0); + m_F5.add_vertex(v0); + } + + int add_v7(1); + + RealVect vIP; + Real distIP; + + // Path 1 + { int cuts(0); + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v1, vIP, distIP)) { + + ++cuts; +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v0--v1: add vIP to F1, F5 and F7 :: " << vIP[0] << "\n"; } + else +#endif + { + m_F1.add_vertex(vIP); + m_F5.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0) ) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v0--v1: intersection ~ 0 :: add vIP to F4\n"; } + else +#endif + { m_F4.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0) ) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v0--v1: intersection ~ 1 :: add vIP to F2\n"; } + else +#endif + { m_F2.add_vertex(vIP); } + } + + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: Add v1 to F1, F2 and F5\n"; } + else +#endif + { + m_F1.add_vertex(v1); + m_F2.add_vertex(v1); + m_F5.add_vertex(v1); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v4, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v1--v4: add vIP to F2, F5 and F7 :: " << vIP[2] << "\n"; } + else +#endif + { + m_F2.add_vertex(vIP); + m_F5.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0) ) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v1--v4: intersection ~ 0 :: add vIP to F1\n"; } + else +#endif + {m_F1.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0) ) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v1--v4: intersection ~ 1 :: add vIP to F3\n"; } + else +#endif + { m_F3.add_vertex(vIP); } + } + + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: Add v4 to F2, F3 and F5\n"; } + else +#endif + { + m_F2.add_vertex(v4); + m_F3.add_vertex(v4); + m_F5.add_vertex(v4); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v4, v7, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v4--v7: add vIP to F2, F3 and F7 :: " << vIP[1] << "\n"; } + else +#endif + { + m_F2.add_vertex(vIP); + m_F3.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v4--v7: intersection ~ 0 :: add vIP to F5\n"; } + else +#endif + { m_F5.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v4--v7: intersection ~ 1 :: add vIP to F6\n"; } + else +#endif + { m_F6.add_vertex(vIP); } + } + } + + if (cuts == 2 && add_v7) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: Add v7 to F2, F3 and F6\n"; } + else +#endif + { + m_F2.add_vertex(v7); + m_F3.add_vertex(v7); + m_F6.add_vertex(v7); + } + + add_v7 = 0; + } + } // end Path 1 + + // Path 4 + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v5, vIP, distIP)) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P4: v1--v5: add vIP to F1, F2 and F7 :: " << vIP[1] << "\n"; } + else +#endif + { + m_F1.add_vertex(vIP); + m_F2.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v1--v5: intersection ~ 0 :: add vIP to F5\n"; } + else +#endif + { m_F5.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: v1--v5: intersection ~ 1 :: add vIP to F6\n"; } + else +#endif + { m_F6.add_vertex(vIP); } + } + } + + + // Path 2 + { int cuts(0); + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v2, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v0--v2: add vIP to F1, F4 and F7 :: " << vIP[1] << "\n"; } + else +#endif + { + m_F1.add_vertex(vIP); + m_F4.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v0--v2: intersection ~ 0 :: add vIP to F5\n"; } + else +#endif + { m_F5.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v0--v2: intersection ~ 1 :: add vIP to F6\n"; } + else +#endif + { m_F6.add_vertex(vIP); } + } + + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: Add v2 to F1, F4, F6\n"; } + else +#endif + { + m_F1.add_vertex(v2); + m_F4.add_vertex(v2); + m_F6.add_vertex(v2); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v5, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v2--v5: add vIP to F1, F6 and F7 :: " << vIP[0] << "\n"; } + else +#endif + { + m_F1.add_vertex(vIP); + m_F6.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v2--v5: intersection ~ 0 :: add vIP to F4\n"; } + else +#endif + { m_F4.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v2--v5: intersection ~ 1 :: add vIP to F2\n"; } + else +#endif + { m_F2.add_vertex(vIP); } + } + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: Add v5 to F1, F2 and F6\n"; } + else +#endif + { + m_F1.add_vertex(v5); + m_F2.add_vertex(v5); + m_F6.add_vertex(v5); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v5, v7, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v5--v7: add vIP to F2, F6 and F7 :: " << vIP[2] << "\n"; } + else +#endif + { + m_F2.add_vertex(vIP); + m_F6.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v5--v7: intersection ~ 0 :: add vIP to F1\n"; } + else +#endif + { m_F1.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: v5--v7: intersection ~ 1 :: add vIP to F3\n"; } + else +#endif + { m_F3.add_vertex(vIP); } + } + } + + if (cuts == 2 && add_v7) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P2: Add v7 to F2, F3 and F6\n"; } + else +#endif + { + m_F2.add_vertex(v7); + m_F3.add_vertex(v7); + m_F6.add_vertex(v7); + } + + add_v7 = 0; + } + + } // end Path 2 + + // Path 5 + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v6, vIP, distIP)) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P5: v2--v6: add vIP to F4, F6 and F7 :: " << vIP[2] << "\n"; } + else +#endif + { + m_F4.add_vertex(vIP); + m_F6.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P5: v2--v6: intersection ~ 0 :: add vIP to F1\n"; } + else +#endif + { m_F1.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P5: v2--v6: intersection ~ 1 :: add vIP to F3\n"; } + else +#endif + { m_F3.add_vertex(vIP); } + } + } + + + // Path 3 + { int cuts(0); + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v3, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v0--v3: add vIP to F4, F5 and F7 :: " << vIP[2] << "\n"; } + else +#endif + { + m_F4.add_vertex(vIP); + m_F5.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v0--v3: intersection ~ 0 :: add vIP to F1\n"; } + else +#endif + { m_F1.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v0--v3: intersection ~ 1 :: add vIP to F3\n"; } + else +#endif + { m_F3.add_vertex(vIP); } + } + + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: Add v3 to F3, F4 and F5\n"; } + else +#endif + { + m_F3.add_vertex(v3); + m_F4.add_vertex(v3); + m_F5.add_vertex(v3); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v6, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v3--v6: add vIP to F3, F4 and F7 :: " << vIP[1] << "\n"; } + else +#endif + { + m_F3.add_vertex(vIP); + m_F4.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v3--v6: intersection ~ 0 :: add vIP to F5\n"; } + else +#endif + { m_F5.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v3--v6: intersection ~ 1 :: add vIP to F6\n"; } + else +#endif + { m_F6.add_vertex(vIP); } + } + + } + + if (cuts%2 == 0) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: Add v6 to F3, F4 and F6\n"; } + else +#endif + { + m_F3.add_vertex(v6); + m_F4.add_vertex(v6); + m_F6.add_vertex(v6); + } + } + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v6, v7, vIP, distIP)) { + + ++cuts; + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v6--v7: add vIP to F3, F6 and F7 :: " << vIP[0] << "\n"; } + else +#endif + { + m_F3.add_vertex(vIP); + m_F6.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v6--v7: intersection ~ 0 :: add vIP to F4\n"; } + else +#endif + { m_F4.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P3: v6--v7: intersection ~ 1 :: add vIP to F2\n"; } + else +#endif + { m_F2.add_vertex(vIP); } + } + } + + if (cuts == 2 && add_v7) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P1: Add v7 to F2, F3 and F6\n"; } + else +#endif + { + m_F2.add_vertex(v7); + m_F3.add_vertex(v7); + m_F6.add_vertex(v7); + } + add_v7 = 0; + } + + } // end Path 3 + + if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v4, vIP, distIP)) { + +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P6: v3--v4: add vIP to F3, F5 and F7 :: " << vIP[0] << "\n"; } + else +#endif + { + m_F3.add_vertex(vIP); + m_F5.add_vertex(vIP); + m_F7.add_vertex(vIP); + } + + if ( almostEqual(distIP, 0.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P6: v3--v4: intersection ~ 0 :: add vIP to F4\n"; } + else +#endif + { m_F4.add_vertex(vIP); } + + } else if ( almostEqual(distIP, 1.0)) { +#ifndef AMREX_USE_GPU + if (a_dry_run) { Print() << "P6: v3--v4: intersection ~ 1 :: add vIP to F2\n"; } + else +#endif + { m_F2.add_vertex(vIP); } + } + } +} + +#endif diff --git a/Source/EB/ERF_EBCutCell.cpp b/Source/EB/ERF_EBCutCell.cpp new file mode 100644 index 000000000..4c1a84881 --- /dev/null +++ b/Source/EB/ERF_EBCutCell.cpp @@ -0,0 +1,43 @@ +#include +#include + +using namespace amrex; + +#ifndef AMREX_USE_GPU +void +eb_cut_cell_:: +debug ( int const a_face ) { + + amrex::Print() << "\n\nDEBUG THIS " + << "--------------------------------------" + << "\nisCovered? " << isCovered() + << "\nisRegular? " << isRegular() + << "\nisSingleValued? " << isSingleValued() + << "\n"; + + if ( isCovered() || isRegular() ) { return; } + + amrex::RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2)); + amrex::RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2)); + + amrex::Print() << "\n" + << "lo: " << v0 << '\n' + << "hi: " << v7 << '\n' + << "p: " << m_eb_point << '\n' + << "n: " << m_eb_normal << '\n' + << "invert? " << m_invert << "\n\n"; + + int const dry_run(1); + amrex::Print() << "Edge intersections:\n"; + calc_edge_intersections(dry_run); + amrex::Print() << '\n'; + + if ( a_face == -1 || a_face == 1 ) { m_F1.report(1, v0); } + if ( a_face == -1 || a_face == 2 ) { m_F2.report(2, v0); } + if ( a_face == -1 || a_face == 3 ) { m_F3.report(3, v0); } + if ( a_face == -1 || a_face == 4 ) { m_F4.report(4, v0); } + if ( a_face == -1 || a_face == 5 ) { m_F5.report(5, v0); } + if ( a_face == -1 || a_face == 6 ) { m_F6.report(6, v0); } + if ( a_face == -1 || a_face == 7 ) { m_F7.report(7, v0); } +} +#endif diff --git a/Source/EB/ERF_EBIF.H b/Source/EB/ERF_EBIF.H index 8746d81fe..dbb7ea2e3 100644 --- a/Source/EB/ERF_EBIF.H +++ b/Source/EB/ERF_EBIF.H @@ -1,216 +1,49 @@ -#ifndef INCFLO_IF_LIST_ -#define INCFLO_IF_LIST_ +#ifndef ERF_EB_IF_H_ +#define ERF_EB_IF_H_ -#include -#include - -#include -#include +using namespace amrex; /******************************************************************************** * * * Union of a list (std::vector) of the same kind of implicit function. * * * ********************************************************************************/ - template -class UnionListIF -{ +class UnionListIF { public: - UnionListIF(const amrex::Vector& a_ifs) - : m_ifs(a_ifs), - empty(!a_ifs.empty()) - {} - - [[nodiscard]] bool is_empty() const - { - return empty; - } - - amrex::Real operator()(const amrex::RealArray& p) const - { - - // NOTE: this assumes that m_ifs is not empty - amrex::Real vmax = m_ifs[0](p); - for(int i = 1; i < m_ifs.size(); i++) - { - amrex::Real vcur = m_ifs[i](p); - if(vmax < vcur) - vmax = vcur; - } - - return vmax; - - // NOTE: this would have been nice, but for some reason it does not work :( - // even though according to https://en.cppreference.com/w/cpp/algorithm/max - // it should ... ? - //F & f_max = std::max( m_ifs, - // [&](const F & f1, const F & f2) { - // return f1(p) < f2(p); - // }); - //return f_max(p); - } - -private: - amrex::Vector m_ifs; - bool empty; -}; - -/******************************************************************************** - * * - * Conditional Implicit Functions => CIF * - * Can be "turned off" based on run-time parameters * - * * - ********************************************************************************/ - -/******************************************************************************** - * * - * Conditional implicit function: a "normal" implicit function is given the * - * additional property `active` (getter: `is_active()`, setter * - * `set_active(bool)`) allowing ConditionalUnion and ConditionalIntersection * - * operations. * - * * - ********************************************************************************/ - -template -class CIF : public F -{ - -public: - CIF(F&& f, bool a_active) - : F(f) - , m_active(a_active) - { - } - - ~CIF() - = default; - - CIF(const CIF& rhs) = default; - CIF(CIF&& rhs) noexcept = default; - CIF& operator=(const CIF& rhs) = default; - CIF& operator=(CIF&& rhs) noexcept = default; - - void set_active(bool a_active) - { - m_active = a_active; - } - [[nodiscard]] bool is_active() const - { - return m_active; - } - -private: - bool m_active; -}; - -/******************************************************************************** - * * - * Conditional Union of two (different) conditional implicit functions. * - * NOTE: at least one of the conditional implicit functions must be active. * - * * - ********************************************************************************/ - -template -class UnionCIF -{ -public: - // Assuming that one always active - UnionCIF(const F1& f1, const F2& f2) - : m_f1(f1) - , m_f2(f2) - , f1_active(f1.is_acive()) - , f2_active(f2.is_active()) + UnionListIF (const Vector & a_ifs) + : m_ifs(a_ifs) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(f1.is_active() || f2.is_active(), - "You must specify at least one active implicit function."); + m_empty = ! a_ifs.empty(); } - ~UnionCIF() - = default; + ~UnionListIF () {} - UnionCIF(const UnionCIF& rhs) = default; - UnionCIF(UnionCIF&& rhs) noexcept = default; - UnionCIF& operator=(const UnionCIF& rhs) = default; - UnionCIF& operator=(UnionCIF&& rhs) noexcept = default; + UnionListIF (const UnionListIF & rhs) = default; + UnionListIF (UnionListIF && rhs) = default; + UnionListIF & operator= (const UnionListIF & rhs) = default; + UnionListIF & operator= (UnionListIF && rhs) = default; - [[nodiscard]] bool is_active() const - { - return true; - } + bool is_empty () const {return m_empty;} - amrex::Real operator()(const amrex::RealArray& p) const - { + Real operator() (const RealArray & p) const { - if(!f1_active) - return m_f2(p); - if(!f2_active) - return m_f1(p); + // NOTE: this assumes that m_ifs is not empty + Real vmax=m_ifs[0](p); + for (int i=1; i< m_ifs.size(); i++) { + Real vcur = m_ifs[i](p); + if ( vmax < vcur ) vmax = vcur; + } - return std::max(m_f1(p), m_f2(p)); + return vmax; } private: - F1 m_f1; - F2 m_f2; - bool f1_active; - bool f2_active; -}; -/******************************************************************************** - * * - * Conditional Intersection of two (different) conditional implicit functions. * - * NOTE: at least one of the conditional implicit functions must be active. * - * * - ********************************************************************************/ - -template -class IntersectionCIF -{ - -public: - // Assuming that one always active - IntersectionCIF(const F1& f1, const F2& f2) - : m_f1(f1) - , m_f2(f2) - , f1_active(f1.is_acive()) - , f2_active(f2.is_active()) - { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(f1.is_active() || f2.is_active(), - "You must specify at least one active implicit function."); - } - - ~IntersectionCIF() - = default; - - IntersectionCIF(const IntersectionCIF& rhs) = default; - IntersectionCIF(IntersectionCIF&& rhs) noexcept = default; - IntersectionCIF& operator=(const IntersectionCIF& rhs) = default; - IntersectionCIF& operator=(IntersectionCIF&& rhs) noexcept = default; - - [[nodiscard]] bool is_active() const - { - return true; - } - - amrex::Real operator()(const amrex::RealArray& p) const - { - - if(!f1_active) - return m_f2(p); - if(!f2_active) - return m_f1(p); - - return std::min(m_f1(p), m_f2(p)); - } - -private: - F1 m_f1; - F2 m_f2; - bool f1_active; - bool f2_active; + Vector m_ifs; + bool m_empty; }; #endif diff --git a/Source/EB/ERF_TerrainIF.H b/Source/EB/ERF_EBIFTerrain.H similarity index 93% rename from Source/EB/ERF_TerrainIF.H rename to Source/EB/ERF_EBIFTerrain.H index 0e0d6ffea..175db9c87 100644 --- a/Source/EB/ERF_TerrainIF.H +++ b/Source/EB/ERF_EBIFTerrain.H @@ -14,7 +14,8 @@ class TerrainIF { public: - TerrainIF (amrex::FArrayBox& a_z_terrain, amrex::Geometry& a_geom, amrex::Gpu::DeviceVector& a_dz_stretched) + TerrainIF (amrex::FArrayBox const& a_z_terrain, amrex::Geometry const& a_geom, + amrex::Gpu::DeviceVector& a_dz_stretched) : terr_arr(a_z_terrain.const_array()), dz_s(a_dz_stretched.data()), dx(a_geom.CellSize(0)), diff --git a/Source/EB/ERF_EBPolygon.H b/Source/EB/ERF_EBPolygon.H new file mode 100644 index 000000000..8bc3c5edb --- /dev/null +++ b/Source/EB/ERF_EBPolygon.H @@ -0,0 +1,196 @@ +#ifndef ERF_EB_POLYGON_H_ +#define ERF_EB_POLYGON_H_ + +#include +#include + +#include "ERF_Constants.H" + +class polygon_ { + + public: + + AMREX_GPU_HOST_DEVICE + polygon_ ( amrex::RealVect a_point, + amrex::RealVect a_normal ) + : m_cell_face(1) + , m_eb_point(a_point) + , m_eb_normal(a_normal) + , m_defined(0) + , m_vertices(0) + , m_area(0.) + , m_sorted(0) + , m_vertex({amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.), + amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.)}) + , m_zdir(0.) + {} + + AMREX_GPU_HOST_DEVICE + polygon_ ( ) + : m_cell_face(0) + , m_eb_point(0.) + , m_eb_normal(0.) + , m_defined(0) + , m_vertices(0) + , m_area(0.) + , m_sorted(0) + , m_vertex({amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.), + amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.)}) + , m_zdir(0.) + {} + + AMREX_GPU_HOST_DEVICE + void add_vertex ( amrex::RealVect const& a_v ) { + AMREX_ASSERT( m_vertices < m_max_vertices ); + m_vertex[m_vertices] = a_v; + m_vertices++; + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void set_area ( amrex::Real const& a_area ) { m_area = a_area; } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void define () { + + AMREX_ALWAYS_ASSERT( m_defined == 0 ); // TODO ---------------------------- remove ALWAYS + + m_defined = 1; + + // We need at least 3 faces. + if (m_vertices < 3) { return; } + + // Check to see if the vertices of the face are inside or outside of the plane. + // If they are outside, this face doesn't belong to the volume. + if ( m_cell_face ) { + for ( int i(0); i= 0.) { + } else { + } + } + } + + // Calculate the centroid of the polygon. + amrex::RealVect centroid = get_centroid(); + + // Shift the vertices relative to the centroid + amrex::Array vertex_cent; + for ( int i(0); i= 0.) ? 0.0 : 2.0*PI); + } + + // Sort counter clockwise based on theta. + for (int i(0); i m_theta[j+1] ) { + amrex::Swap(m_theta[j], m_theta[j+1]); + amrex::Swap(m_vertex[j], m_vertex[j+1]); + amrex::Swap(vertex_cent[j], vertex_cent[j+1]); + } + } // j-loop + } // i-loop + m_sorted = 1; + + // Compute areas of triangles + for (int i(0); i 0. ); // <------------------------------- TODO remove ALWAYS + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + int ok ( ) const noexcept + { return ((m_area > 0. || (m_area == 0. && m_defined == 1)) ? 1 : 0); } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real area ( ) const noexcept { + AMREX_ALWAYS_ASSERT( ok() ); // <-------------------------------------- TODO remove ALWAYS + return m_area; + } + + // Distance from polygon to a_point + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real distance ( amrex::RealVect const& a_point ) const noexcept { + AMREX_ALWAYS_ASSERT( m_defined == 1 ); // <--------------------------- TODO remove ALWAYS + amrex::RealVect x0 = a_point - m_vertex[0]; + return amrex::Math::abs(x0.dotProduct(m_zdir)); + } + + // Calculate the centroid of the polygon. + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::RealVect get_centroid ( ) const noexcept { + amrex::RealVect cent(0.); + for ( int i(0); i(m_vertices); + } + +#ifndef AMREX_USE_GPU + void report ( int const a_id, amrex::RealVect a_v0 ) +#else + void report ( int const /*a_id*/, amrex::RealVect /*a_v0*/ ) +#endif + { + +#ifndef AMREX_USE_GPU + amrex::RealVect centroid = get_centroid(); + amrex::Print() << "Face " << a_id + << " -------------------------------------------\n" + << "\nok? " << (ok() ? "yes" : "no") << "\n\n"; + for (int i(0); i m_vertex; + + amrex::GpuArray m_theta; + + amrex::RealVect m_zdir; +}; + + +#endif diff --git a/Source/EB/ERF_Redistribute.cpp b/Source/EB/ERF_EBRedistribute.cpp similarity index 78% rename from Source/EB/ERF_Redistribute.cpp rename to Source/EB/ERF_EBRedistribute.cpp index 7136ab10d..a2a1f01af 100644 --- a/Source/EB/ERF_Redistribute.cpp +++ b/Source/EB/ERF_EBRedistribute.cpp @@ -6,12 +6,12 @@ using namespace amrex; void -ERF::redistribute_term ( int lev, - MultiFab& result, - MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? - MultiFab const& state, - BCRec const* bc, // this is bc for the state (needed for SRD slopes) - Real const local_dt) +ERF::redistribute_term (int lev, int ncomp, + MultiFab& result, + MultiFab& result_tmp, + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const local_dt) { // ************************************************************************ // Redistribute result_tmp and pass out result @@ -25,17 +25,17 @@ ERF::redistribute_term ( int lev, #endif for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - redistribute_term(mfi, lev, result, result_tmp, state, bc, local_dt); + redistribute_term(mfi, lev, ncomp, result, result_tmp, state, bc, local_dt); } } void -ERF::redistribute_term ( MFIter const& mfi, int lev, - MultiFab& result, - MultiFab& result_tmp, - MultiFab const& state, - BCRec const* bc, // this is bc for the state (needed for SRD slopes) - Real const local_dt) +ERF::redistribute_term (MFIter const& mfi, int lev, int ncomp, + MultiFab& result, + MultiFab& result_tmp, + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const local_dt) { AMREX_ASSERT(result.nComp() == state.nComp()); @@ -50,7 +50,6 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, Array4 out = result.array(mfi); Array4 in = result_tmp.array(mfi); - int ncomp = result.nComp(); if (!regular && !covered) { diff --git a/Source/EB/ERF_EBRegular.cpp b/Source/EB/ERF_EBRegular.cpp deleted file mode 100644 index 710c7cd13..000000000 --- a/Source/EB/ERF_EBRegular.cpp +++ /dev/null @@ -1,14 +0,0 @@ -#include -#include - -#include -#include - -using namespace amrex; - -void ERF::make_eb_regular() -{ - EB2::AllRegularIF my_regular; - auto gshop = EB2::makeShop(my_regular); - EB2::Build(gshop, geom.back(), 0, 100); -} diff --git a/Source/EB/ERF_EBUtils.H b/Source/EB/ERF_EBUtils.H new file mode 100644 index 000000000..509752c09 --- /dev/null +++ b/Source/EB/ERF_EBUtils.H @@ -0,0 +1,60 @@ +#ifndef ERF_EB_UTIL_H_ +#define ERF_EB_UTIL_H_ + +#include + +namespace utils { + +// Harmonic average +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real avg_h ( amrex::Real XXXm, amrex::Real XXXp) { + if ( XXXm + XXXp < std::numeric_limits::min()) { return 0.; } + else { return ((XXXm * XXXp) / (0.5*(XXXm + XXXp))); } +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real trD (int i, int j, int k, + amrex::GpuArray const& dx, + amrex::Array4 const& vel_x, + amrex::Array4 const& vel_y, + amrex::Array4 const& vel_z) +{ + return (vel_x(i+1,j,k)-vel_x(i,j,k))/dx[0] + + (vel_y(i,j+1,k)-vel_y(i,j,k))/dx[1] + + (vel_z(i,j,k+1)-vel_z(i,j,k))/dx[2]; +} + + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int +intersect_plane_edge ( amrex::RealVect const& a_plane_point, + amrex::RealVect const& a_plane_normal, + amrex::RealVect const& a_edge_point0, + amrex::RealVect const& a_edge_point1, + amrex::RealVect& a_intersection_point, + amrex::Real& a_intersection_dist ) +{ + amrex::RealVect const edge(a_edge_point1 - a_edge_point0); + amrex::Real const edge_length = edge.vectorLength(); + + AMREX_ALWAYS_ASSERT(edge_length > 0.); + + amrex::RealVect edge_normal = edge / edge_length; + + amrex::Real np_dot_ne = a_plane_normal.dotProduct(edge_normal); + + if ( amrex::Math::abs(np_dot_ne) < std::numeric_limits::min() ) + { return 0; } + + a_intersection_dist = a_plane_normal.dotProduct(a_plane_point - a_edge_point0); + a_intersection_dist /= np_dot_ne; + + a_intersection_point = a_edge_point0 + a_intersection_dist*edge_normal; + + if (0. <= a_intersection_dist && a_intersection_dist <= edge_length) { return 1;} + + return 0; +} + +} +#endif diff --git a/Source/EB/ERF_WriteEBSurface.cpp b/Source/EB/ERF_EBWriteSurface.cpp similarity index 100% rename from Source/EB/ERF_WriteEBSurface.cpp rename to Source/EB/ERF_EBWriteSurface.cpp diff --git a/Source/EB/ERF_InitEB.cpp b/Source/EB/ERF_InitEB.cpp deleted file mode 100644 index 9230450a4..000000000 --- a/Source/EB/ERF_InitEB.cpp +++ /dev/null @@ -1,63 +0,0 @@ -#include -#include - -#include -#include -#include -#include - -using namespace amrex; - -void ERF::MakeEBGeometry() -{ - /****************************************************************************** - * ERF.geometry= specifies the EB geometry. can be either * - * box or terrain */ - - ParmParse pp("eb2"); - - std::string geom_type = "terrain"; - // pp.query("geometry", geom_type); - - /****************************************************************************** - * * - * CONSTRUCT EB * - * * - ******************************************************************************/ - - int lev = 0; - - int max_coarsening_level; - if (solverChoice.anelastic[0] == 1) { - max_coarsening_level = 100; - } else { - max_coarsening_level = 0; - } - - if (geom_type == "terrain") { - amrex::Print() << "\n Building EB geometry based on idealized terrain." << std::endl; - Real dummy_time = 0.0; - - Box bx(surroundingNodes(Geom(lev).Domain())); bx.grow(3); - FArrayBox terrain_fab(makeSlab(bx,2,0),1); - - prob->init_terrain_surface(Geom(lev), terrain_fab, dummy_time); - - TerrainIF ebterrain(terrain_fab, Geom(lev), stretched_dz_d[lev]); - - auto gshop = EB2::makeShop(ebterrain); - - EB2::Build(gshop, geom.back(), max_level, max_level+max_coarsening_level); - - } else if(geom_type == "box") { - amrex::Print() << "\n Building box geometry." << std::endl; - make_eb_box(); - - } else { - - amrex::Print() << "\n No EB geometry declared in inputs => " - << " Will build all regular geometry." << std::endl; - make_eb_regular(); - } - amrex::Print() << "Done making the geometry ebfactory.\n" << std::endl; -} diff --git a/Source/EB/Make.package b/Source/EB/Make.package index deddc170b..bbdc9dca6 100644 --- a/Source/EB/Make.package +++ b/Source/EB/Make.package @@ -1,9 +1,18 @@ +CEXE_headers += ERF_EB.H +CEXE_sources += ERF_EB.cpp + +CEXE_headers += ERF_EBPolygon.H +CEXE_headers += ERF_EBUtils.H + +CEXE_headers += ERF_EBCutCell.H +CEXE_sources += ERF_EBCutCell.cpp + +CEXE_headers += ERF_EBAux.H +CEXE_sources += ERF_EBAux.cpp -CEXE_sources += main.cpp -CEXE_sources += ERF_InitEB.cpp -CEXE_sources += ERF_EBBox.cpp -CEXE_sources += ERF_EBRegular.cpp -CEXE_sources += ERF_Redistribute.cpp -CEXE_sources += ERF_WriteEBSurface.cpp CEXE_headers += ERF_EBIF.H -CEXE_headers += ERF_TerrainIF.H +CEXE_headers += ERF_EBIFTerrain.H + +CEXE_sources += ERF_EBBox.cpp +CEXE_sources += ERF_EBRedistribute.cpp +CEXE_sources += ERF_EBWriteSurface.cpp diff --git a/Source/ERF.H b/Source/ERF.H index eb693312a..7bf1a1ea0 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -24,6 +24,8 @@ #include #include +#include + #ifdef ERF_USE_FFT #include #endif @@ -474,18 +476,18 @@ public: void MakeEBGeometry (); void make_eb_box (); void make_eb_regular (); - void redistribute_term ( int lev, - amrex::MultiFab& result, - amrex::MultiFab& result_tmp, - amrex::MultiFab const& state, - amrex::BCRec const* bc, - amrex::Real const dt); - void redistribute_term ( amrex::MFIter const& mfi, int lev, - amrex::MultiFab& result, - amrex::MultiFab& result_tmp, - amrex::MultiFab const& state, - amrex::BCRec const* bc, - amrex::Real const dt); + void redistribute_term (int lev, int ncomp, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); + void redistribute_term (amrex::MFIter const& mfi, int lev, int ncomp, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); // more flexible version of AverageDown() that lets you average down across multiple levels void AverageDownTo (int crse_lev, int scomp, int ncomp); // NOLINT diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 623fab5b1..de78590a6 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -360,12 +360,21 @@ ERF::ERF_shared () // We will create each of these in MakeNewLevel.../RemakeLevel m_factory.resize(max_level+1); + // + // Construct the EB data structures and store in a separate class + // // This is needed before initializing level MultiFabs if ( solverChoice.terrain_type == TerrainType::EB || solverChoice.terrain_type == TerrainType::ImmersedForcing) { + int lev = 0; Real dummy_time = 0.0; + Box terrain_bx(surroundingNodes(geom[lev].Domain())); terrain_bx.grow(3); + FArrayBox terrain_fab(makeSlab(terrain_bx,2,0),1); + prob->init_terrain_surface(geom[lev], terrain_fab, dummy_time); + amrex::Print() << "MAKING EB GEOMETRY " << std::endl; - MakeEBGeometry(); + eb_ eb(geom[lev], terrain_fab, stretched_dz_d[lev], solverChoice.anelastic[lev]); + // MakeEBGeometry(); } } diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H index 1ecd4a8a2..67619790e 100644 --- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H @@ -148,14 +148,63 @@ } // not moving terrain } // mfi - if (solverChoice.terrain_type == TerrainType::EB) { - // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF - // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level)); - // dUdt_tmp.setVal(0.); - // dUdt_tmp.FillBoundary(fine_geom.periodicity()); - // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both? - // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp, - // S_sum[IntVars::cons], bc_ptr_h, slow_dt); + if (solverChoice.terrain_type == TerrainType::EB) + { + // Redistribute cons states (cell-centered) + + MultiFab dUdt_tmp(ba, dm, F_slow[IntVars::cons].nComp(), F_slow[IntVars::cons].nGrow(), MFInfo(), Factory(level)); + dUdt_tmp.setVal(0.); + MultiFab::Copy(dUdt_tmp, F_slow[IntVars::cons], 0, 0, F_slow[IntVars::cons].nComp(), F_slow[IntVars::cons].nGrow()); + dUdt_tmp.FillBoundary(fine_geom.periodicity()); + dUdt_tmp.setDomainBndry(1.234e10, 0, ncomp_fast[IntVars::cons], fine_geom); + + BCRec const* bc_ptr_d = domain_bcs_type_d.data(); + S_old[IntVars::cons].FillBoundary(fine_geom.periodicity()); + S_old[IntVars::cons].setDomainBndry(1.234e10, 0, ncomp_fast[IntVars::cons], fine_geom); + + // Update F_slow by Redistribution. + redistribute_term ( level, ncomp_fast[IntVars::cons], F_slow[IntVars::cons], dUdt_tmp, + S_old[IntVars::cons], bc_ptr_d, slow_dt); + + // Update state using the updated F_slow. (NOTE: redistribute_term returns RHS not state variables.) + + for ( MFIter mfi(S_sum[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box bx = mfi.tilebox(); + Vector > ssum_h(ncomp_fast[IntVars::cons]); + Vector > sold_h(ncomp_fast[IntVars::cons]); + Vector > fslow_h(ncomp_fast[IntVars::cons]); + const Array4& dJ_old = detJ_cc[level]->const_array(mfi); + + for (int i = 0; i < ncomp_fast[IntVars::cons]; ++i) { + ssum_h[i] = S_sum[i].array(mfi); + sold_h[i] = S_old[i].array(mfi); + fslow_h[i] = F_slow[i].array(mfi); + } + + Gpu::AsyncVector > sold_d(ncomp_fast[IntVars::cons]); + Gpu::AsyncVector > ssum_d(ncomp_fast[IntVars::cons]); + Gpu::AsyncVector > fslow_d(ncomp_fast[IntVars::cons]); + + Gpu::copy(Gpu::hostToDevice, sold_h.begin(), sold_h.end(), sold_d.begin()); + Gpu::copy(Gpu::hostToDevice, ssum_h.begin(), ssum_h.end(), ssum_d.begin()); + Gpu::copy(Gpu::hostToDevice, fslow_h.begin(), fslow_h.end(), fslow_d.begin()); + + Array4* sold = sold_d.dataPtr(); + Array4* ssum = ssum_d.dataPtr(); + Array4* fslow = fslow_d.dataPtr(); + + ParallelFor(bx, ncomp_fast[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) + { + const int n = scomp_fast[IntVars::cons] + nn; + if (dJ_old(i,j,k) > 0.0) { + ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt * + ( fslow[IntVars::cons](i,j,k,n) ); + } + }); + } + + // Redistribute momentum states (face-centered) will be added here. } // EB } // omp diff --git a/Source/main.cpp b/Source/main.cpp index 4cd2ee0b3..81872ed41 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -4,7 +4,6 @@ #include #include -//#include "IO.H" #include "ERF.H" #ifdef ERF_USE_WW3_COUPLING From 7c969ae071aea6532317aac91dcad9fcc5bf5625 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 18 Jan 2025 15:45:13 -0800 Subject: [PATCH 2/3] fix several errors related to multilevel (#2076) --- Source/ERF.cpp | 14 +------------- Source/ERF_MakeNewLevel.cpp | 25 +++++++------------------ Source/IO/ERF_Plotfile.cpp | 14 ++++++++++++-- 3 files changed, 20 insertions(+), 33 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index de78590a6..d26fb522e 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -715,7 +715,7 @@ ERF::InitData_post () // Create the physbc objects for {cons, u, v, w, base state} // We fill the additional base state ghost cells just in case we have read the old format - for (int lev(0); lev <= max_level; ++lev) { + for (int lev(0); lev <= finest_level; ++lev) { make_physbcs(lev); (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); } @@ -1249,7 +1249,6 @@ ERF::initializeMicrophysics (const int& a_nlevsmax /*!< number of AMR levels */) return; } - #ifdef ERF_USE_WINDFARM void ERF::initializeWindFarm(const int& a_nlevsmax/*!< number of AMR levels */ ) @@ -1261,17 +1260,6 @@ ERF::initializeWindFarm(const int& a_nlevsmax/*!< number of AMR levels */ ) void ERF::restart () { - // TODO: This could be deleted since ba/dm are not created yet? - for (int lev = 0; lev <= finest_level; ++lev) - { - auto& lev_new = vars_new[lev]; - auto& lev_old = vars_old[lev]; - lev_new[Vars::cons].setVal(0.); lev_old[Vars::cons].setVal(0.); - lev_new[Vars::xvel].setVal(0.); lev_old[Vars::xvel].setVal(0.); - lev_new[Vars::yvel].setVal(0.); lev_old[Vars::yvel].setVal(0.); - lev_new[Vars::zvel].setVal(0.); lev_old[Vars::zvel].setVal(0.); - } - #ifdef ERF_USE_NETCDF if (restart_type == "netcdf") { ReadNCCheckpointFile(); diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index 6cf82ee10..6cc8f6ad5 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -263,21 +263,15 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, qmoist[lev][mvar] = micro->Get_Qmoist_Ptr(lev,mvar); } + // ***************************************************************************************************** + // Initialize the boundary conditions (after initializing the terrain but before calling + // initHSE or FillCoarsePatch) + // ***************************************************************************************************** + make_physbcs(lev); + // ******************************************************************************************** - // Update the base state at this level by interpolation from coarser level + // Update the base state at this level by interpolation from coarser level (inside initHSE) // ******************************************************************************************** - // - // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab - // have been pre-filled - this includes ghost cells both inside and outside - // the domain - // - InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(), - IntVect(0,0,0), // do not fill ghost cells outside the domain - base_state[lev-1], 0, 0, base_state[lev].nComp(), - geom[lev-1], geom[lev], - refRatio(lev-1), &cell_cons_interp, - domain_bcs_type, BCVars::cons_bc); - initHSE(lev); // ******************************************************************************************** @@ -285,11 +279,6 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // ******************************************************************************************** update_diffusive_arrays(lev, ba, dm); - // ***************************************************************************************************** - // Initialize the boundary conditions (after initializing the terrain but before calling FillCoarsePatch - // ***************************************************************************************************** - make_physbcs(lev); - // ******************************************************************************************** // Fill data at the new level by interpolation from the coarser level // Note that internal to FillCoarsePatch we will convert velocity to momentum, diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 408143137..e66e69ac5 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1512,14 +1512,24 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data()); } - // Do piecewise interpolation of mf into mf2 + // + // We need to make a temporary that is the size of ncomp_mf + // in order to not get an out of bounds error + // even though the values will not be used + // + Vector temp_domain_bcs_type; + temp_domain_bcs_type.resize(ncomp_mf); + + // + // Do piecewise constant interpolation of mf into mf2 + // for (int lev = 1; lev <= finest_level; ++lev) { Interpolater* mapper_c = &pc_interp; InterpFromCoarseLevel(mf2[lev], t_new[lev], mf[lev], 0, 0, ncomp_mf, geom[lev], g2[lev], null_bc_for_fill, 0, null_bc_for_fill, 0, - r2[lev-1], mapper_c, domain_bcs_type, 0); + r2[lev-1], mapper_c, temp_domain_bcs_type, 0); } // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile From 2d3ed586ec2ff8d32e22987b93ae05ad0d35314e Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 21 Jan 2025 11:04:44 -0800 Subject: [PATCH 3/3] Fix WPS w/ Low Mem & Restart (#2075) * use mask array. * fix compile. * moist flag. * Matched benchmark. --- Source/IO/ERF_ReadFromWRFBdy.cpp | 159 +++++++++++------- .../Initialization/ERF_InitFromWRFInput.cpp | 24 +-- 2 files changed, 112 insertions(+), 71 deletions(-) diff --git a/Source/IO/ERF_ReadFromWRFBdy.cpp b/Source/IO/ERF_ReadFromWRFBdy.cpp index 820c1a23f..63e06b482 100644 --- a/Source/IO/ERF_ReadFromWRFBdy.cpp +++ b/Source/IO/ERF_ReadFromWRFBdy.cpp @@ -518,8 +518,14 @@ convert_wrfbdy_data (const Box& domain, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& cons, + const Geometry& geom, const bool& use_moist) { + // Owner masks for parallel reduce sum + std::unique_ptr mask_c = OwnerMask(cons, geom.periodicity()); + std::unique_ptr mask_u = OwnerMask(xvel, geom.periodicity()); + std::unique_ptr mask_v = OwnerMask(yvel, geom.periodicity()); + int ntimes = bdy_data.size(); for (int nt = 0; nt < ntimes; nt++) { @@ -533,43 +539,66 @@ convert_wrfbdy_data (const Box& domain, for ( MFIter mfi(cons); mfi.isValid(); ++mfi ) { + Box tbx = mfi.tilebox(); + Box xbx = mfi.nodaltilebox(0); + Box ybx = mfi.nodaltilebox(1); + const Box& bx_u = (xbx & bdy_data[nt][WRFBdyVars::U].box()); + const Box& bx_v = (ybx & bdy_data[nt][WRFBdyVars::V].box()); + const Box& bx_t = (tbx & bdy_data[nt][WRFBdyVars::T].box()); + const Box& bx_qv = (tbx & bdy_data[nt][WRFBdyVars::QV].box()); + + // TMP BDY data + Array4 bdy_u_tmp = bdy_data_tmp[WRFBdyVars::U].array(); // This is x-face-centered + Array4 bdy_v_tmp = bdy_data_tmp[WRFBdyVars::V].array(); // This is y-face-centered + Array4 bdy_t_tmp = bdy_data_tmp[WRFBdyVars::T].array(); // This is cell-centered + Array4 bdy_qv_tmp = bdy_data_tmp[WRFBdyVars::QV].array(); // This is cell-centered + + // Mask data + const Array4& mask_c_arr = mask_c->const_array(mfi); + const Array4& mask_u_arr = mask_u->const_array(mfi); + const Array4& mask_v_arr = mask_v->const_array(mfi); + if (nt==0) { - const FArrayBox& xvel_fab = xvel[mfi]; - const FArrayBox& yvel_fab = yvel[mfi]; - const FArrayBox& cons_fab = cons[mfi]; - - bdy_data_tmp[WRFBdyVars::U].template copy(xvel_fab,0,0,1); - bdy_data_tmp[WRFBdyVars::V].template copy(yvel_fab,0,0,1); - bdy_data_tmp[WRFBdyVars::T].template copy(cons_fab,RhoTheta_comp,0,1); - bdy_data_tmp[WRFBdyVars::T].template divide(cons_fab,bdy_data[0][WRFBdyVars::T].box(),Rho_comp,0,1); - - if (use_moist) { - bdy_data_tmp[WRFBdyVars::QV].template copy(cons_fab,RhoQ1_comp,0,1); - bdy_data_tmp[WRFBdyVars::QV].template divide(cons_fab,bdy_data[0][WRFBdyVars::QV].box(),Rho_comp,0,1); - } else { - bdy_data_tmp[WRFBdyVars::QV].template setVal(0.); - } - } else { - Box vbx = mfi.validbox(); - Box xbx = convert(vbx,IntVect(1,0,0)); - Box ybx = convert(vbx,IntVect(0,1,0)); - const Box& bx_u = (xbx & bdy_data[nt][WRFBdyVars::U].box()); - const Box& bx_v = (ybx & bdy_data[nt][WRFBdyVars::V].box()); - const Box& bx_t = (vbx & bdy_data[nt][WRFBdyVars::T].box()); - const Box& bx_qv = (vbx & bdy_data[nt][WRFBdyVars::QV].box()); + const Array4& xvel_fab = xvel.const_array(mfi); + const Array4& yvel_fab = yvel.const_array(mfi); + const Array4& cons_fab = cons.const_array(mfi); - // BDY data - Array4 bdy_u_arr = bdy_data_tmp[WRFBdyVars::U].array(); // This is x-face-centered - Array4 bdy_v_arr = bdy_data_tmp[WRFBdyVars::V].array(); // This is y-face-centered - Array4 bdy_t_arr = bdy_data_tmp[WRFBdyVars::T].array(); // This is cell-centered - Array4 bdy_qv_arr = bdy_data_tmp[WRFBdyVars::QV].array(); // This is cell-centered - Array4 mu_arr = bdy_data[nt][WRFBdyVars::MU].array(); // This is cell-centered + // Define u velocity + ParallelFor(bx_u, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (mask_u_arr(i,j,k)) { bdy_u_tmp(i,j,k) = xvel_fab(i,j,k); } + }); + + // Define v velocity + ParallelFor(bx_v, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (mask_v_arr(i,j,k)) { bdy_v_tmp(i,j,k) = yvel_fab(i,j,k); } + }); + // Define Qv & T + ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (mask_c_arr(i,j,k)) { + bdy_t_tmp(i,j,k) = cons_fab(i,j,k,RhoTheta_comp)/cons_fab(i,j,k,Rho_comp); + bdy_qv_tmp(i,j,k) = (use_moist) ? cons_fab(i,j,k,RhoQ1_comp )/cons_fab(i,j,k,Rho_comp) : + 0.0; + } + }); + + } else { // Populated from read wrfinput Array4 c1h_arr = mf_C1H.const_array(mfi); Array4 c2h_arr = mf_C2H.const_array(mfi); Array4 mub_arr = mf_MUB.const_array(mfi); + // BDY data + Array4 bdy_u_arr = bdy_data[nt][WRFBdyVars::U].array(); // This is x-face-centered + Array4 bdy_v_arr = bdy_data[nt][WRFBdyVars::V].array(); // This is y-face-centered + Array4 bdy_t_arr = bdy_data[nt][WRFBdyVars::T].array(); // This is cell-centered + Array4 bdy_qv_arr = bdy_data[nt][WRFBdyVars::QV].array(); // This is cell-centered + Array4 mu_arr = bdy_data[nt][WRFBdyVars::MU].array(); // This is cell-centered + + // Bounds limiting int ilo = domain.smallEnd()[0]; int ihi = domain.bigEnd()[0]; int jlo = domain.smallEnd()[1]; @@ -578,56 +607,64 @@ convert_wrfbdy_data (const Box& domain, // Define u velocity ParallelFor(bx_u, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu; - if (i == ilo) { - xmu = mu_arr(i,j,0) + mub_arr(i,j,0); - } else if (i > ihi) { - xmu = mu_arr(i-1,j,0) + mub_arr(i-1,j,0); - } else { - xmu = ( mu_arr(i,j,0) + mu_arr(i-1,j,0) - + mub_arr(i,j,0) + mub_arr(i-1,j,0)) * 0.5; + if (mask_u_arr(i,j,k)) { + Real xmu; + if (i == ilo) { + xmu = mu_arr(i,j,0) + mub_arr(i,j,0); + } else if (i > ihi) { + xmu = mu_arr(i-1,j,0) + mub_arr(i-1,j,0); + } else { + xmu = ( mu_arr(i,j,0) + mu_arr(i-1,j,0) + + mub_arr(i,j,0) + mub_arr(i-1,j,0)) * 0.5; + } + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy = bdy_u_arr(i,j,k) / xmu_mult; + bdy_u_tmp(i,j,k) = new_bdy; } - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy = bdy_u_arr(i,j,k) / xmu_mult; - bdy_u_arr(i,j,k) = new_bdy; }); // Define v velocity ParallelFor(bx_v, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu; - if (j == jlo) { - xmu = mu_arr(i,j,0) + mub_arr(i,j,0); - } else if (j > jhi) { - xmu = mu_arr(i,j-1,0) + mub_arr(i,j-1,0); - } else { - xmu = ( mu_arr(i,j,0) + mu_arr(i,j-1,0) - + mub_arr(i,j,0) + mub_arr(i,j-1,0) ) * 0.5; + if (mask_v_arr(i,j,k)) { + Real xmu; + if (j == jlo) { + xmu = mu_arr(i,j,0) + mub_arr(i,j,0); + } else if (j > jhi) { + xmu = mu_arr(i,j-1,0) + mub_arr(i,j-1,0); + } else { + xmu = ( mu_arr(i,j,0) + mu_arr(i,j-1,0) + + mub_arr(i,j,0) + mub_arr(i,j-1,0) ) * 0.5; + } + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy = bdy_v_arr(i,j,k) / xmu_mult; + bdy_v_tmp(i,j,k) = new_bdy; } - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy = bdy_v_arr(i,j,k) / xmu_mult; - bdy_v_arr(i,j,k) = new_bdy; }); // Convert perturbational moist pot. temp. (Th_m) to dry pot. temp. (Th_d) Real theta_ref = 300.; ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy_Th = bdy_t_arr(i,j,k) / xmu_mult + theta_ref; - Real qv_fac = (1. + bdy_qv_arr(i,j,k) / 0.622 / xmu_mult); - new_bdy_Th /= qv_fac; - bdy_t_arr(i,j,k) = new_bdy_Th; + if (mask_c_arr(i,j,k)) { + Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy_Th = bdy_t_arr(i,j,k) / xmu_mult + theta_ref; + Real qv_fac = (1. + bdy_qv_arr(i,j,k) / 0.622 / xmu_mult); + new_bdy_Th /= qv_fac; + bdy_t_tmp(i,j,k) = new_bdy_Th; + } }); // Define Qv ParallelFor(bx_qv, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy_QV = bdy_qv_arr(i,j,k) / xmu_mult; - bdy_qv_arr(i,j,k) = (use_moist) ? new_bdy_QV : 0.; + if (mask_c_arr(i,j,k)) { + Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy_QV = bdy_qv_arr(i,j,k) / xmu_mult; + bdy_qv_tmp(i,j,k) = (use_moist) ? new_bdy_QV : 0.; + } }); } // nt ==0 diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index c530ace0a..f66424f44 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -40,6 +40,7 @@ convert_wrfbdy_data (const Box& domain, const MultiFab& xvel, const MultiFab& yvel, const MultiFab& cons, + const Geometry& geom, const bool& use_moist); void @@ -204,8 +205,8 @@ ERF::init_from_wrfinput (int lev) // Define fabs for holding the initial data int icomp = 0; bool mult_rho = false; + FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; FArrayBox* cur_fab; - if (NC_names[ivar] == "U") { cur_fab = &lev_new[Vars::xvel][mfi]; } else if (NC_names[ivar] == "V") { @@ -233,12 +234,9 @@ ERF::init_from_wrfinput (int lev) if (n_qstate > 3) { icomp = RhoQ4_comp; } } - FArrayBox &cons_fab = lev_new[Vars::cons][mfi]; - cur_fab->template copy(var_fab, 0, icomp, 1); - if (mult_rho) { cur_fab->template mult(cons_fab, Rho_comp, icomp, 1); } - + var_fab.clear(); } // mfi } // valid var (not rho) @@ -426,6 +424,8 @@ ERF::init_from_wrfinput (int lev) msf_fab.template copy(var_fab); } } + + var_fab.clear(); Print() << " DONE\n"; } // ivar Print() << "\n"; @@ -523,16 +523,20 @@ ERF::init_from_wrfinput (int lev) convert_wrfbdy_data(domain,bdy_data_xlo, mf_MUB, mf_C1H, mf_C2H, - lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], use_moist); + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); convert_wrfbdy_data(domain,bdy_data_xhi, mf_MUB, mf_C1H, mf_C2H, - lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], use_moist); + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); convert_wrfbdy_data(domain,bdy_data_ylo, mf_MUB, mf_C1H, mf_C2H, - lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], use_moist); + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); convert_wrfbdy_data(domain,bdy_data_yhi, mf_MUB, mf_C1H, mf_C2H, - lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], use_moist); + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); } // init_type == Real && lev == 0 // Start at the earliest time (read_from_wrfbdy) @@ -570,7 +574,7 @@ init_base_state_from_wrfinput (const Box& domain, const auto& dom_lo = lbound(domain); const auto& dom_hi = ubound(domain); - for ( MFIter mfi(cons, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + for ( MFIter mfi(p_hse, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Box tbx = mfi.tilebox(); Box gtbx = mfi.growntilebox();