diff --git a/src/boundary_conditions/field_boundary_fill/BoundaryPlane.cpp b/src/boundary_conditions/field_boundary_fill/BoundaryPlane.cpp index b66e72b872..1eb13a161c 100644 --- a/src/boundary_conditions/field_boundary_fill/BoundaryPlane.cpp +++ b/src/boundary_conditions/field_boundary_fill/BoundaryPlane.cpp @@ -499,8 +499,7 @@ void BoundaryPlane::initialize_data() if (m_repo.field_exists(fname)) { auto& fld = m_repo.get_field(fname); if (m_io_mode == io_mode::input) { - fld.register_fill_patch_op( - m_mesh, m_time, *this); + fld.add_fill_patch_op(m_mesh, m_time, *this); } m_fields.emplace_back(&fld); } else { diff --git a/src/boundary_conditions/field_boundary_fill/CMakeLists.txt b/src/boundary_conditions/field_boundary_fill/CMakeLists.txt index 29bee7a1da..0bace6af79 100644 --- a/src/boundary_conditions/field_boundary_fill/CMakeLists.txt +++ b/src/boundary_conditions/field_boundary_fill/CMakeLists.txt @@ -4,6 +4,8 @@ target_sources(${kynema_sgf_lib_name} OceanWavesBoundary.cpp PlaneFillInflow.cpp BoundaryPlane.cpp + FillFlather.cpp + Flather.cpp FillMPL.cpp ModulatedPowerLaw.cpp ) \ No newline at end of file diff --git a/src/boundary_conditions/field_boundary_fill/FillFlather.H b/src/boundary_conditions/field_boundary_fill/FillFlather.H new file mode 100644 index 0000000000..cf35264772 --- /dev/null +++ b/src/boundary_conditions/field_boundary_fill/FillFlather.H @@ -0,0 +1,60 @@ +#ifndef FILLFLATHER_H +#define FILLFLATHER_H + +#include "src/core/FieldFillPatchOps.H" +#include "src/core/FieldBCOps.H" +#include "src/boundary_conditions/field_boundary_fill/Flather.H" + +namespace kynema_sgf { + +/** Fill patch operator that applies Flather boundary updates. */ +class FillFlather : public FieldFillPatchOps +{ +public: + FillFlather( + Field& field, + const amrex::AmrCore& mesh, + const SimTime& time, + Flather& flather); + + ~FillFlather() override; + + void fillpatch( + int lev, + amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& nghost, + FieldState fstate = FieldState::New) override; + + void fillpatch_sibling_fields( + int lev, + amrex::Real time, + amrex::Array& mfabs, + amrex::Array& ffabs, + amrex::Array& cfabs, + const amrex::IntVect& nghost, + const amrex::Vector& bcrec, + const amrex::Vector& /* unused */, + FieldState fstate = FieldState::New) override; + + void fillpatch_from_coarse( + int lev, + amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& nghost, + FieldState fstate = FieldState::New) override; + + void fillphysbc( + int lev, + amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& nghost, + FieldState fstate = FieldState::New) override; + +protected: + Flather& m_flather; +}; + +} // namespace kynema_sgf + +#endif /* FILLFLATHER_H */ diff --git a/src/boundary_conditions/field_boundary_fill/FillFlather.cpp b/src/boundary_conditions/field_boundary_fill/FillFlather.cpp new file mode 100644 index 0000000000..e109cad4b3 --- /dev/null +++ b/src/boundary_conditions/field_boundary_fill/FillFlather.cpp @@ -0,0 +1,75 @@ +#include + +#include "src/boundary_conditions/field_boundary_fill/FillFlather.H" + +namespace kynema_sgf { + +FillFlather::FillFlather( + Field& field, + const amrex::AmrCore& mesh, + const SimTime& time, + Flather& flather) + : FieldFillPatchOps( + field, mesh, time, FieldInterpolator::CellConsLinear) + , m_flather(flather) +{} + +FillFlather::~FillFlather() = default; + +void FillFlather::fillpatch( + const int lev, + const amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& /* nghost */, + const FieldState fstate) +{ + if (m_field.base_name() == "velocity") { + m_flather.update_flather_variables(lev, fstate); + m_flather.set_velocity(lev, time, m_field, mfab); + } +} + +void FillFlather::fillpatch_from_coarse( + const int lev, + const amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& /* nghost */, + const FieldState fstate) +{ + if (m_field.base_name() == "velocity") { + m_flather.update_flather_variables(lev, fstate); + m_flather.set_velocity(lev, time, m_field, mfab); + } +} + +void FillFlather::fillphysbc( + const int lev, + const amrex::Real time, + amrex::MultiFab& mfab, + const amrex::IntVect& /* nghost */, + const FieldState fstate) +{ + if (m_field.base_name() == "velocity") { + m_flather.update_flather_variables(lev, fstate); + m_flather.set_velocity(lev, time, m_field, mfab); + } +} + +void FillFlather::fillpatch_sibling_fields( + const int lev, + const amrex::Real time, + amrex::Array& mfabs, + amrex::Array& /* ffabs */, + amrex::Array& /* cfabs */, + const amrex::IntVect& /* nghost */, + const amrex::Vector& /* bcrec */, + const amrex::Vector& /* unused */, + const FieldState fstate) +{ + m_flather.update_flather_variables(lev, fstate, true); + for (int i = 0; std::cmp_less(i, mfabs.size()); ++i) { + m_flather.set_velocity(lev, time, m_field, *mfabs[i], 0, i); + } +} + +} // namespace kynema_sgf diff --git a/src/boundary_conditions/field_boundary_fill/FillMPL.cpp b/src/boundary_conditions/field_boundary_fill/FillMPL.cpp index a324acc19d..874924764c 100644 --- a/src/boundary_conditions/field_boundary_fill/FillMPL.cpp +++ b/src/boundary_conditions/field_boundary_fill/FillMPL.cpp @@ -20,12 +20,9 @@ void FillMPL::fillpatch( const int lev, const amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_abl_mpl.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "temperature") { @@ -37,12 +34,9 @@ void FillMPL::fillpatch_from_coarse( const int lev, const amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch_from_coarse( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_abl_mpl.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "temperature") { @@ -54,12 +48,9 @@ void FillMPL::fillphysbc( const int lev, const amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillphysbc( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_abl_mpl.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "temperature") { @@ -78,39 +69,36 @@ void FillMPL::fillpatch_sibling_fields( const amrex::Vector& /* unused */, const FieldState fstate) { - if (m_field.base_name() == "velocity") { - // For an ABL MPL, we first just foextrap the mac velocities - amrex::Vector lbcrec(m_field.num_comp()); - const auto& ibctype = m_field.bc_type(); - for (amrex::OrientationIter oit; oit != nullptr; ++oit) { - auto ori = oit(); - const auto side = ori.faceDir(); - const auto bct = ibctype[ori]; - const int dir = ori.coordDir(); - for (int i = 0; i < m_field.num_comp(); ++i) { - if ((bct == BC::mass_inflow) || - (bct == BC::mass_inflow_outflow)) { - if (side == amrex::Orientation::low) { - lbcrec[i].setLo(dir, amrex::BCType::foextrap); - } else { - lbcrec[i].setHi(dir, amrex::BCType::foextrap); - } + // For an ABL MPL, we first just foextrap the mac velocities + amrex::Vector lbcrec(m_field.num_comp()); + const auto& ibctype = m_field.bc_type(); + for (amrex::OrientationIter oit; oit != nullptr; ++oit) { + auto ori = oit(); + const auto side = ori.faceDir(); + const auto bct = ibctype[ori]; + const int dir = ori.coordDir(); + for (int i = 0; i < m_field.num_comp(); ++i) { + if ((bct == BC::mass_inflow) || (bct == BC::mass_inflow_outflow)) { + if (side == amrex::Orientation::low) { + lbcrec[i].setLo(dir, amrex::BCType::foextrap); } else { - if (side == amrex::Orientation::low) { - lbcrec[i].setLo(dir, bcrec[i].lo(dir)); - } else { - lbcrec[i].setHi(dir, bcrec[i].hi(dir)); - } + lbcrec[i].setHi(dir, amrex::BCType::foextrap); + } + } else { + if (side == amrex::Orientation::low) { + lbcrec[i].setLo(dir, bcrec[i].lo(dir)); + } else { + lbcrec[i].setHi(dir, bcrec[i].hi(dir)); } } } + } - FieldFillPatchOps::fillpatch_sibling_fields( - lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, lbcrec, fstate); + FieldFillPatchOps::fillpatch_sibling_fields( + lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, lbcrec, fstate); - for (int i = 0; std::cmp_less(i, mfabs.size()); i++) { - m_abl_mpl.set_velocity(lev, time, m_field, *mfabs[i], 0, i); - } + for (int i = 0; std::cmp_less(i, mfabs.size()); i++) { + m_abl_mpl.set_velocity(lev, time, m_field, *mfabs[i], 0, i); } } diff --git a/src/boundary_conditions/field_boundary_fill/Flather.H b/src/boundary_conditions/field_boundary_fill/Flather.H new file mode 100644 index 0000000000..00297c265b --- /dev/null +++ b/src/boundary_conditions/field_boundary_fill/Flather.H @@ -0,0 +1,107 @@ +#ifndef FLATHER_H +#define FLATHER_H + +#include "src/core/Field.H" +#include "src/core/IntField.H" +#include "src/CFDSim.H" +#include "src/boundary_conditions/field_boundary_fill/FieldBoundary.H" +#include "src/utilities/MultiLevelVector.H" +#include "src/utilities/trig_ops.H" +#include "AMReX_REAL.H" +#include + +using namespace amrex::literals; + +namespace kynema_sgf { + +/** Flather-style open boundary for velocity + * + * This boundary reads current boundary ghost-cell values and adjacent interior + * values, computes a blended target, and writes updated boundary values. + */ +class Flather : public FieldBoundary::Register +{ +public: + static std::string identifier() { return "Flather"; } + + explicit Flather(CFDSim& sim); + + void post_init_actions() override; + + void pre_advance_work() override; + + void post_advance_work() override {} + + void set_velocity( + int lev, + amrex::Real time, + const Field& fld, + amrex::MultiFab& mfab, + int dcomp = 0, + int orig_comp = 0) const override; + + void update_flather_variables(const int lev, const FieldState fstate) + { + compute_boundary_z_averages(lev, fstate, false); + } + + void update_flather_variables( + const int lev, const FieldState fstate, const bool use_mac_fields) + { + compute_boundary_z_averages(lev, fstate, use_mac_fields); + } + + void accumulate_boundary( + int current_level, + int idir, + int phase_switch, + bool is_low, + MultiLevelVector& out_uvec, + MultiLevelVector& out_hvec, + bool sample_boundary, + FieldState fstate, + bool use_mac_fields = false) const; + +private: + void compute_internal_z_averages(); + void compute_boundary_z_averages( + int lev, FieldState fstate, bool use_mac_fields = false); + + const CFDSim& m_sim; + const kynema_sgf::SimTime& m_time; + const FieldRepo& m_repo; + const amrex::AmrCore& m_mesh; + Field& m_velocity; + const Field& m_u_mac; + const Field& m_v_mac; + const Field& m_vof; + const IntField* m_terrain_blank{nullptr}; + amrex::Vector m_gravity{0.0_rt, 0.0_rt, -9.81_rt}; + amrex::Real m_rho1{1000.0_rt}; + amrex::Real m_rho2{1.0_rt}; + + MultiLevelVector m_xlo_uliq{FieldLoc::CELL}; + MultiLevelVector m_xhi_uliq{FieldLoc::CELL}; + MultiLevelVector m_ylo_uliq{FieldLoc::CELL}; + MultiLevelVector m_yhi_uliq{FieldLoc::CELL}; + MultiLevelVector m_xlo_umix{FieldLoc::CELL}; + MultiLevelVector m_xhi_umix{FieldLoc::CELL}; + MultiLevelVector m_ylo_umix{FieldLoc::CELL}; + MultiLevelVector m_yhi_umix{FieldLoc::CELL}; + MultiLevelVector m_xlo_bnd_uvof{FieldLoc::CELL}; + MultiLevelVector m_xhi_bnd_uvof{FieldLoc::CELL}; + MultiLevelVector m_ylo_bnd_uvof{FieldLoc::CELL}; + MultiLevelVector m_yhi_bnd_uvof{FieldLoc::CELL}; + MultiLevelVector m_xlo_h_avg{FieldLoc::CELL}; + MultiLevelVector m_xhi_h_avg{FieldLoc::CELL}; + MultiLevelVector m_ylo_h_avg{FieldLoc::CELL}; + MultiLevelVector m_yhi_h_avg{FieldLoc::CELL}; + MultiLevelVector m_xlo_bnd_h_avg{FieldLoc::CELL}; + MultiLevelVector m_xhi_bnd_h_avg{FieldLoc::CELL}; + MultiLevelVector m_ylo_bnd_h_avg{FieldLoc::CELL}; + MultiLevelVector m_yhi_bnd_h_avg{FieldLoc::CELL}; +}; + +} // namespace kynema_sgf + +#endif /* FLATHER_H */ diff --git a/src/boundary_conditions/field_boundary_fill/Flather.cpp b/src/boundary_conditions/field_boundary_fill/Flather.cpp new file mode 100644 index 0000000000..0bb59f3976 --- /dev/null +++ b/src/boundary_conditions/field_boundary_fill/Flather.cpp @@ -0,0 +1,588 @@ +#include "src/CFDSim.H" +#include "src/boundary_conditions/field_boundary_fill/Flather.H" +#include "src/boundary_conditions/field_boundary_fill/FillFlather.H" +#include "src/utilities/index_operations.H" +#include "src/utilities/constants.H" +#include "src/physics/multiphase/MultiPhase.H" +#include "AMReX_MultiFabUtil.H" +#include "AMReX_GpuAtomic.H" +#include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" +#include +#include + +using namespace amrex::literals; + +namespace kynema_sgf { + +Flather::Flather(CFDSim& sim) + : m_sim(sim) + , m_time(m_sim.time()) + , m_repo(m_sim.repo()) + , m_mesh(m_sim.mesh()) + , m_velocity(m_sim.repo().get_field("velocity")) + , m_u_mac(m_sim.repo().get_field("u_mac")) + , m_v_mac(m_sim.repo().get_field("v_mac")) + , m_vof(m_sim.repo().get_field("vof")) +{ + // amrex::ParmParse pp(identifier()); + + if (!m_repo.field_exists("vof")) { + amrex::Abort("Flather BC requires the vof field"); + } + if (m_repo.int_field_exists("terrain_blank")) { + m_terrain_blank = &m_repo.get_int_field("terrain_blank"); + } + + if (m_sim.physics_manager().contains("MultiPhase")) { + m_rho1 = m_sim.physics_manager().get().rho1(); + m_rho2 = m_sim.physics_manager().get().rho2(); + } + + amrex::ParmParse pp("incflo"); + pp.queryarr("gravity", m_gravity); + + // Vector at the x boundary extends in y direction + // Vector at the y boundary extends in x direction + m_xlo_uliq.resize(1, m_mesh.Geom()); + m_xhi_uliq.resize(1, m_mesh.Geom()); + m_ylo_uliq.resize(0, m_mesh.Geom()); + m_yhi_uliq.resize(0, m_mesh.Geom()); + m_xlo_umix.resize(1, m_mesh.Geom()); + m_xhi_umix.resize(1, m_mesh.Geom()); + m_ylo_umix.resize(0, m_mesh.Geom()); + m_yhi_umix.resize(0, m_mesh.Geom()); + m_xlo_bnd_uvof.resize(1, m_mesh.Geom()); + m_xhi_bnd_uvof.resize(1, m_mesh.Geom()); + m_ylo_bnd_uvof.resize(0, m_mesh.Geom()); + m_yhi_bnd_uvof.resize(0, m_mesh.Geom()); + + m_xlo_h_avg.resize(1, m_mesh.Geom()); + m_xhi_h_avg.resize(1, m_mesh.Geom()); + m_ylo_h_avg.resize(0, m_mesh.Geom()); + m_yhi_h_avg.resize(0, m_mesh.Geom()); + m_xlo_bnd_h_avg.resize(1, m_mesh.Geom()); + m_xhi_bnd_h_avg.resize(1, m_mesh.Geom()); + m_ylo_bnd_h_avg.resize(0, m_mesh.Geom()); + m_yhi_bnd_h_avg.resize(0, m_mesh.Geom()); +} + +void Flather::post_init_actions() +{ + m_velocity.add_fill_patch_op(m_mesh, m_time, *this); + compute_internal_z_averages(); +} + +void Flather::pre_advance_work() { compute_internal_z_averages(); } + +void Flather::accumulate_boundary( + int current_level, + int idir, + int phase_switch, + bool is_low, + MultiLevelVector& out_uvec, + MultiLevelVector& out_hvec, + bool sample_boundary, + FieldState fstate, + bool use_mac_fields) const +{ + AMREX_ALWAYS_ASSERT(idir == 0 || idir == 1); + + auto& int_h = out_uvec.host_data(current_level); + auto& dist_h = out_hvec.host_data(current_level); + const int nline = static_cast(int_h.size()); + + amrex::Gpu::DeviceVector uvof_sum_d(nline, 0.0_rt); + amrex::Gpu::DeviceVector vof_sum_d(nline, 0.0_rt); + + const amrex::Real tiny = constants::TIGHT_TOL; + + // Loop through current level and all below + for (int lev = current_level; lev >= 0; --lev) { + + const auto rr = m_mesh.refRatio(lev)[1 - idir]; + + amrex::iMultiFab level_mask; + if (lev < current_level) { + level_mask = makeFineMask( + m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), + m_mesh.boxArray(lev + 1), m_mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + m_mesh.boxArray(lev), m_mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& geom = m_mesh.Geom(lev); + const auto dz = geom.CellSizeArray()[2]; + const auto& dom = geom.Domain(); + const int bidx = is_low ? dom.smallEnd(idir) : dom.bigEnd(idir); + const int shift_to_boundary = sample_boundary ? (is_low ? -1 : 1) : 0; + const auto& src_vel = + use_mac_fields ? ((idir == 0) ? m_u_mac : m_v_mac) : m_velocity; + if (shift_to_boundary != 0) { + AMREX_ALWAYS_ASSERT(src_vel.num_grow()[idir] > 0); + AMREX_ALWAYS_ASSERT(m_vof.num_grow()[idir] > 0); + } + + const auto& vel_mf = + src_vel.state(use_mac_fields ? FieldState::New : fstate)(lev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(vel_mf, amrex::TilingIfNotGPU()); mfi.isValid(); + ++mfi) { + amrex::Box bx = mfi.tilebox() & dom; + if (!bx.ok()) { + continue; + } + if (bidx < bx.smallEnd(idir) || bidx > bx.bigEnd(idir)) { + continue; + } + + // Limit box to along boundary + bx.setSmall(idir, bidx); + bx.setBig(idir, bidx); + + const auto vel_arr = vel_mf.const_array(mfi); + // Due to the order that fillphysbc is called in prepare_boundaries + // (velocity, then scalars), the vof data is not available at + // alternate (nph) states when this is called. + const auto vof_arr = m_vof(lev).const_array(mfi); + const auto mask_arr = level_mask.const_array(mfi); + const bool use_terrain = (m_terrain_blank != nullptr); + const auto terrain_blank_arr = + use_terrain ? (*m_terrain_blank)(lev).const_array(mfi) + : amrex::Array4(); + + amrex::Real* uvof_sum = uvof_sum_d.data(); + amrex::Real* vof_sum = vof_sum_d.data(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const int ii = (idir == 0) ? (i + shift_to_boundary) : i; + const int jj = (idir == 1) ? (j + shift_to_boundary) : j; + const int ii_v = + ii + + static_cast(use_mac_fields && idir == 0 && !is_low); + const int jj_v = + jj + + static_cast(use_mac_fields && idir == 1 && !is_low); + + if (use_terrain && terrain_blank_arr(i, j, k) != 0) { + return; + } + // This index is tangent to the boundary + const int idx_lev = (idir == 0) ? j : i; + // Convert to current level indices + const int idx_min = + idx_lev + idx_lev * (rr - 1) * (current_level - lev); + const int idx_max = + idx_min + + amrex::max(0, rr * (current_level - lev) - 1); + + for (int idx = idx_min; idx <= idx_max; ++idx) { + const auto liquid_height = + vof_arr(ii, jj, k) * dz * mask_arr(i, j, k); + amrex::Gpu::Atomic::Add(&vof_sum[idx], liquid_height); + auto vel_height = liquid_height; + if (phase_switch == 0 && + vof_arr(ii, jj, k) < 1.0_rt - tiny) { + // Needs to be fully liquid to be counted + vel_height = 0.0_rt; + } + if (phase_switch == 1 && + vof_arr(ii, jj, k) >= 1.0_rt - tiny) { + // Needs to be mixture to be counted + vel_height = 0.0_rt; + } + // Fully gas cells never count because of vof multiplier + amrex::Real local_vel = + vel_arr(ii_v, jj_v, k, use_mac_fields ? 0 : idir); + // If summing internal velocities, only allow outflow + if (!sample_boundary) { + local_vel = is_low ? amrex::min(local_vel, 0.0_rt) + : amrex::max(local_vel, 0.0_rt); + } + amrex::Gpu::Atomic::Add( + &uvof_sum[idx], local_vel * vel_height); + } + }); + } + } + + amrex::Gpu::copy( + amrex::Gpu::deviceToHost, uvof_sum_d.begin(), uvof_sum_d.end(), + int_h.begin()); + amrex::Gpu::copy( + amrex::Gpu::deviceToHost, vof_sum_d.begin(), vof_sum_d.end(), + dist_h.begin()); + + amrex::ParallelDescriptor::ReduceRealSum(int_h.data(), nline); + amrex::ParallelDescriptor::ReduceRealSum(dist_h.data(), nline); +} + +void Flather::compute_internal_z_averages() +{ + BL_PROFILE("kynema-sgf::Flather::compute_internal_z_averages"); + + const int nlevels = m_repo.num_active_levels(); + const int nlevels_geom = static_cast(m_mesh.Geom().size()); + if (m_xlo_uliq.size() != nlevels_geom) { + m_xlo_uliq.resize(1, m_mesh.Geom()); + m_xhi_uliq.resize(1, m_mesh.Geom()); + m_ylo_uliq.resize(0, m_mesh.Geom()); + m_yhi_uliq.resize(0, m_mesh.Geom()); + m_xlo_umix.resize(1, m_mesh.Geom()); + m_xhi_umix.resize(1, m_mesh.Geom()); + m_ylo_umix.resize(0, m_mesh.Geom()); + m_yhi_umix.resize(0, m_mesh.Geom()); + m_xlo_bnd_uvof.resize(1, m_mesh.Geom()); + m_xhi_bnd_uvof.resize(1, m_mesh.Geom()); + m_ylo_bnd_uvof.resize(0, m_mesh.Geom()); + m_yhi_bnd_uvof.resize(0, m_mesh.Geom()); + + m_xlo_h_avg.resize(1, m_mesh.Geom()); + m_xhi_h_avg.resize(1, m_mesh.Geom()); + m_ylo_h_avg.resize(0, m_mesh.Geom()); + m_yhi_h_avg.resize(0, m_mesh.Geom()); + m_xlo_bnd_h_avg.resize(1, m_mesh.Geom()); + m_xhi_bnd_h_avg.resize(1, m_mesh.Geom()); + m_ylo_bnd_h_avg.resize(0, m_mesh.Geom()); + m_yhi_bnd_h_avg.resize(0, m_mesh.Geom()); + } + + for (int lev = 0; lev < nlevels; ++lev) { + this->accumulate_boundary( + lev, 0, 0, true, m_xlo_uliq, m_xlo_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 0, 0, false, m_xhi_uliq, m_xhi_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 1, 0, true, m_ylo_uliq, m_ylo_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 1, 0, false, m_yhi_uliq, m_yhi_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 0, 1, true, m_xlo_umix, m_xlo_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 0, 1, false, m_xhi_umix, m_xhi_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 1, 1, true, m_ylo_umix, m_ylo_h_avg, false, FieldState::New); + this->accumulate_boundary( + lev, 1, 1, false, m_yhi_umix, m_yhi_h_avg, false, FieldState::New); + } + + m_xlo_uliq.copy_host_to_device(); + m_xhi_uliq.copy_host_to_device(); + m_ylo_uliq.copy_host_to_device(); + m_yhi_uliq.copy_host_to_device(); + + m_xlo_umix.copy_host_to_device(); + m_xhi_umix.copy_host_to_device(); + m_ylo_umix.copy_host_to_device(); + m_yhi_umix.copy_host_to_device(); + + m_xlo_h_avg.copy_host_to_device(); + m_xhi_h_avg.copy_host_to_device(); + m_ylo_h_avg.copy_host_to_device(); + m_yhi_h_avg.copy_host_to_device(); +} + +void Flather::compute_boundary_z_averages( + int lev, FieldState fstate, bool use_mac_fields) +{ + BL_PROFILE("kynema-sgf::Flather::compute_boundary_z_averages"); + + // accumulating boundaries needs to happen prior to applying fillpatch op + + this->accumulate_boundary( + lev, 0, -1, true, m_xlo_bnd_uvof, m_xlo_bnd_h_avg, true, fstate, + use_mac_fields); + this->accumulate_boundary( + lev, 0, -1, false, m_xhi_bnd_uvof, m_xhi_bnd_h_avg, true, fstate, + use_mac_fields); + this->accumulate_boundary( + lev, 1, -1, true, m_ylo_bnd_uvof, m_ylo_bnd_h_avg, true, fstate, + use_mac_fields); + this->accumulate_boundary( + lev, 1, -1, false, m_yhi_bnd_uvof, m_yhi_bnd_h_avg, true, fstate, + use_mac_fields); + + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_xlo_bnd_uvof.host_data(lev).begin(), + m_xlo_bnd_uvof.host_data(lev).end(), + m_xlo_bnd_uvof.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_xhi_bnd_uvof.host_data(lev).begin(), + m_xhi_bnd_uvof.host_data(lev).end(), + m_xhi_bnd_uvof.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_ylo_bnd_uvof.host_data(lev).begin(), + m_ylo_bnd_uvof.host_data(lev).end(), + m_ylo_bnd_uvof.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_yhi_bnd_uvof.host_data(lev).begin(), + m_yhi_bnd_uvof.host_data(lev).end(), + m_yhi_bnd_uvof.device_data(lev).begin()); + + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_xlo_bnd_h_avg.host_data(lev).begin(), + m_xlo_bnd_h_avg.host_data(lev).end(), + m_xlo_bnd_h_avg.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_xhi_bnd_h_avg.host_data(lev).begin(), + m_xhi_bnd_h_avg.host_data(lev).end(), + m_xhi_bnd_h_avg.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_ylo_bnd_h_avg.host_data(lev).begin(), + m_ylo_bnd_h_avg.host_data(lev).end(), + m_ylo_bnd_h_avg.device_data(lev).begin()); + amrex::Gpu::copyAsync( + amrex::Gpu::hostToDevice, m_yhi_bnd_h_avg.host_data(lev).begin(), + m_yhi_bnd_h_avg.host_data(lev).end(), + m_yhi_bnd_h_avg.device_data(lev).begin()); +} + +void Flather::set_velocity( + const int lev, + const amrex::Real time, + const Field& fld, + amrex::MultiFab& mfab, + const int /* dcomp */, + const int orig_comp) const +{ + BL_PROFILE("kynema-sgf::Flather::set_velocity"); + + const auto& geom = m_mesh.Geom(lev); + const auto& bctype = fld.bc_type(); + const int nghost = 1; + const int numcomp = mfab.nComp(); + const auto& domain = geom.growPeriodicDomain(nghost); + + const auto fstate = time > 0.0_rt ? FieldState::Old : FieldState::New; + const amrex::Real tiny = constants::TIGHT_TOL; + const amrex::Real v_threshold = 1.0e-6_rt; + const auto grav_z = -m_gravity[2]; + const auto rho1 = m_rho1; + const auto rho2 = m_rho2; + + for (amrex::OrientationIter oit; oit != nullptr; ++oit) { + const auto ori = oit(); + if ((bctype[ori] != BC::mass_inflow) && + (bctype[ori] != BC::mass_inflow_outflow)) { + continue; + } + + const int idir = ori.coordDir(); + // Check if orientation aligns with supplied field and choose component + bool skip_fill = false; + int fcomp = 0; + if (numcomp == 1) { + // MAC velocity: only normal field is valid + skip_fill = (orig_comp != idir); + // Only a single component available + } else { + // Cell-centered velocity: field always valid + fcomp = idir; + // Select valid component + } + if (skip_fill) { + continue; + } + + const auto& dbx = ori.isLow() ? amrex::adjCellLo(domain, idir, nghost) + : amrex::adjCellHi(domain, idir, nghost); + const auto shift_to_interior = + amrex::IntVect::TheDimensionVector(idir) * (ori.isLow() ? 1 : -1); + const amrex::Real* xlo_uliq = m_xlo_uliq.device_data(lev).data(); + const amrex::Real* xhi_uliq = m_xhi_uliq.device_data(lev).data(); + const amrex::Real* ylo_uliq = m_ylo_uliq.device_data(lev).data(); + const amrex::Real* yhi_uliq = m_yhi_uliq.device_data(lev).data(); + const amrex::Real* xlo_umix = m_xlo_umix.device_data(lev).data(); + const amrex::Real* xhi_umix = m_xhi_umix.device_data(lev).data(); + const amrex::Real* ylo_umix = m_ylo_umix.device_data(lev).data(); + const amrex::Real* yhi_umix = m_yhi_umix.device_data(lev).data(); + const amrex::Real* xlo_havg = m_xlo_h_avg.device_data(lev).data(); + const amrex::Real* xhi_havg = m_xhi_h_avg.device_data(lev).data(); + const amrex::Real* ylo_havg = m_ylo_h_avg.device_data(lev).data(); + const amrex::Real* yhi_havg = m_yhi_h_avg.device_data(lev).data(); + const amrex::Real* xlo_bnd_u = m_xlo_bnd_uvof.device_data(lev).data(); + const amrex::Real* xhi_bnd_u = m_xhi_bnd_uvof.device_data(lev).data(); + const amrex::Real* ylo_bnd_u = m_ylo_bnd_uvof.device_data(lev).data(); + const amrex::Real* yhi_bnd_u = m_yhi_bnd_uvof.device_data(lev).data(); + const amrex::Real* xlo_bnd_havg = + m_xlo_bnd_h_avg.device_data(lev).data(); + const amrex::Real* xhi_bnd_havg = + m_xhi_bnd_h_avg.device_data(lev).data(); + const amrex::Real* ylo_bnd_havg = + m_ylo_bnd_h_avg.device_data(lev).data(); + const amrex::Real* yhi_bnd_havg = + m_yhi_bnd_h_avg.device_data(lev).data(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (false) +#endif + for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) { + auto gbx = amrex::grow(mfi.validbox(), nghost); + auto shift_to_cc = amrex::IntVect(0); + const auto& bx = utils::face_aware_boundary_box_intersection( + shift_to_cc, gbx, dbx, ori); + if (!bx.ok()) { + continue; + } + + const auto& arr = mfab[mfi].array(); + const auto& ref_arr = fld.state(fstate)(lev)[mfi].const_array(); + + const auto vof_arr = m_vof.state(fstate)(lev).const_array(mfi); + const bool use_terrain = (m_terrain_blank != nullptr); + const auto terrain_blank_arr = + use_terrain ? (*m_terrain_blank)(lev).const_array(mfi) + : amrex::Array4(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const amrex::IntVect iv{i, j, k}; + const amrex::IntVect iv_cc = iv + shift_to_cc; + const amrex::IntVect iv_adj = iv + shift_to_interior; + const amrex::IntVect iv_adj_cc = iv_adj + shift_to_cc; + + amrex::Real boundary_val = arr(iv, fcomp); + amrex::Real interior_val = arr(iv_adj, fcomp); // Unused? + amrex::Real interior_liq = arr(iv_adj, fcomp); + amrex::Real interior_mix = arr(iv_adj, fcomp); + amrex::Real boundary_h = 0.0_rt; + amrex::Real interior_h = 0.0_rt; + // Vectors at x boundaries extend in y direction + // Vectors at y boundaries extend in x direction + if (idir == 0) { + interior_liq = + ori.isLow() ? xlo_uliq[iv_adj[1]] : xhi_uliq[iv_adj[1]]; + interior_mix = + ori.isLow() ? xlo_umix[iv_adj[1]] : xhi_umix[iv_adj[1]]; + interior_val = interior_liq + interior_mix; + interior_h = + ori.isLow() ? xlo_havg[iv_adj[1]] : xhi_havg[iv_adj[1]]; + boundary_val = + ori.isLow() ? xlo_bnd_u[iv[1]] : xhi_bnd_u[iv[1]]; + boundary_h = + ori.isLow() ? xlo_bnd_havg[iv[1]] : xhi_bnd_havg[iv[1]]; + } else { + interior_liq = + ori.isLow() ? ylo_uliq[iv_adj[0]] : yhi_uliq[iv_adj[0]]; + interior_mix = + ori.isLow() ? ylo_umix[iv_adj[0]] : yhi_umix[iv_adj[0]]; + interior_val = interior_liq + interior_mix; + interior_h = + ori.isLow() ? ylo_havg[iv_adj[0]] : yhi_havg[iv_adj[0]]; + boundary_val = + ori.isLow() ? ylo_bnd_u[iv[0]] : yhi_bnd_u[iv[0]]; + boundary_h = + ori.isLow() ? ylo_bnd_havg[iv[0]] : yhi_bnd_havg[iv[0]]; + } + + // Set velocity to zero and skip for terrain + if (use_terrain && terrain_blank_arr(iv_adj_cc) != 0) { + arr(iv, fcomp) = 0.0_rt; + return; + } + // Check interior or boundary vof for liquid + const amrex::Real interior_vof = vof_arr(iv_adj_cc); + const amrex::Real boundary_vof = vof_arr(iv_cc); + // Do nothing here if not liquid + if ((interior_vof < tiny) && (boundary_vof < tiny)) { + return; + } + + // Wave speed + const amrex::Real c = std::sqrt(grav_z * interior_h); + + const auto Flather_val = + boundary_val + (ori.isLow() ? -1.0_rt : 1.0_rt) * c * + (interior_h - boundary_h); + + // Use external (prescribed) velocity if inflow + // Assesses inflow by the whole column, not the local value + // Assumes values are up-to-date from another fillpatch op + const bool prescribed_inflow = + ori.isLow() ? boundary_val > 0.0_rt : boundary_val < 0.0_rt; + + // Clip the velocity in the interior to prevent inflow at + // outflow, this is consistent with the line integral + // calculations + const auto local_internal_vel = + ori.isLow() ? amrex::min(ref_arr(iv_adj_cc, idir), 0.0_rt) + : amrex::max(ref_arr(iv_adj_cc, idir), 0.0_rt); + + // Normal outflow case: + // *) For a given column, apply scale to the interior velocity, + // but only + // to fully liquid cells; leave mixed cells unchanged. + + // Edge cases: + // 1) If the boundary contains no liquid, the flather can be + // very large. + // Use a regular outflow (Neumann) condition. + // 2) If the interior column contains no fully liquid cells, + // there is + // nothing to scale. Instead of scaling the internal profile, + // override the interior velocity with scaled external + // profile. + // 3) If the boundary velocity is 0, then the scaling based on + // external quantities is undefined. Keep the scale_interior + // = 1 + // 4) During initialization, the scaling can be very aggressive. + // Plus, when the internal and external profiles are very + // different, the scaling can lead to rapid acceleration. + // Switch to the external profile when the changes are rapid + + bool override_interior = false; + auto scale_interior = 1.0; + if (std::abs(interior_liq) > v_threshold * interior_h) { + scale_interior = + (Flather_val - interior_mix) / interior_liq; + } else if (std::abs(boundary_val) > v_threshold * boundary_h) { + override_interior = true; + } + + if (scale_interior > 2.0_rt || scale_interior < 0.25_rt) { + override_interior = true; + } + + /* + // Limit the scaling factor to prevent overly aggressive changes + scale_interior = + amrex::max(0.5_rt, amrex::min(scale_interior, 1.2_rt)); + */ + + auto local_vel = 0.0_rt; + auto scaled_vel = 0.0_rt; + if (prescribed_inflow || override_interior) { + local_vel = arr(iv, fcomp); + scaled_vel = local_vel * (Flather_val / boundary_val); + } else { + local_vel = local_internal_vel; + scaled_vel = local_vel * scale_interior; + } + + // Only use if advecting liquid (or zero velocity) + // This is important because the averages used to calculate + // the Flather velocity are only performed on cells containing + // liquid; therefore, the scaling of the velocity is only valid + // on those cells. + const bool outflow = + ori.isLow() ? scaled_vel <= 0.0_rt : scaled_vel >= 0.0_rt; + const bool inflow_any_liq = !outflow && boundary_vof > tiny; + const bool outflow_only_liq = + outflow && interior_vof < 1.0_rt - tiny; + if (boundary_h > tiny && (outflow_only_liq || inflow_any_liq || + (outflow && override_interior))) { + arr(iv, fcomp) = scaled_vel; + } else if (outflow) { + arr(iv, fcomp) = local_vel; + } + }); + } + } +} + +} // namespace kynema_sgf diff --git a/src/boundary_conditions/field_boundary_fill/ModulatedPowerLaw.cpp b/src/boundary_conditions/field_boundary_fill/ModulatedPowerLaw.cpp index 58b3c542ca..f52d6d3924 100644 --- a/src/boundary_conditions/field_boundary_fill/ModulatedPowerLaw.cpp +++ b/src/boundary_conditions/field_boundary_fill/ModulatedPowerLaw.cpp @@ -71,8 +71,8 @@ ModulatedPowerLaw::ModulatedPowerLaw(CFDSim& sim) void ModulatedPowerLaw::post_init_actions() { - m_velocity.register_fill_patch_op(m_mesh, m_time, *this); - m_temperature.register_fill_patch_op(m_mesh, m_time, *this); + m_velocity.add_fill_patch_op(m_mesh, m_time, *this); + m_temperature.add_fill_patch_op(m_mesh, m_time, *this); } void ModulatedPowerLaw::pre_advance_work() diff --git a/src/boundary_conditions/field_boundary_fill/OceanWavesBoundary.cpp b/src/boundary_conditions/field_boundary_fill/OceanWavesBoundary.cpp index da224f969c..636fac9e7b 100644 --- a/src/boundary_conditions/field_boundary_fill/OceanWavesBoundary.cpp +++ b/src/boundary_conditions/field_boundary_fill/OceanWavesBoundary.cpp @@ -30,14 +30,13 @@ void OceanWavesBoundary::post_init_actions() { BL_PROFILE("kynema-sgf::OceanWavesBoundary::post_init_actions"); m_repo.get_field("velocity") - .register_fill_patch_op(m_mesh, m_time, *this); + .add_fill_patch_op(m_mesh, m_time, *this); m_vof_exists = m_repo.field_exists("vof"); if (m_vof_exists) { - m_repo.get_field("vof").register_fill_patch_op( + m_repo.get_field("vof").add_fill_patch_op( + m_mesh, m_time, *this); + m_repo.get_field("density").add_fill_patch_op( m_mesh, m_time, *this); - m_repo.get_field("density") - .register_fill_patch_op( - m_mesh, m_time, *this); } m_terrain_exists = m_repo.int_field_exists("terrain_blank"); diff --git a/src/boundary_conditions/field_boundary_fill/OceanWavesFillInflow.cpp b/src/boundary_conditions/field_boundary_fill/OceanWavesFillInflow.cpp index d2362ecf22..b290804c4c 100644 --- a/src/boundary_conditions/field_boundary_fill/OceanWavesFillInflow.cpp +++ b/src/boundary_conditions/field_boundary_fill/OceanWavesFillInflow.cpp @@ -21,12 +21,9 @@ void OceanWavesFillInflow::fillpatch( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_ow_bndry.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "vof") { @@ -40,12 +37,9 @@ void OceanWavesFillInflow::fillpatch_from_coarse( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch_from_coarse( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_ow_bndry.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "vof") { @@ -59,12 +53,9 @@ void OceanWavesFillInflow::fillphysbc( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillphysbc( - lev, time, mfab, nghost, fstate); - if (m_field.base_name() == "velocity") { m_ow_bndry.set_velocity(lev, time, m_field, mfab); } else if (m_field.base_name() == "vof") { diff --git a/src/boundary_conditions/field_boundary_fill/PlaneFillInflow.cpp b/src/boundary_conditions/field_boundary_fill/PlaneFillInflow.cpp index 36b33cd4ae..32ac242666 100644 --- a/src/boundary_conditions/field_boundary_fill/PlaneFillInflow.cpp +++ b/src/boundary_conditions/field_boundary_fill/PlaneFillInflow.cpp @@ -20,12 +20,9 @@ void PlaneFillInflow::fillpatch( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch( - lev, time, mfab, nghost, fstate); - m_bndry_plane.populate_data(lev, time, m_field, mfab); } @@ -33,12 +30,9 @@ void PlaneFillInflow::fillpatch_from_coarse( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillpatch_from_coarse( - lev, time, mfab, nghost, fstate); - m_bndry_plane.populate_data(lev, time, m_field, mfab); } @@ -46,12 +40,9 @@ void PlaneFillInflow::fillphysbc( int lev, amrex::Real time, amrex::MultiFab& mfab, - const amrex::IntVect& nghost, - const FieldState fstate) + const amrex::IntVect& /* nghost */, + const FieldState /* fstate */) { - FieldFillPatchOps::fillphysbc( - lev, time, mfab, nghost, fstate); - m_bndry_plane.populate_data(lev, time, m_field, mfab); } diff --git a/src/core/Field.H b/src/core/Field.H index dc6ea6230f..3465ad4ee4 100644 --- a/src/core/Field.H +++ b/src/core/Field.H @@ -71,7 +71,7 @@ struct FieldInfo amrex::Vector m_states; //! Function that handles filling patch and physics BC data for this field - std::unique_ptr m_fillpatch_op; + amrex::Vector> m_fillpatch_ops; //! Custom boundary condition actions for this field amrex::Vector> m_bc_func; @@ -222,7 +222,7 @@ public: //! Return a flag indicating whether a fillpatch Op has been registered [[nodiscard]] bool has_fillpatch_op() const { - return static_cast(m_info->m_fillpatch_op); + return !m_info->m_fillpatch_ops.empty(); } /** Setup default BC conditions for fillpatch operations @@ -312,7 +312,25 @@ public: template void register_fill_patch_op(Args&&... args) { - m_info->m_fillpatch_op.reset(new T(*this, std::forward(args)...)); + if (m_info->m_fillpatch_ops.empty()) { + m_info->m_fillpatch_ops.emplace_back( + new T(*this, std::forward(args)...)); + } else if (m_info->m_fillpatch_ops.size() == 1) { + m_info->m_fillpatch_ops[0].reset( + new T(*this, std::forward(args)...)); + } else { + amrex::Abort( + "register_fill_patch_op: More than one fillpatch op already " + "registered. register_fill_patch_op cannot be called after " + "add_fill_patch_op"); + } + } + + template + void add_fill_patch_op(Args&&... args) + { + m_info->m_fillpatch_ops.emplace_back( + new T(*this, std::forward(args)...)); } /** Register a custom boundary conditions class diff --git a/src/core/Field.cpp b/src/core/Field.cpp index 378796c9b5..bc1b5edbd3 100644 --- a/src/core/Field.cpp +++ b/src/core/Field.cpp @@ -173,11 +173,12 @@ void Field::fillpatch( const amrex::IntVect& nghost) { BL_PROFILE("kynema-sgf::Field::fillpatch 2"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); - fop.fillpatch(lev, time, mfab, nghost, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillpatch(lev, time, mfab, nghost, field_state()); + } } void Field::fillpatch_from_coarse( @@ -187,23 +188,24 @@ void Field::fillpatch_from_coarse( const amrex::IntVect& nghost) { BL_PROFILE("kynema-sgf::Field::fillpatch_from_coarse"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); - - fop.fillpatch_from_coarse(lev, time, mfab, nghost, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillpatch_from_coarse(lev, time, mfab, nghost, field_state()); + } } void Field::fillpatch(const amrex::Real time, const amrex::IntVect ng) { BL_PROFILE("kynema-sgf::Field::fillpatch"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); const int nlevels = m_repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { - fop.fillpatch( - lev, time, m_repo.get_multifab(m_id, lev), ng, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillpatch( + lev, time, m_repo.get_multifab(m_id, lev), ng, field_state()); + } } } @@ -215,10 +217,9 @@ void Field::fillpatch_sibling_fields( amrex::Array& fields) const { BL_PROFILE("kynema-sgf::Field::fillpatch array"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); BL_ASSERT(m_info->m_ncomp == static_cast(fields.size())); - auto& fop = *(m_info->m_fillpatch_op); const int nlevels = m_repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { amrex::Array mfabs = {AMREX_D_DECL( @@ -230,9 +231,11 @@ void Field::fillpatch_sibling_fields( } } - fop.fillpatch_sibling_fields( - lev, time, mfabs, mfabs, cfabs, ng, m_info->m_bcrec, - m_info->m_bcrec, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillpatch_sibling_fields( + lev, time, mfabs, mfabs, cfabs, ng, m_info->m_bcrec, + m_info->m_bcrec, field_state()); + } } } @@ -243,22 +246,24 @@ void Field::fillphysbc( const amrex::IntVect& ng) { BL_PROFILE("kynema-sgf::Field::fillphysbc"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); - fop.fillphysbc(lev, time, mfab, ng, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillphysbc(lev, time, mfab, ng, field_state()); + } } void Field::fillphysbc(const amrex::Real time, const amrex::IntVect ng) { BL_PROFILE("kynema-sgf::Field::fillphysbc"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); const int nlevels = m_repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { - fop.fillphysbc( - lev, time, m_repo.get_multifab(m_id, lev), ng, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillphysbc( + lev, time, m_repo.get_multifab(m_id, lev), ng, field_state()); + } } } @@ -268,17 +273,18 @@ void Field::fillphysbc_type( const amrex::Real time, const amrex::BCType::mathematicalBndryTypes bctype) { BL_PROFILE("kynema-sgf::Field::fillphysbc_type"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); // BC does not need to be initialized to fill BCs with a specified type, but // it does need to be copied to device (e.g., if requested type needs data) BL_ASSERT(m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); const int nlevels = m_repo.num_active_levels(); const auto ng = num_grow(); for (int lev = 0; lev < nlevels; ++lev) { - fop.fillphysbc_type( - lev, time, bctype, m_repo.get_multifab(m_id, lev), ng, - field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->fillphysbc_type( + lev, time, bctype, m_repo.get_multifab(m_id, lev), ng, + field_state()); + } } } @@ -297,10 +303,11 @@ void Field::set_inflow( const amrex::IntVect& ng) { BL_PROFILE("kynema-sgf::Field::set_inflow"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); - fop.set_inflow(lev, time, mfab, ng, field_state()); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->set_inflow(lev, time, mfab, ng, field_state()); + } } void Field::set_inflow_sibling_fields( @@ -309,10 +316,11 @@ void Field::set_inflow_sibling_fields( const amrex::Array mfabs) { BL_PROFILE("kynema-sgf::Field::set_inflow_sibling_fields"); - BL_ASSERT(m_info->m_fillpatch_op); + BL_ASSERT(!m_info->m_fillpatch_ops.empty()); BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - auto& fop = *(m_info->m_fillpatch_op); - fop.set_inflow_sibling_fields(lev, time, mfabs); + for (const auto& fop : m_info->m_fillpatch_ops) { + fop->set_inflow_sibling_fields(lev, time, mfabs); + } } void Field::advance_states() @@ -396,7 +404,7 @@ void Field::set_default_fillpatch_bc( bc_op(); } - if (!m_info->m_fillpatch_op) { + if (m_info->m_fillpatch_ops.empty()) { register_fill_patch_op>( repo().mesh(), time); } diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index fb8a0bff80..899f48db97 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -306,7 +306,7 @@ void incflo::ApplyProjection( amrex::Vector vel; for (int lev = 0; lev <= finest_level; ++lev) { vel.push_back(&(velocity(lev))); - vel[lev]->setBndry(0.0_rt); + vel[lev]->setDomainBndry(0.0_rt, geom[lev]); if (!proj_for_small_dt and !incremental) { kynema_sgf::nodal_projection::set_inflow_velocity( m_sim.field_boundaries(), velocity, lev, time, *vel[lev], 1); diff --git a/src/utilities/MultiLevelVector.H b/src/utilities/MultiLevelVector.H index a1d967bd7e..6c7edcd8d9 100644 --- a/src/utilities/MultiLevelVector.H +++ b/src/utilities/MultiLevelVector.H @@ -43,6 +43,11 @@ public: return m_data_d[lev]; }; + amrex::Gpu::DeviceVector& device_data(const int lev) + { + return m_data_d[lev]; + }; + void copy_host_to_device(); void copy_to_field(Field& fld); diff --git a/src/utilities/diagnostics.cpp b/src/utilities/diagnostics.cpp index df26e15039..53d4b0c838 100644 --- a/src/utilities/diagnostics.cpp +++ b/src/utilities/diagnostics.cpp @@ -108,16 +108,23 @@ amrex::Real kynema_sgf::diagnostics::get_vel_max( const int vdir, const amrex::Real factor) { + const int ng = 1; return amrex::ReduceMax( - vel, level_mask, 0, + vel, level_mask, ng, [=] AMREX_GPU_HOST_DEVICE( amrex::Box const& bx, amrex::Array4 const& vel_arr, amrex::Array4 const& mask_arr) -> amrex::Real { amrex::Real max_fab = -1.0e8_rt; amrex::Loop(bx, [=, &max_fab](int i, int j, int k) { + const int ii = amrex::max( + bx.smallEnd(0) + ng, amrex::min(i, bx.bigEnd(0) - ng)); + const int jj = amrex::max( + bx.smallEnd(1) + ng, amrex::min(j, bx.bigEnd(1) - ng)); + const int kk = amrex::max( + bx.smallEnd(2) + ng, amrex::min(k, bx.bigEnd(2) - ng)); max_fab = amrex::max( - max_fab, mask_arr(i, j, k) > 0 + max_fab, mask_arr(ii, jj, kk) > 0 ? factor * vel_arr(i, j, k, vdir) : std::numeric_limits::lowest()); }); @@ -150,8 +157,9 @@ amrex::Real kynema_sgf::diagnostics::get_vel_loc( const amrex::GpuArray problo, const amrex::GpuArray dx) { + const int ng = 1; return amrex::ReduceMax( - vel, level_mask, 0, + vel, level_mask, ng, [=] AMREX_GPU_HOST_DEVICE( amrex::Box const& bx, amrex::Array4 const& vel_arr, @@ -159,9 +167,15 @@ amrex::Real kynema_sgf::diagnostics::get_vel_loc( amrex::Real loc_fab = problo[ldir]; amrex::Loop(bx, [=, &loc_fab](int i, int j, int k) { int idx = (ldir == 0 ? i : (ldir == 1 ? j : k)); + const int ii = amrex::max( + bx.smallEnd(0) + ng, amrex::min(i, bx.bigEnd(0) - ng)); + const int jj = amrex::max( + bx.smallEnd(1) + ng, amrex::min(j, bx.bigEnd(1) - ng)); + const int kk = amrex::max( + bx.smallEnd(2) + ng, amrex::min(k, bx.bigEnd(2) - ng)); amrex::Real offset = 0.5_rt; amrex::Real loc = problo[ldir] + ((idx + offset) * dx[ldir]); - bool mask_check = (mask_arr(i, j, k) > 0); + bool mask_check = (mask_arr(ii, jj, kk) > 0); bool loc_check = (amrex::Math::abs(vel_max - vel_arr(i, j, k, vdir)) < std::numeric_limits::epsilon() * 1.0e6_rt); diff --git a/unit_tests/boundary_conditions/CMakeLists.txt b/unit_tests/boundary_conditions/CMakeLists.txt index 4132aa45d1..85a0e88c3a 100644 --- a/unit_tests/boundary_conditions/CMakeLists.txt +++ b/unit_tests/boundary_conditions/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(${kynema_sgf_unit_test_exe_name} PRIVATE # test cases test_log_law.cpp test_mosd.cpp + test_flather.cpp ) diff --git a/unit_tests/boundary_conditions/test_flather.cpp b/unit_tests/boundary_conditions/test_flather.cpp new file mode 100644 index 0000000000..c2a84c244c --- /dev/null +++ b/unit_tests/boundary_conditions/test_flather.cpp @@ -0,0 +1,225 @@ +#include "gtest/gtest.h" +#include "ks_test_utils/MeshTest.H" +#include "src/boundary_conditions/field_boundary_fill/Flather.H" +#include "AMReX_REAL.H" + +using namespace amrex::literals; + +namespace kynema_sgf_tests { + +namespace { +void initialize_vof( + kynema_sgf::Field& vof, + const amrex::Vector& geom, + amrex::Real wlev) +{ + for (int lev = 0; lev < vof.repo().num_active_levels(); ++lev) { + auto& vof_mfab = vof(lev); + auto vof_arrs = vof_mfab.arrays(); + const auto& zlo = geom[lev].ProbLo(2); + const auto& dz = geom[lev].CellSize(2); + amrex::ParallelFor( + vof_mfab, amrex::IntVect(1), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) { + const amrex::Real z = zlo + (k + 0.5_rt) * dz; + + if (z + 0.5_rt * dz <= wlev) { + vof_arrs[nbx](i, j, k) = 1.0_rt; + } else if (z - 0.5_rt * dz >= wlev) { + vof_arrs[nbx](i, j, k) = 0.0_rt; + } else { + vof_arrs[nbx](i, j, k) = (wlev - (z - 0.5_rt * dz)) / dz; + } + }); + } +} +} // namespace + +class FlatherBoundaryAverageTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + + { + amrex::ParmParse pp("geometry"); + amrex::Vector problo{{0.0_rt, 0.0_rt, 0.0_rt}}; + amrex::Vector probhi{{8.0_rt, 8.0_rt, 8.0_rt}}; + pp.addarr("prob_lo", problo); + pp.addarr("prob_hi", probhi); + amrex::Vector periodic{{0, 0, 0}}; + pp.addarr("is_periodic", periodic); + } + { + amrex::ParmParse pp("amr"); + const amrex::Vector ncell{{m_nx, m_nx, m_nx}}; + pp.add("max_level", 1); + pp.add("max_grid_size", m_nx); + pp.add("blocking_factor", 2); + pp.addarr("n_cell", ncell); + } + + std::stringstream ss; + ss << "1 // Number of levels" << '\n'; + ss << "1 // Number of boxes at this level" << '\n'; + ss << "0 0 2 4 6 6" << '\n'; + + create_mesh_instance(); + std::unique_ptr box_refine( + new kynema_sgf::CartBoxRefinement(sim())); + box_refine->read_inputs(mesh(), ss); + + if (mesh() != nullptr) { + mesh()->refine_criteria_vec().push_back( + std::move(box_refine)); + } + } + + const int m_nx{32}; + const amrex::Real m_wlev{4.0_rt}; +}; + +TEST_F(FlatherBoundaryAverageTest, accumulate_boundary_multilevel) +{ + constexpr amrex::Real u0 = 2.0_rt; + constexpr amrex::Real v0 = 3.0_rt; + constexpr amrex::Real tol = + std::numeric_limits::epsilon() * 1.0e4_rt; + + populate_parameters(); + initialize_mesh(); + + auto& repo = mesh().field_repo(); + auto& velocity = repo.declare_field("velocity", 3, 1); + repo.declare_face_normal_field({"u_mac", "v_mac", "w_mac"}, 1, 1, 1); + auto& vof = repo.declare_field("vof", 1, 1); + + velocity.setVal(u0, 0, 1, 1); + velocity.setVal(v0, 1, 1, 1); + velocity.setVal(0.0_rt, 2, 1, 1); + initialize_vof(vof, mesh().Geom(), m_wlev); + + kynema_sgf::Flather flather(sim()); + + kynema_sgf::MultiLevelVector xlo_uavg; + kynema_sgf::MultiLevelVector xlo_havg; + kynema_sgf::MultiLevelVector xhi_uavg; + kynema_sgf::MultiLevelVector xhi_havg; + kynema_sgf::MultiLevelVector ylo_uavg; + kynema_sgf::MultiLevelVector ylo_havg; + kynema_sgf::MultiLevelVector yhi_uavg; + kynema_sgf::MultiLevelVector yhi_havg; + + xlo_uavg.resize(0, mesh().Geom()); + xlo_havg.resize(0, mesh().Geom()); + xhi_uavg.resize(0, mesh().Geom()); + xhi_havg.resize(0, mesh().Geom()); + ylo_uavg.resize(1, mesh().Geom()); + ylo_havg.resize(1, mesh().Geom()); + yhi_uavg.resize(1, mesh().Geom()); + yhi_havg.resize(1, mesh().Geom()); + + const int nlevels = repo.num_active_levels(); + EXPECT_EQ(nlevels, 2); + + for (int lev = 0; lev < nlevels; ++lev) { + flather.accumulate_boundary( + lev, 0, 0, true, xlo_uavg, xlo_havg, false, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 0, 0, false, xhi_uavg, xhi_havg, false, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 1, 0, true, ylo_uavg, ylo_havg, false, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 1, 0, false, yhi_uavg, yhi_havg, false, + kynema_sgf::FieldState::New); + + const auto xhi_idx = xhi_uavg.ncells(lev) - 1; + const auto yhi_idx = yhi_uavg.ncells(lev) - 1; + + EXPECT_NEAR(xlo_uavg.host_data(lev)[0], u0, tol); + EXPECT_NEAR(xhi_uavg.host_data(lev)[xhi_idx], u0, tol); + EXPECT_NEAR(ylo_uavg.host_data(lev)[0], v0, tol); + EXPECT_NEAR(yhi_uavg.host_data(lev)[yhi_idx], v0, tol); + + EXPECT_NEAR(xlo_havg.host_data(lev)[0], m_wlev, tol); + EXPECT_NEAR(xhi_havg.host_data(lev)[xhi_idx], m_wlev, tol); + EXPECT_NEAR(ylo_havg.host_data(lev)[0], m_wlev, tol); + EXPECT_NEAR(yhi_havg.host_data(lev)[yhi_idx], m_wlev, tol); + } +} + +TEST_F( + FlatherBoundaryAverageTest, accumulate_boundary_multilevel_boundary_cells) +{ + constexpr amrex::Real u0 = 4.0_rt; + constexpr amrex::Real v0 = 1.5_rt; + constexpr amrex::Real tol = + std::numeric_limits::epsilon() * 1.0e4_rt; + + populate_parameters(); + initialize_mesh(); + + auto& repo = mesh().field_repo(); + auto& velocity = repo.declare_field("velocity", 3, 1); + repo.declare_face_normal_field({"u_mac", "v_mac", "w_mac"}, 1, 1, 1); + auto& vof = repo.declare_field("vof", 1, 1); + + velocity.setVal(u0, 0, 1, 1); + velocity.setVal(v0, 1, 1, 1); + velocity.setVal(0.0_rt, 2, 1, 1); + initialize_vof(vof, mesh().Geom(), m_wlev); + + kynema_sgf::Flather flather(sim()); + + kynema_sgf::MultiLevelVector xlo_uavg; + kynema_sgf::MultiLevelVector xlo_havg; + kynema_sgf::MultiLevelVector xhi_uavg; + kynema_sgf::MultiLevelVector xhi_havg; + kynema_sgf::MultiLevelVector ylo_uavg; + kynema_sgf::MultiLevelVector ylo_havg; + kynema_sgf::MultiLevelVector yhi_uavg; + kynema_sgf::MultiLevelVector yhi_havg; + + xlo_uavg.resize(0, mesh().Geom()); + xlo_havg.resize(0, mesh().Geom()); + xhi_uavg.resize(0, mesh().Geom()); + xhi_havg.resize(0, mesh().Geom()); + ylo_uavg.resize(1, mesh().Geom()); + ylo_havg.resize(1, mesh().Geom()); + yhi_uavg.resize(1, mesh().Geom()); + yhi_havg.resize(1, mesh().Geom()); + + const int nlevels = repo.num_active_levels(); + EXPECT_EQ(nlevels, 2); + + for (int lev = 0; lev < nlevels; ++lev) { + flather.accumulate_boundary( + lev, 0, 0, true, xlo_uavg, xlo_havg, true, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 0, 0, false, xhi_uavg, xhi_havg, true, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 1, 0, true, ylo_uavg, ylo_havg, true, + kynema_sgf::FieldState::New); + flather.accumulate_boundary( + lev, 1, 0, false, yhi_uavg, yhi_havg, true, + kynema_sgf::FieldState::New); + + EXPECT_NEAR(xlo_uavg.host_data(lev)[0], u0, tol); + EXPECT_NEAR(xhi_uavg.host_data(lev)[0], u0, tol); + EXPECT_NEAR(ylo_uavg.host_data(lev)[0], v0, tol); + EXPECT_NEAR(yhi_uavg.host_data(lev)[0], v0, tol); + + EXPECT_NEAR(xlo_havg.host_data(lev)[0], m_wlev, tol); + EXPECT_NEAR(xhi_havg.host_data(lev)[0], m_wlev, tol); + EXPECT_NEAR(ylo_havg.host_data(lev)[0], m_wlev, tol); + EXPECT_NEAR(yhi_havg.host_data(lev)[0], m_wlev, tol); + } +} + +} // namespace kynema_sgf_tests