Support spatially variable stream-power-law coefficient (#27)

* update erode_stream_power: k_coef accepts a scalar or (1d, 2d) array * update tests * update doc * update Python bindings * update Python tests * add benchmarks for erode_stream_power * fix benchmark diffusion ADI with variable k_coef * fix test stream power 1-d channel erosion (fix initial profile) * fix unused parameter warning * remove unused includes + add needed ones
fastscape-lem · Aug 21, 2018 · e8177b9 · e8177b9
1 parent d977000
commit e8177b9
Show file tree

Hide file tree

Showing 10 changed files with 512 additions and 47 deletions.
diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
@@ -112,6 +112,8 @@ endif()
 include_directories(${GBENCHMARK_INCLUDE_DIRS})
 
 set(FASTSCAPELIB_BENCHMARK_SRC
+  benchmark_setup.hpp
+  benchmark_bedrock_channel.cpp
   benchmark_hillslope.cpp
   benchmark_sinks.cpp
   main.cpp

diff --git a/benchmark/benchmark_bedrock_channel.cpp b/benchmark/benchmark_bedrock_channel.cpp
@@ -0,0 +1,127 @@
+#include <array>
+#include <cmath>
+#include <type_traits>
+
+#include <benchmark/benchmark.h>
+
+#include "xtensor/xbuilder.hpp"
+#include "xtensor/xtensor.hpp"
+#include "xtensor/xrandom.hpp"
+
+#include "fastscapelib/flow_routing.hpp"
+#include "fastscapelib/bedrock_channel.hpp"
+
+#include "benchmark_setup.hpp"
+
+
+namespace fs = fastscapelib;
+namespace bms = benchmark_setup;
+
+
+namespace fastscapelib
+{
+
+namespace benchmark_bedrock_channel
+{
+
+    template<class K, class S>
+    typename std::enable_if_t<std::is_floating_point<K>::value, K>
+    get_k_coef(S&& shape)
+    {
+        (void) shape;  // do not show warning
+        return 1e-3;
+    }
+
+
+    template<class K, class S>
+    typename std::enable_if_t<xt::is_xexpression<K>::value, K>
+    get_k_coef(S&& shape)
+    {
+        using T = typename K::value_type;
+
+        K k_coef_arr = xt::empty<T>(shape);
+        k_coef_arr.fill(1e-3);
+
+        return k_coef_arr;
+    }
+
+
+    template<class T, class K, int Nd>
+    void bm_erode_stream_power(benchmark::State& state)
+    {
+        auto ns = static_cast<size_t>(state.range(0));
+
+        xt::xtensor<double, Nd> erosion;
+        xt::xtensor<double, Nd> elevation;
+        xt::xtensor<index_t, 1> receivers;
+        xt::xtensor<double, 1> dist2receivers;
+        xt::xtensor<index_t, 1> stack;
+        xt::xtensor<double, Nd> drainage_area;
+
+        if (Nd == 1)
+        {
+            double spacing = 300.;
+            double x0 = 300.;
+            double length = (ns - 1) * spacing;
+
+            xt::xtensor<double, Nd> x = xt::linspace<double>(length + x0, x0, ns);
+            erosion = xt::zeros_like(x);
+            elevation = (length + x0 - x) * 1e-4;
+            receivers = xt::arange<index_t>(-1, ns - 1);
+            dist2receivers = xt::full_like(elevation, spacing);
+            stack = xt::arange<index_t>(0, ns);
+            drainage_area = 6.69 * xt::pow(x, 1.67);   // Hack's law
+        }
+
+        else if (Nd == 2)
+        {
+            auto s = bms::FastscapeSetupBase<bms::SurfaceType::cone, T>(state.range(0));
+
+            fs::compute_receivers_d8(s.receivers, s.dist2receivers,
+                                     s.elevation, s.active_nodes,
+                                     s.dx, s.dy);
+            fs::compute_donors(s.ndonors, s.donors, s.receivers);
+            fs::compute_stack(s.stack, s.ndonors,
+                              s.donors, s.receivers);
+
+            erosion = xt::zeros_like(s.elevation);
+            elevation = s.elevation;
+            receivers = s.receivers;
+            dist2receivers = s.dist2receivers;
+            stack = s.stack;
+            drainage_area = xt::empty_like(s.elevation);
+            fs::compute_drainage_area(drainage_area, stack, receivers, s.dx, s.dy);
+        }
+
+        K k_coef = get_k_coef<K>(elevation.shape());
+        double m_exp = 0.5;
+        double n_exp = 1;
+
+        double dt = 1e4;
+
+        index_t ncorr;
+
+        for (auto _ : state)
+        {
+            ncorr = fs::erode_stream_power(erosion, elevation, stack,
+                                           receivers, dist2receivers, drainage_area,
+                                           k_coef, m_exp, n_exp, dt, 1e-3);
+        }
+    }
+
+
+    BENCHMARK_TEMPLATE(bm_erode_stream_power, double, double, 1)
+    ->Apply(bms::grid_sizes<benchmark::kMicrosecond>);
+
+    BENCHMARK_TEMPLATE(bm_erode_stream_power, double, xt::xtensor<double, 1>, 1)
+    ->Apply(bms::grid_sizes<benchmark::kMicrosecond>);
+
+    BENCHMARK_TEMPLATE(bm_erode_stream_power, double, double, 2)
+    ->Apply(bms::grid_sizes<benchmark::kMillisecond>);
+
+    BENCHMARK_TEMPLATE(bm_erode_stream_power, double, xt::xtensor<double, 2>, 2)
+    ->Apply(bms::grid_sizes<benchmark::kMillisecond>);
+
+}  // namespace benchmark_bedrock_channel
+
+}  // namespace fastscapelib
diff --git a/benchmark/benchmark_hillslope.cpp b/benchmark/benchmark_hillslope.cpp
@@ -1,5 +1,6 @@
 #include <array>
-#include <math.h>
+#include <cmath>
+#include <type_traits>
 
 #include <benchmark/benchmark.h>
 
@@ -9,19 +10,41 @@
 
 #include "fastscapelib/hillslope.hpp"
 
+#include "benchmark_setup.hpp"
+
 
 namespace fs = fastscapelib;
+namespace bms = benchmark_setup;
 
 
 namespace fastscapelib
 {
 
 namespace benchmark_hillslope
 {
-    enum class KCoefType {scalar, array};
 
+    template<class K, class S>
+    typename std::enable_if_t<std::is_floating_point<K>::value, K>
+    get_k_coef(S&& shape)
+    {
+        (void) shape;  // do not show warning
+        return 1e-3;
+    }
+
+
+    template<class K, class S>
+    typename std::enable_if_t<xt::is_xexpression<K>::value, K>
+    get_k_coef(S&& shape)
+    {
+        using T = typename K::value_type;
+
+        K k_coef_arr = xt::empty<T>(shape);
+        k_coef_arr.fill(1e-3);
 
-    template<class T, KCoefType k_coef_type>
+        return k_coef_arr;
+    }
+
+    template<class T, class K>
     void bm_erode_linear_diffusion(benchmark::State& state)
     {
         auto ns = static_cast<size_t>(state.range(0));
@@ -32,17 +55,9 @@ namespace benchmark_hillslope
 
         double dx = 0.4;
         double dy = 0.4;
-        double k_coef = 1e-3;
         double dt = 1e4;
 
-        if (k_coef_type == KCoefType::scalar)
-        {
-            double k_coef = 1e-3;
-        }
-        else if (k_coef_type == KCoefType::array)
-        {
-            auto k_coef = xt::full_like(elevation, 1e-3);
-        }
+        K k_coef = get_k_coef<K>(elevation.shape());
 
         for (auto _ : state)
         {
@@ -52,13 +67,11 @@ namespace benchmark_hillslope
     }
 
 
-    BENCHMARK_TEMPLATE(bm_erode_linear_diffusion, double, KCoefType::scalar)
-    ->Arg(256)->Arg(512)->Arg(1024)->Arg(2048)->Arg(4096)
-    ->Unit(benchmark::kMillisecond);
+    BENCHMARK_TEMPLATE(bm_erode_linear_diffusion, double, double)
+    ->Apply(bms::grid_sizes<benchmark::kMillisecond>);
 
-    BENCHMARK_TEMPLATE(bm_erode_linear_diffusion, double, KCoefType::array)
-    ->Arg(256)->Arg(512)->Arg(1024)->Arg(2048)->Arg(4096)
-    ->Unit(benchmark::kMillisecond);
+    BENCHMARK_TEMPLATE(bm_erode_linear_diffusion, double, xt::xtensor<double, 2>)
+    ->Apply(bms::grid_sizes<benchmark::kMillisecond>);
 
 }  // namespace benchmark_hillslope