from pyro.burgers import burgers_interface



[docs]
def mac_vels(grid,  dt,
             u, v,
             ldelta_ux, ldelta_vx,
             ldelta_uy, ldelta_vy,
             gradp_x, gradp_y,
             source_x=None, source_y=None):
    r"""
    Calculate the MAC velocities in the x and y directions.

    Parameters
    ----------
    grid : Grid2d
        The grid object
    dt : float
        The timestep
    u, v : ndarray
        x-velocity and y-velocity
    ldelta_ux, ldelta_uy: ndarray
        Limited slopes of the x-velocity in the x and y directions
    ldelta_vx, ldelta_vy: ndarray
        Limited slopes of the y-velocity in the x and y directions
    gradp_x, gradp_y : ndarray
        Pressure gradients in the x and y directions
    source_x, source_y : ndarray
        Any other source terms

    Returns
    -------
    out : ndarray, ndarray
        MAC velocities in the x and y directions
    """

    # get the full u and v left and right states (including transverse
    # terms) on both the x- and y-interfaces

    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                burgers_interface.get_interface_states(grid, dt,
                                                       u, v,
                                                       ldelta_ux, ldelta_vx,
                                                       ldelta_uy, ldelta_vy)

    # apply transverse terms on both x- and y- interfaces
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                burgers_interface.apply_transverse_corrections(grid, dt,
                                                               u_xl, u_xr,
                                                               u_yl, u_yr,
                                                               v_xl, v_xr,
                                                               v_yl, v_yr)

    # apply pressure gradient correction terms to the interface state
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                apply_gradp_corrections(dt,
                                        u_xl, u_xr,
                                        u_yl, u_yr,
                                        v_xl, v_xr,
                                        v_yl, v_yr,
                                        gradp_x, gradp_y)

    # apply source terms
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                apply_other_source_terms(dt,
                                        u_xl, u_xr,
                                        u_yl, u_yr,
                                        v_xl, v_xr,
                                        v_yl, v_yr,
                                        source_x, source_y)

    # Riemann problem -- this follows Burger's equation.  We don't use
    # any input velocity for the upwinding.  Also, we only care about
    # the normal states here (u on x and v on y)
    u_MAC = burgers_interface.riemann_and_upwind(grid, u_xl, u_xr)
    v_MAC = burgers_interface.riemann_and_upwind(grid, v_yl, v_yr)

    return u_MAC, v_MAC




[docs]
def states(grid, dt,
           u, v,
           ldelta_ux, ldelta_vx,
           ldelta_uy, ldelta_vy,
           gradp_x, gradp_y,
           u_MAC, v_MAC,
           source_x=None, source_y=None):
    r"""
    This is similar to ``mac_vels``, but it predicts the interface states
    of both u and v on both interfaces, using the MAC velocities to
    do the upwinding.

    Parameters
    ----------
    grid : Grid2d
        The grid object
    dt : float
        The timestep
    u, v : ndarray
        x-velocity and y-velocity
    ldelta_ux, ldelta_uy: ndarray
        Limited slopes of the x-velocity in the x and y directions
    ldelta_vx, ldelta_vy: ndarray
        Limited slopes of the y-velocity in the x and y directions
    gradp_x, gradp_y : ndarray
        Pressure gradients in the x and y directions
    u_MAC, v_MAC : ndarray
        MAC velocities in the x and y directions
    source_x, source_y : ndarray
        Any other source terms

    Returns
    -------
    out : ndarray, ndarray, ndarray, ndarray
        x and y velocities predicted to the interfaces
    """

    # get the full u and v left and right states without transverse terms and gradp
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                burgers_interface.get_interface_states(grid, dt,
                                                       u, v,
                                                       ldelta_ux, ldelta_vx,
                                                       ldelta_uy, ldelta_vy)

    # apply transverse terms on both x- and y- interfaces
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                burgers_interface.apply_transverse_corrections(grid, dt,
                                                               u_xl, u_xr,
                                                               u_yl, u_yr,
                                                               v_xl, v_xr,
                                                               v_yl, v_yr)

    # apply pressure gradient correction terms to the interface state
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                apply_gradp_corrections(dt,
                                        u_xl, u_xr,
                                        u_yl, u_yr,
                                        v_xl, v_xr,
                                        v_yl, v_yr,
                                        gradp_x, gradp_y)

    # apply source terms
    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
                apply_other_source_terms(dt,
                                        u_xl, u_xr,
                                        u_yl, u_yr,
                                        v_xl, v_xr,
                                        v_yl, v_yr,
                                        source_x, source_y)

    # upwind using the MAC velocity to determine which state exists on
    # the interface
    u_xint = burgers_interface.upwind(grid, u_xl, u_xr, u_MAC)
    v_xint = burgers_interface.upwind(grid, v_xl, v_xr, u_MAC)
    u_yint = burgers_interface.upwind(grid, u_yl, u_yr, v_MAC)
    v_yint = burgers_interface.upwind(grid, v_yl, v_yr, v_MAC)

    return u_xint, v_xint, u_yint, v_yint




[docs]
def apply_gradp_corrections(dt,
                            u_xl, u_xr,
                            u_yl, u_yr,
                            v_xl, v_xr,
                            v_yl, v_yr,
                            gradp_x, gradp_y):
    r"""
    Parameters
    ----------
    grid : Grid2d
        The grid object
    dt : float
        The timestep
    u_xl, u_xr : ndarray ndarray
        left and right states of x-velocity in x-interface.
    u_yl, u_yr : ndarray ndarray
        left and right states of x-velocity in y-interface.
    v_xl, v_xr : ndarray ndarray
        left and right states of y-velocity in x-interface.
    v_yl, u_yr : ndarray ndarray
        left and right states of y-velocity in y-interface.

    Returns
    -------
    out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray
        correct the interface states of the left and right states of u and v on
        both the x- and y-interfaces interface states with the pressure gradient
        terms.
    """

    # Apply pressure gradient correction terms

    # transverse term for the u states on x-interfaces
    u_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
    u_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)

    # transverse term for the v states on x-interfaces
    v_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
    v_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)

    # transverse term for the v states on y-interfaces
    v_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
    v_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)

    # transverse term for the u states on y-interfaces
    u_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
    u_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)

    return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr




[docs]
def apply_other_source_terms(dt,
                            u_xl, u_xr,
                            u_yl, u_yr,
                            v_xl, v_xr,
                            v_yl, v_yr,
                            source_x, source_y):
    r"""
    Parameters
    ----------
    grid : Grid2d
        The grid object
    dt : float
        The timestep
    u_xl, u_xr : ndarray ndarray
        left and right states of x-velocity in x-interface.
    u_yl, u_yr : ndarray ndarray
        left and right states of x-velocity in y-interface.
    v_xl, v_xr : ndarray ndarray
        left and right states of y-velocity in x-interface.
    v_yl, u_yr : ndarray ndarray
        left and right states of y-velocity in y-interface.
    source_x, source_y : ndarray
        Any other source terms

    Returns
    -------
    out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray
        correct the interface states of the left and right states of u and v on
        both the x- and y-interfaces interface states with the source terms.
    """

    if source_x is not None:
        u_xl.ip(1, buf=2)[:, :] += 0.5 * dt * source_x.v(buf=2)
        u_xr.v(buf=2)[:, :] += 0.5 * dt * source_x.v(buf=2)
        u_yl.jp(1, buf=2)[:, :] += 0.5 * dt * source_x.v(buf=2)
        u_yr.v(buf=2)[:, :] += 0.5 * dt * source_x.v(buf=2)

    if source_y is not None:
        v_xl.ip(1, buf=2)[:, :] += 0.5 * dt * source_y.v(buf=2)
        v_xr.v(buf=2)[:, :] += 0.5 * dt * source_y.v(buf=2)
        v_yl.jp(1, buf=2)[:, :] += 0.5 * dt * source_y.v(buf=2)
        v_yr.v(buf=2)[:, :] += 0.5 * dt * source_y.v(buf=2)

    return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr