diff --git a/docs/sphinx/user/inputs.rst b/docs/sphinx/user/inputs.rst index 46a631924c..8e35ecdabc 100644 --- a/docs/sphinx/user/inputs.rst +++ b/docs/sphinx/user/inputs.rst @@ -97,6 +97,7 @@ documentation provided here might not work with older releases. inputs_turbulence.rst inputs_Momentum_Sources.rst inputs_Temperature_Sources.rst + inputs_DampingLayer.rst inputs_TKE_Sources.rst inputs_ABL.rst inputs_ABL_meso_forcing.rst diff --git a/docs/sphinx/user/inputs_DampingLayer.rst b/docs/sphinx/user/inputs_DampingLayer.rst new file mode 100644 index 0000000000..0f16b6acf5 --- /dev/null +++ b/docs/sphinx/user/inputs_DampingLayer.rst @@ -0,0 +1,185 @@ +.. _inputs_damping_layer: + +Section: DampingLayer +~~~~~~~~~~~~~~~~~~~~~ + +The DampingLayer physics module creates one or more spatially varying damping +coefficient fields near selected domain boundaries. The DampingLayerSource +source term then uses those coefficients to relax solved fields toward user- +defined targets. + +Activate DampingLayer by including it in :input_param:`incflo.physics`. + +Activate DampingLayerSource for each equation you want to damp, for example: + +- ``ICNS.source_terms = DampingLayerSource`` +- ``temperature.source_terms = DampingLayerSource`` +- ``TKE.source_terms = DampingLayerSource`` + +Additional scalar equations that expose source-term controls can also use +``DampingLayerSource``. + +.. input_param:: DampingLayer.fields + + **type:** List of strings, mandatory when ``DampingLayer`` is active + + Field labels that will have damping-layer coefficient fields created. + Typical values include ``velocity``, ``temperature``, ``density``, ``tke``, + and ``sdr``. + +For each entry in :input_param:`DampingLayer.fields`, parameters are configured +per boundary using: + +``DampingLayer...*`` where ``boundary`` is one of +``xlo``, ``xhi``, ``ylo``, ``yhi``, ``zlo``, or ``zhi``. + +.. input_param:: DampingLayer...thickness + + **type:** Real, optional + + Thickness of the damping region for this boundary. A value ``> 0`` enables + damping on that boundary for this field. + +.. input_param:: DampingLayer...blending_fraction + + **type:** Real, optional, default = 0.0 + + Fraction of ``thickness`` used as a ramp from full damping to zero damping. + The remaining fraction applies full damping. + +.. input_param:: DampingLayer...blending_function_type + + **type:** String, optional, default = cosine + + Shape of the ramp region. Supported values are ``linear``, ``quadratic``, + ``exponential``, and ``cosine``. + +.. input_param:: DampingLayer...minimum_height + + **type:** Real, optional (only for ``x*`` and ``y*`` boundaries) + + Adds a vertical limiter so horizontal boundary damping is only applied above + this height. This parameter is not valid on ``zlo`` or ``zhi``. + +.. input_param:: DampingLayer...vertical_blending_thickness + + **type:** Real, mandatory when ``minimum_height`` is provided + + Vertical blending thickness used by the vertical limiter activated through + ``minimum_height``. This parameter is only valid on ``x*`` and ``y*`` + boundaries. + +.. input_param:: DampingLayer...vertical_blending_function_type + + **type:** String, optional, default = cosine + + Shape for blending in the vertical direction. Supported values are ``linear``, + ``quadratic``, ``exponential``, and ``cosine``. This parameter is only + valid on ``x*`` and ``y*`` boundaries. + +DampingLayerSource target parameters are configured with the same namespace: + +``DampingLayer...*`` + +.. input_param:: DampingLayer...target_type + + **type:** String, mandatory when damping is active on a boundary + + Target mode for the relaxed value. Supported values are ``constant``, + ``profile``, ``function``, and ``field``. + +.. input_param:: DampingLayer...target_value + + **type:** List of Real, mandatory for ``target_type = constant`` + + Constant target value per component. The number of values must match the + number of components in the damped field. + +.. input_param:: DampingLayer...target_profile_heights + + **type:** List of Real, mandatory for ``target_type = profile`` + + Monotone height coordinates used for profile interpolation. + +.. input_param:: DampingLayer...target_profile_values + + **type:** List of Real, optional for ``target_type = profile`` + + Convenience form for scalar or x-component profile values. + +.. input_param:: DampingLayer...target_profile_values_x + + **type:** List of Real, optional for ``target_type = profile`` + + x-component profile values. + +.. input_param:: DampingLayer...target_profile_values_y + + **type:** List of Real, optional for ``target_type = profile`` + + y-component profile values. + +.. input_param:: DampingLayer...target_profile_values_z + + **type:** List of Real, optional for ``target_type = profile`` + + z-component profile values. + +.. input_param:: DampingLayer...target_function + + **type:** String expression, mandatory for ``target_type = function`` + + Expression parsed at runtime. Available variables are ``t``, ``x``, ``y``, + ``z``, and ``n`` (component index). + +.. input_param:: DampingLayer...target_field + + **type:** String, mandatory for ``target_type = field`` + + Name of another field used as the target value. This could work with a field + created by a utility or post-processing routine, such as an averaging field, + or a field created by a dedicated physics module. + +.. input_param:: DampingLayer...damped_components + + **type:** List of Integer, optional + + Per-component on/off mask for damping. The list length must match the + number of components in the damped field, and each entry must be either + ``0`` (off) or ``1`` (on). If omitted, all components default to ``1``. + +Example +^^^^^^^ + +.. code-block:: console + + incflo.physics = ABL DampingLayer + DampingLayer.fields = velocity temperature + + ICNS.source_terms = BoussinesqBuoyancy DampingLayerSource + temperature.source_terms = DampingLayerSource + + DampingLayer.velocity.xlo.thickness = 150.0 + DampingLayer.velocity.xlo.blending_fraction = 0.25 + DampingLayer.velocity.xlo.blending_function_type = linear + DampingLayer.velocity.xlo.minimum_height = 100.0 + DampingLayer.velocity.xlo.vertical_blending_thickness = 80.0 + DampingLayer.velocity.xlo.vertical_blending_function_type = cosine + DampingLayer.velocity.xlo.target_type = constant + DampingLayer.velocity.xlo.target_value = 8.0 0.0 0.0 + DampingLayer.velocity.xlo.damped_components = 1 0 0 + + DampingLayer.velocity.ylo.thickness = 120.0 + DampingLayer.velocity.ylo.blending_fraction = 0.5 + DampingLayer.velocity.ylo.blending_function_type = quadratic + DampingLayer.velocity.ylo.target_type = profile + DampingLayer.velocity.ylo.target_profile_heights = 0.0 500.0 1000.0 + DampingLayer.velocity.ylo.target_profile_values_x = 6.0 7.0 8.0 + DampingLayer.velocity.ylo.target_profile_values_y = 0.0 0.5 1.0 + DampingLayer.velocity.ylo.target_profile_values_z = 0.0 0.0 0.0 + + DampingLayer.temperature.zhi.thickness = 200.0 + DampingLayer.temperature.zhi.blending_fraction = 0.4 + DampingLayer.temperature.zhi.blending_function_type = cosine + DampingLayer.temperature.zhi.target_type = function + DampingLayer.temperature.zhi.target_function = "300.0 + 0.01*z" diff --git a/docs/sphinx/user/inputs_Momentum_Sources.rst b/docs/sphinx/user/inputs_Momentum_Sources.rst index 7aa749bb16..05859f9ad3 100644 --- a/docs/sphinx/user/inputs_Momentum_Sources.rst +++ b/docs/sphinx/user/inputs_Momentum_Sources.rst @@ -307,6 +307,16 @@ Section: Momentum Sources when the form drag is known to be under-resolved. +The following arguments are used when ``DampingLayerSource`` is included in +:input_param:`ICNS.source_terms`. + + See :ref:`inputs_damping_layer` for the full list of + ``DampingLayer.velocity..*`` options, including + ``target_type`` modes (``constant``, ``profile``, ``function``, ``field``), + component filtering with ``damped_components``, and blending controls from + the DampingLayer physics module. + + The following arguments are influential when ``GravityForcing`` is included in :input_param:`ICNS.source_terms`. .. input_param:: ICNS.use_perturb_pressure diff --git a/docs/sphinx/user/inputs_Temperature_Sources.rst b/docs/sphinx/user/inputs_Temperature_Sources.rst index 4d900afa71..0682db5072 100644 --- a/docs/sphinx/user/inputs_Temperature_Sources.rst +++ b/docs/sphinx/user/inputs_Temperature_Sources.rst @@ -30,6 +30,12 @@ Section: Temperature Sources the time step size. +When ``DampingLayerSource`` is included in +:input_param:`temperature.source_terms`, configure it using +``DampingLayer.temperature..*`` parameters documented in +:ref:`inputs_damping_layer`. + + The following list of inputs are used with the `Temperature.source_terms = PerturbationForcing` option to add perturbation to the temperature field to generate flow structures for LES when the inflow data is coarse or uniform flow condition. Not recommended for use with RANS models. diff --git a/docs/sphinx/user/inputs_incflo.rst b/docs/sphinx/user/inputs_incflo.rst index e861006838..cc80451f66 100644 --- a/docs/sphinx/user/inputs_incflo.rst +++ b/docs/sphinx/user/inputs_incflo.rst @@ -12,7 +12,7 @@ as initial conditions and discretization options. Specify a string or a list of strings for each type of physics to initialize and simulate. Physics is additive and more than one type of physics may be used. - Current implemented physics include FreeStream, SyntheticTurbulence, ABL, Actuator, RayleighTaylor, BoussinesqBubble, TaylorGreenVortex, and ScalarAdvection (which is an example of using a passive scalar advection) among others. + Current implemented physics include FreeStream, SyntheticTurbulence, ABL, DampingLayer, Actuator, RayleighTaylor, BoussinesqBubble, TaylorGreenVortex, and ScalarAdvection (which is an example of using a passive scalar advection) among others. For multiphase simulations, the MultiPhase physics must be specified, and for forcing wave profiles into the domain, the OceanWaves physics must be specified as well. For immersed boundary forcing method, :ref:`TerrainDrag ` or :ref:`ChannelBuilder ` can be used. For representing forested regions :ref:`ForestDrag ` must be specified, which has two file-based approaches to parameterize forest drag models. diff --git a/src/equation_systems/CMakeLists.txt b/src/equation_systems/CMakeLists.txt index d7d93f153f..6947bce493 100644 --- a/src/equation_systems/CMakeLists.txt +++ b/src/equation_systems/CMakeLists.txt @@ -11,3 +11,4 @@ add_subdirectory(sdr) add_subdirectory(levelset) add_subdirectory(vof) add_subdirectory(passive_scalar) +add_subdirectory(source_terms) \ No newline at end of file diff --git a/src/equation_systems/source_terms/CMakeLists.txt b/src/equation_systems/source_terms/CMakeLists.txt new file mode 100644 index 0000000000..6924976fe6 --- /dev/null +++ b/src/equation_systems/source_terms/CMakeLists.txt @@ -0,0 +1,3 @@ +target_sources(${kynema_sgf_lib_name} PRIVATE + DampingLayerSource.cpp + ) \ No newline at end of file diff --git a/src/equation_systems/source_terms/DampingLayerSource.H b/src/equation_systems/source_terms/DampingLayerSource.H new file mode 100644 index 0000000000..536fdca1b4 --- /dev/null +++ b/src/equation_systems/source_terms/DampingLayerSource.H @@ -0,0 +1,441 @@ +#ifndef DAMPINGLAYERSOURCE_H +#define DAMPINGLAYERSOURCE_H + +#include + +#include "src/CFDSim.H" +#include "src/core/FieldRepo.H" +#include "src/utilities/linear_interpolation.H" +#include "src/equation_systems/density/DensitySource.H" +#include "src/equation_systems/sdr/SDRSource.H" +#include "src/equation_systems/SourceTerm.H" +#include "src/equation_systems/icns/MomentumSource.H" +#include "src/equation_systems/temperature/TemperatureSource.H" +#include "src/equation_systems/tke/TKESource.H" +#include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" + +using namespace amrex::literals; + +enum class TargetType : std::uint8_t { + Constant = 0, + Profile = 1, + Function = 2, + Field = 3 +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE TargetType +string_to_target_type(const std::string& target_type_str) +{ + std::string target_type_lower = amrex::toLower(target_type_str); + if (target_type_lower == "constant") { + return TargetType::Constant; + } + if (target_type_lower == "profile") { + return TargetType::Profile; + } + if (target_type_lower == "function") { + return TargetType::Function; + } + if (target_type_lower == "field") { + return TargetType::Field; + } + amrex::Abort( + "Invalid target type: " + target_type_str + + ". Must be one of: constant, profile, function, field."); + // Unreachable, but avoids compiler warning + return TargetType::Constant; +} + +namespace kynema_sgf::pde { + +template +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = false; + static constexpr const char* field_name = ""; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = "velocity"; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = "temperature"; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = "density"; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = "tke"; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = "sdr"; +}; + +template <> +struct DampingLayerSourceTraits +{ + static constexpr bool is_specialized = true; + static constexpr const char* field_name = ""; +}; + +template +class DampingLayerSource : public T::template Register> +{ +public: + static_assert( + DampingLayerSourceTraits::is_specialized, + "DampingLayerSourceTraits specialization missing for this source " + "type. Add a specialization with field_name."); + + static std::string identifier() { return "DampingLayerSource"; } + + static std::string alt_identifier() { return "DampingLayer"; } + + static constexpr const char* damped_field_name() + { + return DampingLayerSourceTraits::field_name; + } + + explicit DampingLayerSource(const CFDSim& sim) + : m_mesh(sim.mesh()) + , m_repo(sim.repo()) + , m_time(sim.time()) + , m_damped_field_name(damped_field_name()) + , m_damped_field(sim.repo().get_field(damped_field_name())) + { + + // Insert gets and queries about target field profile + int bc_idx = 0; + for (const auto& name : m_bc_names) { + // Get arguments specific to this boundary + amrex::ParmParse pp( + alt_identifier() + "." + m_damped_field_name + "." + name); + // If thickness > 0, then this boundary has a damping layer + amrex::Real thickness = -1.0_rt; + pp.query("thickness", thickness); + + if (thickness <= 0.0_rt) { + ++bc_idx; + continue; + } + + std::string target_type_string; + std::string target_function_string{"0.0"}; + pp.get("target_type", target_type_string); + auto target_type = string_to_target_type(target_type_string); + m_target_type[bc_idx] = target_type; + if (target_type == TargetType::Constant) { + pp.getarr("target_value", m_target_value[bc_idx]); + } else if (target_type == TargetType::Profile) { + if (pp.contains("target_profile_values")) { + pp.getarr( + "target_profile_values", + m_target_profile_values_x[bc_idx]); + } else { + pp.getarr( + "target_profile_values_x", + m_target_profile_values_x[bc_idx]); + pp.queryarr( + "target_profile_values_y", + m_target_profile_values_y[bc_idx]); + pp.queryarr( + "target_profile_values_z", + m_target_profile_values_z[bc_idx]); + } + pp.getarr( + "target_profile_heights", m_target_profile_heights[bc_idx]); + } else if (target_type == TargetType::Function) { + pp.get("target_function", target_function_string); + } else if (target_type == TargetType::Field) { + pp.get("target_field", m_target_field_name[bc_idx]); + } else { + amrex::Abort( + "DampingLayerSource: unknown target_type: " + + target_type_string); + } + + pp.queryarr("damped_components", m_damped_components[bc_idx]); + if (m_damped_components[bc_idx].empty()) { + // If no components specified, damp all components + m_damped_components[bc_idx].resize(m_damped_field.num_comp()); + for (int comp = 0; comp < m_damped_field.num_comp(); ++comp) { + m_damped_components[bc_idx][comp] = 1; + } + } + // Check values of damped components + if (m_damped_components[bc_idx].size() != + m_damped_field.num_comp()) { + amrex::Abort( + "DampingLayerSource: damped_components for field " + + m_damped_field_name + + " must have the same number of components as the field: " + + std::to_string(m_damped_field.num_comp())); + } + for (const auto comp : m_damped_components[bc_idx]) { + if (comp < 0 || comp > 1) { + amrex::Abort( + "DampingLayerSource: damped_components for field " + + m_damped_field_name + + " must be set as 0 (off) or 1 (on)"); + } + } + + if (target_type == TargetType::Constant) { + if (m_target_value[bc_idx].size() != + m_damped_field.num_comp()) { + amrex::Abort( + "DampingLayerSource: target_value must have the same " + "size as the number of components in the damped " + "field."); + } + const auto num_comp_val = m_target_value[bc_idx].size(); + m_dev_target_value[bc_idx].resize(num_comp_val); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_target_value[bc_idx].begin(), + m_target_value[bc_idx].end(), + m_dev_target_value[bc_idx].begin()); + } + if (target_type == TargetType::Profile) { + if (m_target_profile_values_x[bc_idx].size() != + m_target_profile_heights[bc_idx].size()) { + amrex::Abort( + "DampingLayerSource: target_profile_values and " + "target_profile_heights must have the same size."); + } + const auto num_prof_val = + m_target_profile_values_x[bc_idx].size(); + m_dev_target_profile_values_x[bc_idx].resize(num_prof_val); + m_dev_target_profile_heights[bc_idx].resize(num_prof_val); + + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, + m_target_profile_heights[bc_idx].begin(), + m_target_profile_heights[bc_idx].end(), + m_dev_target_profile_heights[bc_idx].begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, + m_target_profile_values_x[bc_idx].begin(), + m_target_profile_values_x[bc_idx].end(), + m_dev_target_profile_values_x[bc_idx].begin()); + const auto num_prof_val_y = + m_target_profile_values_y[bc_idx].size(); + if (num_prof_val_y > 0) { + if (num_prof_val_y != num_prof_val) { + amrex::Abort( + "DampingLayerSource: target_profile_values_y must " + "have the same size as target_profile_values_x."); + } + m_dev_target_profile_values_y[bc_idx].resize( + num_prof_val_y); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, + m_target_profile_values_y[bc_idx].begin(), + m_target_profile_values_y[bc_idx].end(), + m_dev_target_profile_values_y[bc_idx].begin()); + } + const auto num_prof_val_z = + m_target_profile_values_z[bc_idx].size(); + if (num_prof_val_z > 0) { + if (num_prof_val_z != num_prof_val) { + amrex::Abort( + "DampingLayerSource: target_profile_values_z must " + "have the same size as target_profile_values_x."); + } + m_dev_target_profile_values_z[bc_idx].resize( + num_prof_val_z); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, + m_target_profile_values_z[bc_idx].begin(), + m_target_profile_values_z[bc_idx].end(), + m_dev_target_profile_values_z[bc_idx].begin()); + } + } + if (target_type == TargetType::Function) { + m_function_parser[bc_idx].define(target_function_string); + m_function_parser[bc_idx].registerVariables( + {"t", "x", "y", "z", "n"}); + } + ++bc_idx; + } + } + + ~DampingLayerSource() override = default; + + void operator()( + int lev, FieldState fstate, amrex::MultiFab& src_term) const override + { + Field* damping_layer_ptr{nullptr}; + + auto const& src_arrs = src_term.arrays(); + auto const& damped_field_arrs = + m_damped_field.state(field_impl::dof_state(fstate))(lev) + .const_arrays(); + + const auto& geom = m_mesh.Geom(lev); + const auto& problo = geom.ProbLoArray(); + const auto& dx = geom.CellSizeArray(); + const auto time = m_time.current_time(); + + for (int bc_idx = 0; bc_idx < 6; ++bc_idx) { + const auto damping_coefficient_field_name = + "damping_layer_" + m_damped_field_name + "_" + + m_bc_names[bc_idx]; + // If damping layer exists, coefficient field exists + if (m_repo.field_exists(damping_coefficient_field_name)) { + // Form damping layer field name + damping_layer_ptr = + &m_repo.get_field(damping_coefficient_field_name); + auto& damping_mfab = (*damping_layer_ptr)(lev); + auto const& damping_arrs = damping_mfab.const_arrays(); + + const auto target_type = m_target_type[bc_idx]; + + // Get constant target value if being used + const amrex::Real* const_target_value = + (target_type == TargetType::Constant) + ? m_dev_target_value[bc_idx].data() + : nullptr; + + // Get target profile values, positions if being used + const bool use_profile = (target_type == TargetType::Profile); + const amrex::Real* prof_z = + use_profile ? m_dev_target_profile_heights[bc_idx].data() + : nullptr; + const amrex::Real* prof_z_end = + use_profile ? m_dev_target_profile_heights[bc_idx].end() + : nullptr; + const amrex::Real* prof_v_x = + use_profile ? m_dev_target_profile_values_x[bc_idx].data() + : nullptr; + const amrex::Real* prof_v_y = + use_profile && !m_target_profile_values_y[bc_idx].empty() + ? m_dev_target_profile_values_y[bc_idx].data() + : nullptr; + const amrex::Real* prof_v_z = + use_profile && !m_target_profile_values_z[bc_idx].empty() + ? m_dev_target_profile_values_z[bc_idx].data() + : nullptr; + + // Set up the user-defined function, whether used or not + auto udf_func = m_function_parser[bc_idx] + .template compile(); + + // Get target field arrays if being used + Field* target_field_ptr{nullptr}; + if (target_type == TargetType::Field) { + target_field_ptr = + &m_repo.get_field(m_target_field_name[bc_idx]); + } + auto const& target_field_arrs = + (target_type == TargetType::Field) + ? (*target_field_ptr)(lev).const_arrays() + : amrex::MultiArray4(); + + amrex::GpuArray damp_this_component{ + 0, 0, 0}; + for (int comp = 0; comp < m_damped_field.num_comp(); ++comp) { + damp_this_component[comp] = + m_damped_components[bc_idx][comp]; + } + + amrex::ParallelFor( + damping_mfab, amrex::IntVect(0), m_damped_field.num_comp(), + [=] AMREX_GPU_DEVICE( + int nbx, int i, int j, int k, int n) noexcept { + const amrex::Real x = + problo[0] + ((i + 0.5_rt) * dx[0]); + const amrex::Real y = + problo[1] + ((j + 0.5_rt) * dx[1]); + const amrex::Real z = + problo[2] + ((k + 0.5_rt) * dx[2]); + + // Modify src term based on current value of damped + // field, damping coefficient field, target value + auto target_value = 0.0_rt; + if (target_type == TargetType::Constant) { + target_value = const_target_value[n]; + } else if (target_type == TargetType::Profile) { + const auto idx = + interp::bisection_search(prof_z, prof_z_end, z); + + // Interpolate the target profile values based on + // the cell center height + const amrex::Real* prof = (n == 0) ? prof_v_x + : (n == 1) ? prof_v_y + : prof_v_z; + target_value = + interp::linear_impl(prof_z, prof, z, idx); + } else if (target_type == TargetType::Function) { + // Evaluate the user-defined function + target_value = udf_func( + time, x, y, z, static_cast(n)); + } else if (target_type == TargetType::Field) { + target_value = target_field_arrs[nbx](i, j, k, n); + } + + if (damp_this_component[n] == 1) { + src_arrs[nbx](i, j, k, n) += + damping_arrs[nbx](i, j, k) * + (target_value - + damped_field_arrs[nbx](i, j, k, n)); + } + }); + } + } + } + +private: + const amrex::AmrCore& m_mesh; + const FieldRepo& m_repo; + const SimTime& m_time; + std::string m_damped_field_name; + Field& m_damped_field; + const amrex::Vector m_bc_names = {"xlo", "xhi", "ylo", + "yhi", "zlo", "zhi"}; + + // Can have different targets for each boundary, so use arrays of size 6 + amrex::Array m_target_type; + amrex::Array, 6> m_target_value; + amrex::Array, 6> m_target_profile_values_x; + amrex::Array, 6> m_target_profile_values_y; + amrex::Array, 6> m_target_profile_values_z; + amrex::Array, 6> m_target_profile_heights; + amrex::Array m_function_parser; + amrex::Array m_target_field_name; + amrex::Array, 6> m_damped_components; + + // Device copies of the above vectors + amrex::Array, 6> m_dev_target_value; + amrex::Array, 6> + m_dev_target_profile_values_x; + amrex::Array, 6> + m_dev_target_profile_values_y; + amrex::Array, 6> + m_dev_target_profile_values_z; + amrex::Array, 6> + m_dev_target_profile_heights; +}; + +} // namespace kynema_sgf::pde + +#endif /* DAMPINGLAYERSOURCE_H */ diff --git a/src/equation_systems/source_terms/DampingLayerSource.cpp b/src/equation_systems/source_terms/DampingLayerSource.cpp new file mode 100644 index 0000000000..1f007c1a03 --- /dev/null +++ b/src/equation_systems/source_terms/DampingLayerSource.cpp @@ -0,0 +1,12 @@ +#include "src/equation_systems/source_terms/DampingLayerSource.H" + +namespace kynema_sgf::pde { + +template class DampingLayerSource; +template class DampingLayerSource; +template class DampingLayerSource; +template class DampingLayerSource; +template class DampingLayerSource; +template class DampingLayerSource; + +} // namespace kynema_sgf::pde \ No newline at end of file diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt index d1d4828f73..f7012dc857 100644 --- a/src/physics/CMakeLists.txt +++ b/src/physics/CMakeLists.txt @@ -20,6 +20,7 @@ target_sources(${kynema_sgf_lib_name} ActuatorSourceTagging.cpp Intermittency.cpp ForestDrag.cpp + DampingLayer.cpp ) add_subdirectory(multiphase) diff --git a/src/physics/DampingLayer.H b/src/physics/DampingLayer.H new file mode 100644 index 0000000000..664df51e78 --- /dev/null +++ b/src/physics/DampingLayer.H @@ -0,0 +1,124 @@ +#ifndef DAMPINGLAYER_H +#define DAMPINGLAYER_H + +#include +#include "src/core/Physics.H" +#include "src/core/Field.H" +#include "src/CFDSim.H" +#include "AMReX_Gpu.H" + +namespace kynema_sgf::damping_layer { + +enum class BlendingFunctionType : std::uint8_t { + Linear = 0, + Quadratic = 1, + Exponential = 2, + Cosine = 3 +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BlendingFunctionType +string_to_blending_function_type(const std::string& op_str) +{ + std::string op = amrex::toLower(op_str); + if (op == "linear") { + return BlendingFunctionType::Linear; + } + if (op == "quadratic") { + return BlendingFunctionType::Quadratic; + } + if (op == "exponential") { + return BlendingFunctionType::Exponential; + } + if (op == "cosine") { + return BlendingFunctionType::Cosine; + } + amrex::Abort( + "Invalid blending function type: " + op_str + + ". Must be one of: linear, quadratic, exponential, cosine."); + // Unreachable, but avoids compiler warning + return BlendingFunctionType::Cosine; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real blending_function( + amrex::Real position_nondim, BlendingFunctionType function_type) +{ + switch (function_type) { + case BlendingFunctionType::Linear: + return 1.0_rt - position_nondim; + case BlendingFunctionType::Quadratic: + return 1.0_rt - position_nondim * position_nondim; + case BlendingFunctionType::Exponential: // exponential difference and ratio + return 1.0_rt - (std::expm1(std::pow(position_nondim, 3.5_rt)) / + std::expm1(1.0_rt)); + case BlendingFunctionType::Cosine: + default: + return 0.5_rt + + 0.5_rt * + std::cos(std::numbers::pi_v * position_nondim); + } +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real damping_calc( + amrex::Real distance_from_boundary, + amrex::Real layer_thickness, + amrex::Real blending_fraction, + BlendingFunctionType function_type) +{ + const amrex::Real length_full_damp = + layer_thickness * (1.0_rt - blending_fraction); + const amrex::Real blend_position_nondim = + (distance_from_boundary - length_full_damp) / + (layer_thickness * blending_fraction); + if (blend_position_nondim <= 0.0_rt) { + return 1.0_rt; + } + if (blend_position_nondim <= 1.0_rt) { + return blending_function(blend_position_nondim, function_type); + } + return 0.0_rt; +} + +class DampingLayer : public Physics::Register +{ +public: + static std::string identifier() { return "DampingLayer"; } + + explicit DampingLayer(CFDSim& sim); + + ~DampingLayer() override = default; + + void initialize_fields(int level, const amrex::Geometry& geom) override; + + void pre_init_actions() override {} + + void post_init_actions() override {} + + void post_regrid_actions() override; + + void pre_advance_work() override {} + + void post_advance_work() override {} + +private: + FieldRepo& m_repo; + const amrex::AmrCore& m_mesh; + + const amrex::Vector m_bc_names = {"xlo", "xhi", "ylo", + "yhi", "zlo", "zhi"}; + + //! Quantities for constructing damping layers + amrex::Vector m_field_names; + amrex::Gpu::DeviceVector> m_layers_thickness; + amrex::Gpu::DeviceVector> + m_layers_blending_fraction; + amrex::Gpu::DeviceVector> m_layers_min_height; + amrex::Gpu::DeviceVector> + m_layers_vertical_blending_thickness; + amrex::Gpu::DeviceVector> + m_layers_vertical_blending_function_type; + amrex::Gpu::DeviceVector> + m_layers_blending_function_type; +}; +} // namespace kynema_sgf::damping_layer + +#endif diff --git a/src/physics/DampingLayer.cpp b/src/physics/DampingLayer.cpp new file mode 100644 index 0000000000..4c4fc8d2ff --- /dev/null +++ b/src/physics/DampingLayer.cpp @@ -0,0 +1,275 @@ +#include "src/physics/DampingLayer.H" +#include "src/utilities/math_ops.H" +#include "src/utilities/constants.H" +#include "src/utilities/IOManager.H" +#include "src/CFDSim.H" +#include "AMReX_ParmParse.H" +#include "AMReX_iMultiFab.H" +#include "AMReX_Gpu.H" +#include "AMReX_REAL.H" + +using namespace amrex::literals; + +namespace kynema_sgf::damping_layer { + +DampingLayer::DampingLayer(CFDSim& sim) : m_repo(sim.repo()), m_mesh(sim.mesh()) +{ + amrex::ParmParse pp(identifier()); + pp.getarr("fields", m_field_names); + + amrex::Vector> layers_thickness; + amrex::Vector> layers_blending_fraction; + amrex::Vector> layers_min_height; + amrex::Vector> + layers_vertical_blending_thickness; + amrex::Vector> + layers_vertical_blending_function_type; + amrex::Vector> + layers_blending_function_type; + + for (const auto& lbl : m_field_names) { + const std::string key = identifier() + "." + lbl; + + amrex::Array bc_thickness; + amrex::Array bc_blending_fraction; + amrex::Array bc_min_height; + amrex::Array bc_vertical_blending_thickness; + amrex::Array + bc_vertical_blending_function_type; + amrex::Array bc_blending_function_type; + int bc_index = 0; + for (const auto& name : m_bc_names) { + // Get arguments specific to this boundary + const std::string key_bc = key + "." + name; + amrex::ParmParse pp_bc(key_bc); + amrex::Real thickness = -1.0_rt; + pp_bc.query("thickness", thickness); + amrex::Real blending_fraction = 0.0_rt; + pp_bc.query("blending_fraction", blending_fraction); + amrex::Real min_height = constants::LOW_NUM; + pp_bc.query("minimum_height", min_height); + std::string blending_function_str = "cosine"; + pp_bc.query("blending_function_type", blending_function_str); + BlendingFunctionType blending_function_type = + string_to_blending_function_type(blending_function_str); + amrex::Real vert_blend_thickness = -1.0_rt; + if (pp_bc.contains("minimum_height")) { + pp_bc.get("vertical_blending_thickness", vert_blend_thickness); + } + std::string vert_blend_function_str = "cosine"; + pp_bc.query( + "vertical_blending_function_type", vert_blend_function_str); + BlendingFunctionType vert_blend_function_type = + string_to_blending_function_type(vert_blend_function_str); + // Abort statement if height is specified for a z boundary + if (name == "zlo" || name == "zhi") { + if (pp_bc.contains("minimum_height")) { + amrex::Abort( + "DampingLayer: minimum_height is not supported for z " + "boundaries."); + } + if (pp_bc.contains("vertical_blending_thickness")) { + amrex::Abort( + "DampingLayer: vertical_blending_thickness is not " + "supported for z boundaries."); + } + if (pp_bc.contains("vertical_blending_function_type")) { + amrex::Abort( + "DampingLayer: vertical_blending_function_type is not " + "supported for z boundaries."); + } + } + // Create field to go with this boundary damping layer + bc_thickness[bc_index] = thickness; + if (thickness > 0.0_rt) { + const auto field_name = "damping_layer_" + lbl + "_" + name; + m_repo.declare_field(field_name, 1, 1, 1); + sim.io_manager().register_io_var(field_name); + // Record damping layer parameters for field creation + bc_blending_fraction[bc_index] = blending_fraction; + bc_blending_function_type[bc_index] = blending_function_type; + if (name != "zlo" && name != "zhi") { + bc_min_height[bc_index] = min_height; + bc_vertical_blending_thickness[bc_index] = + vert_blend_thickness; + bc_vertical_blending_function_type[bc_index] = + vert_blend_function_type; + } + } + ++bc_index; + } + + layers_thickness.emplace_back(bc_thickness); + layers_blending_fraction.emplace_back(bc_blending_fraction); + layers_min_height.emplace_back(bc_min_height); + layers_vertical_blending_thickness.emplace_back( + bc_vertical_blending_thickness); + layers_vertical_blending_function_type.emplace_back( + bc_vertical_blending_function_type); + layers_blending_function_type.emplace_back(bc_blending_function_type); + } + + const int nfields = static_cast(layers_thickness.size()); + m_layers_thickness.resize(nfields); + m_layers_blending_fraction.resize(nfields); + m_layers_min_height.resize(nfields); + m_layers_blending_function_type.resize(nfields); + m_layers_vertical_blending_thickness.resize(nfields); + m_layers_vertical_blending_function_type.resize(nfields); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, layers_thickness.begin(), + layers_thickness.end(), m_layers_thickness.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, layers_blending_fraction.begin(), + layers_blending_fraction.end(), m_layers_blending_fraction.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, layers_min_height.begin(), + layers_min_height.end(), m_layers_min_height.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, layers_vertical_blending_thickness.begin(), + layers_vertical_blending_thickness.end(), + m_layers_vertical_blending_thickness.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, + layers_vertical_blending_function_type.begin(), + layers_vertical_blending_function_type.end(), + m_layers_vertical_blending_function_type.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, layers_blending_function_type.begin(), + layers_blending_function_type.end(), + m_layers_blending_function_type.begin()); +} + +void DampingLayer::initialize_fields(int level, const amrex::Geometry& geom) +{ + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + + const int nfields = static_cast(m_field_names.size()); + const std::string* field_names_ptr = m_field_names.data(); + const amrex::Array* bc_thickness_ptr = + m_layers_thickness.data(); + const amrex::Array* bc_blending_fraction_ptr = + m_layers_blending_fraction.data(); + const amrex::Array* bc_min_height_ptr = + m_layers_min_height.data(); + const amrex::Array* bc_blending_function_type_ptr = + m_layers_blending_function_type.data(); + + for (int field_idx = 0; field_idx < nfields; ++field_idx) { + const std::string lbl = field_names_ptr[field_idx]; + + const amrex::Array& bc_thickness = + bc_thickness_ptr[field_idx]; + const amrex::Array& bc_blending_fraction = + bc_blending_fraction_ptr[field_idx]; + const amrex::Array& bc_min_height = + bc_min_height_ptr[field_idx]; + const amrex::Array& bc_vertical_blending_thickness = + m_layers_vertical_blending_thickness[field_idx]; + const amrex::Array& + bc_vertical_blending_function_type = + m_layers_vertical_blending_function_type[field_idx]; + const amrex::Array& bc_blending_function_type = + bc_blending_function_type_ptr[field_idx]; + + // Field pointer to the damping layer field for this boundary condition + Field* damping_layer_ptr{nullptr}; + + for (int bc_idx = 0; bc_idx < 6; ++bc_idx) { + if (bc_thickness[bc_idx] > 0.0_rt) { + // Form damping layer field name + const auto field_name = + "damping_layer_" + lbl + "_" + m_bc_names[bc_idx]; + damping_layer_ptr = &m_repo.get_field(field_name); + + // Get other damping layer parameters for this boundary + const amrex::Real blending_fraction = + bc_blending_fraction[bc_idx]; + const amrex::Real min_height = + (bc_idx < 4) ? bc_min_height[bc_idx] : constants::LOW_NUM; + const amrex::Real vertical_blending_thickness = + (bc_idx < 4) ? bc_vertical_blending_thickness[bc_idx] + : -1.0_rt; + const BlendingFunctionType vertical_blending_function_type = + (bc_idx < 4) ? bc_vertical_blending_function_type[bc_idx] + : BlendingFunctionType::Cosine; + const BlendingFunctionType blending_function_type = + bc_blending_function_type[bc_idx]; + + auto& damping_layer_mfab = (*damping_layer_ptr)(level); + auto damping_layer_arrs = damping_layer_mfab.arrays(); + + amrex::ParallelFor( + damping_layer_mfab, + [=] AMREX_GPU_DEVICE( + int nbx, int i, int j, int k) noexcept { + const amrex::Real x = prob_lo[0] + (i + 0.5_rt) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5_rt) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5_rt) * dx[2]; + + // Determine the distance from the boundary + amrex::Real distance_from_boundary = 0.0_rt; + switch (bc_idx) { + case 0: // xlo + distance_from_boundary = x - prob_lo[0]; + break; + case 1: // xhi + distance_from_boundary = prob_hi[0] - x; + break; + case 2: // ylo + distance_from_boundary = y - prob_lo[1]; + break; + case 3: // yhi + distance_from_boundary = prob_hi[1] - y; + break; + case 4: // zlo + distance_from_boundary = z - prob_lo[2]; + break; + case 5: // zhi + distance_from_boundary = prob_hi[2] - z; + break; + default: + amrex::Abort( + "Invalid boundary index in " + "DampingLayer::initialize_fields"); + } + + amrex::Real damping_coeff = damping_calc( + distance_from_boundary, bc_thickness[bc_idx], + blending_fraction, blending_function_type); + + if (min_height > constants::LOW_NUM) { + const amrex::Real distance_from_zhi = + prob_hi[2] - z; + const amrex::Real vertical_thickness = + prob_hi[2] - min_height; + const amrex::Real vertical_blending_fraction = + vertical_blending_thickness / distance_from_zhi; + const amrex::Real vertical_damping_coeff = + damping_calc( + distance_from_zhi, vertical_thickness, + vertical_blending_fraction, + vertical_blending_function_type); + damping_coeff = + std::min(damping_coeff, vertical_damping_coeff); + } + // Set the damping coefficient in the field + damping_layer_arrs[nbx](i, j, k, 0) = damping_coeff; + }); + amrex::Gpu::streamSynchronize(); + } + } + } +} + +void DampingLayer::post_regrid_actions() +{ + const int nlevels = m_repo.num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + initialize_fields(lev, m_mesh.Geom(lev)); + } +} + +} // namespace kynema_sgf::damping_layer diff --git a/src/physics/TerrainDrag.H b/src/physics/TerrainDrag.H index af7d566594..1b776178d9 100644 --- a/src/physics/TerrainDrag.H +++ b/src/physics/TerrainDrag.H @@ -66,7 +66,7 @@ private: amrex::Real m_uniform_z0{0.1_rt}; Field& m_terrainz0; Field& m_terrain_height; - Field& m_terrain_damping; + Field* m_terrain_damping{nullptr}; //! Lateral Damping - X amrex::Real m_damp_east_slope{0.0_rt}; diff --git a/src/physics/TerrainDrag.cpp b/src/physics/TerrainDrag.cpp index a07481ee23..b20ddc03b9 100644 --- a/src/physics/TerrainDrag.cpp +++ b/src/physics/TerrainDrag.cpp @@ -9,6 +9,7 @@ #include "src/utilities/io_utils.H" #include "src/utilities/linear_interpolation.H" #include "AMReX_REAL.H" +#include "src/physics/TerrainDragDamping.H" using namespace amrex::literals; @@ -24,7 +25,6 @@ TerrainDrag::TerrainDrag(CFDSim& sim) , m_terrain_drag(sim.repo().declare_int_field("terrain_drag", 1, 1, 1)) , m_terrainz0(sim.repo().declare_field("terrainz0", 1, 1, 1)) , m_terrain_height(sim.repo().declare_field("terrain_height", 1, 1, 1)) - , m_terrain_damping(sim.repo().declare_field("terrain_damping", 1, 1, 1)) { m_terrain_is_waves = sim.physics_manager().contains("OceanWaves") && @@ -51,28 +51,26 @@ TerrainDrag::TerrainDrag(CFDSim& sim) m_sim.io_manager().register_output_int_var("terrain_blank"); m_sim.io_manager().register_io_var("terrainz0"); m_sim.io_manager().register_io_var("terrain_height"); - m_sim.io_manager().register_io_var("terrain_damping"); m_terrain_blank.setVal(0); m_terrain_drag.setVal(0); - m_terrain_damping.setVal(0); m_terrainz0.set_default_fillpatch_bc(m_sim.time()); m_terrain_height.set_default_fillpatch_bc(m_sim.time()); - m_terrain_damping.set_default_fillpatch_bc(m_sim.time()); + amrex::ParmParse pp(identifier()); - pp.query("damp_east_slope", m_damp_east_slope); - pp.query("damp_east_full", m_damp_east_full); - pp.query("damp_west_slope", m_damp_west_slope); - pp.query("damp_west_full", m_damp_west_full); - pp.query("damp_north_slope", m_damp_north_slope); - pp.query("damp_north_full", m_damp_north_full); - pp.query("damp_south_slope", m_damp_south_slope); - pp.query("damp_south_full", m_damp_south_full); - pp.query("horizontal_time_scale", m_horizontal_tau); - pp.query("horizontal_abl_height", m_horizontal_abl_height); - pp.query("horizontal_slope_end", m_horizontal_slope_end); - pp.query("vertical_slope", m_vertical_slope); - pp.query("vertical_full", m_vertical_full); + if (damping_legacy::has_damping_inputs(pp)) { + m_terrain_damping = + &sim.repo().declare_field("terrain_damping", 1, 1, 1); + m_sim.io_manager().register_io_var("terrain_damping"); + m_terrain_damping->setVal(0); + m_terrain_damping->set_default_fillpatch_bc(m_sim.time()); + damping_legacy::query_damping_inputs( + pp, m_damp_east_slope, m_damp_east_full, m_damp_west_slope, + m_damp_west_full, m_damp_north_slope, m_damp_north_full, + m_damp_south_slope, m_damp_south_full, m_horizontal_tau, + m_horizontal_abl_height, m_horizontal_slope_end, m_vertical_slope, + m_vertical_full); + } } void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) @@ -101,12 +99,10 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); - const auto& prob_hi = geom.ProbHiArray(); auto& blanking = m_terrain_blank(level); auto& terrainz0 = m_terrainz0(level); auto& terrain_height = m_terrain_height(level); auto& drag = m_terrain_drag(level); - auto& damping = m_terrain_damping(level); const auto xterrain_size = xterrain.size(); const auto yterrain_size = yterrain.size(); const auto zterrain_size = zterrain.size(); @@ -148,7 +144,6 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto levelDrag = drag.arrays(); auto levelz0 = terrainz0.arrays(); auto levelheight = terrain_height.arrays(); - auto levelDamping = damping.arrays(); const auto uniform_z0 = m_uniform_z0; amrex::ParallelFor( @@ -188,113 +183,14 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) }); amrex::Gpu::streamSynchronize(); - // Lateral East - const amrex::Real horizontal_tau = m_horizontal_tau; - const amrex::Real horizontal_abl_height = m_horizontal_abl_height; - const amrex::Real z_sloped = m_horizontal_slope_end; - const amrex::Real vertical_slope = m_vertical_slope; - const amrex::Real vertical_full = m_vertical_full; - const amrex::Real damping_east_start = - prob_hi[0] - (m_damp_east_full + m_damp_east_slope); - const amrex::Real damping_east_end = prob_hi[0] - m_damp_east_full; - // West - const amrex::Real damping_west_start = - prob_lo[0] + (m_damp_west_full + m_damp_west_slope); - const amrex::Real damping_west_end = prob_lo[0] + m_damp_west_full; - // North - const amrex::Real damping_north_start = - prob_hi[1] - (m_damp_north_full + m_damp_north_slope); - const amrex::Real damping_north_end = prob_hi[1] - m_damp_north_full; - // South - const amrex::Real damping_south_start = - prob_lo[1] + (m_damp_south_full + m_damp_south_slope); - const amrex::Real damping_south_end = prob_lo[1] + m_damp_south_full; - amrex::ParallelFor( - damping, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - amrex::Real horizontal_coeff_east = 0.0_rt; - amrex::Real horizontal_coeff_north = 0.0_rt; - amrex::Real horizontal_coeff_west = 0.0_rt; - amrex::Real horizontal_coeff_south = 0.0_rt; - amrex::Real vertical_coeff = 0.0_rt; - const amrex::Real x = prob_lo[0] + (i + 0.5_rt) * dx[0]; - const amrex::Real y = prob_lo[1] + (j + 0.5_rt) * dx[1]; - const amrex::Real z = prob_lo[2] + (k + 0.5_rt) * dx[2]; - if (x < damping_east_start) { - horizontal_coeff_east = 0.0_rt; - } else if (x >= damping_east_end) { - horizontal_coeff_east = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (x - damping_east_start) / - (damping_east_end - damping_east_start)); - horizontal_coeff_east = term * term; - } - if (x > damping_west_start) { - horizontal_coeff_west = 0.0_rt; - } else if (x <= damping_west_end) { - horizontal_coeff_west = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (x - damping_west_start) / - (damping_west_end - damping_west_start)); - horizontal_coeff_west = term * term; - } - if (y < damping_north_start) { - horizontal_coeff_north = 0.0_rt; - } else if (y >= damping_north_end) { - horizontal_coeff_north = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (y - damping_north_start) / - (damping_north_end - damping_north_start)); - horizontal_coeff_north = term * term; - } - if (y > damping_south_start) { - horizontal_coeff_south = 0.0_rt; - } else if (y <= damping_south_end) { - horizontal_coeff_south = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (y - damping_south_start) / - (damping_south_end - damping_south_start)); - horizontal_coeff_south = term * term; - } - if (z <= horizontal_abl_height) { - vertical_coeff = 0.0_rt; - } else if (z > z_sloped) { - vertical_coeff = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (z - horizontal_abl_height) / - (z_sloped - horizontal_abl_height)); - vertical_coeff = term * term; - } - levelDamping[nbx](i, j, k, 0) = - vertical_coeff * - (horizontal_coeff_east + horizontal_coeff_north + - horizontal_coeff_west + horizontal_coeff_south); - //! Add the full vertical - if (z <= vertical_slope) { - vertical_coeff = 0.0_rt; - } else if (z > vertical_full) { - vertical_coeff = 1.0_rt; - } else { - const amrex::Real term = std::sin( - std::numbers::pi_v * 0.5_rt * - (z - vertical_slope) / - (vertical_full - vertical_slope + 1e-15_rt)); - vertical_coeff = term * term; - } - levelDamping[nbx](i, j, k, 0) = - amrex::min( - vertical_coeff + levelDamping[nbx](i, j, k, 0), 1.0_rt) / - horizontal_tau; - }); + if (m_terrain_damping != nullptr) { + damping_legacy::compute_terrain_damping( + *m_terrain_damping, level, geom, m_damp_east_slope, + m_damp_east_full, m_damp_west_slope, m_damp_west_full, + m_damp_north_slope, m_damp_north_full, m_damp_south_slope, + m_damp_south_full, m_horizontal_tau, m_horizontal_abl_height, + m_horizontal_slope_end, m_vertical_slope, m_vertical_full); + } } void TerrainDrag::post_init_actions() diff --git a/src/physics/TerrainDragDamping.H b/src/physics/TerrainDragDamping.H new file mode 100644 index 0000000000..c744c8c3d6 --- /dev/null +++ b/src/physics/TerrainDragDamping.H @@ -0,0 +1,194 @@ +#ifndef TERRAINDRAGDAMPING_H +#define TERRAINDRAGDAMPING_H + +#include +#include + +#include "AMReX_Geometry.H" +#include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" +#include "src/core/Field.H" + +namespace kynema_sgf::terraindrag::damping_legacy { + +inline bool has_damping_inputs(const amrex::ParmParse& pp) +{ + return pp.contains("damp_east_slope") || pp.contains("damp_east_full") || + pp.contains("damp_west_slope") || pp.contains("damp_west_full") || + pp.contains("damp_north_slope") || pp.contains("damp_north_full") || + pp.contains("damp_south_slope") || pp.contains("damp_south_full") || + pp.contains("horizontal_time_scale") || + pp.contains("horizontal_abl_height") || + pp.contains("horizontal_slope_end") || + pp.contains("vertical_slope") || pp.contains("vertical_full"); +} + +inline void query_damping_inputs( + amrex::ParmParse& pp, + amrex::Real& damp_east_slope, + amrex::Real& damp_east_full, + amrex::Real& damp_west_slope, + amrex::Real& damp_west_full, + amrex::Real& damp_north_slope, + amrex::Real& damp_north_full, + amrex::Real& damp_south_slope, + amrex::Real& damp_south_full, + amrex::Real& horizontal_tau, + amrex::Real& horizontal_abl_height, + amrex::Real& horizontal_slope_end, + amrex::Real& vertical_slope, + amrex::Real& vertical_full) +{ + pp.query("damp_east_slope", damp_east_slope); + pp.query("damp_east_full", damp_east_full); + pp.query("damp_west_slope", damp_west_slope); + pp.query("damp_west_full", damp_west_full); + pp.query("damp_north_slope", damp_north_slope); + pp.query("damp_north_full", damp_north_full); + pp.query("damp_south_slope", damp_south_slope); + pp.query("damp_south_full", damp_south_full); + pp.query("horizontal_time_scale", horizontal_tau); + pp.query("horizontal_abl_height", horizontal_abl_height); + pp.query("horizontal_slope_end", horizontal_slope_end); + pp.query("vertical_slope", vertical_slope); + pp.query("vertical_full", vertical_full); +} + +inline void compute_terrain_damping( + Field& terrain_damping, + int level, + const amrex::Geometry& geom, + amrex::Real damp_east_slope, + amrex::Real damp_east_full, + amrex::Real damp_west_slope, + amrex::Real damp_west_full, + amrex::Real damp_north_slope, + amrex::Real damp_north_full, + amrex::Real damp_south_slope, + amrex::Real damp_south_full, + amrex::Real horizontal_tau, + amrex::Real horizontal_abl_height, + amrex::Real horizontal_slope_end, + amrex::Real vertical_slope, + amrex::Real vertical_full) +{ + using namespace amrex::literals; + + auto& damping = terrain_damping(level); + auto level_damping = damping.arrays(); + + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + const auto& prob_hi = geom.ProbHiArray(); + + const amrex::Real damping_east_start = + prob_hi[0] - (damp_east_full + damp_east_slope); + const amrex::Real damping_east_end = prob_hi[0] - damp_east_full; + const amrex::Real damping_west_start = + prob_lo[0] + (damp_west_full + damp_west_slope); + const amrex::Real damping_west_end = prob_lo[0] + damp_west_full; + const amrex::Real damping_north_start = + prob_hi[1] - (damp_north_full + damp_north_slope); + const amrex::Real damping_north_end = prob_hi[1] - damp_north_full; + const amrex::Real damping_south_start = + prob_lo[1] + (damp_south_full + damp_south_slope); + const amrex::Real damping_south_end = prob_lo[1] + damp_south_full; + + amrex::ParallelFor( + damping, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + amrex::Real horizontal_coeff_east = 0.0_rt; + amrex::Real horizontal_coeff_north = 0.0_rt; + amrex::Real horizontal_coeff_west = 0.0_rt; + amrex::Real horizontal_coeff_south = 0.0_rt; + amrex::Real vertical_coeff = 0.0_rt; + const amrex::Real x = prob_lo[0] + (i + 0.5_rt) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5_rt) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5_rt) * dx[2]; + + if (x < damping_east_start) { + horizontal_coeff_east = 0.0_rt; + } else if (x >= damping_east_end) { + horizontal_coeff_east = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (x - damping_east_start) / + (damping_east_end - damping_east_start)); + horizontal_coeff_east = term * term; + } + + if (x > damping_west_start) { + horizontal_coeff_west = 0.0_rt; + } else if (x <= damping_west_end) { + horizontal_coeff_west = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (x - damping_west_start) / + (damping_west_end - damping_west_start)); + horizontal_coeff_west = term * term; + } + + if (y < damping_north_start) { + horizontal_coeff_north = 0.0_rt; + } else if (y >= damping_north_end) { + horizontal_coeff_north = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (y - damping_north_start) / + (damping_north_end - damping_north_start)); + horizontal_coeff_north = term * term; + } + + if (y > damping_south_start) { + horizontal_coeff_south = 0.0_rt; + } else if (y <= damping_south_end) { + horizontal_coeff_south = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (y - damping_south_start) / + (damping_south_end - damping_south_start)); + horizontal_coeff_south = term * term; + } + + if (z <= horizontal_abl_height) { + vertical_coeff = 0.0_rt; + } else if (z > horizontal_slope_end) { + vertical_coeff = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (z - horizontal_abl_height) / + (horizontal_slope_end - horizontal_abl_height)); + vertical_coeff = term * term; + } + + level_damping[nbx](i, j, k, 0) = + vertical_coeff * + (horizontal_coeff_east + horizontal_coeff_north + + horizontal_coeff_west + horizontal_coeff_south); + + if (z <= vertical_slope) { + vertical_coeff = 0.0_rt; + } else if (z > vertical_full) { + vertical_coeff = 1.0_rt; + } else { + const amrex::Real term = std::sin( + std::numbers::pi_v * 0.5_rt * + (z - vertical_slope) / + (vertical_full - vertical_slope + 1e-15_rt)); + vertical_coeff = term * term; + } + + level_damping[nbx](i, j, k, 0) = + amrex::min( + vertical_coeff + level_damping[nbx](i, j, k, 0), 1.0_rt) / + horizontal_tau; + }); +} + +} // namespace kynema_sgf::terraindrag::damping_legacy + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7c7b005e2b..87e38ae19e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -217,6 +217,7 @@ add_test_r(abl_amd_wenoz) add_test_r(abl_anelastic) add_test_r(abl_bndry_output_amr_native) add_test_r(abl_bndry_output_native) +add_test_r(abl_damping_layer_source) add_test_r(abl_godunov_geostrophic_timetable) add_test_r(abl_godunov_mpl) add_test_r(abl_godunov_mpl_amr) diff --git a/test/test_files/abl_damping_layer_source/abl_damping_layer_source.inp b/test/test_files/abl_damping_layer_source/abl_damping_layer_source.inp new file mode 100644 index 0000000000..f694b1b751 --- /dev/null +++ b/test/test_files/abl_damping_layer_source/abl_damping_layer_source.inp @@ -0,0 +1,131 @@ +# Regression test for DampingLayer + DampingLayerSource feature combinations + +# Simulation control +time.stop_time = 22000.0 +time.max_step = 10 + +time.fixed_dt = 0.5 +time.cfl = 0.95 + +# Output +# test/CMakeLists overrides plot/checkpoint controls for CI runs + +# Physics and numerics +incflo.gravity = 0.0 0.0 -9.81 +incflo.density = 1.0 + +incflo.use_godunov = 1 +incflo.diffusion_type = 1 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +transport.reference_temperature = 300.0 + +incflo.physics = ABL DampingLayer +ICNS.source_terms = BoussinesqBuoyancy DampingLayerSource +Temperature.source_terms = DampingLayerSource + +incflo.post_processing = avg1 +avg1.type = TimeAveraging +avg1.labels = means +avg1.averaging_window = 10.0 +avg1.means.fields = velocity +avg1.means.averaging_type = ReAveraging + +incflo.velocity = 6.0 2.0 0.0 +ABL.temperature_heights = 0.0 1000.0 +ABL.temperature_values = 300.0 304.0 +ABL.kappa = 0.41 +ABL.surface_roughness_z0 = 0.15 + +# DampingLayer physics configuration +DampingLayer.fields = velocity temperature + +DampingLayer.velocity.xlo.thickness = 200.0 +DampingLayer.velocity.xlo.blending_fraction = 0.25 +DampingLayer.velocity.xlo.blending_function_type = linear +DampingLayer.velocity.xlo.minimum_height = 150.0 +DampingLayer.velocity.xlo.vertical_blending_thickness = 120.0 +DampingLayer.velocity.xlo.vertical_blending_function_type = cosine + +DampingLayer.velocity.xhi.thickness = 180.0 +DampingLayer.velocity.xhi.blending_fraction = 0.40 +DampingLayer.velocity.xhi.blending_function_type = quadratic + +DampingLayer.velocity.ylo.thickness = 150.0 +DampingLayer.velocity.ylo.blending_fraction = 0.50 +DampingLayer.velocity.ylo.blending_function_type = exponential + +DampingLayer.velocity.yhi.thickness = 150.0 +DampingLayer.velocity.yhi.blending_fraction = 0.35 +DampingLayer.velocity.yhi.blending_function_type = cosine +DampingLayer.velocity.yhi.minimum_height = 120.0 +DampingLayer.velocity.yhi.vertical_blending_thickness = 100.0 +DampingLayer.velocity.yhi.vertical_blending_function_type = linear + +DampingLayer.velocity.zlo.thickness = 80.0 +DampingLayer.velocity.zlo.blending_fraction = 0.20 +DampingLayer.velocity.zlo.blending_function_type = cosine + +DampingLayer.velocity.zhi.thickness = 120.0 +DampingLayer.velocity.zhi.blending_fraction = 0.30 +DampingLayer.velocity.zhi.blending_function_type = cosine + +DampingLayer.temperature.xlo.thickness = 120.0 +DampingLayer.temperature.xlo.blending_fraction = 0.5 +DampingLayer.temperature.xlo.blending_function_type = linear + +DampingLayer.temperature.zhi.thickness = 140.0 +DampingLayer.temperature.zhi.blending_fraction = 0.4 +DampingLayer.temperature.zhi.blending_function_type = cosine + +# DampingLayerSource targets for velocity +DampingLayer.velocity.xlo.target_type = constant +DampingLayer.velocity.xlo.target_value = 8.0 1.0 0.0 +DampingLayer.velocity.xlo.damped_components = 1 1 0 + +DampingLayer.velocity.xhi.target_type = function +DampingLayer.velocity.xhi.target_function = "6.0 + 0.001*z + 0.2*n" + +DampingLayer.velocity.ylo.target_type = profile +DampingLayer.velocity.ylo.target_profile_heights = 0.0 500.0 1000.0 +DampingLayer.velocity.ylo.target_profile_values_x = 5.5 6.0 6.5 +DampingLayer.velocity.ylo.target_profile_values_y = 0.5 1.0 1.5 +DampingLayer.velocity.ylo.target_profile_values_z = 0.0 0.0 0.0 + +DampingLayer.velocity.yhi.target_type = field +DampingLayer.velocity.yhi.target_field = velocity_mean_avg1 + +DampingLayer.velocity.zlo.target_type = constant +DampingLayer.velocity.zlo.target_value = 6.0 0.0 0.0 +DampingLayer.velocity.zlo.damped_components = 1 0 0 + +DampingLayer.velocity.zhi.target_type = constant +DampingLayer.velocity.zhi.target_value = 6.0 0.0 0.0 +DampingLayer.velocity.zhi.damped_components = 1 0 1 + +# DampingLayerSource targets for temperature +# Use legacy profile key to exercise compatibility path +DampingLayer.temperature.xlo.target_type = profile +DampingLayer.temperature.xlo.target_profile_heights = 0.0 500.0 1000.0 +DampingLayer.temperature.xlo.target_profile_values = 300.0 301.5 303.0 + +DampingLayer.temperature.zhi.target_type = function +DampingLayer.temperature.zhi.target_function = "300.0 + 0.003*z + 0.1*t" + +# Mesh +amr.n_cell = 32 32 32 +amr.max_level = 0 + +# Domain +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1000.0 1000.0 1000.0 +geometry.is_periodic = 1 1 0 + +# Boundary conditions +zlo.type = wall_model +zhi.type = slip_wall +zhi.temperature_type = fixed_gradient +zhi.temperature = 0.003 + +incflo.verbose = 0 diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index af7d2ae46f..fb3aa3e473 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -17,6 +17,7 @@ add_subdirectory(ocean_waves) add_subdirectory(projection) add_subdirectory(boundary_conditions) add_subdirectory(eos_models) +add_subdirectory(physics) if(KYNEMA_SGF_ENABLE_MASA) add_subdirectory(mms) diff --git a/unit_tests/equation_systems/CMakeLists.txt b/unit_tests/equation_systems/CMakeLists.txt index e76ca68794..913d39b89b 100644 --- a/unit_tests/equation_systems/CMakeLists.txt +++ b/unit_tests/equation_systems/CMakeLists.txt @@ -5,4 +5,5 @@ target_sources(${kynema_sgf_unit_test_exe_name} test_icns_gravityforcing.cpp test_icns_init.cpp test_explicit_diffusion_rk2.cpp + test_damping_layer_source.cpp ) diff --git a/unit_tests/equation_systems/test_damping_layer_source.cpp b/unit_tests/equation_systems/test_damping_layer_source.cpp new file mode 100644 index 0000000000..c85f644541 --- /dev/null +++ b/unit_tests/equation_systems/test_damping_layer_source.cpp @@ -0,0 +1,28 @@ +#include "gtest/gtest.h" +#include "src/equation_systems/source_terms/DampingLayerSource.H" + +namespace kynema_sgf_tests { + +TEST(DampingLayerSource, target_type_parse) +{ + EXPECT_EQ(string_to_target_type("constant"), TargetType::Constant); + EXPECT_EQ(string_to_target_type("PROFILE"), TargetType::Profile); + EXPECT_EQ(string_to_target_type("Function"), TargetType::Function); + EXPECT_EQ(string_to_target_type("field"), TargetType::Field); +} + +TEST(DampingLayerSource, trait_field_names) +{ + using namespace kynema_sgf::pde; + + EXPECT_STREQ( + DampingLayerSourceTraits::field_name, "velocity"); + EXPECT_STREQ( + DampingLayerSourceTraits::field_name, "temperature"); + EXPECT_STREQ( + DampingLayerSourceTraits::field_name, "density"); + EXPECT_STREQ(DampingLayerSourceTraits::field_name, "tke"); + EXPECT_STREQ(DampingLayerSourceTraits::field_name, "sdr"); +} + +} // namespace kynema_sgf_tests diff --git a/unit_tests/physics/CMakeLists.txt b/unit_tests/physics/CMakeLists.txt new file mode 100644 index 0000000000..a005822ead --- /dev/null +++ b/unit_tests/physics/CMakeLists.txt @@ -0,0 +1,4 @@ +target_sources( + ${kynema_sgf_unit_test_exe_name} PRIVATE + test_damping_layer.cpp + ) \ No newline at end of file diff --git a/unit_tests/physics/test_damping_layer.cpp b/unit_tests/physics/test_damping_layer.cpp new file mode 100644 index 0000000000..a6ef6487f4 --- /dev/null +++ b/unit_tests/physics/test_damping_layer.cpp @@ -0,0 +1,92 @@ +#include "gtest/gtest.h" +#include "src/physics/DampingLayer.H" +#include "src/utilities/constants.H" +#include "AMReX_REAL.H" + +using namespace amrex::literals; + +namespace kynema_sgf_tests { + +TEST(DampingLayerMath, blending_function_type_parse) +{ + using namespace kynema_sgf::damping_layer; + + EXPECT_EQ( + string_to_blending_function_type("linear"), + BlendingFunctionType::Linear); + EXPECT_EQ( + string_to_blending_function_type("QUADRATIC"), + BlendingFunctionType::Quadratic); + EXPECT_EQ( + string_to_blending_function_type("Exponential"), + BlendingFunctionType::Exponential); + EXPECT_EQ( + string_to_blending_function_type("cosine"), + BlendingFunctionType::Cosine); +} + +TEST(DampingLayerMath, blending_function_values) +{ + using namespace kynema_sgf::damping_layer; + + constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL; + + EXPECT_NEAR( + blending_function(0.0_rt, BlendingFunctionType::Linear), 1.0_rt, tol); + EXPECT_NEAR( + blending_function(1.0_rt, BlendingFunctionType::Linear), 0.0_rt, tol); + + EXPECT_NEAR( + blending_function(0.5_rt, BlendingFunctionType::Quadratic), 0.75_rt, + tol); + + const amrex::Real exp_mid = + blending_function(0.5_rt, BlendingFunctionType::Exponential); + EXPECT_GT(exp_mid, 0.0_rt); + EXPECT_LT(exp_mid, 1.0_rt); + + EXPECT_NEAR( + blending_function(0.0_rt, BlendingFunctionType::Cosine), 1.0_rt, tol); + EXPECT_NEAR( + blending_function(0.5_rt, BlendingFunctionType::Cosine), 0.5_rt, tol); + EXPECT_NEAR( + blending_function(1.0_rt, BlendingFunctionType::Cosine), 0.0_rt, tol); +} + +TEST(DampingLayerMath, damping_calc_piecewise_behavior) +{ + using namespace kynema_sgf::damping_layer; + + constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL; + constexpr amrex::Real thickness = 10.0_rt; + constexpr amrex::Real blend_frac = 0.3_rt; + constexpr amrex::Real full_damp_len = thickness * (1.0_rt - blend_frac); + + EXPECT_NEAR( + damping_calc( + 2.0_rt, thickness, blend_frac, BlendingFunctionType::Cosine), + 1.0_rt, tol); + + const amrex::Real blend_pos = + (8.0_rt - full_damp_len) / (thickness * blend_frac); + const amrex::Real expected_mid = + blending_function(blend_pos, BlendingFunctionType::Cosine); + EXPECT_NEAR( + damping_calc( + 8.0_rt, thickness, blend_frac, BlendingFunctionType::Cosine), + expected_mid, tol); + + EXPECT_NEAR( + damping_calc( + 11.0_rt, thickness, blend_frac, BlendingFunctionType::Cosine), + 0.0_rt, tol); + + EXPECT_NEAR( + damping_calc(5.0_rt, thickness, 0.0_rt, BlendingFunctionType::Linear), + 1.0_rt, tol); + EXPECT_NEAR( + damping_calc(10.01_rt, thickness, 0.0_rt, BlendingFunctionType::Linear), + 0.0_rt, tol); +} + +} // namespace kynema_sgf_tests