bug: failure to load mf6 test models

[wpbonelli]

Describe the bug
FloPy fails to load several MF6 models from the test models repository:

• test001a_Tharmonic_tabs

flopy.mf6.mfbase.MFDataException: An error occurred in data element "perioddata" package "flow15.tdis". The error occurred while loading data list from package file in the "read_list_data_from_file" method.
Additional Information:
(1) Data "perioddata" with value "1.00000000,1,1.00000000" can not be converted to float.
(2) Unable to process line 1 of data list: "                        1.00000000,1,1.00000000                                         Items:  PERLEN  NSTP    TSMULT

• test004_bcfss

flopy.mf6.mfbase.MFDataException: An error occurred in data element "icelltype" model "gwf_1" package "npf". The error occurred while reading data from file <truncated path>/bcf2ss.npf in the "_load_layer" method.

• test014_NWTP3Low_dev

flopy.mf6.mfbase.MFDataException: An error occurred in data element "dev_minimum_saturated_thickness" model "gwf_1" package "npf". The error occurred while converting data to string in the "load_from_package" method.
Additional Information:
(1) Data "dev_minimum_saturated_thickness" with value "1D-9" can not be converted to float.
(2) Could not convert "1D-9" of type "DatumType.double_precision" to a string

• test041_flowdivert_nwt_dev

flopy.mf6.mfbase.MFDataException: An error occurred in data element "dev_minimum_saturated_thickness" model "gwf_1" package "npf". The error occurred while converting data to string in the "load_from_package" method.
Additional Information:
(1) Data "dev_minimum_saturated_thickness" with value "1D-9" can not be converted to float.
(2) Could not convert "1D-9" of type "DatumType.double_precision" to a string.

To Reproduce

E.g. MFSimulation.load(sim_ws="path/to/simulation/directory")

[christianlangevin]

Nice find, @wpbonelli. Maybe @spaulins-usgs can have a look at this?

[spaulins-usgs]

@langevin-uses, of the four test failures above, "test001a_Tharmonic_tabs", "test014_NWTP3Low_dev", and "test041_flowdivert_nwt_dev" are fairly easy fixes.

However, with "test004_bcfss" there are comments at the end of the line which make things more difficult. Data like this:

internal
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  row 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  row 2

Is a challenge for flopy to deal with while trying to do an efficient fast load. Flopy currently splits the line by the comment marker (line.split("#")) to quickly separate the comments from the data. In the case above there are no comment characters separating the data from the comment. The only way I can think of to correctly process cases like this is to split each line by the delimiter (already done) and then loop through resulting list looking for non-numeric elements. Using a python loop to access and examine each individual data entry like this will most likely significantly slow down the loading process for large data files. My guess is it will be at least twice as slow.

One possible way around this is to load as we currently do, and if the load fails load it the slow way looking for inferred comments. There are however some technical issues with doing it this way. Trying to load the data a second time requires getting back to the original file location, but tell() and seek() are not 100% reliable for some of our data files. To get around this flopy could restart the load process of the entire file, but I believe this will have its own set of challenges.

The easiest and cleanest solution I can think of is to require that all input file comments must start with a comment character (#). The above example would change to:

internal
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  # row 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  # row 2

This would make it much easier for flopy to quickly determine what is data and what is comments and would also make the input files more human readable.

[christianlangevin]

Hey @spaulins-usgs, the following is the fortran code for mf6 for reading a 2d integer array. For a variable like icelltype, it is called multiple times for a variable like icelltype if the layered option is specified and the grid is structured.

  subroutine read_ascii(this)
    class(Integer2dReaderType) :: this
    integer(I4B) :: i, j
    integer(I4B) :: istat
    do i = 1, size(this%int2d, dim=2)
      read (this%input_unit, *, iostat=istat, iomsg=errmsg) &
        (this%int2d(j, i), j=1, size(this%int2d, dim=1))
    end do
    if (istat /= 0) then
      errmsg = 'Error reading data for array '//trim(this%array_name)// &
               '.  '//trim(errmsg)
      call store_error(errmsg)
      call store_error_unit(this%input_unit)
    end if
  end subroutine read_ascii

I guess the point here is that MODFLOW reads a row at a time, which allows comments to be placed at the end of a row as long as the number of data columns is equal to the number of model columns. I'm wondering if there is a way to duplicate this concept in FloPy? Perhaps the innards of FloPy couldn't be easily refactored for this. I could be wrong, but it seems to me like the older flopy code could handle this type of array input with comments at the end of a line. @mjreno implemented this in the IDM refactoring for mf6 and might have thoughts.

If it's too complicated to refactor flopy, then we can look at other options. But I guess that means that users may have valid mf6 models that run fine but cannot be loaded by flopy.

[spaulins-usgs]

@langevin-usgs, I think I understand now. I can change flopy's behavior to be like MODFLOW. Though, in order to reproduce MODFLOW's behavior flopy will always need to know the dimensions of array data at the time of loading, which currently is not always the case. One case I know of where flopy does not know the dimensions of array data at load time is the time array series "tas-array" data. The reason for this is that flopy loads referenced packages immediately after reading the "FILEIN" line, meaning that flopy loads the tas package before finishing loading the package that references it. Therefore at the time of loading flopy does not know what data the "tas_array" is for. I can refactor this so that the reference package loads before tas and then add some code to resolve the data shape for "tas-array". Note that this will be the reverse order that ts loads in, since the referencing package needs to know the ts attribute names before loading.

Anyway, this is definitely doable but will require a little work.

[christianlangevin]

The definition files have a "layered" attribute as well as "shape", which is what we use in the IDM loading. Seems like you should be able to do something similar so loading behavior between flopy and mf6 is the same.

[spaulins-usgs]

@wpbonelli, each of the test model load failures have been address and flopy can now load the models described above. The fixes allow flopy to better handle comma delimited text with space delimited comments, translate floating point notation that uses a "D" instead of an "E" (1D-4), and read in array data based on the dimensions of the array (only read data entries up to the "number of columns" from each line).

[wpbonelli]

The 4 original models load successfully. There is a different one mf6_test019_VilhelmsenLGR[_nr] failing. Also test014_NWTP3Low_dev loads but fails to converge

[christianlangevin]

I think I see what is wrong with mf6_test019_VilhelmsenLGR[_nr] and am working on a fix for flopy. Looks like this commit broke the idomain loading for that particular problem. Will also check on test014_NWTP3Low_dev.