modflowapi interface

[jlarsen-usgs]

Interface code and unit tests for modflowapi python package

Includes the addition of:

• Simulation: container for wrapping modflow simulations and allows users to access underlying data
• Model: model container that allows user to access model data and is a container that loads and provides access to packages
• ArrayPackage, ListPackage, AdvancedPackage: Containers for storing package information and loading/editing package data
• ArrayInput, ArrayPointer, ListInput: Containers for package data that allow the user to dynamically edit and set modflow stress data.
• run_model: a simulation runner for the modflowapi
• Callbacks: an enumeration scheme for creating custom callback functions to adjust model data
• unit tests: for interface methods and for running modflow6-examples
• example notebooks

[jlarsen-usgs]

@jdhughes-usgs @langevin-usgs @spaulins-usgs

Draft pull request for the modflowapi interface methods. The previous one had a sticky CI issue that required the repository to be recreated.

[christianlangevin]

@jlarsen-usgs, this is a great contribution. You've clearly put a lot of work and thought into this design, and I like where it's headed. I have only one small comment, and that is I wonder if it would be beneficial to rename some of the classes. For example, you have Simulation, Model, and PackageBase classes. Would it be worth considering ApiSimulation, ApiModel, and ApiPackageBase? I'm not completely sold on this renaming, but I can see people (mostly me) getting confused looking at a script that has a mix of flopy and modflowapi objects. To complicate things further for mf6 developers, we have yet another Simulation class/object, which has been confusing for us, and is on our list for refactoring.

[jlarsen-usgs]

@mjr-deltares and @w-bonelli

[mjr-deltares]

Hey @jlarsen-usgs , nice job! I am pretty sure that people have been waiting for this capability.

A couple of things came to mind when looking at the current implementation, first more at the design level:

• I agree with Chris to think about a different 'domain language' for the Model and Simulation concepts to avoid confusion. Also the name 'interface' for this framework on top of the API (=interface) might not be optimal. Maybe 'extensions'?

• MODFLOW API currently is mostly XMI + our get_var_address utility function, but it was not intended to be just that. Have you considered leveraging the information available in the Core by exposing it in convenient units. For example, I see you do an analysis of input var names to figure out solution names, model names etc. That works fine at the moment, but it is rather implicit and could be sensitive to future modifications (for parallel, the IDM, or even when a subcomponent is added to an exchange). We can introduce API functions to expose the required information in a stable way, say, a get_nr_models, get_models, get_solutions, get_exchanges, etc. And then it would also be easier to test/constrain its behavior in the MODFLOW 6 autotests.

• I was kind of thinking about Exchanges, which do not seem to be part of the current object model. I guess that might be a work in progress, but then I thought more general, how you might want to replicate the full MODFLOW internal object model as a tree of containers with certain input types (your ListInput, ArrayInput, etc.). In that case you would eventually have a simulation, which holds multiple solution(groups), a time discretization, models, exchanges. You can then drop the general XMI terminology such as 'subcomponent' and just start calling it solution, because that's what it is in case of the MODFLOW API.

• And finally (and I haven't thought this entirely through to be honest), could flopy also provide that structure/object model such that, after constructing the simulation with flopy, the API models and packages can be accessed as

  flopy.mf6.ModflowGwfchd.API.stress_period_data = ...

Some minor things:

• It should probably be 'run_simulation' instead of 'run_model'?
• I think I noticed a few occasions where you create the variable address string assuming a separator '/'. We have added get_var_address on the Fortran side to help with that operation without depending on the actual implementation in the memory manager. I can imagine you might sometimes need the reverse operation, this is available in the MemoryHelper in MODFLOW but has not yet been exposed through the API. I guess there would be nothing holding us back from adding that.

If it's more convenient to iterate on these points in a meeting, I am available.

[jlarsen-usgs]

@mjr-deltares

I've looked through your comments and made some updates to the draft PR.

• I've refactored the code and changed .interface to .extensions
• Updated class names Simulation to ApiSimulation, Model to ApiModel, and package names ex. WelPackage is now ApiWelPackage
• run_model updated to run_simulation()
• replaced all instances of manual variable address construction with mf6.get_var_address() calls.

I think it'd be worth having a meeting to discuss the other comments on the design level review and where/how to best implement some of this.

[jdhughes-dev]

@jlarsen-usgs wondering if

model = sim.get_model()

could be made to work so that it uses the first model name/id if a model name or id is not specified. This is what we have done for the solution groups in xmpy,

[jlarsen-usgs]

@jdhughes-usgs

I could update:

model = sim.get_model()

to return the first model name/id if a model isn't specified.

[jdhughes-dev]

@jlarsen-usgs
when I execute the following:

sim = ApiSimulation.load(mf6)
model = sim.get_model(sim.model_names[0])
model

I get:

MODEL_DFMF, 2 Layer, 10 Row, 11, Column model
Packages accessible include: 
  ArrayPackage objects:
    dis: <class 'modflowapi.extensions.apimodel.ApiDisPackage'>
    npf: <class 'modflowapi.extensions.apimodel.ApiNpfPackage'>
    sto: <class 'modflowapi.extensions.apimodel.ApiStoPackage'>
    ic: <class 'modflowapi.extensions.apimodel.ApiIcPackage'>
  ListPackage objects:
    drn_0: <class 'modflowapi.extensions.apimodel.ApiDrnPackage'>
    ghb_0: <class 'modflowapi.extensions.apimodel.ApiGhbPackage'>
    rch_0: <class 'modflowapi.extensions.apimodel.ApiRchPackage'>
  AdvancedPackage objects:
    buy: <class 'modflowapi.extensions.apimodel.ApiBuyPackage'>
    vsc: <class 'modflowapi.extensions.apimodel.ApiVscPackage'>
    gnc: <class 'modflowapi.extensions.apimodel.ApiGncPackage'>
    hfb: <class 'modflowapi.extensions.apimodel.ApiHfbPackage'>
    csub: <class 'modflowapi.extensions.apimodel.ApiCsubPackage'>
    mvr: <class 'modflowapi.extensions.apimodel.ApiMvrPackage'>

The model does not have any of the advanced packages. Is it intended that these objects are created even if they are not in a model?

[jdhughes-dev]

@jlarsen-usgs could we add a @property to get totim like the kstp, kper, etc. properties?

[jdhughes-dev]

@jlarsen-usgs I think there is also a need for a Callbacks.finalize. I could use it to finalize D-FLOW FM.

[jdhughes-dev]

Might be better if this string were Normal termination of simulation. in order to be consistent with the MODFLOW 6 normal termination message.

[jdhughes-dev]

@jlarsen-usgs How can I get and set auxvar from a boundary package?

[jlarsen-usgs]

@jdhughes-usgs auxiliary variables are now included in the stress period data recarray when they exist (naux > 0) within a modflow package.

[jdhughes-dev]

@jlarsen-usgs The auxvar modification works great. I see the dtype uses a default aux name

[('nodelist', 'O'), ('elev', '<f8'), ('cond', '<f8'), ('aux_0', '<f8')]

An enhancement would be to use the actual auxname.

[jlarsen-usgs]

@jdhughes-usgs I agree that the actual auxname in the numpy dtype would be a nice enhancement. At the moment it looks like either xmipy or the bmi doesn't support returning lists of strings and the call to auxname crashes.

[christianlangevin]

Hey @jlarsen-usgs, you should be able to get a copy of AUXNAME_CST through the xmipy, I think.

[jdhughes-dev]

@jlarsen-usgs Would be nice to be able to get [(nodelist), simvals] numpy recarray from the boundary packages. This would be useful in cases where simulated MODFLOW flows need to be passed to another model. In the past I have used custom code like

    nbound = mf6.get_value(
        mf6.get_var_address("NBOUND", modelname, drn_packagename)
    )[0]
    node_list = get_nodelist_ptr(modelname, mf6, drn_packagename)[:nbound]
    q = mf6.get_value(
        mf6.get_var_address("SIMVALS", modelname, drn_packagename)
    )[:nbound]
    drn_q[node_list - 1] = q[:]

to get this data from the MODFLOW model. Seems something similar to the way stress period data is accessed would work here as well.

[mjr-deltares]

Hey @jlarsen-usgs, you should be able to get a copy of AUXNAME_CST through the xmipy, I think.

Hey @jlarsen-usgs, here is a test of the new functionality for getting a string array (in this case it is BOUNDNAME_CST):

https://github.com/Deltares/xmipy/blob/96976220712f4118629f9ccaf5e4b669cd8d6128/tests/test_mf6_dis_bmi.py#L223

[jlarsen-usgs]

@mjr-deltares

It looks like xmipy is working correctly for getting a string array. I tested it with BOUNDNAME_CST and it works well.
The issue looks like nothing is being stored in AUXNAME_CST on the modflow side. When I get the variable shape xmipy returns 0. To futher investigate I replaced the strtype = "<S" + str(ilen + 1) with strtype = "O" in my local version of xmipy. This ends up returning an empty array for the AUXNAME_CST variable when get_value() is called even though there are auxilary variables defined in the WEL package file.

[mjr-deltares]

Hi @jlarsen-usgs, I just ran a minimal example with a GWF model + WEL with 2 aux variables. It looks like the get_value call returns AUXNAME_CST correctly for this case. If you want, we can set up a call so that I can show you what i did.

[jlarsen-usgs]

@mjr-deltares updated my libmf6.dll and I am now able to access AUXNAME_CST.

[mjr-deltares]

excellent!

[christianlangevin]

sweet