Add more tests (#19)

* add tests initial conditions

* update FastscapelibContext

* include newline in test file

* add tests boundary

* add tests tectonics

* remove xsimlab import in tests initial conditions

* apply review comments

---------

Co-authored-by: Benoit Bovy <[email protected]>

fastscape/processes/initial.py

@@ -56,7 +56,7 @@ class Escarpment:
     elevation = xs.foreign(SurfaceTopography, "elevation", intent="out")
 
     def initialize(self):
-        self.elevation = np.full(self.shape, self.elevation_left)
+        self.elevation = np.full(self.shape, self.elevation_left, dtype=np.double)
 
         # align scarp limit locations
         idx_left = np.argmax(self.x > self.x_left)

fastscape/tests/test_processes_boundary.py

@@ -0,0 +1,61 @@
+import numpy as np
+import pytest
+
+from fastscape.processes import BorderBoundary
+
+
+def test_border_boundary_broadcast():
+    p = BorderBoundary(status="fixed_value")
+
+    p.initialize()
+    np.testing.assert_equal(p.border, np.array(["left", "right", "top", "bottom"]))
+    np.testing.assert_equal(p.border_status, ["fixed_value"] * 4)
+
+
+@pytest.mark.parametrize(
+    "status, expected_ibc",
+    [
+        (["fixed_value", "fixed_value", "fixed_value", "fixed_value"], 1111),
+        (["core", "fixed_value", "fixed_value", "fixed_value"], 1110),
+        (["fixed_value", "core", "fixed_value", "fixed_value"], 1011),
+        (["fixed_value", "fixed_value", "core", "fixed_value"], 111),
+        (["fixed_value", "fixed_value", "fixed_value", "core"], 1101),
+        (["looped", "looped", "fixed_value", "fixed_value"], 1010),
+        (["fixed_value", "fixed_value", "looped", "looped"], 101),
+    ],
+)
+def test_border_boundary_ibc(status, expected_ibc):
+    p = BorderBoundary(status=status)
+
+    p.initialize()
+    np.testing.assert_equal(p.ibc, expected_ibc)
+
+
+@pytest.mark.parametrize(
+    "status, error_msg",
+    [
+        (["fixed_value", "fixed_value", "core"], "Border status should be defined for all borders"),
+        ("invalid_status", "Invalid border status"),
+        (["looped", "looped", "core", "core"], "There must be at least one border with status"),
+        (
+            ["looped", "fixed_value", "looped", "fixed_value"],
+            "Periodic boundary conditions must be symmetric",
+        ),
+    ],
+)
+def test_border_boundary_error(status, error_msg):
+    with pytest.raises(ValueError, match=error_msg):
+        BorderBoundary(status=status)
+
+
+@pytest.mark.parametrize(
+    "status, warning_msg",
+    [
+        (["core", "core", "fixed_value", "fixed_value"], "Left and right"),
+        (["fixed_value", "fixed_value", "core", "core"], "Top and bottom"),
+    ],
+)
+def test_border_boundary_warning(status, warning_msg):
+    p = BorderBoundary(status=status)
+    with pytest.warns(UserWarning, match=warning_msg):
+        p.initialize()

fastscape/tests/test_processes_initial.py

@@ -0,0 +1,83 @@
+import numpy as np
+import pytest
+
+from fastscape.processes import (
+    BareRockSurface,
+    Escarpment,
+    FlatSurface,
+    NoErosionHistory,
+)
+
+
+def test_flat_surface():
+    shape = (3, 2)
+    rs = np.random.RandomState(seed=1234)
+    elevation = rs.rand(*shape)
+
+    np.random.default_rng(1234)
+    p = FlatSurface(shape=shape, seed=1234)
+    p.initialize()
+
+    np.testing.assert_equal(shape, p.shape)
+    np.testing.assert_allclose(elevation, p.elevation)
+
+    p2 = FlatSurface(shape=shape, seed=None)
+    p2.initialize()
+    assert np.all(p2.elevation > 0.0)
+    assert np.all(p2.elevation <= 1.0)
+
+
+@pytest.mark.parametrize(
+    "inputs",
+    [
+        (
+            {
+                "x_left": 10,
+                "x_right": 20,
+                "elevation_left": 0.0,
+                "elevation_right": 100.0,
+                "shape": (11, 31),
+                "x": np.linspace(0, 30, 31),
+            }
+        ),
+        (
+            {
+                "x_left": 15,
+                "x_right": 15,
+                "elevation_left": 0.0,
+                "elevation_right": 100.0,
+                "shape": (11, 31),
+                "x": np.linspace(0, 30, 31),
+            }
+        ),
+    ],
+)
+def test_escarpment(inputs):
+    p = Escarpment(**inputs)
+    p.initialize()
+
+    # test invariant elevation along the rows (y-axis) up to random values
+    assert np.all(np.abs(p.elevation - p.elevation[0, :]) < 1.0)
+
+    # shape and x-coordinate values chosen so that the escaprement limits
+    # match the grid
+    assert abs(p.elevation[0, int(p.x_left)] - p.elevation_left) < 1.0
+    assert abs(p.elevation[0, int(p.x_right) + 1] - p.elevation_right) < 1.0
+
+
+def test_bare_rock_surface():
+    elevation = np.array([[2, 3], [4, 1]])
+
+    p = BareRockSurface(surf_elevation=elevation)
+    p.initialize()
+
+    np.testing.assert_equal(elevation, p.bedrock_elevation)
+    # bedrock_elevation must be a copy of surf_elevation
+    assert p.bedrock_elevation.base is not p.surf_elevation
+
+
+def test_no_erosion_history():
+    p = NoErosionHistory()
+    p.initialize()
+
+    assert p.height == 0

fastscape/tests/test_processes_tectonics.py

@@ -0,0 +1,109 @@
+import numpy as np
+import pytest
+
+from fastscape.processes import (
+    BlockUplift,
+    SurfaceAfterTectonics,
+    TectonicForcing,
+    TwoBlocksUplift,
+)
+from fastscape.processes.context import FastscapelibContext
+
+
+def test_tectonic_forcing():
+    grid_shape = (3, 2)
+    uplift = np.full(grid_shape, 1.0)
+    isostasy = np.full(grid_shape, 2.0)
+    bedrock_advect = np.full(grid_shape, 3.0)
+    surf_advect = np.full(grid_shape, 4.0)
+    area = 300.0
+    dt = 10.0
+
+    p = TectonicForcing(
+        bedrock_forcing_vars=[uplift, isostasy, bedrock_advect],
+        surface_forcing_vars=[uplift, isostasy, surf_advect],
+        grid_area=area,
+    )
+
+    p.run_step(dt)
+
+    # uplift + isostasy + bedrock_advect
+    np.testing.assert_equal(p.bedrock_upward, np.full(grid_shape, 1.0 + 2.0 + 3.0))
+
+    # uplift + isostasy + surf_advect
+    np.testing.assert_equal(p.surface_upward, np.full(grid_shape, 1.0 + 2.0 + 4.0))
+
+    # test scalar values
+    p2 = TectonicForcing(surface_forcing_vars=[1.0, 2.0, 3.0], grid_area=area)
+    p2.run_step(dt)
+    assert p2.surface_upward == 6.0
+    assert p2.bedrock_upward == 0.0  # no variables given
+    assert p2._domain_rate() == 6.0 * area / dt
+
+
+def test_surface_after_tectonics():
+    grid_shape = (3, 2)
+    topo_elevation = np.full(grid_shape, 2.0)
+    forced_motion = np.full(grid_shape, 3.0)
+
+    p = SurfaceAfterTectonics(topo_elevation=topo_elevation, forced_motion=forced_motion)
+
+    p.run_step()
+
+    expected = topo_elevation + forced_motion
+    np.testing.assert_equal(p.elevation, expected)
+
+
+@pytest.mark.parametrize(
+    "b_status, expected_uplift",
+    [
+        (
+            np.array(["fixed_value", "fixed_value", "fixed_value", "fixed_value"]),
+            np.array([[0.0, 0.0, 0.0, 0.0], [0.0, 50.0, 50.0, 0.0], [0.0, 0.0, 0.0, 0.0]]),
+        ),
+        (
+            np.array(["fixed_value", "core", "core", "fixed_value"]),
+            np.array([[0.0, 50.0, 50.0, 50.0], [0.0, 50.0, 50.0, 50.0], [0.0, 0.0, 0.0, 0.0]]),
+        ),
+    ],
+)
+def test_block_uplift(b_status, expected_uplift):
+    rate = 5
+    shape = (3, 4)
+    dt = 10.0
+
+    # dummy context
+    f = FastscapelibContext(shape=shape, length=(10.0, 30.0), ibc=1010)
+    f.initialize()
+    f.run_step(dt)
+
+    p = BlockUplift(rate=rate, shape=shape, status=b_status, fs_context=f)
+
+    p.initialize()
+    p.run_step(dt)
+    np.testing.assert_equal(p.uplift, expected_uplift)
+
+    # test variable rate
+    p2 = BlockUplift(rate=np.full(shape, 5.0), shape=shape, status=b_status, fs_context=f)
+    p2.initialize()
+    p2.run_step(dt)
+    np.testing.assert_equal(p2.uplift, expected_uplift)
+
+
+def test_two_blocks_uplift():
+    x = np.array([1, 2, 3])
+    x_pos = 1
+    rate_l = 2
+    rate_r = 3
+    grid = (3, 4)
+    dt = 10.0
+
+    p = TwoBlocksUplift(x_position=x_pos, rate_left=rate_l, rate_right=rate_r, shape=grid, x=x)
+
+    p.initialize()
+    p.run_step(dt)
+
+    expected = np.array(
+        [[20.0, 30.0, 30.0, 30.0], [20.0, 30.0, 30.0, 30.0], [20.0, 30.0, 30.0, 30.0]]
+    )
+    np.testing.assert_equal(p.uplift, expected)