Skip to content

Commit

Permalink
make diffusive routing capable of running on refactored hydrofabric o…
Browse files Browse the repository at this point in the history 
…f Lower Colorado River, TX (#561)

- Enabled diffusive wave routing for refactored hydrofabric with ~10x compute speed increase compared to original hydrofabric

- Assimilated USGS observed streamflow data to diffusive streamflow at stream segments where the observed data is available
  • Loading branch information
@kumdonoaa
kumdonoaa authored May 24, 2022
1 parent a01b1be commit c6d8c9d
Show file tree
Hide file tree
Showing 15 changed files with 2,619 additions and 330 deletions.
435 changes: 299 additions & 136 deletions src/kernel/diffusive/diffusive.f90
View file

Large diffs are not rendered by default.

46 changes: 27 additions & 19 deletions src/kernel/diffusive/pydiffusive.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,43 +5,51 @@ module diffusive_interface

implicit none
contains
subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, &
nts_qtrib_g, mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, &
tw_ar_g, twcc_ar_g, mann_ar_g, manncc_ar_g, so_ar_g, dx_ar_g, &
iniq, frnw_col, frnw_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, &
mann_bathy_g, size_bathy_g, q_ev_g, elv_ev_g) bind(c)
subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qtrib_g, nts_da_g, &
mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, mann_ar_g, &
manncc_ar_g, so_ar_g, dx_ar_g, &
iniq, frnw_col, frnw_ar_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, mann_bathy_g, size_bathy_g, &
usgs_da_g, usgs_da_reach_g, rdx_ar_g, cwnrow_g, cwncol_g, crosswalk_g, &
q_ev_g, elv_ev_g) bind(c)

integer(c_int), intent(in) :: nts_ql_g, nts_ub_g, nts_db_g, nts_qtrib_g
integer(c_int), intent(in) :: nts_ql_g, nts_ub_g, nts_db_g, nts_qtrib_g, nts_da_g
integer(c_int), intent(in) :: ntss_ev_g
integer(c_int), intent(in) :: mxncomp_g, nrch_g
integer(c_int), intent(in) :: frnw_col
integer(c_int), intent(in) :: paradim
integer(c_int), intent(in) :: mxnbathy_g
integer(c_int), dimension(nrch_g, frnw_col), intent(in) :: frnw_g
integer(c_int), intent(in) :: mxnbathy_g
integer(c_int), intent(in) :: cwnrow_g
integer(c_int), intent(in) :: cwncol_g
integer(c_int), dimension(nrch_g), intent(in) :: usgs_da_reach_g
integer(c_int), dimension(nrch_g, frnw_col), intent(in) :: frnw_ar_g
integer(c_int), dimension(mxncomp_g, nrch_g), intent(in) :: size_bathy_g
real(c_double), dimension(nts_db_g), intent(in) :: dbcd_g
real(c_double), dimension(paradim), intent(in) :: para_ar_g
real(c_double), dimension(:), intent(in) :: timestep_ar_g(8)
real(c_double), dimension(:), intent(in) :: timestep_ar_g(9)
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: z_ar_g, bo_ar_g, traps_ar_g
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: tw_ar_g, twcc_ar_g
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: mann_ar_g, manncc_ar_g
real(c_double), dimension(mxncomp_g, nrch_g), intent(inout) :: so_ar_g
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g, iniq
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: dx_ar_g, rdx_ar_g
real(c_double), dimension(mxncomp_g, nrch_g), intent(in) :: iniq
real(c_double), dimension(nts_ub_g, nrch_g), intent(in) :: ubcd_g
real(c_double), dimension(nts_qtrib_g, nrch_g), intent(in) :: qtrib_g
real(c_double), dimension(nts_da_g, nrch_g), intent(in) :: usgs_da_g
real(c_double), dimension(nts_ql_g, mxncomp_g, nrch_g), intent(in) :: qlat_g
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g
real(c_double), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g

call diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, &
nts_qtrib_g, mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, &
tw_ar_g, twcc_ar_g, mann_ar_g, manncc_ar_g, so_ar_g, dx_ar_g, &
iniq, frnw_col, frnw_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, &
mann_bathy_g, size_bathy_g, q_ev_g, elv_ev_g)
real(c_double), dimension(cwnrow_g, cwncol_g), intent(in ) :: crosswalk_g
real(c_double), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g

call diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qtrib_g, nts_da_g, &
mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, mann_ar_g, &
manncc_ar_g, so_ar_g, dx_ar_g, &
iniq, frnw_col, frnw_ar_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, mann_bathy_g, size_bathy_g, &
usgs_da_g, usgs_da_reach_g, rdx_ar_g, cwnrow_g, cwncol_g, crosswalk_g, &
q_ev_g, elv_ev_g)

end subroutine c_diffnw
end module diffusive_interface
85 changes: 81 additions & 4 deletions src/troute-network/troute/nhd_network_utilities_v02.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -511,7 +511,7 @@ def build_da_sets(da_params, run_sets, t0):
"usace_timeslices_folder",
None
)

# User-specified DA ON/OFF preferences
usace_da = False
usgs_da = False
Expand Down Expand Up @@ -552,7 +552,7 @@ def build_da_sets(da_params, run_sets, t0):
dt_timeslice = timedelta(minutes = 15)

da_sets = [] # initialize list to store TimeSlice set lists

# Loop through each run set and build lists of available TimeSlice files
for (i, set_dict) in enumerate(run_sets):

Expand Down Expand Up @@ -753,8 +753,10 @@ def build_data_assimilation_lastobs(data_assimilation_parameters):
da_parameter_dict = {}
da_parameter_dict["da_decay_coefficient"] = data_assimilation_parameters.get(
"da_decay_coefficient",
120
)
120)
da_parameter_dict["diffusive_streamflow_nudging"] = streamflow_da_parameters.get(
"diffusive_streamflow_nudging", False)

# TODO: Add parameters here for interpolation length (14/59), QC threshold (1.0)

return lastobs_df, da_parameter_dict
Expand Down Expand Up @@ -798,3 +800,78 @@ def build_data_assimilation_folder(data_assimilation_parameters, run_parameters)
usgs_df = pd.DataFrame()

return usgs_df


def build_refac_connections(diff_network_parameters):
'''
Construct network connections network for refacored dataset.
Arguments
---------
diff_network_parameters (dict): User input network parameters
Returns:
--------
connections (dict int: [int]): Network connections
'''

# crosswalking dictionary between variables names in input dataset and
# variable names recognized by troute.routing module.
cols = diff_network_parameters.get(
'columns',
{
'key' : 'link',
'downstream': 'to',
'dx' : 'Length',
'n' : 'n',
'waterbody' : 'NHDWaterbodyComID',
'gages' : 'gages',
'alt' : 'z',
'line_d' : 'xid_d',
}
)

# read parameter dataframe
param_df = nhd_io.read(pathlib.Path(diff_network_parameters["geo_file_path"]))

# numeric code used to indicate network terminal segments
terminal_code = set(param_df.to.unique()) - set(param_df.link.unique())

# select the column names specified in the values in the cols dict variable
param_df = param_df[list(cols.values())]

# rename dataframe columns to keys in the cols dict variable
param_df = param_df.rename(columns=nhd_network.reverse_dict(cols))

# set parameter dataframe index as segment id number, sort
param_df = param_df.set_index("key").sort_index()

# get and apply domain mask
if "mask_file_path" in diff_network_parameters:
data_mask = nhd_io.read_mask(
pathlib.Path(diff_network_parameters["mask_file_path"]),
layer_string=diff_network_parameters.get("mask_layer_string", None),
)
data_mask = data_mask.set_index(data_mask.columns[0])
param_df = param_df.filter(data_mask.index, axis=0)

# There can be an externally determined terminal code -- that's this first value
terminal_codes = set()
terminal_codes.update(terminal_code)
# ... but there may also be off-domain nodes that are not explicitly identified
# but which are terminal (i.e., off-domain) as a result of a mask or some other
# an interior domain truncation that results in a
# otherwise valid node value being pointed to, but which is masked out or
# being intentionally separated into another domain.
terminal_codes = terminal_codes | set(
param_df[~param_df["downstream"].isin(param_df.index)]["downstream"].values
)

param_df_unique = param_df.drop_duplicates("downstream")
# build connections dictionary
connections = nhd_network.extract_connections(
param_df_unique, "downstream", terminal_codes=terminal_codes
)

return connections

11 changes: 10 additions & 1 deletion src/troute-nwm/src/nwm_routing/__main__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,8 @@ def nwm_route(
diffusive_parameters,
diffusive_network_data,
topobathy_data,
refactored_diffusive_domain,
refactored_reaches,
subnetwork_list,
):

Expand Down Expand Up @@ -133,13 +135,14 @@ def nwm_route(

LOG.debug("MC computation complete in %s seconds." % (time.time() - start_time_mc))
start_time_diff = time.time()

# call diffusive wave simulation and append results to MC results
results.extend(
compute_diffusive_routing(
results,
diffusive_network_data,
cpu_pool,
t0,
dt,
nts,
q0,
Expand All @@ -151,6 +154,8 @@ def nwm_route(
diffusive_parameters,
waterbodies_df,
topobathy_data,
refactored_diffusive_domain,
refactored_reaches,
)
)
LOG.debug("Diffusive computation complete in %s seconds." % (time.time() - start_time_diff))
Expand Down Expand Up @@ -365,6 +370,8 @@ def main_v03(argv):
usace_lake_gage_crosswalk,
diffusive_network_data,
topobathy_data,
refactored_diffusive_domain,
refactored_reaches,
) = nwm_network_preprocess(
supernetwork_parameters,
waterbody_parameters,
Expand Down Expand Up @@ -499,6 +506,8 @@ def main_v03(argv):
diffusive_parameters,
diffusive_network_data,
topobathy_data,
refactored_diffusive_domain,
refactored_reaches,
subnetwork_list,
)

Expand Down
Loading

0 comments on commit c6d8c9d

Please sign in to comment.