Improve terrain search and use bilinear interpolation. (#1328)

amr-wind/physics/TerrainDrag.cpp

@@ -6,6 +6,8 @@
 #include "AMReX_ParReduce.H"
 #include "amr-wind/utilities/trig_ops.H"
 #include "amr-wind/utilities/IOManager.H"
+#include "amr-wind/utilities/io_utils.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 namespace amr_wind::terraindrag {
 
@@ -21,27 +23,23 @@ TerrainDrag::TerrainDrag(CFDSim& sim)
     , m_terrainz0(sim.repo().declare_field("terrainz0", 1, 1, 1))
     , m_terrain_height(sim.repo().declare_field("terrain_height", 1, 1, 1))
 {
-    std::string terrainfile("terrain.amrwind");
-    std::ifstream file(terrainfile, std::ios::in);
-    if (!file.good()) {
-        amrex::Abort("Cannot find terrain.amrwind file");
-    }
-    amrex::Real value1, value2, value3;
-    while (file >> value1 >> value2 >> value3) {
-        m_xterrain.push_back(value1);
-        m_yterrain.push_back(value2);
-        m_zterrain.push_back(value3);
-    }
-    file.close();
+    std::string terrain_file("terrain.amrwind");
+    std::string roughness_file("terrain.roughness");
+    amrex::ParmParse pp(identifier());
+    pp.query("terrain_file", terrain_file);
+    pp.query("roughness_file", roughness_file);
+
+    ioutils::read_flat_grid_file(
+        terrain_file, m_xterrain, m_yterrain, m_zterrain);
+
     // No checks for the file as it is optional currently
-    std::string roughnessfile("terrain.roughness");
-    std::ifstream file1(roughnessfile, std::ios::in);
-    while (file1 >> value1 >> value2 >> value3) {
-        m_xrough.push_back(value1);
-        m_yrough.push_back(value2);
-        m_z0rough.push_back(value3);
+    std::ifstream file(roughness_file, std::ios::in);
+    if (file.good()) {
+        ioutils::read_flat_grid_file(
+            roughness_file, m_xrough, m_yrough, m_z0rough);
     }
-    file1.close();
+    file.close();
+
     m_sim.io_manager().register_output_int_var("terrain_drag");
     m_sim.io_manager().register_output_int_var("terrain_blank");
     m_sim.io_manager().register_io_var("terrainz0");
@@ -63,12 +61,12 @@ void TerrainDrag::post_init_actions()
         auto& terrain_height = m_terrain_height(level);
         auto& drag = m_terrain_drag(level);
         // copy terrain data to gpu
-        amrex::Gpu::DeviceVector<amrex::Real> device_xterrain(
-            m_xterrain.size());
-        amrex::Gpu::DeviceVector<amrex::Real> device_yterrain(
-            m_xterrain.size());
-        amrex::Gpu::DeviceVector<amrex::Real> device_zterrain(
-            m_xterrain.size());
+        const auto xterrain_size = m_xterrain.size();
+        const auto yterrain_size = m_yterrain.size();
+        const auto zterrain_size = m_zterrain.size();
+        amrex::Gpu::DeviceVector<amrex::Real> device_xterrain(xterrain_size);
+        amrex::Gpu::DeviceVector<amrex::Real> device_yterrain(yterrain_size);
+        amrex::Gpu::DeviceVector<amrex::Real> device_zterrain(zterrain_size);
         amrex::Gpu::copy(
             amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(),
             device_xterrain.begin());
@@ -82,9 +80,12 @@ void TerrainDrag::post_init_actions()
         const auto* yterrain_ptr = device_yterrain.data();
         const auto* zterrain_ptr = device_zterrain.data();
         // Copy Roughness to gpu
-        amrex::Gpu::DeviceVector<amrex::Real> device_xrough(m_xrough.size());
-        amrex::Gpu::DeviceVector<amrex::Real> device_yrough(m_xrough.size());
-        amrex::Gpu::DeviceVector<amrex::Real> device_z0rough(m_xrough.size());
+        const auto xrough_size = m_xrough.size();
+        const auto yrough_size = m_yrough.size();
+        const auto z0rough_size = m_z0rough.size();
+        amrex::Gpu::DeviceVector<amrex::Real> device_xrough(xrough_size);
+        amrex::Gpu::DeviceVector<amrex::Real> device_yrough(yrough_size);
+        amrex::Gpu::DeviceVector<amrex::Real> device_z0rough(z0rough_size);
         amrex::Gpu::copy(
             amrex::Gpu::hostToDevice, m_xrough.begin(), m_xrough.end(),
             device_xrough.begin());
@@ -103,67 +104,36 @@ void TerrainDrag::post_init_actions()
             auto levelDrag = drag.array(mfi);
             auto levelz0 = terrainz0.array(mfi);
             auto levelheight = terrain_height.array(mfi);
-            const unsigned terrainSize = m_xterrain.size();
-            const unsigned roughnessSize = m_xrough.size();
             amrex::ParallelFor(
                 vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                    // compute the source term
                     const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
                     const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
                     const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
-                    // Terrain Height
-                    amrex::Real residual = 10000;
-                    amrex::Real terrainHt = 0.0;
-                    for (unsigned ii = 0; ii < terrainSize; ++ii) {
-                        const amrex::Real radius = std::sqrt(
-                            std::pow(x - xterrain_ptr[ii], 2) +
-                            std::pow(y - yterrain_ptr[ii], 2));
-                        if (radius < residual) {
-                            residual = radius;
-                            terrainHt = zterrain_ptr[ii];
-                        }
-                    }
+                    const amrex::Real terrainHt = interp::bilinear(
+                        xterrain_ptr, xterrain_ptr + xterrain_size,
+                        yterrain_ptr, yterrain_ptr + yterrain_size,
+                        zterrain_ptr, x, y);
                     levelBlanking(i, j, k, 0) =
                         static_cast<int>(z <= terrainHt);
                     levelheight(i, j, k, 0) =
                         std::max(std::abs(z - terrainHt), 0.5 * dx[2]);
-                    residual = 10000;
+
                     amrex::Real roughz0 = 0.1;
-                    for (unsigned ii = 0; ii < roughnessSize; ++ii) {
-                        const amrex::Real radius = std::sqrt(
-                            std::pow(x - xrough_ptr[ii], 2) +
-                            std::pow(y - yrough_ptr[ii], 2));
-                        if (radius < residual) {
-                            residual = radius;
-                            roughz0 = z0rough_ptr[ii];
-                        }
-                        if (radius < dx[0]) {
-                            break;
-                        }
+                    if (xrough_size > 0) {
+                        roughz0 = interp::bilinear(
+                            xrough_ptr, xrough_ptr + xrough_size, yrough_ptr,
+                            yrough_ptr + yrough_size, z0rough_ptr, x, y);
                     }
                     levelz0(i, j, k, 0) = roughz0;
                 });
+
             amrex::ParallelFor(
                 vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                    // Terrain Height
-                    amrex::Real residual = 10000;
-                    amrex::Real terrainHt = 0.0;
-                    const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
-                    const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
-                    const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
-                    for (unsigned ii = 0; ii < terrainSize; ++ii) {
-                        const amrex::Real radius = std::sqrt(
-                            std::pow(x - xterrain_ptr[ii], 2) +
-                            std::pow(y - yterrain_ptr[ii], 2));
-                        if (radius < residual) {
-                            residual = radius;
-                            terrainHt = zterrain_ptr[ii];
-                        }
-                    }
-                    levelDrag(i, j, k, 0) = 0;
-                    if (z > terrainHt && k > 0 &&
-                        levelBlanking(i, j, k - 1, 0) == 1) {
+                    if ((levelBlanking(i, j, k, 0) == 0) && (k > 0) &&
+                        (levelBlanking(i, j, k - 1, 0) == 1)) {
                         levelDrag(i, j, k, 0) = 1;
+                    } else {
+                        levelDrag(i, j, k, 0) = 0;
                     }
                 });
         }

amr-wind/utilities/CMakeLists.txt

@@ -6,6 +6,7 @@ target_sources(${amr_wind_lib_name}
       bc_ops.cpp
       console_io.cpp
       index_operations.cpp
+      io_utils.cpp
       IOManager.cpp
       FieldPlaneAveraging.cpp
       FieldPlaneAveragingFine.cpp

amr-wind/utilities/io_utils.H

@@ -74,6 +74,18 @@ inline void assert_with_message(const bool val, const std::string& msg)
         amrex::Abort(msg);
     }
 }
+
+/** Read a flattened 2D grid file
+ *
+ *  File is assumed to contain a single column with the following data:
+ *  [nx, ny, x values (len = nx), y values (len = ny), z values (len = nx * ny)]
+ */
+void read_flat_grid_file(
+    const std::string& fname,
+    amrex::Vector<amrex::Real>& xs,
+    amrex::Vector<amrex::Real>& ys,
+    amrex::Vector<amrex::Real>& zs);
+
 } // namespace amr_wind::ioutils
 
 #endif /* IO_UTILS_H */

amr-wind/utilities/io_utils.cpp

@@ -0,0 +1,39 @@
+#include <fstream>
+
+#include "amr-wind/utilities/io_utils.H"
+
+namespace amr_wind::ioutils {
+
+void read_flat_grid_file(
+    const std::string& fname,
+    amrex::Vector<amrex::Real>& xs,
+    amrex::Vector<amrex::Real>& ys,
+    amrex::Vector<amrex::Real>& zs)
+{
+    std::ifstream file(fname, std::ios::in);
+
+    if (!file.good()) {
+        amrex::Abort("Cannot find file");
+    }
+
+    size_t nx = 0;
+    size_t ny = 0;
+    file >> nx;
+    file >> ny;
+    AMREX_ALWAYS_ASSERT(nx > 0);
+    AMREX_ALWAYS_ASSERT(ny > 0);
+    xs.resize(nx);
+    ys.resize(ny);
+    zs.resize(nx * ny);
+    for (size_t n = 0; n < nx; n++) {
+        file >> xs[n];
+    }
+    for (size_t n = 0; n < ny; n++) {
+        file >> ys[n];
+    }
+    for (size_t n = 0; n < nx * ny; n++) {
+        file >> zs[n];
+    }
+    file.close();
+}
+} // namespace amr_wind::ioutils

amr-wind/utilities/linear_interpolation.H

@@ -270,6 +270,79 @@ inline void linear_angle(
     }
 }
 
+template <typename C1, typename C2>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    typename std::iterator_traits<C2>::value_type
+    bilinear_impl(
+        const C1 xbegin,
+        const C1 ybegin,
+        const int ny,
+        const C2 zinp,
+        const typename std::iterator_traits<C1>::value_type& xout,
+        const typename std::iterator_traits<C1>::value_type& yout,
+        const Index& xidx,
+        const Index& yidx)
+{
+    using DType1 = typename std::iterator_traits<C1>::value_type;
+
+    if ((xidx.lim == Limits::LOWLIM) || (xidx.lim == Limits::UPLIM) ||
+        (yidx.lim == Limits::LOWLIM) || (yidx.lim == Limits::UPLIM)) {
+        return zinp[xidx.idx * ny + yidx.idx];
+    }
+
+    const int i = xidx.idx;
+    const int j = yidx.idx;
+
+    static constexpr DType1 eps = 1.0e-8;
+    const auto xdenom = (xbegin[i + 1] - xbegin[i]);
+    const auto xfacR = (xdenom > eps) ? ((xout - xbegin[i]) / xdenom) : 1.0;
+    const auto xfacL = static_cast<DType1>(1.0) - xfacR;
+    const auto ydenom = (ybegin[j + 1] - ybegin[j]);
+    const auto yfacR = (ydenom > eps) ? ((yout - ybegin[j]) / ydenom) : 1.0;
+    const auto yfacL = static_cast<DType1>(1.0) - yfacR;
+
+    const auto z00 = zinp[i * ny + j];
+    const auto z10 = zinp[(i + 1) * ny + j];
+    const auto z01 = zinp[i * ny + (j + 1)];
+    const auto z11 = zinp[(i + 1) * ny + (j + 1)];
+
+    return z00 * xfacL * yfacL + z10 * xfacR * yfacL + z01 * xfacL * yfacR +
+           z11 * xfacR * yfacR;
+}
+
+template <typename C1, typename C2>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    typename std::iterator_traits<C2>::value_type
+    bilinear(
+        const C1 xbegin,
+        const C1 xend,
+        const C1 ybegin,
+        const C1 yend,
+        const C2 zinp,
+        const typename std::iterator_traits<C1>::value_type& xout,
+        const typename std::iterator_traits<C1>::value_type& yout)
+{
+    const auto xidx = bisection_search(xbegin, xend, xout);
+    const auto yidx = bisection_search(ybegin, yend, yout);
+    return bilinear_impl(
+        xbegin, ybegin, static_cast<int>(yend - ybegin), zinp, xout, yout, xidx,
+        yidx);
+}
+
+template <typename C1, typename C2>
+inline typename C2::value_type bilinear(
+    const C1& xinp,
+    const C1& yinp,
+    const C2& zinp,
+    const typename C1::value_type& xout,
+    const typename C1::value_type& yout)
+{
+    AMREX_ALWAYS_ASSERT((xinp.size() * yinp.size()) == zinp.size());
+    return bilinear(
+        xinp.data(), (xinp.data() + xinp.size()), yinp.data(),
+        (yinp.data() + yinp.size()), zinp.data(), xout, yout);
+}
+
 } // namespace amr_wind::interp
 
 #endif /* LINEAR_INTERPOLATION_H */

amr-wind/utilities/sampling/VolumeSampler.cpp

@@ -60,7 +60,7 @@ void VolumeSampler::check_bounds()
         }
     }
     if (!all_ok) {
-        amrex::Print() << "WARNING: LineSampler: Out of domain line was "
+        amrex::Print() << "WARNING: VolumeSampler: Out of domain line was "
                           "truncated to match domain"
                        << std::endl;
     }

docs/sphinx/user/inputs_incflo.rst

@@ -14,7 +14,7 @@ as initial conditions and discretization options.
    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).
    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 TerrainDrag must be specified and the folder should include a `terrain.amrwind` file. 
+   For immersed boundary forcing method TerrainDrag must be specified and the folder should include a terrain file (default name: `terrain.amrwind`, user control is `TerrainDrag.terrain_file`) file.
 
 .. input_param:: incflo.density
 

docs/sphinx/walkthrough/terrain.rst

@@ -15,10 +15,15 @@ The setup for the terrain follows the typical simulation of the atmospheric boun
 large eddy simulation or Reynolds-averaged Navier Stokes turbulence models. The IBFM can be used with
 periodic or inflow-outflow boundary conditions with few modifications. 
 
-The first step in including the terrain is to set the terrain variables. This is accomplished 
-by modifying the ABL physics to include the ``TerrainDrag`` flow physics: ``incflo.physics = ABL TerrainDrag``.
-This looks for the ``terrain.amrwind`` text file in the case folder. The file contains the terrain height 
-in the format ``x y terrainHt``. 
+The first step in including the terrain is to set the terrain
+variables. This is accomplished by modifying the ABL physics to
+include the ``TerrainDrag`` flow physics: ``incflo.physics = ABL
+TerrainDrag``.  This looks for the ``terrain.amrwind`` text file in
+the case folder (this is the default name, the user can modify the
+file it searches for by specifying the ``TerrainDrag.terrain_file``
+input parameter). The file contains the terrain height as a single
+column organized as: ``nx, ny, x values (of length nx), y values (of
+length ny), terrain height values (of length nx x ny)``.
 
 The second step is the inclusion of the terrain forcing in the momentum and energy equations. This is 
 accomplished by adding ``DragForcing`` and ``DragTempForcing`` terms to ``ICNS.source_terms`` and 

test/test_files/abl_kosovic_neutral_ib/terrain.amrwind

@@ -1,8 +1,17 @@
-1024 1024 0
-1040 1024 0
-1056 1024 0
-1072 1024 0
-1088 1024 0
-1104 1024 0
-1120 1024 0
-
+7
+1
+1024.0
+1040.0
+1056.0
+1072.0
+1088.0
+1104.0
+1120.0
+1024.0
+0.0
+0.0
+0.0
+0.0
+0.0
+0.0
+0.0