Very slow MF6 external file writing

[briochh]

The MODFLOW simulation write step (or specifically, the .set_all_data_external() step) seems to be prohibitively slow, especially for models with a large number of rows and columns, or stress periods. e.g.:

ar = dis.top.array.copy()
ar.shape
(302, 501)
%timeit np.savetxt(os.path.join(gwf.model_ws, 'top_ar.txt'), ar, delimiter=' ', fmt="%15.8f")
45.5 ms ± 3.96 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit dis.top.store_as_external_file(os.path.join(gwf.model_ws, 'top.txt'))
929 ms ± 77.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

1 sec write time for a single array quickly causes issues when the model contains many stress periods!

This seems to relate to the amount of time spent on overhead for file formatting (esp mffileaccess.get_data_string()). Is there any way to speed this up? Alternatively, would it be possible to support passing no data or empty arrays to Flopy package constructers, set-up the open/close lines in the packages, and then allow the user to write the external arrays etc in a more optimised way (e.g. numpy.savetxt(), and pandas.to_csv() etc.)?

[spaulins-usgs]

Yes you can tell Flopy that you already have external data set up in a file by specifying the "filename" key and omitting the "data" key in a dictionary defining the data. For example, the following points the top array to an existing "top.txt", but creates a botm array with constant values of 50.0.

top = {"filename": "top.txt"}
botm = {"filename": "botm.txt", "data": 50.0}
dis_package = flopy.mf6.ModflowGwfdis(
    model,
    length_units="FEET",
    nlay=1,
    nrow=1,
    ncol=10,
    delr=500.0,
    delc=500.0,
    top=top,
    botm=botm,
    filename="{}.dis".format(model_name),
    pname="mydispkg",
)

[briochh]

Thanks @spaulins-usgs, that might work for me.

What is the approach if we wanted to define pre-written layer-based files (e.g. for botm)?

Do we pass a dict with layer as keys (botm={i: {"filename": f"botm_{i}.txt"} for i in range(nlay)}) to the package constructor? or a list of filenames (botm={"filename": [ f"botm_{i}.txt" for i in range(nlay)]})? Sorry if I am missing something in the Flopy docs.

[briochh]

Ok think I have got it... botm=[{"filename": f"botm_{i}.txt"} for i in range(nlay)]?

[spaulins-usgs]

Yes that will work.

[briochh]

Hi @spaulins-usgs, A continuation of this. I am using the above approach for the GWT sources stress period data.
E.g.

# where df is say t,k,i,j,c
spd = {}
for t in range(per):
    fname = f"src_stress_period_data_{t}.txt"
    df.loc[df.t==t, ['k', 'i', 'j']].to_csv(os.path.join(gwt.model_ws, fname),
              sep=' ', index=False, header=False)
    spd[t] = {'filename': fname} 
src = flopy.mf6.ModflowGwtsrc( gwt, stress_period_data=spd)

I have 735 stress periods, and am specifying a source in every surface cell -- so a fair bit to set up. However, the overheads in instantiating that object are pretty killer. Writing the stress period data to external files in the loop takes about 1 minute. But initialising that flopy.mf6.ModflowGwtsrc package object by passing a dict of filenames takes almost 1 hour. It seems very slow to read in the external data files.

Looking at the profiling it looks that there are (2*number of files) calls to modelgrid (and subsequent functions) in the check_valid_cellids() method in mfdatalist.py, at over 1 second per call this takes up about half of the package initialisation time.

Cacheing what we need from modelgrid might be an opportunity for a quick speed gain? Alternatively, can we pass an arg to the flopy.mf6.ModflowGwtsrc( ) that sets check_data = False (maybe in that spd dictionary)?

The other expensive operation appears to relate to reading in the data, specifically in the split_data_line() method in datautil.py. Not sure what can be done with this, as I understand the desire for flexibility and being able split out comments etc., but I think there should be some gains available.

I can look into it but there is always a danger that I break more than I fix!

[spaulins-usgs]

I think it is reasonable to have an option for flopy to skip the data check for external data. Try changing line 427 of mfdatalist.py from:

    if check_data:

to:

    if check_data and data_check is not None:

If that solves your problem I will add a simulation-level flag to turn off data checking, possibly one flag for just external data and one flag for all data.

[briochh]

@spaulins-usgs, yep that speeds things up. It now takes a split second to instantiate that package object.

I wonder if we could add flexibility to whether checks are carried out or not. I have hacked in specifying check: False in the stress_period_data dictionary, e.g. from above:

# where df is say t,k,i,j,c
spd = {}
for t in range(per):
    fname = f"src_stress_period_data_{t}.txt"
    df.loc[df.t==t, ['k', 'i', 'j']].to_csv(os.path.join(gwt.model_ws, fname),
              sep=' ', index=False, header=False)
    spd[t] = {'filename': fname, 'check': False} 
src = flopy.mf6.ModflowGwtsrc( gwt, stress_period_data=spd)

This together with a very quick lookup around line 1517 of mfdatalist.py and propagating through the set_data() and _set_data() function args, would support turning off checking file by file. I can PR this if you want.

Generally though, this might me more treating the symptom rather than the cause. I think we need to be careful with calls to modelgrid. It looks like maybe the use of super() in set_data() (around 1518 of mfdatalist.py) or some other function passing means that when _check_valid_cellids() runs hasattr(self._model_or_sim, "modelgrid") it calls a new instance of _model_or_sim which requires a fresh setup of modelgrid which is 'very' expensive (when hit a number of times). The only part of modelgrid which is used in _check_valid_cellids() seems to be idomain so there is a lot of unnecessary overhead loading and building the entire modelgrid object. A hack would be to explicitly pass idomain through the set_data() methods (calling modelgrid only outside the stress period loop). This seems to reduce package object initialisation time by about 50%. The data loading is still pretty slow but the gains in there might be harder to come by/beyond my skills.

UPDATE: Unfortunately is seems not running the checks at initialisation time just shifts the time burden of reading in all the external datasets to package write time -- no free lunch (_update_size_defs()?). Do we need to be read the externals at write time? What are the implications of setting simulation_data.auto_set_sizes = False? Possibly max benefits might be to speed up that list reading though...

[briochh]

@spaulins-usgs, appreciate your help on this. It seems like I can cheat my way to efficiency by switching check_data off (through whichever hack works) and setting sim.simulation_data.auto_set_sizes = False on my parent simulation object.

I have put together a toy model and script here that can be used to explore the scalability of the methods. Also my horrible hacks are on https://github.com/briochh/flopy/tree/hotfix_listwrite . Although they are horrible the runtime for building a writing 70 stress periods of data for a package is reduced from 108 s to 3 s. (it just might break other stuff.)

[spaulins-usgs]

@briochh Your code that adds support for a file-by-file check looks good. Your code to minimize the number of modelgrids constructed will work, but I am going to implement my own solution for that. I plan to implement a "_get_modelgrid" method which supports modelgrid caching in certain cases. That way the code will be reusable in case other methods need to be optimized in the same way. I will try to put together a PR today.

[spaulins-usgs]

Regarding _update_size_defs, turning that off with simulation_data.auto_set_sizes = False is fine. _update_size_defs is just a convenience method that adjusts MODFLOW size variables to agree with the size of the arrays they are referring to. For example, if you added wells to your WEL stress period data but didn't update MAXBOUND to reflect the new wells, _update_size_defs would automatically update MAXBOUND for you. Setting auto_set_sizes to False just means that you can not count of flopy to update size variables, like MAXBOUND, for you.

[spaulins-usgs]

@briochh A fork with performance improvement changes including the ones you have suggested can be found here:

https://github.com/scottrp/flopy

I am still working to further optimize the code and will merge all my changes next week.

I modified how sim.simulation_data.verify_data works, when set to false it now disables all data verification checks for the entire simulation. For maximum performance I recommend turning off both auto_set_sizes and verify_data.

[briochh]

Looks good @spaulins-usgs. Thanks for digging into this. I think it'll help with the practical usability of Flopy for these larger problems.

[spaulins-usgs]

I added one more optimization today and merged all the changes into the main fork.

On my laptop your optimization test code initially took 119 seconds to run. After adding your code optimization and mine, and switching auto_set_sizes and verify_data to false, it took 65 seconds to run. You can reduce the run time further by not changing the simulation path, I got the time down to 58 seconds by doing so. Note that when you change the simulation path flopy may copy all the external files to the new path depending on what "ext_file_action" is set to when you call "write_simulation".