Skip to content

Boundary condition change conditional for setting tempslope#24

Open
dvalters wants to merge 3 commits intomasterfrom
tempslope-x-minus-one
Open

Boundary condition change conditional for setting tempslope#24
dvalters wants to merge 3 commits intomasterfrom
tempslope-x-minus-one

Conversation

@dvalters
Copy link
Owner

@dvalters dvalters commented Jul 1, 2018

Changes part of flow_route where tempslope is set to the edgeslope parameter. Currently on master tempslope is set where y and x are <=2. Opposite boundary uses imax and jmax.

Now replacing:

            if (x == imax) tempslope = edgeslope;	             
            if (x <= 2) tempslope = 0 - edgeslope;	

with:

            if (x == imax) tempslope = edgeslope;
            if (x == 0) tempslope = edgeslope;

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.

@dvalters
Copy link
Owner Author

dvalters commented Jul 1, 2018

Regression test results

Ok, strictly speaking this breaks the tests as the results are not exactly identical:

chart_compare_x_lt_1

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.

@dvalters dvalters self-assigned this Jul 1, 2018
@dvalters dvalters added the eCSE project Issue arising from the eCSE project label Jul 1, 2018
(elev[x][y] + water_depth[x][y])) / DX;

if (x == imax) tempslope = edgeslope;
if (x <= 2) tempslope = 0 - edgeslope;
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is tempslope negative at the Western domain boundary, and positive at the Eastern domain boundary. (And similarly for N-S)?

This seems odd...

@dvalters
Copy link
Owner Author

dvalters commented Jul 2, 2018

@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.

equal_edgeslope

i.e.

            if (x == imax) tempslope = edgeslope;
            if (x == 0) tempslope = edgeslope;

and

            if (y == jmax) tempslope = edgeslope;
            if (y == 0) tempslope = edgeslope;

@dvalters
Copy link
Owner Author

dvalters commented Jul 2, 2018

@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:

equal_edgeslope_flippeddems

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...

@dvalters dvalters changed the title Boundary condition change conditional for setting temslope Boundary condition change conditional for setting tempslope Jul 2, 2018
@dvalters
Copy link
Owner Author

dvalters commented Jul 4, 2018

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!

@dvalters
Copy link
Owner Author

dvalters commented Jul 4, 2018

Rough sketch of how this looks (at least in my head... 🤔)

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)

img_20180704_174158924_hdr

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

eCSE project Issue arising from the eCSE project

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant