performance of loading tabular mf6 datasets

consider this model, which has a well package with 1M fluxes:

import time
import cProfile
import numpy as np
import pandas as pd
import flopy

model_name = 'junk'
workspace = 'tmp'
nper = 100

sim = flopy.mf6.MFSimulation(sim_name=model_name, version='mf6', exe_name='mf6', 
                             sim_ws=workspace)
tdis_rc = [(1.0, 1, 1.0)] * nper
tdis = flopy.mf6.ModflowTdis(sim, 
                             nper=nper, perioddata=tdis_rc)
gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name,
                           model_nam_file='{}.nam'.format(model_name))
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=1, nrow=2, ncol=2)
nwells = 10000
spd = {}
for kper in range(100):
    kspd = flopy.mf6.ModflowGwfwel.stress_period_data.empty(gwf,
                                                      maxbound=nwells,
                                                      boundnames=True)[0]
    kspd['cellid'] = [(0, 0, 0)] * nwells
    kspd['q'] = np.random.randn(nwells)
    kspd['boundnames'] = 'well'
    spd[kper] = kspd
wel = flopy.mf6.ModflowGwfwel(gwf, pname='wel', #maxbound=nwells,
                              stress_period_data=spd,
                              save_flows=True)
sim.write_simulation()

if you profile loading the model back in, you can see that almost the entire time (50 sec on my machine) is spent on the _load_list_line (1M calls) and _append_data_list (2M calls) methods in mffileaccess.py. (the load actually errors out, but the time is similar to what I got for a real-life example, suggesting that it was almost done).

pr = cProfile.Profile()
pr.enable()
flopy.mf6.MFSimulation.load('junk', sim_ws='tmp')
pr.disable()
pr.print_stats(sort='cumtime')

50 sec might seem reasonable for 1 million lines, but I can load the same data with pandas in about a half second.

import time
import pandas as pd
t0 = time.time()
df = pd.read_csv('tmp/junk.wel', skiprows=9,
            header=None,
            delim_whitespace=True, error_bad_lines=False).dropna(axis=0)
print("finished in {:.2f}s\n".format(time.time() - t0))
finished in 0.44s

In the real model that I am working with, the majority of the load time seems to be consumed by the well package and perhaps SFR, suggesting that speeding up the loading of list-like data could substantially cut the overall load time for transient regional models with lots of wells.

I realize that my pandas example is an overly simple kludge (and that there may not be an easy fix to this issue for all use cases), but I think the greater issue of performance is something to think about. _load_list_line and _append_data_list have quite a bit of complexity (they are respectively about 300 and 100 lines long, not including calls to other methods), and perhaps for good reason. But is this really needed for most uses? Isn't one of the paradigms of MODFLOW-6 supposed to be simple, clean input structures that are easily readable?

[spaulins-usgs]

I did some limited optimization of this code about a year or two ago. I ran into several bottlenecks when optimizing the code. The main problem was variable amounts of data and types of data on each line. While the well package is one of the simpler cases, it can still have 1, 2, or 3 integers for the cellid and possible optional information on the end of each line (aux variables and/or boundname). There is therefore some logic required, even for the simple packages, to determine how to store the data being read from the file. In the case of the simple packages there are some ways I can further optimize the code in cases where I can assume a consistent amount and type of data on every line. Though note that it is unlikely to be near the speed of pandas since pandas processes the data using mostly compiled code while the flopy code is still likely to be bottlenecked by interpreted code. I will try a few things to optimize some special cases like your example case and post a followup message once I have some results.

For the advanced packages I don't have any good ideas on how to improve performance. For example, MAW's period data has variable length and has variable types depending on each line's mawsetting. In this case data is converted and stored differently depending on the mawsetting keyword. Last time I profiled this code the main bottleneck was due to having to re-run the interpreting and processing code on each line because flopy could not assume every line had the same amount or type of data.

[aleaf]

With the Well package at least, we should know how many items to expect on each line from the options block (where aux variables and the boundnames option are specified), and from the type of DIS package (determining the number of ints for cellid). With that information, we should be able to do something like this:

t0 = time.time()
data = []
# define the expected columns up front
columns = ['per', 'k', 'i', 'j', 'q']
cellid_cols = ['k', 'i', 'j']

with open('tmp/junk.wel') as src:
    for line in src:
        if 'begin period' in line.lower():
            per = int(line.strip().split()[-1])
            for line in src:
                if 'end period' in line.lower():
                    break
                items = line.strip().split()
                data.append([per] + items)
df = pd.DataFrame(data, columns=columns)
for c in cellid_cols:
    df[c] = df[c].astype(int)
df['q'] = df['q'].astype(float)

print("finished in {:.2f}s\n".format(time.time() - t0))
finished in 3.17s

Quite a bit slower than read_csv, but still more than an order of magnitude faster than the current parser. For cases that need more specialization, maybe we should dispatch to different parser functions that work as efficiently as practical for their problems. Simple parser functions that focus on a particular problem might also be easier to test and maintain.

[christianlangevin]

@aleaf-usgs, I think it's great that you are always making sure that the code works well for real models! Having slow loads is really a bummer, and I agree that we want to be efficient. One of our biggest issues is that MODFLOW will read and run fine with certain files that cannot be loaded with flopy unless we do extra processing line by line. For example, the following wel package is perfectly valid for MODFLOW 6 (notice the comma, #, blank lines, and 1.d0):

BEGIN options
  auxiliary  CONCENTRATION
END options

BEGIN dimensions
  MAXBOUND  40
END dimensions

BEGIN period  1
  1 1, 1       0.14256000       0.00000000
#  2 1 1       0.14256000       0.00000000


  3 1 1       0.14256000       1.d0
  4 1 1       0.14256000       0.00000000
END period

I don't think this would load with your version here, which means that we need a line parser to handle these special cases. Not sure what the answer is here. I suppose we could try a fast load first wrapped in a try/except and then the slower more rigorous load in the except? In my opinion we want to do the following in this order of priority:

1. Load any valid file that works with MODFLOW
2. Do (1) as fast as possible

My 2 cents.

[jwhite-usgs]

This will probably get me in trouble, but maybe the MODFLOW 6 input format could be modified to be more amenable to using industry-standard ASCII parsing tools like pandas? Im not at all familiar with MAW but if MF6 could use the same number of entries for each row in list-directed input, that would make things a lot easier for everyone...

[aleaf]

@langevin-usgs @spaulins-usgs I suspect some of the special cases like the comment and blank lines above could be handled with at speed about as good as the one above. In lieu of what @jwhite-usgs mentioned, I would suggest dispatching to optimized parsers first based on the expected input, and then falling back to the "god" parser for the most difficult situations. I think this would cover the majority of use cases, and we could say explicitly in the documentation that if you put certain things in your input files, you can expect a slower load time.

[spaulins-usgs]

@aleaf-usgs I have already written a simpler optimized version of code but it apparently is not getting called any more. The code is not as simple as your code since it needs to still be somewhat robust as @langevin-usgs mentioned and it also needs to do things like decrement cellid integers to make them zero based (by flopy's convention), but from what I remember it is still faster than the more complex code. Once I'm done fixing the deepcopy issue I will start troubleshooting why this code isn't getting called for the well package. In case you are curious the simpler code starts on line 939 of mffileaccess.py.

[spaulins-usgs]

@aleaf-usgs I modified the optimized version of the code to work with an expanded set of cases including your well example above. Additionally I added a verify_data flag you can set when calling the MFSimulation.load method. With verify_data off flopy will not do any verification checks on the data which further speeds up the load. On my computer your example loaded over three times faster when using the optimized code with verify_data set to False. The optimized code has been checked in and will automatically run for simple cases. In general the optimized code will run for the basic packages that do not reference time-series files.

[aleaf]

thanks @spaulins-usgs. When I run the example above with these changes (assuming they are all incorporated in 06343a7), it takes about 20 sec with verify_data=True and 14 sec with verify_data=False. Certainly an improvement over 50 sec., but to be honest still a little unsatisfying considering the speed of the basic python example above. I realize though, that there is probably only so much that can practically be done at the moment without major changes to the flopy code structure or MODFLOW-6 input.

[jdhughes-dev]

@spaulins-usgs I think this can be closed reopen if that is not the case.