diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.H b/amr-wind/equation_systems/icns/source_terms/DragForcing.H index 283ccd7bdc..2561640a86 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.H @@ -34,6 +34,7 @@ private: const CFDSim& m_sim; const amrex::AmrCore& m_mesh; const Field& m_velocity; + const Field* m_target_vel{nullptr}; amrex::Gpu::DeviceVector m_device_vel_ht; amrex::Gpu::DeviceVector m_device_vel_vals; amrex::Real m_drag_coefficient{10.0}; @@ -48,6 +49,7 @@ private: int m_sponge_south{0}; int m_sponge_north{1}; bool m_is_laminar{false}; + bool m_terrain_is_waves{false}; }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 2a6c6654c1..0d81284482 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -4,6 +4,7 @@ #include "AMReX_Gpu.H" #include "AMReX_Random.H" #include "amr-wind/wind_energy/ABL.H" +#include "amr-wind/physics/TerrainDrag.H" #include "amr-wind/utilities/linear_interpolation.H" namespace amr_wind::pde::icns { @@ -42,6 +43,13 @@ DragForcing::DragForcing(const CFDSim& sim) } else { m_sponge_strength = 0.0; } + if (phy_mgr.contains("OceanWaves")) { + const auto terrain_phys = + m_sim.physics_manager().get(); + const auto target_vel_name = terrain_phys.wave_velocity_field_name(); + m_target_vel = &sim.repo().get_field(target_vel_name); + m_terrain_is_waves = true; + } } DragForcing::~DragForcing() = default; @@ -68,6 +76,12 @@ void DragForcing::operator()( const auto& drag = (*m_terrain_drag)(lev).const_array(mfi); auto* const m_terrainz0 = &this->m_sim.repo().get_field("terrainz0"); const auto& terrainz0 = (*m_terrainz0)(lev).const_array(mfi); + + const bool is_waves = m_terrain_is_waves; + const auto& target_vel_arr = is_waves + ? (*m_target_vel)(lev).const_array(mfi) + : amrex::Array4(); + const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); @@ -158,20 +172,29 @@ void DragForcing::operator()( Dyz = -ustar * ustar * uy1 / (tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2]; } + amrex::Real target_u = 0.; + amrex::Real target_v = 0.; + amrex::Real target_w = 0.; + if (is_waves) { + target_u = target_vel_arr(i, j, k, 0); + target_v = target_vel_arr(i, j, k, 1); + target_w = target_vel_arr(i, j, k, 2); + } + const amrex::Real CdM = std::min(Cd / (m + tiny), cd_max / scale_factor); src_term(i, j, k, 0) -= - (CdM * m * ux1 * blank(i, j, k) + Dxz * drag(i, j, k) + + (CdM * m * (ux1 - target_u) * blank(i, j, k) + Dxz * drag(i, j, k) + bc_forcing_x * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (ux1 - sponge_density * spongeVelX)); src_term(i, j, k, 1) -= - (CdM * m * uy1 * blank(i, j, k) + Dyz * drag(i, j, k) + + (CdM * m * (uy1 - target_v) * blank(i, j, k) + Dyz * drag(i, j, k) + bc_forcing_y * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uy1 - sponge_density * spongeVelY)); src_term(i, j, k, 2) -= - (CdM * m * uz1 * blank(i, j, k) + + (CdM * m * (uz1 - target_w) * blank(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uz1 - sponge_density * spongeVelZ)); }); diff --git a/amr-wind/ocean_waves/OceanWaves.H b/amr-wind/ocean_waves/OceanWaves.H index a3850589e7..35669a75b0 100644 --- a/amr-wind/ocean_waves/OceanWaves.H +++ b/amr-wind/ocean_waves/OceanWaves.H @@ -52,6 +52,8 @@ protected: private: void relaxation_zones(); + void update_target_volume_fraction(); + CFDSim& m_sim; //! Unique pointer to the ocean waves model @@ -68,6 +70,11 @@ private: //! Ocean waves pressure gradient // Field& m_ow_pressure; + + //! Multiphase mode forces the velocity and vof in the domain + //! When off, single-phase mode simply provides the target variables, and + //! forcing is handled in other parts of the code + bool m_multiphase_mode{true}; }; } // namespace ocean_waves diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 86372b2c7f..364d67e126 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -18,7 +18,7 @@ OceanWaves::OceanWaves(CFDSim& sim) sim.repo().declare_field("ow_velocity", AMREX_SPACEDIM, 3, 1)) { if (!sim.physics_manager().contains("MultiPhase")) { - amrex::Abort("OceanWaves requires Multiphase physics to be active"); + m_multiphase_mode = false; } m_ow_levelset.set_default_fillpatch_bc(sim.time()); m_ow_vof.set_default_fillpatch_bc(sim.time()); @@ -32,6 +32,13 @@ void OceanWaves::pre_init_actions() BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_init_actions"); amrex::ParmParse pp(identifier()); + if (!(m_multiphase_mode || + m_sim.physics_manager().contains("TerrainDrag"))) { + amrex::Abort( + "OceanWaves requires MultiPhase or TerrainDrag physics to be " + "active"); + } + std::string label; pp.query("label", label); const std::string& tname = label; @@ -54,13 +61,18 @@ void OceanWaves::pre_init_actions() void OceanWaves::initialize_fields(int level, const amrex::Geometry& geom) { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::initialize_fields"); - m_owm->init_waves(level, geom); + m_owm->init_waves(level, geom, m_multiphase_mode); } void OceanWaves::post_init_actions() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions"); - relaxation_zones(); + if (m_multiphase_mode) { + relaxation_zones(); + } else { + m_owm->update_target_fields(true); + m_owm->update_target_volume_fraction(); + } } void OceanWaves::post_regrid_actions() @@ -72,12 +84,18 @@ void OceanWaves::post_regrid_actions() void OceanWaves::pre_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work"); + if (!m_multiphase_mode) { + m_owm->update_target_fields(true); + m_owm->update_target_volume_fraction(); + } } void OceanWaves::post_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions"); - relaxation_zones(); + if (m_multiphase_mode) { + relaxation_zones(); + } } /** Update ocean waves relaxation zones @@ -85,12 +103,19 @@ void OceanWaves::post_advance_work() */ void OceanWaves::relaxation_zones() { - BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_relaxation_zones"); - m_owm->update_relax_zones(); + BL_PROFILE("amr-wind::ocean_waves::OceanWaves::relaxation_zones"); + m_owm->update_target_fields(false); m_owm->apply_relax_zones(); m_owm->reset_regrid_flag(); } +void OceanWaves::update_target_volume_fraction() +{ + BL_PROFILE( + "amr-wind::ocean_waves::OceanWaves::update_target_volume_fraction"); + m_owm->update_target_volume_fraction(); +} + void OceanWaves::prepare_outputs() { const std::string post_dir = m_sim.io_manager().post_processing_directory(); diff --git a/amr-wind/ocean_waves/OceanWavesModel.H b/amr-wind/ocean_waves/OceanWavesModel.H index 86fa87bdbf..0a07fce81c 100644 --- a/amr-wind/ocean_waves/OceanWavesModel.H +++ b/amr-wind/ocean_waves/OceanWavesModel.H @@ -32,12 +32,14 @@ public: virtual void read_inputs(const ::amr_wind::utils::MultiParser&) = 0; - virtual void init_waves(int, const amrex::Geometry&) = 0; + virtual void init_waves(int, const amrex::Geometry&, bool) = 0; - virtual void update_relax_zones() = 0; + virtual void update_target_fields(bool) = 0; virtual void apply_relax_zones() = 0; + virtual void update_target_volume_fraction() = 0; + virtual void prepare_outputs(const std::string&) = 0; virtual void write_outputs() = 0; @@ -91,9 +93,9 @@ public: m_out_op.read_io_options(pp); } - void update_relax_zones() override + void update_target_fields(bool for_forcing_term) override { - ops::UpdateRelaxZonesOp()(m_data); + ops::UpdateTargetFieldsOp()(m_data, for_forcing_term); } void apply_relax_zones() override @@ -101,6 +103,11 @@ public: ops::ApplyRelaxZonesOp()(m_data); } + void update_target_volume_fraction() override + { + ops::UpdateTargetVolumeFractionOp()(m_data); + } + void prepare_outputs(const std::string& out_dir) override { m_out_op.prepare_outputs(out_dir); @@ -112,9 +119,11 @@ public: void write_outputs() override { m_out_op.write_outputs(); } - void init_waves(int level, const amrex::Geometry& geom) override + void init_waves( + int level, const amrex::Geometry& geom, bool multiphase_mode) override { - ops::InitDataOp()(m_data, level, geom); + ops::InitDataOp()( + m_data, level, geom, multiphase_mode); } }; diff --git a/amr-wind/ocean_waves/OceanWavesOps.H b/amr-wind/ocean_waves/OceanWavesOps.H index bb48cc8d3f..529d508465 100644 --- a/amr-wind/ocean_waves/OceanWavesOps.H +++ b/amr-wind/ocean_waves/OceanWavesOps.H @@ -32,11 +32,14 @@ template struct InitDataOp; template -struct UpdateRelaxZonesOp; +struct UpdateTargetFieldsOp; template struct ApplyRelaxZonesOp; +template +struct UpdateTargetVolumeFractionOp; + template struct ProcessOutputsOp; diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index 8b624d3a5b..afa16dc584 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -1,7 +1,6 @@ #ifndef LINEAR_WAVES_OPS_H #define LINEAR_WAVES_OPS_H -#include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/ocean_waves/relaxation_zones/LinearWaves.H" #include "amr-wind/ocean_waves/OceanWavesOps.H" #include "amr-wind/ocean_waves/OceanWaves.H" @@ -19,8 +18,11 @@ struct ReadInputsOp auto& wdata = data.meta(); auto& info = data.info(); relaxation_zones::read_inputs(wdata, info, pp); - // Get gravity from MultiPhase physics, assume negative z - wdata.g = -data.sim().physics_manager().get().gravity()[2]; + // Get gravity, assume negative z + amrex::Vector gravity{0.0, 0.0, -9.81}; + amrex::ParmParse pp_incflo("incflo"); + pp_incflo.queryarr("gravity", gravity); + wdata.g = -gravity[2]; pp.get("wave_length", wdata.wave_length); pp.get("wave_height", wdata.wave_height); @@ -31,36 +33,43 @@ template <> struct InitDataOp { void operator()( - LinearWaves::DataType& data, int level, const amrex::Geometry& geom) + LinearWaves::DataType& data, + int level, + const amrex::Geometry& geom, + bool multiphase_mode) { const auto& wdata = data.meta(); auto& sim = data.sim(); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } // cppcheck-suppress constVariableReference - auto& m_levelset = sim.repo().get_field("levelset"); - // cppcheck-suppress constVariableReference - auto& m_velocity = sim.repo().get_field("velocity"); + auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - const auto& phi = m_levelset(level).arrays(); - const auto& vel = m_velocity(level).arrays(); + const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real zero_sea_level = wdata.zsl; const amrex::Real gen_length = wdata.gen_length; const amrex::Real beach_length = wdata.beach_length; const amrex::Real g = wdata.g; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; const amrex::Real water_depth = wdata.water_depth; amrex::ParallelFor( - m_velocity(level), amrex::IntVect(3), + velocity(level), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; @@ -98,32 +107,39 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); } -}; +}; // namespace amr_wind::ocean_waves::ops template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { - void operator()(LinearWaves::DataType& data) + void operator()(LinearWaves::DataType& data, bool for_forcing_term) { const auto& wdata = data.meta(); auto& sim = data.sim(); - const auto& time = sim.time().new_time(); + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // cppcheck-suppress constVariableReference - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); // cppcheck-suppress constVariableReference - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); auto geom = sim.mesh().Geom(); @@ -132,8 +148,8 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - const auto& phi = m_ow_levelset(lev).arrays(); - const auto& vel = m_ow_velocity(lev).arrays(); + const auto& phi = ow_levelset(lev).arrays(); + const auto& vel = ow_velocity(lev).arrays(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -142,7 +158,7 @@ struct UpdateRelaxZonesOp const amrex::Real g = wdata.g; amrex::ParallelFor( - m_ow_velocity(lev), amrex::IntVect(3), + ow_velocity(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H index e015c61057..c939ff43ab 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H @@ -28,6 +28,8 @@ void init_data_structures(CFDSim&); */ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata); +void update_target_vof(CFDSim& sim); + void prepare_netcdf_file( const std::string& /*ncfile*/, const RelaxZonesBaseData& /*meta*/, @@ -51,7 +53,7 @@ struct UseDefaultOp }; template -struct UpdateRelaxZonesOp< +struct UpdateTargetFieldsOp< WaveTheoryTrait, typename std::enable_if_t< std::is_base_of_v>> @@ -76,6 +78,22 @@ struct ApplyRelaxZonesOp< } }; +template +struct UpdateTargetVolumeFractionOp< + WaveTheoryTrait, + typename std::enable_if_t< + std::is_base_of_v>> +{ + void operator()(typename WaveTheoryTrait::DataType& data) + { + BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_target_vof"); + + auto& sim = data.sim(); + + relaxation_zones::update_target_vof(sim); + } +}; + template struct ProcessOutputsOp< WaveTheoryTrait, diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp index 4c98d83a0b..82f84a9728 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp @@ -228,6 +228,27 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) density.fillpatch(time); } +void update_target_vof(CFDSim& sim) +{ + const int nlevels = sim.repo().num_active_levels(); + auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& m_ow_vof = sim.repo().get_field("ow_vof"); + const auto& geom = sim.mesh().Geom(); + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = geom[lev].CellSizeArray(); + const auto target_phi = m_ow_levelset(lev).arrays(); + auto target_volfrac = m_ow_vof(lev).arrays(); + const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]); + amrex::ParallelFor( + m_ow_vof(lev), amrex::IntVect(2), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + target_volfrac[nbx](i, j, k) = + multiphase::levelset_to_vof(i, j, k, eps, target_phi[nbx]); + }); + } +} + void prepare_netcdf_file( const std::string& ncfile, const RelaxZonesBaseData& meta, diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 0ff477842c..706b3f5cdc 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -1,7 +1,6 @@ #ifndef STOKES_WAVES_OPS_H #define STOKES_WAVES_OPS_H -#include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/ocean_waves/relaxation_zones/StokesWaves.H" #include "amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H" #include "amr-wind/ocean_waves/OceanWavesOps.H" @@ -21,8 +20,12 @@ struct ReadInputsOp auto& info = data.info(); relaxation_zones::read_inputs(wdata, info, pp); - // Get gravity from MultiPhase physics, assume negative z - wdata.g = -data.sim().physics_manager().get().gravity()[2]; + // Get gravity, assume negative z + amrex::Vector gravity{0.0, 0.0, -9.81}; + amrex::ParmParse pp_incflo("incflo"); + pp_incflo.queryarr("gravity", gravity); + wdata.g = -gravity[2]; + // Get wave attributes pp.get("wave_height", wdata.wave_height); pp.get("order", wdata.order); @@ -72,21 +75,28 @@ template <> struct InitDataOp { void operator()( - StokesWaves::DataType& data, int level, const amrex::Geometry& geom) + StokesWaves::DataType& data, + int level, + const amrex::Geometry& geom, + bool multiphase_mode) { const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariableReference - auto& m_levelset = sim.repo().get_field("levelset"); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } // cppcheck-suppress constVariableReference auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - const auto& phi = m_levelset(level).arrays(); const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -96,8 +106,8 @@ struct InitDataOp const amrex::Real beach_length = wdata.beach_length; const amrex::Real g = wdata.g; const int order = wdata.order; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; amrex::ParallelFor( velocity(level), amrex::IntVect(3), @@ -123,32 +133,39 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); } }; template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { - void operator()(StokesWaves::DataType& data) + void operator()(StokesWaves::DataType& data, bool for_forcing_term) { const auto& wdata = data.meta(); auto& sim = data.sim(); - const auto& time = sim.time().new_time(); + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // cppcheck-suppress constVariableReference - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); // cppcheck-suppress constVariableReference - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); auto geom = sim.mesh().Geom(); @@ -157,8 +174,8 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - const auto& phi = m_ow_levelset(lev).arrays(); - const auto& vel = m_ow_velocity(lev).arrays(); + const auto& phi = ow_levelset(lev).arrays(); + const auto& vel = ow_velocity(lev).arrays(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -168,7 +185,7 @@ struct UpdateRelaxZonesOp const int order = wdata.order; amrex::ParallelFor( - m_ow_velocity(lev), amrex::IntVect(3), + ow_velocity(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index b7c53a4c9d..ace4308cb9 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -23,7 +23,7 @@ int evaluate_read_resize( const int new_ntime, const amrex::Real wtinit, const amrex::Real wdt, - const amrex::Real newtime) + const amrex::Real time) { // Flag to indicate that data must be read twice for interpolation int double_data = 0; @@ -41,7 +41,7 @@ int evaluate_read_resize( // Sim time to go with recorded data wtime = new_ntime * wdt + wtinit; // If double reading is deemed necessary, check for convenience - if (double_data == 1 && std::abs(wtime - newtime) <= 1e-10) { + if (double_data == 1 && std::abs(wtime - time) <= 1e-10) { // Reading can be done just once, w2a fields replace ow fields double_data = 2; } @@ -54,7 +54,7 @@ int evaluate_read_resize( resize_flag = true; // Confirm that new time is not coincident with modes time step - if (std::abs(wtime - newtime) > 1e-10) { + if (std::abs(wtime - time) > 1e-10) { // Data must be read before and after wtime double_data = 1; } else { @@ -72,7 +72,7 @@ int evaluate_read_resize( // interpolation is ready to go } // Record latest time as 'last' for next timestep - t_last = newtime; + t_last = time; // Return flag regarding double reading return double_data; } @@ -294,7 +294,8 @@ struct ReadInputsOp if (std::abs(depth - (wdata.zsl - prob_lo_input[2])) > 1e-3 * depth) { amrex::Print() << "WARNING: Mismatch between water depths from AMR-Wind " - "domain and HOS data interpreted by Waves2AMR"; + "domain and HOS data interpreted by Waves2AMR. \n ^This " + "warning is not a concern when using waves as terrain.\n"; } // Allocate pointers for FFTW @@ -332,59 +333,51 @@ struct ReadInputsOp std::to_string(flag)); } - // If init_wave_field is activated and initialization will be done, get - // modes on every processor - if (wdata.init_wave_field && data.sim().time().time_index() == 0) { - bool no_EOF = wdata.rmodes.get_data( + // Get modes on every processor (unless requested to skip?) + bool no_EOF = wdata.rmodes.get_data( + wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, + wdata.mFS); + if (!no_EOF) { + // End of file detected, reset reading + wdata.n_offset = update_offset_timestep(wdata.ntime, wdata.n_winit); + // Print warning to screen + amrex::Print() + << "WARNING (waves2amr_ops): end of mode data file " + "detected, resetting to beginning of mode data.\n"; + // Read data again, now from a valid timestep + no_EOF = wdata.rmodes.get_data( wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, wdata.mFS); + // If no valid data is detected at this point, abort if (!no_EOF) { - // End of file detected, reset reading - wdata.n_offset = - update_offset_timestep(wdata.ntime, wdata.n_winit); - // Print warning to screen - amrex::Print() - << "WARNING (waves2amr_ops): end of mode data file " - "detected, resetting to beginning of mode data.\n"; - // Read data again, now from a valid timestep - no_EOF = wdata.rmodes.get_data( - wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, - wdata.mFS); - // If no valid data is detected at this point, abort - if (!no_EOF) { - amrex::Abort( - "waves2amr_ops: end of mode data file detected after " - "resetting to beginning; please evaluate HOS_init_time " - "or HOS_init_timestep and check the length of the mode " - "file."); - } + amrex::Abort( + "waves2amr_ops: end of mode data file detected after " + "resetting to beginning; please evaluate HOS_init_time " + "or HOS_init_timestep and check the length of the mode " + "file."); } + } - // Convert modes to spatial data - modes_hosgrid::copy_complex( - wdata.n0, wdata.n1, wdata.mFS, wdata.eta_mptr); - wdata.sp_eta_vec.resize( - static_cast(wdata.n0) * static_cast(wdata.n1), - 0.0); - modes_hosgrid::populate_hos_eta( - wdata.rmodes, wdata.plan, wdata.eta_mptr, wdata.sp_eta_vec); - // Mesh is not yet created, so get data at every height - const auto n_hts = wdata.hvec.size(); - wdata.sp_u_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - wdata.sp_v_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - wdata.sp_w_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - for (int iht = 0; iht < static_cast(n_hts); ++iht) { - // Get sample height - amrex::Real ht = wdata.hvec[iht]; - // Sample velocity - modes_hosgrid::populate_hos_vel( - wdata.rmodes, ht, wdata.mX, wdata.mY, wdata.mZ, wdata.plan, - wdata.u_mptr, wdata.v_mptr, wdata.w_mptr, wdata.sp_u_vec, - wdata.sp_v_vec, wdata.sp_w_vec, iht * wdata.n0 * wdata.n1); - } + // Convert modes to spatial data + modes_hosgrid::copy_complex( + wdata.n0, wdata.n1, wdata.mFS, wdata.eta_mptr); + wdata.sp_eta_vec.resize( + static_cast(wdata.n0) * static_cast(wdata.n1), 0.0); + modes_hosgrid::populate_hos_eta( + wdata.rmodes, wdata.plan, wdata.eta_mptr, wdata.sp_eta_vec); + // Mesh is not yet created, so get data at every height + const auto n_hts = wdata.hvec.size(); + wdata.sp_u_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + wdata.sp_v_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + wdata.sp_w_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + for (int iht = 0; iht < static_cast(n_hts); ++iht) { + // Get sample height + amrex::Real ht = wdata.hvec[iht]; + // Sample velocity + modes_hosgrid::populate_hos_vel( + wdata.rmodes, ht, wdata.mX, wdata.mY, wdata.mZ, wdata.plan, + wdata.u_mptr, wdata.v_mptr, wdata.w_mptr, wdata.sp_u_vec, + wdata.sp_v_vec, wdata.sp_w_vec, iht * wdata.n0 * wdata.n1); } // Declare fields for HOS @@ -406,7 +399,8 @@ struct InitDataOp void // cppcheck-suppress constParameterReference operator()( - W2AWaves::DataType & data, int level, const amrex::Geometry & geom) + W2AWaves::DataType & + data, int level, const amrex::Geometry & geom, bool multiphase_mode) { #ifdef AMR_WIND_USE_W2A @@ -417,9 +411,14 @@ struct InitDataOp // Fill ow fields, then populate flow fields according to setup auto& ow_levelset = sim.repo().get_field("ow_levelset"); auto& ow_velocity = sim.repo().get_field("ow_velocity"); - auto& levelset = sim.repo().get_field("levelset"); + auto& velocity = sim.repo().get_field("velocity"); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } + const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); @@ -448,15 +447,17 @@ struct InitDataOp // Populate flow fields according to intended forcing and init setup const auto& ow_phi = ow_levelset(level).const_arrays(); const auto& ow_vel = ow_velocity(level).const_arrays(); - const auto& phi = levelset(level).arrays(); const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real gen_length = wdata.gen_length; const amrex::Real beach_length = wdata.beach_length; const amrex::Real zero_sea_level = wdata.zsl; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; amrex::ParallelFor( velocity(level), amrex::IntVect(3), @@ -479,14 +480,17 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); // Start w2a fields at 0, some areas will not be modified @@ -495,28 +499,32 @@ struct InitDataOp w2a_levelset.setVal(0.0); w2a_velocity.setVal(0.0); #else - amrex::ignore_unused(data, level, geom); + amrex::ignore_unused(data, level, geom, multiphase_mode); #endif } }; // namespace ocean_waves template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { // cppcheck-suppress constParameterReference - void operator()(W2AWaves::DataType& data) + void operator()(W2AWaves::DataType& data, bool for_forcing_term) { #ifdef AMR_WIND_USE_W2A auto& wdata = data.meta(); auto& sim = data.sim(); - // Nudge the solution toward where it should be - const amrex::Real newtime = sim.time().new_time(); + // Provide target solution based on intended part of the algorithm + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // Update ow fields every time - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); // Update HOS fields when necessary auto& w2a_levelset = sim.repo().get_field("w2a_levelset"); auto& w2a_velocity = sim.repo().get_field("w2a_velocity"); @@ -531,10 +539,10 @@ struct UpdateRelaxZonesOp bool read_flag = false; // Check if time indicates reading must take place int new_ntime = - wdata.rmodes.time2step(newtime + wdata.t_winit, wdata.ntime); + wdata.rmodes.time2step(time + wdata.t_winit, wdata.ntime); int double_data = evaluate_read_resize( wdata.ntime, read_flag, wdata.resize_flag, wdata.t, wdata.t_last, - new_ntime, wdata.t_winit, wdata.dt_modes, newtime); + new_ntime, wdata.t_winit, wdata.dt_modes, time); // Check if stop time is exceeded, introduce offset to ntime if (read_flag) { // Need to only check when reading is occurring @@ -568,7 +576,7 @@ struct UpdateRelaxZonesOp // Get heights for this processor, check overlap in z flag_z = (interp_to_mfab::get_local_height_indices( - wdata.indvec, wdata.hvec, m_ow_velocity.vec_ptrs(), + wdata.indvec, wdata.hvec, ow_velocity.vec_ptrs(), geom) == 1); // No overlap from heights definitely means no interp @@ -576,17 +584,24 @@ struct UpdateRelaxZonesOp const int dir = 0; flag_xlo = (interp_to_mfab::check_lateral_overlap_lo( - wdata.gen_length, dir, m_ow_velocity.vec_ptrs(), - geom) == 1); + wdata.gen_length, dir, ow_velocity.vec_ptrs(), geom) == + 1); // No overlap with gen region means no interp, unless ... if (wdata.has_outprofile) { // ... if overlap exists here, needing interp flag_xhi = (interp_to_mfab::check_lateral_overlap_hi( - wdata.beach_length, dir, m_ow_velocity.vec_ptrs(), + wdata.beach_length, dir, ow_velocity.vec_ptrs(), geom) == 1); } + // Wave generation zones are irrelevant for single-phase mode, + // indicated by the use for a forcing term + if (for_forcing_term) { + flag_xlo = true; + flag_xhi = true; + } + if (flag_z && (flag_xlo || flag_xhi)) { // Interpolation is needed wdata.do_interp = true; @@ -620,30 +635,30 @@ struct UpdateRelaxZonesOp if (double_data == 1) { if (wdata.do_interp) { populate_fields_all_levels( - wdata, geom, m_ow_levelset, m_ow_velocity, -1); + wdata, geom, ow_levelset, ow_velocity, -1); } // Average down to get fine information on coarse grid where // possible (may be unnecessary) for (int lev = nlevels - 1; lev > 0; --lev) { amrex::average_down( - m_ow_velocity(lev), m_ow_velocity(lev - 1), 0, + ow_velocity(lev), ow_velocity(lev - 1), 0, AMREX_SPACEDIM, sim.mesh().refRatio(lev - 1)); amrex::average_down( - m_ow_levelset(lev), m_ow_levelset(lev - 1), 0, 1, + ow_levelset(lev), ow_levelset(lev - 1), 0, 1, sim.mesh().refRatio(lev - 1)); } // Fill patch to get correct ghost cells after average down - m_ow_velocity.fillpatch(sim.time().new_time()); - m_ow_levelset.fillpatch(sim.time().new_time()); + ow_velocity.fillpatch(sim.time().new_time()); + ow_levelset.fillpatch(sim.time().new_time()); // Prior t_last (corresponding to ow fields) t_last = (wdata.ntime - 1) * wdata.dt_modes; } else if (double_data == 2) { // Restarting simulation or taking a big step, new time at ntime // Initialize ow fields to 0 for time interp, will be replaced - m_ow_levelset.setVal(0.0); - m_ow_velocity.setVal(0.0); + ow_levelset.setVal(0.0); + ow_velocity.setVal(0.0); // No modification needed for t_last, leads to interp factor = 1 } @@ -670,32 +685,32 @@ struct UpdateRelaxZonesOp // Temporally interpolate at every timestep to get target solution for (int lev = 0; lev < nlevels; ++lev) { - auto phi = m_ow_levelset(lev).arrays(); - auto vel = m_ow_velocity(lev).arrays(); - auto W2A_phi = w2a_levelset(lev).arrays(); - auto W2A_vel = w2a_velocity(lev).arrays(); + auto phi = ow_levelset(lev).arrays(); + auto vel = ow_velocity(lev).arrays(); + const auto W2A_phi = w2a_levelset(lev).const_arrays(); + const auto W2A_vel = w2a_velocity(lev).const_arrays(); const amrex::Real W2A_t = wdata.t; amrex::ParallelFor( - m_ow_levelset(lev), amrex::IntVect(3), + ow_levelset(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { // Interpolate temporally every time phi[nbx](i, j, k) += (W2A_phi[nbx](i, j, k) - phi[nbx](i, j, k)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 0) += (W2A_vel[nbx](i, j, k, 0) - vel[nbx](i, j, k, 0)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 1) += (W2A_vel[nbx](i, j, k, 1) - vel[nbx](i, j, k, 1)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 2) += (W2A_vel[nbx](i, j, k, 2) - vel[nbx](i, j, k, 2)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); }); } #else - amrex::ignore_unused(data); + amrex::ignore_unused(data, for_forcing_term); #endif } }; diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 6497be1d53..7e9e485cd5 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -26,14 +26,23 @@ public: void pre_init_actions() override {} - void post_init_actions() override {} + void post_init_actions() override; void post_regrid_actions() override; - void pre_advance_work() override {} + void pre_advance_work() override; void post_advance_work() override {} + //! Terrain Drag for waves + + void convert_waves_to_terrain_fields(); + + std::string wave_velocity_field_name() const + { + return m_wave_velocity_name; + } + private: CFDSim& m_sim; const FieldRepo& m_repo; @@ -54,6 +63,14 @@ private: amrex::Vector m_xrough; amrex::Vector m_yrough; amrex::Vector m_z0rough; + + //! Terrain Drag for waves + bool m_terrain_is_waves{false}; + Field* m_wave_volume_fraction{nullptr}; + Field* m_wave_negative_elevation{nullptr}; + const std::string m_wave_volume_fraction_name{"ow_vof"}; + const std::string m_wave_negative_elevation_name{"ow_levelset"}; + const std::string m_wave_velocity_name{"ow_velocity"}; }; } // namespace amr_wind::terraindrag diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 2b6b661757..17a0cece59 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -22,22 +22,31 @@ TerrainDrag::TerrainDrag(CFDSim& sim) , m_terrainz0(sim.repo().declare_field("terrainz0", 1, 1, 1)) , m_terrain_height(sim.repo().declare_field("terrain_height", 1, 1, 1)) { - std::string terrain_file("terrain.amrwind"); - std::string roughness_file("terrain.roughness"); - amrex::ParmParse pp(identifier()); - pp.query("terrain_file", terrain_file); - pp.query("roughness_file", roughness_file); - - ioutils::read_flat_grid_file( - terrain_file, m_xterrain, m_yterrain, m_zterrain); - - // No checks for the file as it is optional currently - std::ifstream file(roughness_file, std::ios::in); - if (file.good()) { + + m_terrain_is_waves = sim.physics_manager().contains("OceanWaves"); + + if (!m_terrain_is_waves) { + std::string terrain_file("terrain.amrwind"); + std::string roughness_file("terrain.roughness"); + amrex::ParmParse pp(identifier()); + pp.query("terrain_file", terrain_file); + pp.query("roughness_file", roughness_file); + ioutils::read_flat_grid_file( - roughness_file, m_xrough, m_yrough, m_z0rough); + terrain_file, m_xterrain, m_yterrain, m_zterrain); + + // No checks for the file as it is optional currently + std::ifstream file(roughness_file, std::ios::in); + if (file.good()) { + ioutils::read_flat_grid_file( + roughness_file, m_xrough, m_yrough, m_z0rough); + } + file.close(); + } else { + m_wave_volume_fraction = &m_repo.get_field(m_wave_volume_fraction_name); + m_wave_negative_elevation = + &m_repo.get_field(m_wave_negative_elevation_name); } - file.close(); m_sim.io_manager().register_output_int_var("terrain_drag"); m_sim.io_manager().register_output_int_var("terrain_blank"); @@ -52,6 +61,10 @@ TerrainDrag::TerrainDrag(CFDSim& sim) void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) { + if (m_terrain_is_waves) { + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::initialize_fields"); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); @@ -59,6 +72,12 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto& terrainz0 = m_terrainz0(level); auto& terrain_height = m_terrain_height(level); auto& drag = m_terrain_drag(level); + + auto levelBlanking = blanking.arrays(); + auto levelDrag = drag.arrays(); + auto levelz0 = terrainz0.arrays(); + auto levelHeight = terrain_height.arrays(); + const auto xterrain_size = m_xterrain.size(); const auto yterrain_size = m_yterrain.size(); const auto zterrain_size = m_zterrain.size(); @@ -96,10 +115,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) const auto* xrough_ptr = device_xrough.data(); const auto* yrough_ptr = device_yrough.data(); const auto* z0rough_ptr = device_z0rough.data(); - auto levelBlanking = blanking.arrays(); - auto levelDrag = drag.arrays(); - auto levelz0 = terrainz0.arrays(); - auto levelheight = terrain_height.arrays(); + amrex::ParallelFor( blanking, m_terrain_blank.num_grow(), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { @@ -111,7 +127,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) yterrain_ptr + yterrain_size, zterrain_ptr, x, y); levelBlanking[nbx](i, j, k, 0) = static_cast((z <= terrainHt) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = + levelHeight[nbx](i, j, k, 0) = std::max(std::abs(z - terrainHt), 0.5 * dx[2]); amrex::Real roughz0 = 0.1; @@ -122,6 +138,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) } levelz0[nbx](i, j, k, 0) = roughz0; }); + amrex::ParallelFor( blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && @@ -132,11 +149,79 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) } }); } + +void TerrainDrag::post_init_actions() +{ + if (!m_terrain_is_waves) { + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::post_init_actions"); + convert_waves_to_terrain_fields(); +} + +void TerrainDrag::pre_advance_work() +{ + if (!m_terrain_is_waves) { + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::pre_advance_work"); + convert_waves_to_terrain_fields(); +} + void TerrainDrag::post_regrid_actions() +{ + if (m_terrain_is_waves) { + convert_waves_to_terrain_fields(); + } else { + const int nlevels = m_sim.repo().num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + initialize_fields(lev, m_sim.repo().mesh().Geom(lev)); + } + } +} + +void TerrainDrag::convert_waves_to_terrain_fields() { const int nlevels = m_sim.repo().num_active_levels(); - for (int lev = 0; lev < nlevels; ++lev) { - initialize_fields(lev, m_sim.repo().mesh().Geom(lev)); + // Uniform, low roughness for waves + m_terrainz0.setVal(1e-4); + for (int level = 0; level < nlevels; ++level) { + const auto geom = m_sim.repo().mesh().Geom(level); + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + auto& blanking = m_terrain_blank(level); + auto& terrain_height = m_terrain_height(level); + auto& drag = m_terrain_drag(level); + + auto levelBlanking = blanking.arrays(); + auto levelDrag = drag.arrays(); + auto levelHeight = terrain_height.arrays(); + + const auto negative_wave_elevation = + (*m_wave_negative_elevation)(level).const_arrays(); + const auto wave_vol_frac = + (*m_wave_volume_fraction)(level).const_arrays(); + + // Get terrain blanking from ocean waves fields + amrex::ParallelFor( + blanking, m_terrain_blank.num_grow(), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + levelBlanking[nbx](i, j, k, 0) = static_cast( + (wave_vol_frac[nbx](i, j, k) >= 0.5) && (z > prob_lo[2])); + levelHeight[nbx](i, j, k, 0) = + -negative_wave_elevation[nbx](i, j, k); + }); + amrex::ParallelFor( + blanking, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && + (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { + levelDrag[nbx](i, j, k, 0) = 1; + } else { + levelDrag[nbx](i, j, k, 0) = 0; + } + }); } } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8799e0eeee..3829504144 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -206,6 +206,7 @@ add_test_re(abl_unstable_local_wall_model) add_test_re(abl_unstable_schumann_wall_model) add_test_re(abl_anelastic) add_test_re(abl_multiphase_laminar) +add_test_re(abl_waves_terrain) add_test_re(act_abl_joukowskydisk) add_test_re(act_abl_uniformctdisk) add_test_re(act_fixed_wing) @@ -319,6 +320,7 @@ endif() if(AMR_WIND_ENABLE_W2A) add_test_re(ow_w2a) add_test_re(abl_multiphase_w2a) + add_test_red(abl_w2a_terrain abl_bndry_output_native) endif() #============================================================================= diff --git a/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp new file mode 100644 index 0000000000..9bbc0de2a1 --- /dev/null +++ b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp @@ -0,0 +1,98 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 12.25 # Max (simulated) time to evolve +time.max_step = 10 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.2 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +io.restart_file = "../abl_bndry_output_native/chk00005" +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.outputs = density velocity p ow_levelset ow_velocity + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0.0 0.0 -9.81 # Gravitational force (3D) +incflo.density = 1.225 +transport.model = ConstTransport +transport.viscosity = 1e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + +incflo.physics = ABL OceanWaves TerrainDrag +ICNS.source_terms = CoriolisForcing DragForcing + +ABL.reference_temperature = 290.0 +CoriolisForcing.east_vector = 1.0 0.0 0.0 +CoriolisForcing.north_vector = 0.0 1.0 0.0 +CoriolisForcing.latitude = 90.0 +CoriolisForcing.rotational_time_period = 125663.706143592 +GeostrophicForcing.geostrophic_wind = 10.0 0.0 0.0 +DragForcing.sponge_strength = 0 + +temperature.source_terms = DragTempForcing +DragTempForcing.reference_temperature = 265.0 + +incflo.velocity = 10.0 0.0 0.0 +ABL.temperature_heights = 0.0 2000.0 +ABL.temperature_values = 290.0 290.0 +ABL.perturb_temperature = false +ABL.cutoff_height = 50.0 +ABL.perturb_velocity = true +ABL.perturb_ref_height = 50.0 +ABL.Uperiods = 4.0 +ABL.Vperiods = 4.0 +ABL.deltaU = 1.0 +ABL.deltaV = 1.0 +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.01 +ABL.bndry_file = "../abl_bndry_output_native/bndry_files" +ABL.bndry_io_mode = 1 +ABL.bndry_var_names = velocity temperature +ABL.bndry_output_format = native + +OceanWaves.label = W2A1 +OceanWaves.W2A1.type = W2AWaves +OceanWaves.W2A1.HOS_modes_filename = ../ow_w2a/modes_HOS_SWENSE.dat +# Offset to include wave troughs in domain +OceanWaves.W2A1.zero_sea_level = 10 +# These variables should change with resolution in z +OceanWaves.W2A1.number_interp_points_in_z = 6 +OceanWaves.W2A1.interp_spacing_at_surface = 2. +OceanWaves.W2A1.number_interp_above_surface = 3 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 +amr.max_level = 2 +tagging.labels = sr +tagging.sr.type = CartBoxRefinement +tagging.sr.static_refinement_def = static_box.refine +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1) +# Boundary conditions +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.temperature = 0.0 +xhi.type = "pressure_outflow" +ylo.type = "mass_inflow" +ylo.density = 1.0 +ylo.temperature = 0.0 +yhi.type = "pressure_outflow" +zlo.type = "slip_wall" +zhi.type = "slip_wall" diff --git a/test/test_files/abl_w2a_terrain/static_box.refine b/test/test_files/abl_w2a_terrain/static_box.refine new file mode 100644 index 0000000000..a08fe9a30a --- /dev/null +++ b/test/test_files/abl_w2a_terrain/static_box.refine @@ -0,0 +1,5 @@ +2 +1 +150 150 0 1000 1000 25 +1 +170 170 0 1000 1000 20 \ No newline at end of file diff --git a/test/test_files/abl_waves_terrain/abl_waves_terrain.inp b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp new file mode 100644 index 0000000000..f6476b7aae --- /dev/null +++ b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp @@ -0,0 +1,78 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 28800.0 # Max (simulated) time to evolve +time.max_step = 10 # Max number of time steps +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.05 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.outputs = velocity_mueff temperature_mueff ow_vof ow_levelset ow_velocity +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0.0 0.0 -9.81 # Gravitational force (3D) +incflo.density = 1.225 +transport.model = ConstTransport +transport.viscosity = 1e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +turbulence.model = Kosovic +Kosovic.refMOL = -1e30 +incflo.physics = ABL OceanWaves TerrainDrag +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing NonLinearSGSTerm DragForcing +BoussinesqBuoyancy.reference_temperature = 263.5 +CoriolisForcing.east_vector = 1.0 0.0 0.0 +CoriolisForcing.north_vector = 0.0 1.0 0.0 +CoriolisForcing.latitude = 90.0 +CoriolisForcing.rotational_time_period = 90405.5439881955 +GeostrophicForcing.geostrophic_wind = 8.0 0.0 0.0 +temperature.source_terms = DragTempForcing +DragTempForcing.reference_temperature = 265.0 + +incflo.velocity = 8.0 0.0 0.0 +ABL.reference_temperature = 263.5 +ABL.temperature_heights = 0.0 100 400.0 +ABL.temperature_values = 265.0 265.0 268.0 +ABL.perturb_temperature = false +ABL.cutoff_height = 50.0 +ABL.perturb_velocity = true +ABL.perturb_ref_height = 50.0 +ABL.Uperiods = 4.0 +ABL.Vperiods = 4.0 +ABL.deltaU = 1e-5 +ABL.deltaV = 1e-5 +ABL.normal_direction = 2 + +OceanWaves.label = Wave1 +OceanWaves.Wave1.type = LinearWaves +OceanWaves.Wave1.wave_height=2.0 +OceanWaves.Wave1.wave_length=40.0 +OceanWaves.Wave1.water_depth=400 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel +amr.max_level = 2 # Max AMR level in hierarchy +tagging.labels = static +tagging.static.type = CartBoxRefinement +tagging.static.static_refinement_def = static_box.txt +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. -20. # Lo corner coordinates +geometry.prob_hi = 400. 400. 380. # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) +# Boundary conditions +zlo.type = "slip_wall" +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.01 +incflo.verbose = 0 diff --git a/test/test_files/abl_waves_terrain/static_box.txt b/test/test_files/abl_waves_terrain/static_box.txt new file mode 100644 index 0000000000..fa81648cb5 --- /dev/null +++ b/test/test_files/abl_waves_terrain/static_box.txt @@ -0,0 +1,5 @@ +2 +1 +0.0 0.0 -20. 400. 400. 20. +1 +0.0 0.0 -8. 400. 400. 8.