Test failures for many derived models in terrainbento.

[alexmitchell]

Describe the bug
When running terrainbento tests, 30 tests for derived models fail. They all seem to be array mismatch issues. I have not looked into it at all, but I suspect that maybe updates or bug fixes in landlab broke these tests. In particular, I suspect the flow routing bug fix in landlab. See attached file for complete test output.

To Reproduce
Steps to reproduce the behavior:

1. Go to the tests directory
2. Run pytest. (the notebooks are fine, so you can skip those excruciatingly slow tests.)

Expected behavior
I would expect no failures.

Environment (please complete the following information):

• OS: Fedora 32
• Python Version 3.7.9
• Landlab Version (developer) 2.3.0.dev0+56.gaf924396
• terrainbento version (developer) v2.0.0b1+22.g50b1730

Additional context
Here is the output from pytest running over all tests. While my version of terrainbento includes the new output writer code, the failures predated the new code.
pytest_all_15_02_2021.txt

[kbarnhart]

@alexmitchell thanks for pointing this out. I'm on it.

[kbarnhart]

@gregtucker I may need some help with you on this one as I think many of these errors has something to do with the DepressionFinderAndRouter.

Namely, the node located at the top of a ridge is getting labeled as a lake. It thus does not erode... And thus get too high, breaking the tests. I suspect that @alexmitchell is correct that it is related to this fix landlab/landlab/pull/1248

# modification of the tests for basicDdHy

import sys
import numpy as np
from numpy.testing import assert_array_almost_equal

import matplotlib.pyplot as plt

from landlab import RasterModelGrid
from terrainbento import Clock
from terrainbento import BasicDdHy, NotCoreNodeBaselevelHandler


U = 0.0001
K = 0.01

grid = RasterModelGrid((3, 21), xy_spacing=100.0)
grid.set_closed_boundaries_at_grid_edges(False, True, False, True)
grid.add_zeros("node", "topographic__elevation")
grid.add_ones("node", "soil__depth")
grid.add_zeros("node", "lithology_contact__elevation")
grid_1 = grid

clock_simple = Clock(step=1000.0, stop=5.1e6)

solver = "basic"
threshold = 0
m_sp = 0.5
n_sp = 1
depression_finder = "DepressionFinderAndRouter"

ncnblh = NotCoreNodeBaselevelHandler(
    grid_1, modify_core_nodes=True, lowering_rate=-U
)
F_f = 0.0
v_sc = 0.001
thresh_change_per_depth = 0
# construct dictionary. note that D is turned off here
params = {
    "grid": grid_1,
    "clock": clock_simple,
    "regolith_transport_parameter": 0.0,
    "water_erodibility": K,
    "settling_velocity": v_sc,
    "fraction_fines": F_f,
    "water_erosion_rule__threshold": threshold,
    "water_erosion_rule__thresh_depth_derivative": thresh_change_per_depth,
    "solver": solver,
    "m_sp": m_sp,
    "n_sp": n_sp,
    "depression_finder": depression_finder,
    "boundary_handlers": {"NotCoreNodeBaselevelHandler": ncnblh},
}

# construct and run model
model = BasicDdHy(**params)
for i in range(500):
  
    model.run_one_step(10)
    lakes = model.flow_accumulator.depression_finder.lake_codes
    if len(lakes)>0:
        print(i, lakes) # lakes at all timesteps. 
        # after few tens, its always the middle node. 
        
model.flow_accumulator.run_one_step()
# construct actual and predicted slopes
actual_slopes = model.grid.at_node["topographic__steepest_slope"]
actual_areas = model.grid.at_node["surface_water__discharge"]
predicted_slopes = np.power(
    ((U * v_sc) / (K * np.power(actual_areas, m_sp)))
    + (U / (K * np.power(actual_areas, m_sp))),
    1.0 / n_sp,
)


sel = model.grid.core_nodes[1:-1]
plt.plot(model.grid.x_of_node[sel], model.z[sel],  '.')
plt.plot(model.grid.x_of_node[lakes], model.z[lakes],  'r.')

# assert actual and predicted slopes are the same.
assert_array_almost_equal(
    actual_slopes[model.grid.core_nodes[1:-1]],
    predicted_slopes[model.grid.core_nodes[1:-1]],
    decimal=4,
)