From e8177b9ed0e56560a308f9c4318c508140f9c688 Mon Sep 17 00:00:00 2001 From: Benoit Bovy Date: Tue, 21 Aug 2018 13:51:33 +0200 Subject: [PATCH] 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 --- benchmark/CMakeLists.txt | 2 + benchmark/benchmark_bedrock_channel.cpp | 127 ++++++++++++ benchmark/benchmark_hillslope.cpp | 49 +++-- benchmark/benchmark_setup.hpp | 180 ++++++++++++++++++ doc/source/api_cpp.rst | 5 +- include/fastscapelib/bedrock_channel.hpp | 92 ++++++++- .../tests/test_bedrock_channel.py | 18 +- python/src/main.cpp | 11 +- test/test_bedrock_channel.cpp | 73 +++++-- test/test_hillslope.cpp | 2 - 10 files changed, 512 insertions(+), 47 deletions(-) create mode 100644 benchmark/benchmark_bedrock_channel.cpp create mode 100644 benchmark/benchmark_setup.hpp diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 5d667fa5..31a87cbd 100644 --- 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 new file mode 100644 index 00000000..f0491279 --- /dev/null +++ b/benchmark/benchmark_bedrock_channel.cpp @@ -0,0 +1,127 @@ +#include +#include +#include + +#include + +#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 + typename std::enable_if_t::value, K> + get_k_coef(S&& shape) + { + (void) shape; // do not show warning + return 1e-3; + } + + + template + typename std::enable_if_t::value, K> + get_k_coef(S&& shape) + { + using T = typename K::value_type; + + K k_coef_arr = xt::empty(shape); + k_coef_arr.fill(1e-3); + + return k_coef_arr; + } + + + template + void bm_erode_stream_power(benchmark::State& state) + { + auto ns = static_cast(state.range(0)); + + xt::xtensor erosion; + xt::xtensor elevation; + xt::xtensor receivers; + xt::xtensor dist2receivers; + xt::xtensor stack; + xt::xtensor drainage_area; + + if (Nd == 1) + { + double spacing = 300.; + double x0 = 300.; + double length = (ns - 1) * spacing; + + xt::xtensor x = xt::linspace(length + x0, x0, ns); + erosion = xt::zeros_like(x); + elevation = (length + x0 - x) * 1e-4; + receivers = xt::arange(-1, ns - 1); + dist2receivers = xt::full_like(elevation, spacing); + stack = xt::arange(0, ns); + drainage_area = 6.69 * xt::pow(x, 1.67); // Hack's law + } + + else if (Nd == 2) + { + auto s = bms::FastscapeSetupBase(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(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_TEMPLATE(bm_erode_stream_power, double, xt::xtensor, 1) + ->Apply(bms::grid_sizes); + + BENCHMARK_TEMPLATE(bm_erode_stream_power, double, double, 2) + ->Apply(bms::grid_sizes); + + BENCHMARK_TEMPLATE(bm_erode_stream_power, double, xt::xtensor, 2) + ->Apply(bms::grid_sizes); + +} // namespace benchmark_bedrock_channel + +} // namespace fastscapelib diff --git a/benchmark/benchmark_hillslope.cpp b/benchmark/benchmark_hillslope.cpp index 76c96571..0e438000 100644 --- a/benchmark/benchmark_hillslope.cpp +++ b/benchmark/benchmark_hillslope.cpp @@ -1,5 +1,6 @@ #include -#include +#include +#include #include @@ -9,8 +10,11 @@ #include "fastscapelib/hillslope.hpp" +#include "benchmark_setup.hpp" + namespace fs = fastscapelib; +namespace bms = benchmark_setup; namespace fastscapelib @@ -18,10 +22,29 @@ namespace fastscapelib namespace benchmark_hillslope { - enum class KCoefType {scalar, array}; + template + typename std::enable_if_t::value, K> + get_k_coef(S&& shape) + { + (void) shape; // do not show warning + return 1e-3; + } + + + template + typename std::enable_if_t::value, K> + get_k_coef(S&& shape) + { + using T = typename K::value_type; + + K k_coef_arr = xt::empty(shape); + k_coef_arr.fill(1e-3); - template + return k_coef_arr; + } + + template void bm_erode_linear_diffusion(benchmark::State& state) { auto ns = static_cast(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(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_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) + ->Apply(bms::grid_sizes); } // namespace benchmark_hillslope diff --git a/benchmark/benchmark_setup.hpp b/benchmark/benchmark_setup.hpp new file mode 100644 index 00000000..f0da4499 --- /dev/null +++ b/benchmark/benchmark_setup.hpp @@ -0,0 +1,180 @@ +#pragma once + +#include +#include +#include + +#include "xtensor/xmath.hpp" +#include "xtensor/xtensor.hpp" +#include "xtensor/xrandom.hpp" +#include "xtensor/xview.hpp" + + +namespace benchmark_setup +{ + + +/* + * Run benchmarks for different grid sizes. + * + * Assumes a square grid, i.e., the total number of grid points is s^2. + * + * Use this function with benchmark macros, e.g., + * BENCHMARK(...)->Apply(grid_sizes) + */ +template +static void grid_sizes(benchmark::internal::Benchmark* bench) { + std::array sizes {256, 512, 1024, 2048, 4096}; + + for (int s : sizes) + { + bench->Args({s})->Unit(time_unit); + } +} + + +enum class SurfaceType {cone, cone_inv, cone_noise, flat_noise, gauss, custom}; + + +/* + * A class for generating one of the following synthetic topography on + * a ``n`` x ``n`` square grid: + * + * cone + * Cone surface (smooth, regular slope, no depression) + * cone_inv + * Inverted cone surface (a single, big depression) + * cone_noise + * Cone surface with random perturbations so that the surface + * has many depressions of different sizes + * flat_noise + * Nearly flat surface with random perturbations + * (many small depressions) + * gauss + * Gaussian surface (smooth, no depression) + */ +template +class SyntheticTopography +{ +public: + int seed = 0; + + SyntheticTopography(int n) + { + nnodes_side = n; + auto n_ = static_cast(n); + shape = {n_, n_}; + + grid = xt::meshgrid(xt::linspace(-1, 1, n), + xt::linspace(-1, 1, n)); + } + + xt::xtensor get_elevation() + { + xt::random::seed(seed); + + switch (surf_type) { + case SurfaceType::cone: + elevation = get_elevation_cone(); + break; + + case SurfaceType::cone_inv: + elevation = -get_elevation_cone(); + break; + + case SurfaceType::cone_noise: + elevation = (get_elevation_cone() + + xt::random::rand(shape) * 5. / nnodes_side); + break; + + case SurfaceType::flat_noise: + elevation = xt::random::rand(shape); + break; + + case SurfaceType::gauss: + elevation = xt::exp(-(xt::pow(std::get<0>(grid), 2) / 2. + + xt::pow(std::get<1>(grid), 2) / 2.)); + break; + + default: + elevation = xt::zeros(shape); + break; + } + + return elevation; + } + +private: + int nnodes_side; + std::array shape; + std::tuple, xt::xtensor> grid; + xt::xtensor elevation; + + auto get_elevation_cone() + { + return std::sqrt(2.) - xt::sqrt(xt::pow(std::get<0>(grid), 2) + + xt::pow(std::get<1>(grid), 2)); + } + +}; + +/* + * Set fixed boundary conditions for each of the 4 sides of the grid. + */ +template +void set_fixed_boundary_faces(A&& active_nodes) +{ + auto active_nodes_ = xt::view(active_nodes, xt::all(), xt::all()); + active_nodes_ = true; + + const auto shape = active_nodes.shape(); + const std::array rows_idx {0, shape[0] - 1}; + const std::array cols_idx {0, shape[1] - 1}; + + auto row_bounds = xt::view(active_nodes, xt::keep(rows_idx), xt::all()); + row_bounds = false; + + auto col_bounds = xt::view(active_nodes, xt::all(), xt::keep(cols_idx)); + col_bounds = false; +} + + +/* + * A base setup common to various benchmarks. + */ +template +struct FastscapeSetupBase +{ + int nnodes; + double dx = 100.; + double dy = 100.; + xt::xtensor elevation; + xt::xtensor active_nodes; + xt::xtensor receivers; + xt::xtensor dist2receivers; + xt::xtensor ndonors; + xt::xtensor donors; + xt::xtensor stack; + xt::xtensor basins; + + FastscapeSetupBase(int n) + { + auto topo = SyntheticTopography(n); + + elevation = topo.get_elevation(); + + active_nodes = xt::ones(elevation.shape()); + set_fixed_boundary_faces(active_nodes); + + nnodes = n * n; + + receivers = xt::empty({nnodes}); + dist2receivers = xt::empty({nnodes}); + ndonors = xt::empty({nnodes}); + donors = xt::empty({nnodes, 8}); + stack = xt::empty({nnodes}); + basins = xt::empty({nnodes}); + } +}; + +} // namespace benchmark_setup diff --git a/doc/source/api_cpp.rst b/doc/source/api_cpp.rst index 3b3e3688..b591bc52 100644 --- a/doc/source/api_cpp.rst +++ b/doc/source/api_cpp.rst @@ -94,7 +94,10 @@ Functions used to drive the evolution of bedrock channels. Defined in ``fastscapelib/bedrock_channel.hpp``. -.. doxygenfunction:: fastscapelib::erode_stream_power +.. doxygenfunction:: fastscapelib::erode_stream_power(xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, double, double, double, double, double) + :project: fastscapelib + +.. doxygenfunction:: fastscapelib::erode_stream_power(xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, const xtensor_t&, double, double, double, double) :project: fastscapelib Hillslope diff --git a/include/fastscapelib/bedrock_channel.hpp b/include/fastscapelib/bedrock_channel.hpp index cab8394b..24949b7e 100644 --- a/include/fastscapelib/bedrock_channel.hpp +++ b/include/fastscapelib/bedrock_channel.hpp @@ -3,8 +3,10 @@ */ #pragma once +#include #include #include +#include #include "xtensor/xtensor.hpp" #include "xtensor/xstrided_view.hpp" @@ -19,6 +21,29 @@ namespace detail { +template +auto k_coef_as_array(T k_coef, + S&& shape, + typename std::enable_if_t::value>* = 0) +{ + return xt::broadcast(std::forward(k_coef), shape); +} + + +template +auto k_coef_as_array(K&& k_coef, + S&& shape, + typename std::enable_if_t::value>* = 0) +{ + auto k_coef_arr = xt::flatten(k_coef); + + assert(k_coef_arr.shape() == shape); + (void) shape; // TODO: still unused parameter warning despite assert? + + return k_coef_arr; +} + + /** * erode_stream_power implementation. * @@ -28,14 +53,14 @@ namespace detail * between a node and its receiver, rather than on the node's * elevation itself. This allows saving some operations. */ -template +template index_t erode_stream_power_impl(Er&& erosion, El&& elevation, S&& stack, R&& receivers, Di&& dist2receivers, Dr&& drainage_area, - double k_coef, + K&& k_coef, double m_exp, double n_exp, double dt, @@ -48,6 +73,8 @@ index_t erode_stream_power_impl(Er&& erosion, const auto elevation_flat = xt::flatten(elevation); const auto drainage_area_flat = xt::flatten(drainage_area); + const auto k_coef_arr = k_coef_as_array(k_coef, elevation_flat.shape()); + index_t n_corr = 0; for (const auto& istack : stack) @@ -71,7 +98,8 @@ index_t erode_stream_power_impl(Er&& erosion, continue; } - auto factor = k_coef * dt * std::pow(drainage_area_flat(istack), m_exp); + auto factor = (k_coef_arr(istack) * dt * + std::pow(drainage_area_flat(istack), m_exp)); T delta_0 = istack_elevation - irec_elevation; T delta_k; @@ -200,4 +228,62 @@ index_t erode_stream_power(xtensor_t& erosion, dt, tolerance); } + +/** + * Compute bedrock channel erosion during a single time step using the + * Stream Power Law. + * + * This version accepts a spatially variable stream-power law coefficient. + * + * @param erosion : ``[intent=out, shape=(nrows, ncols)||(nnodes)]`` + * Erosion at grid node. + * @param elevation : ``[intent=in, shape=(nrows, ncols)||(nnodes)]`` + * Elevation at grid node. + * @param stack :``[intent=in, shape=(nnodes)]`` + * Stack position at grid node. + * @param receivers : ``[intent=in, shape=(nnodes)]`` + * Index of flow receiver at grid node. + * @param dist2receivers : ``[intent=out, shape=(nnodes)]`` + * Distance to receiver at grid node. + * @param drainage_area : ``[intent=out, shape=(nrows, ncols)||(nnodes)]`` + * Drainage area at grid node. + * @param k_coef : ``[intent=in, shape=(nrows, ncols)||(nnodes)]`` + * Stream Power Law coefficient. + * @param m_exp : ``[intent=in]`` + * Stream Power Law drainage area exponent. + * @param n_exp : ``[intent=in]`` + * Stream Power Law slope exponent. + * @param dt : ``[intent=in]`` + * Time step duration. + * @param tolerance : ``[intent=in]`` + * Tolerance used for Newton's iterations (``n_exp`` != 1). + * + * @returns + * Total number of nodes for which erosion has been + * arbitrarily limited to ensure consistency. + */ +template +index_t erode_stream_power(xtensor_t& erosion, + const xtensor_t& elevation, + const xtensor_t& stack, + const xtensor_t& receivers, + const xtensor_t& dist2receivers, + const xtensor_t& drainage_area, + const xtensor_t& k_coef, + double m_exp, + double n_exp, + double dt, + double tolerance) +{ + return detail::erode_stream_power_impl(erosion.derived_cast(), + elevation.derived_cast(), + stack.derived_cast(), + receivers.derived_cast(), + dist2receivers.derived_cast(), + drainage_area.derived_cast(), + k_coef.derived_cast(), + m_exp, n_exp, + dt, tolerance); +} + } // namespace fastscapelib diff --git a/python/fastscapelib/tests/test_bedrock_channel.py b/python/fastscapelib/tests/test_bedrock_channel.py index 74734574..e9adfa17 100644 --- a/python/fastscapelib/tests/test_bedrock_channel.py +++ b/python/fastscapelib/tests/test_bedrock_channel.py @@ -1,9 +1,11 @@ +import pytest import numpy as np import fastscapelib -def test_erode_stream_power(): +@pytest.mark.parametrize("k_coef_type", ["constant", "variable"]) +def test_erode_stream_power(k_coef_type): # Test on a tiny (2x2) 2-d square grid with a planar surface # tilted in y (rows) and with all outlets on the 1st row. spacing = 300. @@ -27,10 +29,16 @@ def test_erode_stream_power(): dt = 1 # use small time step (compare with explicit scheme) tolerance = 1e-3 - n_corr = fastscapelib.erode_stream_power_d( - erosion, elevation, stack, - receivers, dist2receivers, drainage_area, - k_coef, m_exp, n_exp, dt, tolerance) + if k_coef_type == 'constant': + func = fastscapelib.erode_stream_power_d + k = k_coef + elif k_coef_type == 'variable': + func = fastscapelib.erode_stream_power_var_d + k = np.full_like(elevation, k_coef) + + n_corr = func(erosion, elevation, stack, + receivers, dist2receivers, drainage_area, + k, m_exp, n_exp, dt, tolerance) slope = h / spacing err = dt * k_coef * a**m_exp * slope**n_exp diff --git a/python/src/main.cpp b/python/src/main.cpp index 427d1851..2d308efc 100644 --- a/python/src/main.cpp +++ b/python/src/main.cpp @@ -119,14 +119,14 @@ void compute_drainage_area_grid_py(xt::pytensor& drainage_area, } -template +template index_t erode_stream_power_py(xt::pyarray& erosion, const xt::pyarray& elevation, const xt::pytensor& stack, const xt::pytensor& receivers, const xt::pytensor& dist2receivers, const xt::pyarray& drainage_area, - double k_coef, + const K k_coef, double m_exp, double n_exp, double dt, @@ -194,10 +194,15 @@ PYBIND11_MODULE(_fastscapelib_py, m) m.def("fill_sinks_sloped_d", &fill_sinks_sloped_py, "Fill depressions in elevation data (slightly sloped surfaces)."); - m.def("erode_stream_power_d", &erode_stream_power_py, + m.def("erode_stream_power_d", &erode_stream_power_py, "Compute bedrock channel erosion during a single time step " "using the Stream Power Law."); + m.def("erode_stream_power_var_d", &erode_stream_power_py&, double>, + "Compute bedrock channel erosion during a single time step " + "using the Stream Power Law.\n\n" + "Version with spatially variable stream power coefficient."); + m.def("erode_linear_diffusion_d", &erode_linear_diffusion_py, "Compute hillslope erosion by linear diffusion on a 2-d regular " "grid using finite differences with an Alternating Direction" diff --git a/test/test_bedrock_channel.cpp b/test/test_bedrock_channel.cpp index ab2d2f05..6a92fdf3 100644 --- a/test/test_bedrock_channel.cpp +++ b/test/test_bedrock_channel.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "gtest/gtest.h" #include "xtensor/xtensor.hpp" @@ -21,19 +22,23 @@ namespace fs = fastscapelib; * result of uplift vs. erosion using the Stream-Power law. * * Drainage area along the profile is estimated using Hack's law. + * + * The template parameter is used to set if the stream-power + * coefficient ``k_coef`` is a scalar or an array. */ +template struct ChannelProfile1D { - double k_coef = 1e-3; + int nnodes = 101; + double spacing = 300.; + double x0 = 300.; + + K k_coef; double m_exp = 0.5; double n_exp = 1; double u_rate = 1e-3; - int nnodes = 101; - double spacing = 300.; - double x0 = 300.; - double hack_coef = 6.69; double hack_exp = 1.67; @@ -41,6 +46,8 @@ struct ChannelProfile1D ChannelProfile1D() { + init_k_coef(); + // fixed boundary (outlet) receivers(0) = 0; dist2receivers(0) = 0.; @@ -51,12 +58,14 @@ struct ChannelProfile1D xt::xtensor get_steady_slope_analytical(void); private: + double k_coef_scalar = 1e-3; + double length = (nnodes - 1) * spacing; xt::xtensor x = xt::linspace(length + x0, x0, static_cast(nnodes)); xt::xtensor drainage_area = hack_coef * xt::pow(x, hack_exp); - xt::xtensor elevation = (x - (length + x0)) * 0.001; + xt::xtensor elevation = (length + x0 - x) * 1e-4; xt::xtensor receivers = xt::arange(-1, nnodes - 1); xt::xtensor dist2receivers = xt::ones({nnodes}) * spacing; @@ -64,6 +73,18 @@ struct ChannelProfile1D xt::xtensor uplift = xt::ones({nnodes}) * u_rate; xt::xtensor erosion = xt::zeros({nnodes}); + + template + void init_k_coef(typename std::enable_if_t::value>* = 0) + { + k_coef = xt::full_like(elevation, k_coef_scalar); + } + + template + void init_k_coef(typename std::enable_if_t::value>* = 0) + { + k_coef = k_coef_scalar; + } }; @@ -71,8 +92,9 @@ struct ChannelProfile1D * Get slope along the channel profile from the numerical solution at * steady state. */ -xt::xtensor ChannelProfile1D::get_steady_slope_numerical(double dt, - double tolerance) +template +xt::xtensor ChannelProfile1D::get_steady_slope_numerical(double dt, + double tolerance) { // run model until steady state is reached double elevation_diff = std::numeric_limits::max(); @@ -103,14 +125,18 @@ xt::xtensor ChannelProfile1D::get_steady_slope_numerical(double dt, /** * Get slope along the channel profile from the analytical solution at - * steady state. + * steady state + * + * Note: this should be used to compare with the numerical solution only if + * ``k_coef`` is invariant in space! */ -xt::xtensor ChannelProfile1D::get_steady_slope_analytical(void) +template +xt::xtensor ChannelProfile1D::get_steady_slope_analytical(void) { // exclude 1st node for proper comparison with the numerical solution auto drainage_area_ = xt::view(drainage_area, xt::range(1, _)); - xt::xtensor slope = (std::pow(u_rate / k_coef, 1. / n_exp) * + xt::xtensor slope = (std::pow(u_rate / k_coef_scalar, 1. / n_exp) * xt::pow(drainage_area_, -m_exp / n_exp)); return slope; @@ -127,7 +153,7 @@ TEST(bedrock_channel, erode_stream_power) "along a 1-d profile at steady-state for n_exp = " + std::to_string(n_exp)); - auto profile_1d = ChannelProfile1D(); + auto profile_1d = ChannelProfile1D(); profile_1d.n_exp = n_exp; auto slope_n = profile_1d.get_steady_slope_numerical(1e4, 1e-3); @@ -137,10 +163,23 @@ TEST(bedrock_channel, erode_stream_power) EXPECT_TRUE(profile_1d.n_corr == 0); } + { + SCOPED_TRACE("test against analytical solution of slope " + "using a 1-d array for k_coef"); + + auto profile_1d = ChannelProfile1D>(); + + auto slope_n = profile_1d.get_steady_slope_numerical(1e4, 1e-3); + auto slope_a = profile_1d.get_steady_slope_analytical(); + + EXPECT_TRUE(xt::allclose(slope_n, slope_a, 1e-5, 1e-5)); + EXPECT_TRUE(profile_1d.n_corr == 0); + } + { SCOPED_TRACE("test arbitrary limitation of erosion"); - auto profile_1d = ChannelProfile1D(); + auto profile_1d = ChannelProfile1D(); profile_1d.n_exp = 0.8; auto slope_n = profile_1d.get_steady_slope_numerical(1e4, 1e-3); @@ -150,7 +189,8 @@ TEST(bedrock_channel, erode_stream_power) { SCOPED_TRACE("test on a tiny (2x2) 2-d square grid " - "with a planar surface tilted in y (rows)"); + "with a planar surface tilted in y (rows) " + "and set k_coef with a 2-d array"); double k_coef = 1e-3; double m_exp = 0.5; @@ -170,18 +210,21 @@ TEST(bedrock_channel, erode_stream_power) xt::xtensor erosion = xt::empty({2, 2}); + auto k_coef_arr = xt::full_like(elevation, k_coef); + double dt = 1.; // use small time step (compare with explicit scheme) double tolerance = 1e-3; index_t n_corr = fs::erode_stream_power( erosion, elevation, stack, receivers, dist2receivers, drainage_area, - k_coef, m_exp, n_exp, dt, tolerance); + k_coef_arr, m_exp, n_exp, dt, tolerance); double slope = h / dy; double err = dt * k_coef * std::pow(a, m_exp) * std::pow(slope, n_exp); xt::xtensor expected_erosion = {{0., 0.}, {err, err}}; EXPECT_TRUE(xt::allclose(erosion, expected_erosion, 1e-5, 1e-5)); + EXPECT_TRUE(n_corr == 0); } } diff --git a/test/test_hillslope.cpp b/test/test_hillslope.cpp index 10559c37..214d4053 100644 --- a/test/test_hillslope.cpp +++ b/test/test_hillslope.cpp @@ -7,8 +7,6 @@ #include "fastscapelib/hillslope.hpp" -#include - namespace fs = fastscapelib;