modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
517 stars 313 forks source link

loading a MODFLOW-NWT model, resulting flopy version gives slightly different results #89

Closed hwreeves-USGS closed 8 years ago

hwreeves-USGS commented 8 years ago

I loaded a MODFLOW-NWT model using flopy from the 'develop' branch. Mike and Andy found some problems and got it to where the model loads. I copied the .nwt file from my original location to the new model working directory and then ran the modflow executable.

two issues: (1) Mike is already changing, needed to add the FREE option to the .bas file written by flopy because the CHD and WEL files created were not appropriate for fixed-format reading.

(2) the flopy version generates a CHD file with values that are not the same as the original CHD file due to rounding errors. I haven't checked all the other input files.

The list file shows them as exactly the same, but the CHD differs, for example:

Original:

MODFLOW2000 Constant-Head Boundary Package (CHD)

  6551
  6551         0
     1       536       266  1040.000  1040.000
     1       537       266  1040.000  1040.000
     1       538       267  1039.900  1039.900
     1       539       267  1040.000  1040.000
     1       540       268  1040.000  1040.000
     1       541       268  1040.000  1040.000
     1       542       268  1040.400  1040.400
     1       543       268  1040.200  1040.200
     1       544       268  1040.400  1040.400
     1       545       268  1040.800  1040.800

flopy version that is created after loading the original model and then writing to the working directory:

CHD for MODFLOW, generated by Flopy.

  6551
  6551         0 # stress period 0
     1       536       266 1040.000000 1040.000000
     1       537       266 1040.000000 1040.000000
     1       538       267 1039.900024 1039.900024
     1       539       267 1040.000000 1040.000000
     1       540       268 1040.000000 1040.000000
     1       541       268 1040.000000 1040.000000
     1       542       268 1040.400024 1040.400024
     1       543       268 1040.199951 1040.199951
     1       544       268 1040.400024 1040.400024
     1       545       268 1040.800049 1040.800049

so when you run the model, the small differences resulted in a slightly different mass balance in the end. Maybe not a big deal overall, but frustrating for archiving or reproducing a model using flopy- we'd like to load a model and then use the flopy tools to explore it and not have to worry about these differences. I didn't try to run the model after loading within the python notebook as that was failing with the earlier versions, these results were running MODFLOW-NWT_64 from a shell, so it may be you can load the original version and run it without a problem, the only issue will be if the user decides to write out the loaded model.

hwreeves-USGS commented 8 years ago

flopy can't run a model that hasn't been written out - so this is an issue for loading and running in a new workspace.

jtwhite79 commented 8 years ago

Issue 1.) is a bug - when an old style control record model (e.g. not using internal, external, open/close) is loaded, the original unit numbers are retained in the Util2d.locat values, even though flopy is resetting the package unit numbers arbitrarily. So if the array is written in the package file (internally), the wrong (original) unit number is being written to the new file. The easy fix is to set the .array_free_format flag after the load.

I'll work on this...probably going to be painful to fix it the right way...we really need a better system for handling unit numbers/name file entries.

jtwhite79 commented 8 years ago

I just pushed several mods and bug fixes to develop - I'm now able to run the LPR model with just standard flopy workflow steps. I also attempted to address the float point garbage being written the BC packages - let me know if that is mo' better.

hwreeves-USGS commented 8 years ago

I ran the latest 'develop' branch and the output is not exactly the same as the original files I posted for you - the problem appears to be in the CHD and WEL files: they are rounding 3-digit input to 1 or 2-digits in the new versions of the files generated with ml.write_files()

I'll send the two versions to you Jeremy.

HWR

jtwhite79 commented 8 years ago

Alright Howard, I've made one last attempt to get this working for you - please pull the latest develop branch commits. If you load your model with fixed format lists, set the ModflowBas.ifrefm flag to free, then write the model input files, the resulting list file water budget exactly matches the list water budget of the original input files - at least for me.

There is still a lingering issue with loading fixed format then writing fixed format because of the way numpy represents the floats as strings with the format specification minilangage...I'm not inclined to do more on that front. So for now, the best option to reproduce the results of your model exactly is to load as fixed then reset to free..plus this has the added benefit of making the input files more human readable.

hwreeves-USGS commented 8 years ago

Jeremy, this set of changes to develop worked for me and the steady-state Little Plover model. I added ml.set_ifrefm(True) to my notebook. I did not need ml.free_format=True. The load, change working directory, then write and run MODFLOW-NWT gave the same results as the original model.