Open dvalters opened 6 years ago
Ok, strictly speaking this breaks the tests as the results are not exactly identical:
However the results are broadly comparable. There is some slight timing difference at certain points, but this might be expected since we have modified where tempslope is set for 2 grid cell on the West and Northern edges of the model domain.
@aproeme Setting edgeslope/tempslope at the boundary edges in addition to the fix above does not seem to produce wildly different results, and seems more logically consistent.
i.e.
if (x == imax) tempslope = edgeslope;
if (x == 0) tempslope = edgeslope;
and
if (y == jmax) tempslope = edgeslope;
if (y == 0) tempslope = edgeslope;
@aproeme I flipped the DEM in this one so that the catchment is on the Eastern edge. Tempslope is positive and equal to edgeslope now on all the boundaries:
There are slightly higher peak discharges when edgeslope is positive on the E edge, but I'm wondering if this is actually more 'correct' since having a negative edgeslope on the boundary edge would act like a little 'speedbump' slowing the water down slightly...
Will hold off on merging this into master this until we reconcile comments from Tom:
...is there a reason for it being x<=2 rather than x==0? (It seems like should apply only to the leftmost column of cells, mirroring the opposite side of the domain (x == xmax)
Its because the flow routing is calculated via looking at the +ve or -ve water surface slope between (x,y) and (x-1,y) and between (x,y) and (x,y-1). The DEM runs from 1 (not 0) to Xmax so to work out what goes onto row 1 (the edge of the DEM where water is effectively removed) you need to work on row2. For Xmax you look at xmax as it looks at fluxes to the left (x-1) and above (y-1).
Secondly, is there a reason for tempslope being negative where x<=2 (and y<=2) yet positive on the opposing edges?
See above. As it looks at the fluxes both ways between (eg) (x,y) and (x-1,y) then on the LH edge the slope is -ve and on the RH side its +ve.
Maybe there is a way to get around this for the libgeodecomp port... will think about it!
Rough sketch of how this looks (at least in my head... :thinking:)
Drew this (confusingly) using matrix row/column notation even though x/y implies cartesian...sorry! - i.e. [0][0] is the top left corner here (NW corner)
Changes part of
flow_route
wheretempslope
is set to theedgeslope
parameter. Currently onmaster
tempslope is set wherey
andx
are<=2
. Opposite boundary usesimax
andjmax
.Now replacing:
with:
And same for the
y
directions. Also removing negative edgeslope at opposite ends of the model domain.Corresponds to jmax/imax at other boundary conditions of model domain.