Modifying time series data
pya opened this issue · 13 comments
I run some tests modifying values of boundary conditions. Changing values of variables such as <model_name>/WEL/BOUND works and impacts the model results. But there are no effects on model runs if the data is read from time series entries. While the time series values "show up" in variables such as <model_nam>/WEL/BOUND, their modification via xmipy has no impact on the model results. Is this a known problem? If so, is there a workaround? I can create some reproducible examples. This will take some time. So far I was working with ex-gwf-advtidal from the MODFLOW 6 examples. It has some time series data and some directly specified boundary condition values mixed in files such as ex-gwf-advtidal.wel.
Hello @pya, quick check: how did you run the time loop, does it include using prepare_time_step ? It could be that the values you set through the API are subsequently overwritten with the usual Read-Prepare step inside MODFLOW. The following picture from the MODFLOW API paper (https://doi.org/10.1016/j.envsoft.2021.105257) shows the construct to use:
and then you should be able to set the BOUND values after the prepare_time_step call. Please let us know if this solves your issue or not.
Yes, I use prepare_time_step before modifying values and calling do_time_step. This is needed to see the impact of modifying BOUND values that are defined as non-time-series input, i.e. in files like *.wel. So it works for BOUND in general but some part of the MODFLOW 6 code must read the time series data (defined in *.tsfiles) again after modification of the values and before the time step is finalized. This seems to be the only explaination that makes sense to me.
BTW, thanks for the pointer to the paper. I was not aware of it. Maybe I just overlook some important part of the documentation. I would like to suggest to link to it from the xmipy README. I used this source file as documentation.
Sorry for the delay, I've been out of office last week. Maybe you can produce a script that runs on one of the MODFLOW examples and reproduce this behavior? Then I will give it a spin to see where this behavior comes from.
Yes, creating a re-produceable example is on my task list. It will take me a few days until I get to it.
I took longer anticipated, but I put together a (hopefully) reproducible test. Here are the details.
Hi @pya, thanks for taking the effort to put the example together! I am working on other projects for the coming two weeks, but I will have a look at it when I get the chance.
Hi @mjr-deltares, take your time. Everything should be re-producible. Let me know if not.
Hi @mjr-deltares, any news regarding my test? Do you need any further information?
Hi @pya, I've had a look at the working example you provided and found the issue. When running with timeseries, other than with stress period data, the values are read from file in the Advance (AD) stage, which is part of the do_time_step call. The value you were trying to set programmatically through Python is then subsequently overwritten by the one from file. Now, you can avoid that by interacting with MODFLOW one level deeper and explicitly call the solve methods. I have modified and tested your example to see that it works. Replace the time loop in run_model from the file run_xmipy_model.py with the following:
while current_time < end_time:
dt = mf6.get_time_step()
mf6.prepare_time_step(dt)
mf6.prepare_solve(1)
if get_value is not None:
value = get_value(current_time)
wel[:] = value
if verbosity_level > 0:
print(f'{current_time = :6.2f}, {value = :6.2f}')
mxit_tag = mf6.get_var_address("MXITER", "SLN_1")
max_iter = mf6.get_value_ptr(mxit_tag)[0]
kiter = 0
while kiter < max_iter:
has_converged = mf6.solve(1)
kiter += 1
if has_converged:
break
mf6.finalize_solve(1)
mf6.finalize_time_step()
current_time = mf6.get_current_time()
(NB: for clarity I have hardcoded for a single numerical solution, i.e. the "1" in the calls, instead of a loop over multiple solutions, which would be a more general approach)
Let me know if this works as intended.
Works like charm. Thanks a lot. This is my version:
def run_model(nam_file, get_value=None, verbosity_level=0, dll_path=DLL_PATH):
"""Run a Model via `xmipy`."""
model_path = Path(nam_file).parent
model_name = model_path.name.upper()
with cd(model_path):
mf6 = XmiWrapper(dll_path)
mf6.initialize(nam_file)
wel_tag = mf6.get_var_address('BOUND', model_name, 'WEL_0')
wel = mf6.get_value_ptr(wel_tag)
mxit_tag = mf6.get_var_address("MXITER", "SLN_1")
max_iter = mf6.get_value_ptr(mxit_tag)[0]
current_time = mf6.get_current_time()
end_time = mf6.get_end_time()
while current_time < end_time:
dt = mf6.get_time_step()
mf6.prepare_time_step(dt)
mf6.prepare_solve(1)
if get_value is not None:
value = get_value(current_time)
wel[:] = value
if verbosity_level > 0:
print(f'{current_time = :6.2f}, {value = :6.2f}')
for _ in range(max_iter):
has_converged = mf6.solve(1)
if has_converged:
break
mf6.finalize_solve(1)
mf6.finalize_time_step()
current_time = mf6.get_current_time()
mf6.finalize()I made two changes:
- I moved the definition of
max_iterout of the loop. This assumes it is constant and does NOT change over time. Is this assumption correct? - I replaced the
whilewith aforloop. Should be equivalent. Is this correct or did I overlook something?
I will make the solution work with multiple solutions. I will post my attempt here to make sure I understand the multi-solution business correctly.
This is my version for multiple solutions:
def run_model(nam_file, get_value=None, verbosity_level=0, dll_path=DLL_PATH):
"""Run a Model via `xmipy`."""
model_path = Path(nam_file).parent
model_name = model_path.name.upper()
with cd(model_path):
mf6 = XmiWrapper(dll_path)
mf6.initialize(nam_file)
wel_tag = mf6.get_var_address('BOUND', model_name, 'WEL_0')
wel = mf6.get_value_ptr(wel_tag)
nsol = sum(1 for name in mf6.get_input_var_names()
if name.startswith('SLN_') and name.endswith('MXITER'))
sol_numbers = list(range(1, nsol + 1))
max_iters = []
for sol_number in sol_numbers:
mxit_tag = mf6.get_var_address('MXITER', f'SLN_{sol_number}')
max_iters.append(mf6.get_value_ptr(mxit_tag)[0])
current_time = mf6.get_current_time()
end_time = mf6.get_end_time()
while current_time < end_time:
dt = mf6.get_time_step()
mf6.prepare_time_step(dt)
for sol_number in sol_numbers:
mf6.prepare_solve(sol_number)
if get_value is not None:
value = get_value(current_time)
wel[:] = value
if verbosity_level > 0:
print(f'{current_time = :6.2f}, {value = :6.2f}')
for sol_number, max_iter in zip(sol_numbers, max_iters):
for _ in range(max_iter):
has_converged = mf6.solve(sol_number)
if has_converged:
break
for sol_number in sol_numbers:
mf6.finalize_solve(sol_number)
mf6.finalize_time_step()
current_time = mf6.get_current_time()
mf6.finalize()Works for my example. Is this correct?
The method for finding the number of solutions seems a bit crude, but works. ;)
Hi @pya, nice to see you have everything working now, and yes, it looks good to me. I will close this issue now.
A shortcut to getting the number of solutions would be the function mf6.get_subcomponent_count(). The reason we don't have a cleaner way to get to the solution names/variables here is that that is specific to MODFLOW where xmipy is meant to be more general. There are quite a few recent developments in the modflowapi Python project that might interest you. If you are exclusively dealing with MODFLOW 6, that would be a great place to have a look: https://github.com/MODFLOW-USGS/modflowapi. All your code should work with minimal effort since we have modflowapi extending xmipy.