alternative implementations of flooding_correction()

[amoodie]

description

I have examined some alternative implementations of the flooding_correction() method.

The motivation of this method, as best as I can tell, is to ensure there are no "windows" (geology reference that arise by chance of not being in the path of any water parcels, but would be flooded due to "water runs downhill" in a real-world situation.

At present, the method implementation uses a few other helping methods to build neighbor-stage arrays and updates the stage at the center cell with the maximum stage of its neighbors. At face value, this seems reasonable to me.

1. The current implementation is pretty convoluted and difficult to follow. I believe that my implementation below (with _flag set instead to max reimplements the same thing as is currently implemented, but is simpler. This needs to be verified by a consistency check between two runs, or simultaneous computations through a run.
2. The original matlab implementation actually does something subtly different. Due to the order of looping through the neighbors used in this implementation, the stage value updated is not the maximum of the neighbors, but rather the last cell in iteration that is wet and has a stage higher than the center dry cell. This iteration follows an E --> clockwise iteration. I have implemented (I think) the original matlab version below, when the _flag is set to legacy.

Prelim tests

I have run a set of simulations with a) the current implementation, and b) my implementation of the matlab version.

old	new
errored

and if I'm being honest, they all kinda don't look excellent, and I can't really discern any qualitative difference.

As a result, I've decided not to make any change to the flooding_correction in #133.

def flooding_correction(self):
        """Flood dry cells along the shore if necessary.

        Check the neighbors of all dry cells. If any dry cells have wet
        neighbors, check that their stage is not higher than the bed elevation
        of the center cell.
        If it is, flood the dry cell.
        """
        _msg = 'Computing flooding correction'
        self.log_info(_msg, verbosity=2)

        pad_wet_mask = np.copy(np.pad(self.wet_mask, 1, 'edge'))
        pad_stage = np.copy(np.pad(self.stage, 1, 'edge'))

        _ord = np.array([5, 8, 7, 6, 3, 0, 1, 2])
        for i in np.arange(self.L):
            for j in np.arange(self.W):
                if self.wet_mask[i, j] == 0: # locate dry nodes

                    # slice the padded array neighbors
                    wm_nbrs = pad_wet_mask[
                        i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]
                    stage_nbrs = pad_stage[
                        i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]

                    wm_nbrs_flat = wm_nbrs.ravel()
                    stage_nbrs_flat = stage_nbrs.ravel()

                    _flag = 'legacy'
                    if _flag == 'max':
                        wm_nbrs_flat[4] = 0

                        wm_nbrs_stage_flat = stage_nbrs_flat[wm_nbrs_flat]
                        if wm_nbrs_stage_flat.size > 0:
                            max_wm_nbrs_stage = np.max(wm_nbrs_stage_flat)

                            if max_wm_nbrs_stage > self.eta[i, j]:
                                self.stage[i, j] = max_wm_nbrs_stage
                    else:
                        wm_nbrs_ord = wm_nbrs_flat[_ord]
                        stage_nbrs_ord = stage_nbrs_flat[_ord]
                        stage_above_eta = (stage_nbrs_ord > self.eta[i, j])
                        both_conds = np.logical_and(wm_nbrs_ord, stage_above_eta)
                        if np.any(both_conds):
                            idx = np.where(both_conds)[0][-1]
                            self.stage[i, j] = stage_nbrs_ord[idx]

[amoodie]

also of note in this regard, is that the wet_mask is calculated before water iteration in the matlab implementation.

so, add near the top of water iteration:

self.wet_mask = self.depth >= self.dry_depth