Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ten diffusive domains w schism boundary #573

Merged

Conversation

kumdonoaa
Copy link
Contributor

@kumdonoaa kumdonoaa commented Jul 7, 2022

The changes made here enabled for diffusive wave channel routing to run in serial mode on all ten available mainstem domains where channel cross section survey data are available. Also, schism's water depth output was integrated with t-route so that the water depth values along shorelines were used as downstream boundary conditions for the diffusive wave.

Additions

  • diffusive.f90: 1) Using schism water depth values as diffusive wave's downstream boundary condition is activated with minimum depth enforced for numerical stability; 2) Commented out maxCourant variable as it is solely used for bookkeeping with the initial value undefined
  • nhd_io.py: 1) In Lines 505~510, if statement was added just in case when flowveldepth.index doesn't have an index of link_gage_df, 2) in line 1582, a new equation to compute water depth from a channel bottom using Schism ouputs was implemented.
  • input.py: diffusive_parameters dictionary was simply replaced by hybrid_parameters that include diffusive_parameters.
  • compute.py: In lines 1158 to1159 and 1172 and 1173, initialization of related numpy arrays was introduced for running multiple diffusive domains
  • diffusive_utils.py: Lines 1055~1060 prevents dbcd_df (dataframe of schism water depth for diffusive boundary condition) from having any zero or negative water depth.
  • HurricaneLaura/domain/coastal_10diffusive.txt: mask files for Hurricane Laura (or Atl/Gulf covering all ten diffusive domains) domain
  • HurricaneLaura/domain/coastal_domain_subset.yaml: lists of links of diffusive domains on original hydrofabric
  • HurricaneLaura/domain/refactored_coastal_domain_subset.yaml: lists of links of diffusive domains on refactored hydrofabric (refactored so as to remove too short stream segments by converging to adjacent segments)
  • HurricaneLaura/domain/final_diffusive_natural_xs.nc: natural channel cross section data on original hydrofabric including (x, z) vertices data and manning's N values at each point
  • HurricaneLaura/domain/refac_final_diffusive_natural_xs.nc: natural channel cross section data on refactored hydrofabric including (x, z) vertices data and manning's N values at each point

Removals

Changes

Testing

  1. In order to run diffusive/MC on Hurricane Laura domain, that is, diffusive on ten mainstems and MC on the rest tributary networks), run test_AnA.yaml in test/HurricaneLaura. The command is "python3 -m nwm_routing -f test_AnA.yaml" after compiling using ./complier.sh in t-route directory.

Screenshots

Notes

Todos

Checklist

  • PR has an informative and human-readable title
  • Changes are limited to a single goal (no scope creep)
  • Code can be automatically merged (no conflicts)
  • Code follows project standards (link if applicable)
  • Passes all existing automated tests
  • Any change in functionality is tested
  • New functions are documented (with a description, list of inputs, and expected output)
  • Placeholder code is flagged / future todos are captured in comments
  • Visually tested in supported browsers and devices (see checklist below 👇)
  • Project documentation has been updated (including the "Unreleased" section of the CHANGELOG)
  • Reviewers requested with the Reviewers tool ➡️

Testing checklist

Target Environment support

  • Windows
  • Linux
  • Browser

Accessibility

  • Keyboard friendly
  • Screen reader friendly

Other

  • Is useable without CSS
  • Is useable without JS
  • Flexible from small to large screens
  • No linting errors or warnings
  • JavaScript tests are passing

@kumdonoaa kumdonoaa requested a review from hellkite500 July 7, 2022 22:39
Copy link
Contributor

@hellkite500 hellkite500 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we need to include any of the test/HurricaneLaura in the repository itself. We can make an archive of the relevant domain data and (relative) configurations, but they shouldn't get commited to the repo.

@@ -485,12 +502,20 @@ def write_chanobs(
-------------

'''
# TODO: for and if statements are improvised just in case when link of gage location is not in flowveldepth, which should be fixed
link_na = []
for link in list(link_gage_df.index):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think you need to convert this to a list to simply iterate it, same with the next .index...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

list() removed. Thanks.

for link in list(link_gage_df.index):
if link not in list(flowveldepth.index):
link_na.append(link)
link_gage_df_nona = link_gage_df.drop(link_na)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could probably just drop this inplace...link_gaga_df.drop(link_na, inplace=True)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

doesn't look like this change got committed/pushed???

start_date = ds.variables['time'].units

start_date = dparser.parse(start_date,fuzzy=True)
dt_schism = 3600 # [sec]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could this be dynamically inferred from the timesteps array:

if len(timesteps) > 1:
    dt_schism = timesteps[1]-timesteps[0]
else:
    raise RuntimeError("schism provided less than 2 time steps")

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

changed. Thanks for your code.

start_date = dparser.parse(start_date,fuzzy=True)
dt_schism = 3600 # [sec]
dt_timeslice = timedelta(minutes=dt_schism/60.0)
start_date = start_date # + dt_timeslice
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove this line

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done.

@@ -1547,7 +1612,8 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds)
double precision :: x1, x2, df, dxx, dxold, f, fh, fl, temp, xh, xl
double precision :: y_norm, y_ulm_multi, y_llm_multi, elv_norm, y_old
double precision :: rtsafe

open(unit=301, file="./output/rtsafe_check.txt")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove open file unit

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed.

else ! newton step acceptable. take it.
dxold = dxx
dxx = f / df
temp = rtsafe
rtsafe = rtsafe - dxx
if (temp == rtsafe) return
if (temp == rtsafe) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove this whole block since it is a noop now anyways.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed. thanks.

end if

if (abs(dxx) < xacc) return ! convergence criterion.
if (abs(dxx) < xacc) then
!write(301, *) "fnal:", i, j, iter, rtsafe
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

another noop block here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done. Thanks.


allocate(el1(nel), a1(nel), peri1(nel), redi1(nel), redi1All(nel))
allocate(equiv_mann(nel), conv1(nel), tpW1(nel))
allocate(newdKdA(nel))
allocate(compoundSKK(nel), elev(nel))
allocate(i_start(nel), i_end(nel))

!open(unit=201, file="./output/xsec_natural.txt")
open(unit=301, file="./output/xsec_natural_check.txt")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove this file unit

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

!iv = 0
!do ii2 = iel, nel
! iv = iv + 1
! dmyarr(iv) = newdKdA(ii2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove any unused/unwanted code from this funnction

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

para_ar_g[2] = 50.0 # lower limit of diffusivity (default: 50)
para_ar_g[3] = 5000.0 # upper limit of diffusivity (default: 1000)
para_ar_g[2] = 10.0 # lower limit of diffusivity (default: 50)
para_ar_g[3] = 40000.0 # upper limit of diffusivity (default: 1000)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If these are the default values, then the comments should be updated appropriately...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed. Thanks.

@hellkite500
Copy link
Contributor

what is this file for? coastal_10diffusive.txt and is it the same in both lower colorado and hurricane laura?

@kumdonoaa
Copy link
Contributor Author

what is this file for? coastal_10diffusive.txt and is it the same in both lower colorado and hurricane laura?

It is a mask file for diffusive domains. Even it can be used for both cases, I can slim it down for Lower Colorado. I will adjust it shortly.

@hellkite500
Copy link
Contributor

Is the one in the hurricane laura directory still needed? test/HurricaneLaura/domain/coastal_10diffusive.txt

Copy link
Contributor

@hellkite500 hellkite500 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this can be merged for now. I think some additional refactoring and cleanup would be useful in places, but we can address that once the functionality is merged and can changes can be tested against a "stable" implementation/

for link in list(link_gage_df.index):
if link not in list(flowveldepth.index):
link_na.append(link)
link_gage_df_nona = link_gage_df.drop(link_na)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

doesn't look like this change got committed/pushed???

@kumdonoaa kumdonoaa merged commit 0c9ec38 into NOAA-OWP:master Jul 18, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants