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

Improve time discretization #257

Merged
merged 27 commits into from
Sep 4, 2023
Merged

Improve time discretization #257

merged 27 commits into from
Sep 4, 2023

Conversation

dbrakenhoff
Copy link
Collaborator

@dbrakenhoff dbrakenhoff commented Aug 28, 2023

This PR aims to simplify and improve the handling of time discretization in nlmod.

To add time index to your model dataset use nlmod.time.set_ds_time()*.

*NOTE: this function replaces the old function with the same name, but works differently, so this is a minor breaking change. The old function can still be called through nlmod.time.set_ds_time_deprecated() for now.

The arguments to nlmod.time.set_ds_time() are:

  • time is an array-like with the timestamps at the end of each stress period (note that this does not include the starting time of the model)
  • start is the start datetime of the model
  • steady is an array-like indicating whether a stress period is steady-state (True or 1) or transient (False or 0). The default is steady-state for all timesteps.
  • nstp and tsmult are optional arguments to indicate the number of time steps per stress period and the timestep multiplier. These are used when writing the Modflow6 TDIS file.
# basic call signature
ds = nlmod.time.set_ds_time(ds, time, start, steady)

# set ds time with array of elapsed time periods
ds = nlmod.time.set_ds_time(ds, [100, 200], "2020-01-01", steady=[1, 0])

# set ds time with perlen
perlen = [100, 100]
ds = nlmod.time.set_ds_time(ds, np.cumsum(perlen), "2020-01-01", steady=[1, 0])  # note the cumsum

# set ds time with dates
ds = nlmod.time.set_ds_time(ds, ["2020-04-10", "2020-07-19"], "2020-01-01", steady=[1, 0])

A function has been added to convert flopy TDIS period data settings (perlen, nstp, tsmult) into a time index:

# for getting the time index for all time steps
ds = nlmod.time.ds_time_idx_from_tdis_settings(start, perlen, nstp, tsmult)

# for getting the time index at the end of each stress period, omit the nstp and tsmult keyword arguments
ds = nlmod.time.ds_time_idx_from_tdis_settings(start, perlen)

EDIT: Some additional comments on this PR after @OnnoEbbens review:

The variables steady, nstp and tsmult are now stored as data variables in the model dataset with dimension time. This should make them more visible and easier to modify for users as well.

A check is included whether the first entry of the time array (after conversion to a datetime object) is equal to or lies before the start datetime. This should help users figure out if they've made a mistake when assigning these variables.

- add new simplified set_ds_time() function
- deprecate old set_ds_time() by renaming to set_ds_time_deprecated()
- rename functions ds_time_idx_from_model and ds_time_idx_from_modeltime to make clear they return a time index
- handle new dataset time discretization
- add message if no STO pkg is created because model is steady state
- allow overriding of stored nstp and tsmult.
- use utilility methods to obtain nstp and tsmult from ds
Copy link
Collaborator

@OnnoEbbens OnnoEbbens left a comment

Choose a reason for hiding this comment

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

Good job! some comments from me nothing big.

nlmod/dims/time.py Outdated Show resolved Hide resolved
nlmod/dims/time.py Show resolved Hide resolved
nlmod/sim/sim.py Outdated Show resolved Hide resolved
nlmod/sim/sim.py Outdated Show resolved Hide resolved
tests/test_001_model.py Outdated Show resolved Hide resolved
tests/test_001_model.py Outdated Show resolved Hide resolved
- steady, nstp, tsmult become data variables with dim time
- add docstring to ds_time_idx_from_tdis_settings
- add TODO comments for next hydropandas release
- steady, nstp, tsmult become data variables with dim time
- add docstring to ds_time_idx_from_tdis_settings
- add TODO comments for next hydropandas release
@rubencalje
Copy link
Collaborator

rubencalje commented Aug 31, 2023

I liked the original method because of the following default values:

  • the first period would be steady state, following periods would be transient
  • the first steady state period had a length of 3652 days (10 years)

Multiple steady state stress periods is probably only useful in transport calculations, so I think a transient run for multiple stress periods should be the default. With the above default values the starting conditions of the transient run are determined by a steady-state run of 10 years. To get the same behavior of the new method would require a few lines of code from the user.

Also, the method set_ds_time_deprecated is not that useful, as it does not store the time settings in the right way for the rest of nlmod anymore.

Therefore, I would be more in favor of fixing/improving the original method, in which we:

  • deprecate some input-parameters like perlen and transient_timesteps (with a deprecation message, or just an appropriate exception message)
  • keep 'steady_start' with a default value of True (like it is now)
  • replace steady_state by steady, with a default value of False, and use the proposed method in this PR to store the steady_state stress periods in ds['steady'].

@dbrakenhoff
Copy link
Collaborator Author

dbrakenhoff commented Aug 31, 2023

Thanks for the review! I'm a bit biased, but i find the old code almost unreadable, so I'd be in favor of modifying this PR further to modify the defaults to what you suggest.

I agree we most often want to default to one steady state period and then subsequent transient periods. That's something we can implement with an extra keyword argument? And then explicitly entering the 10 year period isn't that much effort on the part of the user, just subtract 10 years from time[0] or just entering a year. This has the added advantage of being much more explicit :)?

So to reproduce current behavior it would become something like this (with these values being default for the steady keyword args):

nlmod.time.set_ds_time(np.cumsum([3652, 100, 100]), start="2010", steady=False, steady_start=True)

or in the pandas case:

time = pd.date_range("2020-01-01", "2025-01-01", "MS")
nlmod.time.set_ds_time(time, start=time[0] - pd.Timedelta(3652, "D"), steady=False, steady_start=True)

- add steady_start kwarg to set_ds_time
- modify set_ds_time_deprecated, also set steady, nstp, and tsmult in ds
@rubencalje
Copy link
Collaborator

The old code was complicated, because it needed to support perlen (which is mentioned as a positive in issue #124) and transient_timesteps. Let's forget the old code, and work on the lines in this PR.

It is a good idea to set the length of the steady-state start period more explicit. I can work with your proposed code:
time = pd.date_range("2020-01-01", "2025-01-01", "MS")
nlmod.time.set_ds_time(time, start="2010", steady_start=True)

@dbrakenhoff
Copy link
Collaborator Author

Alright, sounds good!

As for perlen, its still supported but you have to take the cumulative sum in this new function.

# set ds time with perlen
perlen = [100, 100]
ds = nlmod.time.set_ds_time(ds, np.cumsum(perlen), start="2020")  # note the cumsum

@rubencalje
Copy link
Collaborator

I added perlen, also because it is requested in issue #124, and set a default value for start.

I would remove the log-message at the start of the method, as this will get annoying, I think.

@dbrakenhoff
Copy link
Collaborator Author

Nice! My thoughts on these additions:

  • Adding a perlen option is nice! I do think it should be treated exactly the same as the time option (so type should always be array-like. Or we modify time to also support single values (but I'm more for option 1).
  • I'm not a big fan of the steady default change though. I wanted users to explicitly provide this value because it makes it very obvious what the time period you're modeling is going to be. I do like the option of allowing a number and interpreting that as a Timedelta relative to the first entry in time. This does add a layer of difficulty if time is passed as elapsed time, or if perlen is provided instead of time, in which case there is no reference point, but I saw you already included a check for that.

So in short, I'm for perlen always array-like, and getting rid of the steady defaults and forcing the user to pick something.

Perhaps we can get @OnnoEbbens to weigh in and help us pick a way to go :).

@rubencalje
Copy link
Collaborator

@OnnoEbbens agrees with you, and I agree that removing the default value for start is clearer to the user.

@rubencalje
Copy link
Collaborator

I fixed some bugs relating to previous changes and the new pandas version. They are not really related to this pull request, but were needed to test all the notebooks for the time-related changes. I think this PR is ready to be merged.

We use a newer version of modpath, so this is not a problem anymore
Copy link
Collaborator Author

@dbrakenhoff dbrakenhoff left a comment

Choose a reason for hiding this comment

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

I have one fix for nb 11, for the rest all good to go!

nlmod/sim/sim.py Outdated Show resolved Hide resolved
tests/test_001_model.py Outdated Show resolved Hide resolved
docs/examples/00_model_from_scratch.ipynb Outdated Show resolved Hide resolved
@dbrakenhoff dbrakenhoff merged commit 89fe3ea into dev Sep 4, 2023
1 check passed
@dbrakenhoff dbrakenhoff deleted the time_discretization_issues branch September 4, 2023 15:09
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
3 participants