NOAA-OWP / t-route

Tree based hydrologic and hydraulic routing
Other
43 stars 48 forks source link

Enable Courant output in v02 framework #225

Closed awlostowski-noaa closed 2 months ago

awlostowski-noaa commented 3 years ago

Courant metrics are not immediately available in the v02 framework and we need to change this. While the Fortran Muskingum-Cunge routing subroutine muskingcungenwm in MCsingleSegStime_f2py_NOLOOP.f90 calculates and outputs Courant metrics (cn, ck, and X), these variables are silently dropped in the subroutine c_muskingcungenwm in pyMCsingleSegStime_NoLoop.f90.

Proposed solution steps

  1. Edit pyMCsingleSegStime_NoLoop.f90 such that c_muskingcungenwm returns ck, cn, and X. Specifically, designate these variables as intent out and expand the subroutine argument list accordingly.
  2. Edit fortran_wrappers.pyx such that c_muskingcungenwm arguments include ck, cn, and X.
  3. Edit reach.pyx as above. not sure what *rv means, but seems important in packaging return variables. Should we add Courant metrics to this grouping?
  4. Edit mc_reach.compute_reach_kernel to handle additional returns. Perhaps output_buf can carry additional Courant metric returns?
  5. Edit mc_reach.compute_network. Expand size of out_view to accept Courant metrics? Create a new array, like flowveldepth, to store network-wide Courant results.
  6. Add argument to mc_reach.compute_network to optionally return Courant metrics.

ping @jameshalgren

awlostowski-noaa commented 3 years ago

Step 1

subroutine c_muskingcungenwm(dt, qup, quc, qdp, ql, dx, bw, tw, twcc,&
    n, ncc, cs, s0, velp, depthp, qdc, velc, depthc, ck, cn, X) bind(c) # adjust args

    real(c_float), intent(in) :: dt
    real(c_float), intent(in) :: qup, quc, qdp, ql
    real(c_float), intent(in) :: dx, bw, tw, twcc, n, ncc, cs, s0
    real(c_float), intent(in) :: velp, depthp
    real(c_float), intent(out) :: qdc, velc, depthc
    real(c_float), intent(out) :: ck, cn, X # change these to intent out

    call muskingcungenwm(dt, qup, quc, qdp, ql, dx, bw, tw, twcc,&
    n, ncc, cs, s0, velp, depthp, qdc, velc, depthc, ck, cn, X)
awlostowski-noaa commented 3 years ago
@cython.boundscheck(False)
cdef void muskingcunge(float dt,
        float qup,
        float quc,
        float qdp,
        float ql,
        float dx,
        float bw,
        float tw,
        float twcc,
        float n,
        float ncc,
        float cs,
        float s0,
        float velp,
        float depthp,
        QVD *rv) nogil:
    cdef:
        float qdc = 0.0
        float depthc = 0.0
        float velc = 0.0

    c_muskingcungenwm(
        &dt,
        &qup,
        &quc,
        &qdp,
        &ql,
        &dx,
        &bw,
        &tw,
        &twcc,
        &n,
        &ncc,
        &cs,
        &s0,
        &velp,
        &depthp,
        &qdc,
        &velc,
        &depthc) # three additional returns in here. 
    rv.qdc = qdc
    rv.depthc = depthc
    rv.velc = velc
groutr commented 3 years ago

@awlostowski-noaa to help clarify some of your questions:

  1. *rv is a C level pointer. It points to a struct to hold output values. It is there to avoid allocating memory for every call to muskingcunge (this was a critical optimization early on). We allocate the struct once, and reuse it for all computations, copying the values when needed.
  2. To keep consistent with paradigm here, output_buf should be expanded to contain the new return values. This is an attempt to keep the return values contiguous in memory so we can hopefully quickly access them.
JurgenZach-NOAA commented 2 months ago

Addressed in PRs 226, 380