Problem with bc

[jeanbraun]

Hi @benbovy. I just noticed that there might be an issue with the boundary conditions in FastScape. I am using 0.1.0beta3

When I run a simple model with a uniform uplift and this type of bc: 'boundary__status': ['core', 'fixed_value','looped', 'looped'], I get a solution that looks like it should. But if I flip the bc on their side: 'boundary__status': ['looped', 'looped','core', 'fixed_value'], I get something that does not make any sense...

It looks like there is a bug in where the uplift function is applied in two of the corners: The uplift is imposed but no erosion can take place.

Good solution :

Bad solution:

[benbovy]

Thanks for reporting @jeanbraun, I could reproduce this.

The issue is because grid-y boundaries are incorrectly mapped from fastscape to fastscapelib_fortran (top/bottom boundaries are inverted). This could be easily fixed here:

https://github.com/fastscape-lem/fastscape/blob/29c181aafadfc01dc0047aee0b890b6f47bbe915/fastscape/processes/boundary.py#L104-L105

with:

        # different border order
+        self.ibc = sum(arr_bc * np.array([1, 100, 1000, 10]))
-        self.ibc = sum(arr_bc * np.array([1, 100, 10, 1000]))

In the meantime, a quick and dirty workaround would be to set periodic boundaries on the x-axis and then transpose the output dataset to get the desired result:

in_ds = xs.create_setup(
    ...
    input_vars={
        ...
        'boundary__status': ['fixed_value', 'core', 'looped', 'looped'],
    },
)

out_ds = in_ds.xsimlab.run(...).transpose()

On a related note, we should also probably clarify to what exactly correspond the top and bottom borders. Currently in fastscape, the top border corresponds to y = 0 in both array indices and grid (x,y) coordinate labels. However, xarray (and maybe other) plotting functions show y = 0 at the bottom, which is confusing.

[benbovy]

(Note: to have consistent "x" and "y" axes in the workaround, we need to rename the axes after transposing it)

out_ds.transpose().swap_dims({"y": "x_", "x": "y_"}).rename({"x_": "x", "y_": "y"})