alternative implementation of updating water discharge

[amoodie]

I began work on an alternative implementation for the water discharge updating. My thought was that it would be helpful to have the flexibility to only update water discharge at the end of water routing. This would enable particles that were unsatisfactory to be re-routed, without having to mess with un-doing their impact on the water discharge array. I think that this approach of rerouting might be useful for the ocassional issue where 1) the delta digs itself a hole, 2) no new water parcels reach the domain edge, 3) thus no new parcels are counted in the free surface, 4) the free surface decays to the sea level, 5) everything blows up.

Ultimately, I decided not to pursue this any further, but I wanted to record some of my efforts, and have made a branch that begins to implement this ability, in case we want to go down that road in the future. HERE is a link to the feature branch.

@njit
def _accumulate_Qs(free_surf_walk_inds, directions, cell_type,
                   Qp_water, dx, iwalk_flat, jwalk_flat):
    """Update discharge field values after one set of water parcel steps."""

    # nparcels = free_surf_walk_inds.shape[0]
    domain_shape = cell_type.shape
    L, W = domain_shape
    cell_type_flat = cell_type.ravel()

    qxn = np.zeros((L * W))
    qyn = np.zeros((L * W))
    qwn = np.zeros((L * W))

    idx = 0
    current_inds = free_surf_walk_inds[:, idx]

    qxn[current_inds] += 1
    qyn[current_inds] += 0  # this could be omitted...
    qwn[current_inds] += Qp_water / dx / 2

    while np.sum(current_inds) > 0:

        # extract next inds *as recorded*
        next_inds = free_surf_walk_inds[:, idx + 1]

        # get the direction of the step between current_inds and next_inds
        steps_direction = directions[:, idx + 1]

        # get the stepped-ness
        dist, istep, jstep, astep = shared_tools.get_steps(
            steps_direction,
            iwalk_flat,
            jwalk_flat)

        # get the potential new indices
        #     (if no step correction were applied)
        potential_inds = _calculate_new_ind(
            current_inds,
            steps_direction,
            iwalk_flat,
            jwalk_flat,
            domain_shape)

        # any cells that will end the timestep at the 0 index,
        #     should not be counted
        astep[next_inds == 0] = 0

        # update the current and potential steps
        qxn = _update_dirQfield(
            qxn[:], dist, current_inds,
            astep, jstep)
        qyn = _update_dirQfield(
            qyn[:], dist, current_inds,
            astep, istep)
        qwn = _update_absQfield(
            qwn[:], dist, current_inds,
            astep, Qp_water, dx)

        qxn = _update_dirQfield(
            qxn[:], dist, potential_inds,
            astep, jstep)
        qyn = _update_dirQfield(
            qyn[:], dist, potential_inds,
            astep, istep)
        qwn = _update_absQfield(
            qwn[:], dist, potential_inds,
            astep, Qp_water, dx)

        # if potential was edge, balance flux at edge
        whr_edge = cell_type_flat[potential_inds] == -1
        qxn = _update_dirQfield(
            qxn[:], dist[whr_edge], potential_inds[whr_edge],
            astep[whr_edge], jstep[whr_edge])
        qyn = _update_dirQfield(
            qyn[:], dist[whr_edge], potential_inds[whr_edge],
            astep[whr_edge], istep[whr_edge])
        qwn = _update_absQfield(
            qwn[:], dist[whr_edge], potential_inds[whr_edge],
            astep[whr_edge], Qp_water, dx)

        # step to the next location
        current_inds = free_surf_walk_inds[:, idx + 1]
        idx += 1

    return qxn, qyn, qwn