From 8ad05ae57e8db712c0b3edad99734ca08bde8675 Mon Sep 17 00:00:00 2001 From: Jon Rood Date: Wed, 15 Jan 2025 09:59:47 -0700 Subject: [PATCH 1/5] Increase test timeout. (#1448) --- test/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 11e6aecedc..c36d15eec2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -87,7 +87,7 @@ function(add_test_r TEST_NAME) add_test(${TEST_NAME} bash -c "set -o pipefail && ${MPI_COMMANDS} ${CMAKE_BINARY_DIR}/${amr_wind_exe_name} ${MPIEXEC_POSTFLAGS} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.inp ${RUNTIME_OPTIONS} 2>&1 | tee ${TEST_NAME}.log ${SAVE_GOLDS_COMMAND} ${FCOMPARE_COMMAND}") # Set properties for test set_tests_properties(${TEST_NAME} PROPERTIES - TIMEOUT 7200 + TIMEOUT 10800 PROCESSORS ${TEST_NP} WORKING_DIRECTORY "${CURRENT_TEST_BINARY_DIR}/" LABELS "regression" @@ -157,7 +157,7 @@ function(add_test_pps TEST_NAME FORMAT NINT NSTEPS) add_test("${TEST_NAME}_sampling_${FORMAT}" bash -c "set -o pipefail && ${MPI_COMMANDS} ${CMAKE_BINARY_DIR}/${amr_wind_exe_name} ${MPIEXEC_POSTFLAGS} ${CURRENT_TEST_BINARY_DIR}/${TEST_NAME}.inp ${RUNTIME_OPTIONS} ${ADDL_RUNTIME_OPTIONS} 2>&1 | tee ${TEST_NAME}_sampling_${FORMAT}.log && ${Python3_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tools/sampling_${TEST_NAME}_${FORMAT}.py") # Set properties for test set_tests_properties("${TEST_NAME}_sampling_${FORMAT}" PROPERTIES - TIMEOUT 7200 + TIMEOUT 10800 PROCESSORS ${TEST_NP} WORKING_DIRECTORY "${CURRENT_TEST_BINARY_DIR}/" LABELS "post_processing;no_ci" From d0a56a426a0d15f0b838f5d86ea77b1ebf512c3c Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Wed, 15 Jan 2025 14:19:02 -0700 Subject: [PATCH 2/5] Fix sampling when more than 2 labels (#1449) --- .../utilities/sampling/SamplingContainer.cpp | 5 +++-- test/test_files/abl_sampling/abl_sampling.inp | 10 +++++++++- .../abl_sampling_netcdf.inp | 18 +++++++++++++++++- 3 files changed, 29 insertions(+), 4 deletions(-) diff --git a/amr-wind/utilities/sampling/SamplingContainer.cpp b/amr-wind/utilities/sampling/SamplingContainer.cpp index ae5df0a0ef..00b19df1d5 100644 --- a/amr-wind/utilities/sampling/SamplingContainer.cpp +++ b/amr-wind/utilities/sampling/SamplingContainer.cpp @@ -80,8 +80,9 @@ void SamplingContainer::initialize_particles( continue; } - const auto uid_offset = - (iprobe == 0) ? 0 : samplers[iprobe - 1]->num_points(); + const int uid_offset = std::accumulate( + samplers.begin(), samplers.begin() + iprobe, 0, + [&](int sum, const auto& s) { return sum + s->num_points(); }); const auto probe_id = probe->id(); amrex::Gpu::DeviceVector dlocs(npts); amrex::Gpu::copy( diff --git a/test/test_files/abl_sampling/abl_sampling.inp b/test/test_files/abl_sampling/abl_sampling.inp index 04de7dc6b9..ef4741c0bb 100644 --- a/test/test_files/abl_sampling/abl_sampling.inp +++ b/test/test_files/abl_sampling/abl_sampling.inp @@ -128,7 +128,7 @@ probe_sampling.probe1.offsets = 0.0 200.0 400.0 600.0 plane_sampling.output_frequency = 1 plane_sampling.output_format = native plane_sampling.fields = velocity -plane_sampling.labels = plane1 plane2 +plane_sampling.labels = plane1 plane2 plane3 plane_sampling.plane1.type = PlaneSampler plane_sampling.plane1.axis1 = 400.0 0.0 0.0 plane_sampling.plane1.axis2 = 0.0 0.0 400.0 @@ -145,6 +145,14 @@ plane_sampling.plane2.num_points = 10 10 plane_sampling.plane2.offset_vector = 0.0 0.0 1.0 plane_sampling.plane2.offsets = 0.0 499.0 +plane_sampling.plane3.type = PlaneSampler +plane_sampling.plane3.axis1 = 400.0 0.0 0.0 +plane_sampling.plane3.axis2 = 0.0 0.0 400.0 +plane_sampling.plane3.origin = 400.0 200.0 200.0 +plane_sampling.plane3.num_points = 10 10 +plane_sampling.plane3.offset_vector = 0.0 1.0 1.0 +plane_sampling.plane3.offsets = 0.0 20.0 30.0 200.0 + line_sampling.output_frequency = 1 line_sampling.output_format = native line_sampling.fields = velocity diff --git a/test/test_files/abl_sampling_netcdf/abl_sampling_netcdf.inp b/test/test_files/abl_sampling_netcdf/abl_sampling_netcdf.inp index 4b0df21292..f96df35670 100644 --- a/test/test_files/abl_sampling_netcdf/abl_sampling_netcdf.inp +++ b/test/test_files/abl_sampling_netcdf/abl_sampling_netcdf.inp @@ -128,7 +128,7 @@ probe_sampling.probe1.offsets = 0.0 200.0 400.0 600.0 plane_sampling.output_frequency = 1 plane_sampling.output_format = netcdf plane_sampling.fields = velocity -plane_sampling.labels = plane1 +plane_sampling.labels = plane1 plane2 plane3 plane_sampling.plane1.type = PlaneSampler plane_sampling.plane1.axis1 = 400.0 0.0 0.0 plane_sampling.plane1.axis2 = 0.0 0.0 400.0 @@ -137,6 +137,22 @@ plane_sampling.plane1.num_points = 10 10 plane_sampling.plane1.offset_vector = 0.0 1.0 1.0 plane_sampling.plane1.offsets = 0.0 20.0 30.0 200.0 +plane_sampling.plane2.type = PlaneSampler +plane_sampling.plane2.axis1 = 1000.0 0.0 0.0 +plane_sampling.plane2.axis2 = 0.0 1000.0 0.0 +plane_sampling.plane2.origin = 0.0 0.0 500.0 +plane_sampling.plane2.num_points = 10 10 +plane_sampling.plane2.offset_vector = 0.0 0.0 1.0 +plane_sampling.plane2.offsets = 0.0 499.0 + +plane_sampling.plane3.type = PlaneSampler +plane_sampling.plane3.axis1 = 400.0 0.0 0.0 +plane_sampling.plane3.axis2 = 0.0 0.0 400.0 +plane_sampling.plane3.origin = 400.0 200.0 200.0 +plane_sampling.plane3.num_points = 10 10 +plane_sampling.plane3.offset_vector = 0.0 1.0 1.0 +plane_sampling.plane3.offsets = 0.0 20.0 30.0 200.0 + line_sampling.output_frequency = 1 line_sampling.output_format = netcdf line_sampling.fields = velocity From 2dd88056cc5a71eab3f4839496bd619398e1632d Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Thu, 16 Jan 2025 07:19:36 -0700 Subject: [PATCH 3/5] Update anelastic to define and use a reference temperature field (#1446) --- .../icns/source_terms/ABLMeanBoussinesq.H | 3 - .../icns/source_terms/ABLMeanBoussinesq.cpp | 9 +-- .../icns/source_terms/BoussinesqBuoyancy.H | 3 - .../icns/source_terms/BoussinesqBuoyancy.cpp | 8 +-- .../source_terms/DragTempForcing.H | 3 - .../source_terms/DragTempForcing.cpp | 7 +- .../tke/source_terms/KransAxell.H | 3 - .../tke/source_terms/KransAxell.cpp | 5 +- amr-wind/transport_models/ConstTransport.H | 65 ++++++++++++++++--- amr-wind/transport_models/TransportModel.H | 9 +++ amr-wind/transport_models/TwoPhaseTransport.H | 65 ++++++++++++++++--- amr-wind/utilities/constants.H | 10 +++ amr-wind/wind_energy/ABLAnelastic.H | 1 + amr-wind/wind_energy/ABLAnelastic.cpp | 20 ++++++ .../abl_anelastic/abl_anelastic.inp | 3 +- 15 files changed, 170 insertions(+), 44 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H index d4da314be0..767e108e31 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H @@ -46,9 +46,6 @@ private: //! Transport model const transport::TransportModel& m_transport; - //! Reference temperature - std::unique_ptr m_ref_theta; - int m_axis{2}; bool m_const_profile{false}; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp index 9d6d4031b7..d15391f8d7 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp @@ -16,8 +16,6 @@ ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim) const auto& abl = sim.physics_manager().get(); abl.register_mean_boussinesq_term(this); - m_ref_theta = m_transport.ref_theta(); - // gravity in `incflo` namespace amrex::ParmParse pp_incflo("incflo"); pp_incflo.queryarr("gravity", m_gravity); @@ -96,7 +94,10 @@ void ABLMeanBoussinesq::operator()( amrex::Array4 const& beta_arr = beta_fab.array(); m_transport.beta_impl(lev, mfi, bx, beta_arr); - const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); + amrex::FArrayBox ref_theta_fab(bx, 1, amrex::The_Async_Arena()); + amrex::Array4 const& ref_theta_arr = ref_theta_fab.array(); + m_transport.ref_theta_impl(lev, mfi, bx, ref_theta_arr); + const amrex::GpuArray gravity{ m_gravity[0], m_gravity[1], m_gravity[2]}; @@ -110,7 +111,7 @@ void ABLMeanBoussinesq::operator()( amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::IntVect iv(i, j, k); const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir]; - const amrex::Real T0 = ref_theta(i, j, k); + const amrex::Real T0 = ref_theta_arr(i, j, k); const amrex::Real temp = amr_wind::interp::linear(theights, theights_end, tvals, ht); const amrex::Real fac = beta_arr(i, j, k) * (temp - T0); diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H index 4774fd79ed..2466479f67 100644 --- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H +++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H @@ -37,9 +37,6 @@ private: //! Transport model const transport::TransportModel& m_transport; - - //! Reference temperature - std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp index ab46706467..eb2d9b4621 100644 --- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp +++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp @@ -13,8 +13,6 @@ BoussinesqBuoyancy::BoussinesqBuoyancy(const CFDSim& sim) : m_temperature(sim.repo().get_field("temperature")) , m_transport(sim.transport_model()) { - m_ref_theta = m_transport.ref_theta(); - // gravity in `incflo` namespace amrex::ParmParse pp_incflo("incflo"); pp_incflo.queryarr("gravity", m_gravity); @@ -38,10 +36,12 @@ void BoussinesqBuoyancy::operator()( amrex::FArrayBox beta_fab(bx, 1, amrex::The_Async_Arena()); amrex::Array4 const& beta_arr = beta_fab.array(); m_transport.beta_impl(lev, mfi, bx, beta_arr); - const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); + amrex::FArrayBox ref_theta_fab(bx, 1, amrex::The_Async_Arena()); + amrex::Array4 const& ref_theta_arr = ref_theta_fab.array(); + m_transport.ref_theta_impl(lev, mfi, bx, ref_theta_arr); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real T = temp(i, j, k, 0); - const amrex::Real T0 = ref_theta(i, j, k); + const amrex::Real T0 = ref_theta_arr(i, j, k); const amrex::Real fac = beta_arr(i, j, k) * (T0 - T); src_term(i, j, k, 0) += gravity[0] * fac; diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H index 6c2c27aafd..5403b403d0 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H @@ -33,9 +33,6 @@ private: //! Transport model const transport::TransportModel& m_transport; - - //! Reference temperature - std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::temperature diff --git a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp index f699820d0c..1c7f89da8c 100644 --- a/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp @@ -17,7 +17,6 @@ DragTempForcing::DragTempForcing(const CFDSim& sim) { amrex::ParmParse pp("DragTempForcing"); pp.query("drag_coefficient", m_drag_coefficient); - m_ref_theta = m_transport.ref_theta(); if (pp.contains("reference_temperature")) { amrex::Abort( "DragTempForcing.reference_temperature has been deprecated. Please " @@ -50,7 +49,9 @@ void DragTempForcing::operator()( const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const amrex::Real drag_coefficient = m_drag_coefficient / dx[2]; - const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi); + amrex::FArrayBox ref_theta_fab(bx, 1, amrex::The_Async_Arena()); + amrex::Array4 const& ref_theta_arr = ref_theta_fab.array(); + m_transport.ref_theta_impl(lev, mfi, bx, ref_theta_arr); const auto tiny = std::numeric_limits::epsilon(); const amrex::Real cd_max = 10.0; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -60,7 +61,7 @@ void DragTempForcing::operator()( const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1); const amrex::Real Cd = std::min(drag_coefficient / (m + tiny), cd_max / dx[2]); - const amrex::Real T0 = ref_theta(i, j, k); + const amrex::Real T0 = ref_theta_arr(i, j, k); src_term(i, j, k, 0) -= (Cd * (temperature(i, j, k, 0) - T0) * blank(i, j, k, 0)); }); diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.H b/amr-wind/equation_systems/tke/source_terms/KransAxell.H index a8e0f0d80a..2089a8cfc5 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.H +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.H @@ -48,9 +48,6 @@ private: //! Transport model const transport::TransportModel& m_transport; - - //! Reference temperature - std::unique_ptr m_ref_theta; }; } // namespace amr_wind::pde::tke diff --git a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp index df24297bbe..2d81bf5dc5 100644 --- a/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KransAxell.cpp @@ -30,7 +30,6 @@ KransAxell::KransAxell(const CFDSim& sim) amrex::ParmParse pp_incflow("incflo"); pp_incflow.queryarr("gravity", m_gravity); } - m_ref_theta = m_transport.ref_theta(); } KransAxell::~KransAxell() = default; @@ -49,7 +48,9 @@ void KransAxell::operator()( const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi); const auto& dissip_arr = (this->m_dissip)(lev).array(mfi); const auto& tke_arr = m_tke(lev).array(mfi); - const auto& ref_theta_arr = (*m_ref_theta)(lev).const_array(mfi); + amrex::FArrayBox ref_theta_fab(bx, 1, amrex::The_Async_Arena()); + amrex::Array4 const& ref_theta_arr = ref_theta_fab.array(); + m_transport.ref_theta_impl(lev, mfi, bx, ref_theta_arr); const auto& geom = m_mesh.Geom(lev); const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& probhi = m_mesh.Geom(lev).ProbHiArray(); diff --git a/amr-wind/transport_models/ConstTransport.H b/amr-wind/transport_models/ConstTransport.H index f7d8780906..f1a8857590 100644 --- a/amr-wind/transport_models/ConstTransport.H +++ b/amr-wind/transport_models/ConstTransport.H @@ -2,7 +2,6 @@ #define CONSTTRANSPORT_H #include "amr-wind/transport_models/TransportModel.H" -#include "amr-wind/utilities/constants.H" #include "AMReX_ParmParse.H" namespace amr_wind::transport { @@ -176,14 +175,26 @@ public: const amrex::Array4& beta) const override { - const amrex::Real beta_val = (m_constant_beta > 0.0) - ? m_constant_beta - : 1.0 / m_reference_temperature; - - amrex::ParallelFor( - bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - beta(i, j, k) = beta_val; - }); + if (m_constant_beta > 0.0) { + const amrex::Real beta_val = m_constant_beta; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = beta_val; + }); + } else if (m_repo.field_exists("reference_temperature")) { + const auto& temp0 = m_repo.get_field("reference_temperature"); + const auto& temp0_arr = temp0(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = 1.0 / temp0_arr(i, j, k); + }); + } else { + const amrex::Real beta_val = 1.0 / m_reference_temperature; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = beta_val; + }); + } if (m_repo.field_exists("vof")) { const auto& vof = m_repo.get_field("vof"); @@ -211,11 +222,45 @@ public: auto ref_theta = m_repo.create_scratch_field(1, m_ngrow); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - (*ref_theta)(lev).setVal(m_reference_temperature); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi((*ref_theta)(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& ref_theta_arr = (*ref_theta)(lev).array(mfi); + ref_theta_impl(lev, mfi, bx, ref_theta_arr); + } } return ref_theta; } + //! Compute the reference temperature + inline void ref_theta_impl( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const amrex::Array4& ref_theta) const override + { + if (m_reference_temperature < 0.0) { + amrex::Abort("Reference temperature was not set"); + } + + if (m_repo.field_exists("reference_temperature")) { + auto& temp0 = m_repo.get_field("reference_temperature"); + const auto& temp0_arr = temp0(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ref_theta(i, j, k) = temp0_arr(i, j, k); + }); + } else { + const amrex::Real ref_theta_val = m_reference_temperature; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ref_theta(i, j, k) = ref_theta_val; + }); + } + } + private: //! Reference to the field repository (for creating scratch fields) FieldRepo& m_repo; diff --git a/amr-wind/transport_models/TransportModel.H b/amr-wind/transport_models/TransportModel.H index 2a93636b78..9d6f1a650a 100644 --- a/amr-wind/transport_models/TransportModel.H +++ b/amr-wind/transport_models/TransportModel.H @@ -4,6 +4,8 @@ #include "amr-wind/core/Factory.H" #include "amr-wind/CFDSim.H" #include "amr-wind/core/FieldRepo.H" +#include "amr-wind/utilities/constants.H" +#include "amr-wind/core/field_ops.H" namespace amr_wind::transport { @@ -56,6 +58,13 @@ public: //! Reference temperature virtual amrex::Real reference_temperature() const = 0; + //! Reference temperature + virtual void ref_theta_impl( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const amrex::Array4& ref_theta) const = 0; + //! Reference temperature virtual std::unique_ptr ref_theta() const = 0; diff --git a/amr-wind/transport_models/TwoPhaseTransport.H b/amr-wind/transport_models/TwoPhaseTransport.H index 6e1a390d04..11fdaff90f 100644 --- a/amr-wind/transport_models/TwoPhaseTransport.H +++ b/amr-wind/transport_models/TwoPhaseTransport.H @@ -3,7 +3,6 @@ #include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/transport_models/TransportModel.H" -#include "amr-wind/utilities/constants.H" #include "AMReX_ParmParse.H" namespace amr_wind::transport { @@ -309,9 +308,27 @@ public: const amrex::Box& bx, const amrex::Array4& beta) const override { - const amrex::Real beta_val = (m_constant_beta > 0.0) - ? m_constant_beta - : 1.0 / m_reference_temperature; + + if (m_constant_beta > 0.0) { + const amrex::Real beta_val = m_constant_beta; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = beta_val; + }); + } else if (m_repo.field_exists("reference_temperature")) { + const auto& temp0 = m_repo.get_field("reference_temperature"); + const auto& temp0_arr = temp0(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = 1.0 / temp0_arr(i, j, k); + }); + } else { + const amrex::Real beta_val = 1.0 / m_reference_temperature; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + beta(i, j, k) = beta_val; + }); + } if (m_ifacetype == InterfaceCapturingMethod::VOF) { const auto& vof = m_repo.get_field("vof"); @@ -320,8 +337,6 @@ public: bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { if (vof_arr(i, j, k) > constants::TIGHT_TOL) { beta(i, j, k) = 0.0; - } else { - beta(i, j, k) = beta_val; } }); } else if (m_ifacetype == InterfaceCapturingMethod::LS) { @@ -342,7 +357,7 @@ public: 1.0 / M_PI * std::sin(phi_arr(i, j, k) * M_PI / eps)); } - beta(i, j, k) = beta_val * (1 - smooth_heaviside); + beta(i, j, k) *= (1 - smooth_heaviside); }); } } @@ -361,11 +376,45 @@ public: auto ref_theta = m_repo.create_scratch_field(1, m_ngrow); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { - (*ref_theta)(lev).setVal(m_reference_temperature); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi((*ref_theta)(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& ref_theta_arr = (*ref_theta)(lev).array(mfi); + ref_theta_impl(lev, mfi, bx, ref_theta_arr); + } } return ref_theta; } + //! Compute the reference temperature + inline void ref_theta_impl( + const int lev, + const amrex::MFIter& mfi, + const amrex::Box& bx, + const amrex::Array4& ref_theta) const override + { + if (m_reference_temperature < 0.0) { + amrex::Abort("Reference temperature was not set"); + } + + if (m_repo.field_exists("reference_temperature")) { + auto& temp0 = m_repo.get_field("reference_temperature"); + const auto& temp0_arr = temp0(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ref_theta(i, j, k) = temp0_arr(i, j, k); + }); + } else { + const amrex::Real ref_theta_val = m_reference_temperature; + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + ref_theta(i, j, k) = ref_theta_val; + }); + } + } + private: //! Reference to the CFD sim const CFDSim& m_sim; diff --git a/amr-wind/utilities/constants.H b/amr-wind/utilities/constants.H index 3d164671c0..4c2bf732a6 100644 --- a/amr-wind/utilities/constants.H +++ b/amr-wind/utilities/constants.H @@ -36,5 +36,15 @@ is_close(const amrex::Real a, const amrex::Real b) return amrex::Math::abs(a - b) < EPS; } +//! Boltzmann constant (J/K) +static constexpr amrex::Real BOLTZMANN_CONSTANT = 1.380649 * 1e-23; + +//! Avogadro's constant (mol^-1) +static constexpr amrex::Real AVOGADRO_CONSTANT = 6.02214076 * 1e23; + +//! Universal gas constant (J/K mol) +static constexpr amrex::Real UNIVERSAL_GAS_CONSTANT = + AVOGADRO_CONSTANT * BOLTZMANN_CONSTANT; + } // namespace amr_wind::constants #endif diff --git a/amr-wind/wind_energy/ABLAnelastic.H b/amr-wind/wind_energy/ABLAnelastic.H index f81a20401c..67377ca4d9 100644 --- a/amr-wind/wind_energy/ABLAnelastic.H +++ b/amr-wind/wind_energy/ABLAnelastic.H @@ -1,6 +1,7 @@ #ifndef ABLANELASTIC_H #define ABLANELASTIC_H #include "amr-wind/CFDSim.H" +#include "amr-wind/utilities/constants.H" #include "amr-wind/utilities/MultiLevelVector.H" #include "AMReX_ParmParse.H" diff --git a/amr-wind/wind_energy/ABLAnelastic.cpp b/amr-wind/wind_energy/ABLAnelastic.cpp index 82291ab182..387f5a23fd 100644 --- a/amr-wind/wind_energy/ABLAnelastic.cpp +++ b/amr-wind/wind_energy/ABLAnelastic.cpp @@ -35,6 +35,9 @@ ABLAnelastic::ABLAnelastic(CFDSim& sim) : m_sim(sim) "reference_density", 1, density.num_grow()[0], 1); ref_density.set_default_fillpatch_bc(m_sim.time()); m_sim.repo().declare_nd_field("reference_pressure", 1, 0, 1); + auto& ref_theta = m_sim.repo().declare_field( + "reference_temperature", 1, density.num_grow()[0], 1); + ref_theta.set_default_fillpatch_bc(m_sim.time()); } } @@ -82,5 +85,22 @@ void ABLAnelastic::initialize_data() m_pressure.copy_to_field(p0); rho0.fillpatch(m_sim.time().current_time()); + + auto& temp0 = m_sim.repo().get_field("reference_temperature"); + const amrex::Real air_molar_mass = 0.02896492; // kg/mol + const amrex::Real Rair = constants::UNIVERSAL_GAS_CONSTANT / air_molar_mass; + for (int lev = 0; lev < m_sim.repo().num_active_levels(); ++lev) { + const auto& rho0_arrs = rho0(lev).const_arrays(); + const auto& p0_arrs = p0(lev).const_arrays(); + const auto& temp0_arrs = temp0(lev).arrays(); + amrex::ParallelFor( + temp0(lev), amrex::IntVect(0), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + temp0_arrs[nbx](i, j, k) = + p0_arrs[nbx](i, j, k) / (Rair * rho0_arrs[nbx](i, j, k)); + }); + } + amrex::Gpu::synchronize(); + temp0.fillpatch(m_sim.time().current_time()); } } // namespace amr_wind diff --git a/test/test_files/abl_anelastic/abl_anelastic.inp b/test/test_files/abl_anelastic/abl_anelastic.inp index 91b55ac3f1..7372941a5f 100644 --- a/test/test_files/abl_anelastic/abl_anelastic.inp +++ b/test/test_files/abl_anelastic/abl_anelastic.inp @@ -38,12 +38,13 @@ transport.reference_temperature = 300.0 turbulence.model = Smagorinsky Smagorinsky_coeffs.Cs = 0.135 -io.outputs = reference_density reference_pressure +io.outputs = reference_density reference_pressure reference_temperature incflo.physics = ABL ICNS.source_terms = GravityForcing ABL.perturb_temperature = false ABL.perturb_velocity = false +ABL.stats_output_frequency = 1 incflo.velocity = 0.0 0.0 0.0 From 658384220af8ed2061555ec9c2a1507292b0f890 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Thu, 16 Jan 2025 09:36:15 -0700 Subject: [PATCH 4/5] add tolerance to boundary file write time check (#1447) --- amr-wind/wind_energy/ABLBoundaryPlane.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.cpp b/amr-wind/wind_energy/ABLBoundaryPlane.cpp index 4b789a2384..e98d4ba3f5 100644 --- a/amr-wind/wind_energy/ABLBoundaryPlane.cpp +++ b/amr-wind/wind_energy/ABLBoundaryPlane.cpp @@ -581,7 +581,7 @@ void ABLBoundaryPlane::write_file() // Only output data if at the desired timestep if ((t_step % m_write_frequency != 0) || ((m_io_mode != io_mode::output)) || - (time < m_out_start_time)) { + (time < m_out_start_time - constants::LOOSE_TOL)) { return; } From 9fee431aa3571ae0340e67949be3218fd410b63a Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Thu, 16 Jan 2025 14:32:29 -0700 Subject: [PATCH 5/5] Enable optional coupling to openfast 4.0 (#1441) --- CMakeLists.txt | 14 +++ .../actuator/turbine/fast/FastIface.H | 4 + .../actuator/turbine/fast/FastIface.cpp | 28 ++++- .../actuator/turbine/fast/fast_types.H | 8 +- .../actuator/turbine/fast/fast_wrapper.H | 100 ++++++++++++++++-- docs/sphinx/user/build.rst | 8 ++ test/CMakeLists.txt | 2 +- 7 files changed, 147 insertions(+), 17 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1b62d7cb15..1594ddd4db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -173,6 +173,20 @@ if(AMR_WIND_ENABLE_OPENFAST) target_compile_definitions(${amr_wind_lib_name} PUBLIC AMR_WIND_USE_OPENFAST) target_include_directories(${amr_wind_lib_name} PUBLIC ${OpenFAST_INCLUDE_DIRS}) target_link_libraries(${amr_wind_lib_name} PUBLIC ${OpenFAST_LIBRARIES}) + + set(AMR_WIND_OPENFAST_VERSION "3" CACHE STRING "OpenFAST major version") + if((AMR_WIND_OPENFAST_VERSION STREQUAL "master") OR (AMR_WIND_OPENFAST_VERSION STREQUAL "main") OR (AMR_WIND_OPENFAST_VERSION STREQUAL "develop")) + set(AMR_WIND_OPENFAST_VERSION "4.0") + endif() + if(NOT AMR_WIND_OPENFAST_VERSION VERSION_LESS "5") + message(FATAL_ERROR "AMR_WIND_OPENFAST_VERSION must be less than 5.") + endif() + string(REPLACE "." ";" OPENFAST_VERSION_LIST ${AMR_WIND_OPENFAST_VERSION}) + list(LENGTH OPENFAST_VERSION_LIST OPENFAST_VERSION_LIST_LENGTH) + if(OPENFAST_VERSION_LIST_LENGTH GREATER_EQUAL 1) + list(GET OPENFAST_VERSION_LIST 0 OPENFAST_VERSION_MAJOR) + endif() + target_compile_definitions(${amr_wind_lib_name} PUBLIC OPENFAST_VERSION_MAJOR=${OPENFAST_VERSION_MAJOR}) endif() if(AMR_WIND_ENABLE_ASCENT) diff --git a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.H b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.H index 4a2048350e..1ae0d521a9 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.H +++ b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.H @@ -81,6 +81,10 @@ protected: int m_num_sc_inputs{0}; int m_num_sc_outputs{0}; +#if OPENFAST_VERSION_MAJOR == 4 + int m_inflow_type{2}; +#endif + int m_num_sc_inputs_glob{0}; float m_init_sc_inputs_glob{0.0}; float m_init_sc_inputs_turbine{0.0}; diff --git a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp index 60995dbac1..b59dd27703 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp +++ b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp @@ -130,7 +130,7 @@ void FastIface::init_solution(const int local_id) AMREX_ALWAYS_ASSERT(m_is_initialized); auto& fi = *m_turbine_data[local_id]; - fast_func(FAST_OpFM_Solution0, &fi.tid_local); + fast_func(FAST_Solution0, &fi.tid_local); fi.is_solution0 = false; } @@ -168,7 +168,7 @@ void FastIface::advance_turbine(const int local_id) write_velocity_data(fi); for (int i = 0; i < fi.num_substeps; ++i, ++fi.time_index) { - fast_func(FAST_OpFM_Step, &fi.tid_local); + fast_func(FAST_Step, &fi.tid_local); } if (fi.chkpt_interval > 0 && @@ -220,6 +220,17 @@ void FastIface::fast_init_turbine(FastTurbine& fi) amrex::Array inp_file; copy_filename(fi.input_file, inp_file.begin()); +#if OPENFAST_VERSION_MAJOR == 4 + char out_file[fast_strlen()]; + fast_func( + FAST_ExtInfw_Init, &fi.tid_local, &fi.stop_time, inp_file.begin(), + &fi.tid_global, out_file, &m_num_sc_inputs_glob, &m_num_sc_inputs, + &m_num_sc_outputs, &m_init_sc_inputs_glob, &m_init_sc_inputs_turbine, + &fi.num_pts_blade, &fi.num_pts_tower, fi.base_pos.begin(), &abort_lev, + &fi.dt_cfd, &fi.dt_fast, &m_inflow_type, &fi.num_blades, + &fi.num_blade_elem, &fi.num_tower_elem, &fi.chord_cluster_type, + &fi.to_cfd, &fi.from_cfd, &fi.to_sc, &fi.from_sc); +#else fast_func( FAST_OpFM_Init, &fi.tid_local, &fi.stop_time, inp_file.begin(), &fi.tid_global, &m_num_sc_inputs_glob, &m_num_sc_inputs, @@ -227,6 +238,7 @@ void FastIface::fast_init_turbine(FastTurbine& fi) &fi.num_pts_blade, &fi.num_pts_tower, fi.base_pos.begin(), &abort_lev, &fi.dt_fast, &fi.num_blades, &fi.num_blade_elem, &fi.chord_cluster_type, &fi.to_cfd, &fi.from_cfd, &fi.to_sc, &fi.from_sc); +#endif { #ifdef AMR_WIND_USE_OPENFAST @@ -289,13 +301,13 @@ void FastIface::fast_replay_turbine(FastTurbine& fi) // restart fi.time_index = 0; read_velocity_data(fi, ncf, 0); - fast_func(FAST_OpFM_Solution0, &fi.tid_local); + fast_func(FAST_Solution0, &fi.tid_local); fi.is_solution0 = false; for (int ic = 0; ic < num_cfd_steps; ++ic) { read_velocity_data(fi, ncf, ic); for (int i = 0; i < fi.num_substeps; ++i, ++fi.time_index) { - fast_func(FAST_OpFM_Step, &fi.tid_local); + fast_func(FAST_Step, &fi.tid_local); } } AMREX_ALWAYS_ASSERT(fi.time_index == num_steps); @@ -321,10 +333,18 @@ void FastIface::fast_restart_turbine(FastTurbine& fi) amrex::Array chkpt_file; copy_filename(fi.checkpoint_file, chkpt_file.begin()); +#if OPENFAST_VERSION_MAJOR == 4 + fast_func( + FAST_ExtInfw_Restart, &fi.tid_local, chkpt_file.begin(), &abort_lev, + &fi.dt_fast, &m_inflow_type, &fi.num_blades, &fi.num_blade_elem, + &fi.num_tower_elem, &fi.time_index, &fi.to_cfd, &fi.from_cfd, &fi.to_sc, + &fi.from_sc); +#else fast_func( FAST_OpFM_Restart, &fi.tid_local, chkpt_file.begin(), &abort_lev, &fi.dt_fast, &fi.num_blades, &fi.num_blade_elem, &fi.time_index, &fi.to_cfd, &fi.from_cfd, &fi.to_sc, &fi.from_sc); +#endif { #ifdef AMR_WIND_USE_OPENFAST diff --git a/amr-wind/wind_energy/actuator/turbine/fast/fast_types.H b/amr-wind/wind_energy/actuator/turbine/fast/fast_types.H index ddd271ff30..ec08d0dc07 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/fast_types.H +++ b/amr-wind/wind_energy/actuator/turbine/fast/fast_types.H @@ -57,6 +57,10 @@ struct FastTurbine //! Total number of elements along the blade int num_blade_elem; +#if OPENFAST_VERSION_MAJOR == 4 + //! Total number of elements along the tower + int num_tower_elem; +#endif //! Node cluster type for the chord int chord_cluster_type{0}; @@ -87,8 +91,8 @@ struct FastTurbine // Data structures that are used to exchange between fast/cfd - exw_fast::OpFM_InputType to_cfd; - exw_fast::OpFM_OutputType from_cfd; + exw_fast::OfInputType to_cfd; + exw_fast::OfOutputType from_cfd; exw_fast::SC_DX_InputType to_sc; exw_fast::SC_DX_OutputType from_sc; diff --git a/amr-wind/wind_energy/actuator/turbine/fast/fast_wrapper.H b/amr-wind/wind_energy/actuator/turbine/fast/fast_wrapper.H index ccff5d36eb..12733d5990 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/fast_wrapper.H +++ b/amr-wind/wind_energy/actuator/turbine/fast/fast_wrapper.H @@ -2,6 +2,7 @@ #define FAST_WRAPPER_H namespace exw_fast { + #ifdef AMR_WIND_USE_OPENFAST extern "C" { #include "FAST_Library.h" @@ -19,11 +20,6 @@ inline constexpr int fast_strlen() { return INTERFACE_STRING_LENGTH; } #define ErrID_Severe 3 #define ErrID_Fatal 4 -struct OpFM_InputType -{}; -struct OpFM_OutputType -{}; - struct SC_DX_InputType {}; struct SC_DX_OutputType @@ -37,11 +33,7 @@ inline void FAST_AllocateTurbines(int* /*unused*/, int* /*unused*/, char* /*unused*/) {} inline void FAST_DeallocateTurbines(int* /*unused*/, char* /*unused*/) {} -inline void -FAST_OpFM_Solution0(int* /*unused*/, int* /*unused*/, char* /*unused*/) -{} -inline void FAST_OpFM_Step(int* /*unused*/, int* /*unused*/, char* /*unused*/) -{} + inline void FAST_CreateCheckpoint( int* /*unused*/, char* /*unused*/, int* /*unused*/, char* /*unused*/) {} @@ -55,6 +47,79 @@ inline void FAST_HubPosition( char* /*unused*/) {} +#if OPENFAST_VERSION_MAJOR == 4 +struct ExtInfw_InputType_t +{}; + +struct ExtInfw_OutputType_t +{}; + +inline void +FAST_CFD_Solution0(int* /*unused*/, int* /*unused*/, char* /*unused*/) +{} + +inline void FAST_CFD_Step(int* /*unused*/, int* /*unused*/, char* /*unused*/) {} + +inline void FAST_ExtInfw_Init( + int* /*unused*/, + double* /*unused*/, + const char* /*unused*/, + int* /*unused*/, + const char* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + float* /*unused*/, + float* /*unused*/, + int* /*unused*/, + int* /*unused*/, + float* /*unused*/, + int* /*unused*/, + double* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + OpFM_InputType* /*unused*/, + OpFM_OutputType* /*unused*/, + SC_DX_InputType* /*unused*/, + SC_DX_OutputType* /*unused*/, + int* /*unused*/, + char* /*unused*/) +{} + +inline void FAST_ExtInfw_Restart( + int* /*unused*/, + char* /*unused*/, + int* /*unused*/, + double* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + int* /*unused*/, + OpFM_InputType* /*unused*/, + OpFM_OutputType* /*unused*/, + SC_DX_InputType* /*unused*/, + SC_DX_OutputType* /*unused*/, + int* /*unused*/, + char* /*unused*/) +{} + +#else +struct OpFM_InputType +{}; + +struct OpFM_OutputType +{}; + +inline void +FAST_OpFM_Solution0(int* /*unused*/, int* /*unused*/, char* /*unused*/) +{} + +inline void FAST_OpFM_Step(int* /*unused*/, int* /*unused*/, char* /*unused*/) +{} + inline void FAST_OpFM_Init( int* /*unused*/, double* /*unused*/, @@ -97,6 +162,21 @@ inline void FAST_OpFM_Restart( char* /*unused*/) {} #endif + +#endif + +#if OPENFAST_VERSION_MAJOR == 4 +static constexpr auto& FAST_Solution0 = FAST_CFD_Solution0; +static constexpr auto& FAST_Step = FAST_CFD_Step; +using OfInputType = ExtInfw_InputType_t; +using OfOutputType = ExtInfw_OutputType_t; +#else +static constexpr auto& FAST_Solution0 = FAST_OpFM_Solution0; +static constexpr auto& FAST_Step = FAST_OpFM_Step; +using OfInputType = OpFM_InputType; +using OfOutputType = OpFM_OutputType; +#endif + } // namespace exw_fast #endif /* FAST_WRAPPER_H */ diff --git a/docs/sphinx/user/build.rst b/docs/sphinx/user/build.rst index 7b4f270adc..4cdc052fbb 100644 --- a/docs/sphinx/user/build.rst +++ b/docs/sphinx/user/build.rst @@ -103,6 +103,14 @@ Dependencies Enable NetCDF outputs. Default: OFF +.. cmakeval:: AMR_WIND_ENABLE_OPENFAST + + Enable OpenFAST coupling. Default: OFF + +.. cmakeval:: AMR_WIND_OPENFAST_VERSION + + OpenFAST version. Default: 3 + Other AMR-Wind specific options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c36d15eec2..e3cb8d7795 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -323,7 +323,7 @@ if(AMR_WIND_ENABLE_W2A) add_test_re(abl_multiphase_w2a) endif() -if(AMR_WIND_ENABLE_OPENFAST) +if(AMR_WIND_ENABLE_OPENFAST AND AMR_WIND_OPENFAST_VERSION VERSION_LESS "4") set(ACT_UNIFORM_ALM_TEST_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/test_files/act_uniform_alm) set(ACT_UNIFORM_ALM_TEST_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/test_files/act_uniform_alm) set(RTEST_ACT_UNIFORM_ALM_TEST_BINARY_DIR ${ACT_UNIFORM_ALM_TEST_BINARY_DIR}/r-test)