examples/Notebooks/flopy3_mf6_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flopy MODFLOW 6 (MF6) Support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Flopy library contains classes for creating, saving, running, loading, and modifying MF6 simulations.  The MF6 portion of the flopy library is located in:\n",
    "\n",
    "*flopy.mf6*\n",
    "\n",
    "While there are a number of classes in flopy.mf6, to get started you only need to use the main classes summarized below:\n",
    "\n",
    "flopy.mf6.MFSimulation  \n",
    "* MODFLOW Simulation Class.  Entry point into any MODFLOW simulation.\n",
    "\n",
    "flopy.mf6.ModflowGwf\n",
    "* MODFLOW Groundwater Flow Model Class.  Represents a single model in a simulation.\n",
    "\n",
    "flopy.mf6.Modflow[pc]\n",
    "    * MODFLOW package classes where [pc] is the abbreviation of the package name.  Each package is a separate class.  \n",
    "\n",
    "For packages that are part of a groundwater flow model, the abbreviation begins with \"Gwf\".  For example, \"flopy.mf6.ModflowGwfdis\" is the Discretization package.\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "from shutil import copyfile\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "try:\n",
    "    import flopy\n",
    "except:\n",
    "    fpth = os.path.abspath(os.path.join('..', '..'))\n",
    "    sys.path.append(fpth)\n",
    "    import flopy\n",
    "    \n",
    "print(sys.version)\n",
    "print('numpy version: {}'.format(np.__version__))\n",
    "print('matplotlib version: {}'.format(mpl.__version__))\n",
    "print('flopy version: {}'.format(flopy.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a MF6 Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A MF6 simulation is created by first creating a simulation object \"MFSimulation\".  When you create the simulation object you can define the simulation's name, version, executable name, workspace path, and the name of the tdis file.  All of these are optional parameters, and if not defined each one will default to the following:\n",
    "\n",
    "sim_name='modflowtest'\n",
    "\n",
    "version='mf6'\n",
    "\n",
    "exe_name='mf6.exe'\n",
    "\n",
    "sim_ws='.'\n",
    "\n",
    "sim_tdis_file='modflow6.tdis'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "from shutil import copyfile\n",
    "try:\n",
    "    import flopy\n",
    "except:\n",
    "    fpth = os.path.abspath(os.path.join('..', '..'))\n",
    "    sys.path.append(fpth)\n",
    "    import flopy\n",
    "\n",
    "sim_name = 'example_sim'\n",
    "sim_path = os.path.join('data', 'example_project')\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, version='mf6', exe_name='mf6', \n",
    "                             sim_ws=sim_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create a tdis package object \"ModflowTdis\".  The first parameter of the ModflowTdis class is a simulation object, which ties a ModflowTdis object to a specific simulation.  The other parameters and their definitions can be found in the docstrings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis = flopy.mf6.ModflowTdis(sim, pname='tdis', time_units='DAYS', nper=2, \n",
    "                             perioddata=[(1.0, 1, 1.0), (10.0, 5, 1.0)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next one or more models are created using the ModflowGwf class.  The first parameter of the ModflowGwf class is the simulation object that the model will be a part of."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_name = 'example_model'\n",
    "model = flopy.mf6.ModflowGwf(sim, modelname=model_name,\n",
    "                             model_nam_file='{}.nam'.format(model_name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next create one or more Iterative Model Solution (IMS) files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ims_package = flopy.mf6.ModflowIms(sim, pname='ims', print_option='ALL',\n",
    "                                   complexity='SIMPLE', outer_hclose=0.00001,\n",
    "                                   outer_maximum=50, under_relaxation='NONE',\n",
    "                                   inner_maximum=30, inner_hclose=0.00001,\n",
    "                                   linear_acceleration='CG',\n",
    "                                   preconditioner_levels=7,\n",
    "                                   preconditioner_drop_tolerance=0.01,\n",
    "                                   number_orthogonalizations=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each ModflowGwf object needs to be associated with an ModflowIms object.  This is done by calling the MFSimulation object's \"register_ims_package\" method.  The first parameter in this method is the ModflowIms object and the second parameter is a list of model names (strings) for the models to be associated with the ModflowIms object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.register_ims_package(ims_package, [model_name])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next add packages to each model.  The first package added needs to be a spatial discretization package since flopy uses  information from the spatial discretization package to help you build other packages.  There are three spatial discretization packages to choose from:\n",
    "\n",
    "DIS (ModflowGwfDis) - Structured discretization\n",
    "DISV (ModflowGwfdisv) - Discretization with vertices\n",
    "DISU (ModflowGwfdisu) - Unstructured discretization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dis_package = flopy.mf6.ModflowGwfdis(model, pname='dis', length_units='FEET', nlay=2,\n",
    "                                      nrow=2, ncol=5, delr=500.0,\n",
    "                                      delc=500.0,\n",
    "                                      top=100.0, botm=[50.0, 20.0],\n",
    "                                      filename='{}.dis'.format(model_name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing Namefiles\n",
    "\n",
    "Namefiles are automatically built for you by flopy.  However, there are some options contained in the namefiles that you may want to set. To get the namefile object access the name_file attribute in either a simulation or model object to get the simulation or model namefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the nocheck property in the simulation namefile\n",
    "sim.name_file.nocheck = True\n",
    "# set the print_input option in the model namefile\n",
    "model.name_file.print_input = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying Options\n",
    "\n",
    "Option that appear alone are assigned a boolean value, like the print_input option above.  Options that have additional optional parameters are assigned using a tuple, with the entries containing the names of the optional parameters to turn on.  Use a tuple with an empty string to indicate no optional parameters and use a tuple with None to turn the option off. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn Newton option on with under relaxation \n",
    "model.name_file.newtonoptions = ('UNDER_RELAXATION')\n",
    "# Turn Newton option on without under relaxation\n",
    "model.name_file.newtonoptions = ('')\n",
    "# Turn off Newton option \n",
    "model.name_file.newtonoptions = (None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MFArray Templates\n",
    "\n",
    "Lastly define all other packages needed.  \n",
    "\n",
    "Note that flopy supports a number of ways to specify data for a package.  A template, which defines the data array shape for you, can be used to specify the data.  Templates are built by calling the empty of the data type you are building.  For example, to build a template for k in the npf package you would call:\n",
    "\n",
    "ModflowGwfnpf.k.empty()\n",
    "\n",
    "The empty method for \"MFArray\" data templates (data templates whose size is based on the structure of the model grid)  take up to four parameters:\n",
    "\n",
    "* model - The model object that the data is a part of.  A valid model object with a discretization package is required in order to build the proper array dimensions.  This parameter is required.\n",
    "\n",
    "* layered - True or false whether the data is layered or not.\n",
    "\n",
    "* data_storage_type_list - List of data storage types, one for each model layer.  If the template is not layered, only one data storage type needs to be specified.  There are three data storage types supported, internal_array, internal_constant, and external_file. \n",
    "\n",
    "* default_value - The initial value for the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build a data template for k that stores the first layer as an internal array and the second\n",
    "# layer as a constant with the default value of k for all layers set to 100.0 \n",
    "layer_storage_types = [flopy.mf6.data.mfdatastorage.DataStorageType.internal_array, \n",
    "                       flopy.mf6.data.mfdatastorage.DataStorageType.internal_constant]\n",
    "k_template = flopy.mf6.ModflowGwfnpf.k.empty(model, True, layer_storage_types, 100.0)\n",
    "# change the value of the second layer to 50.0\n",
    "k_template[0]['data'] = [65.0, 60.0, 55.0, 50.0, 45.0, 40.0, 35.0, 30.0, 25.0, 20.0]\n",
    "k_template[0]['factor'] = 1.5\n",
    "print(k_template)\n",
    "# create npf package using the k template to define k\n",
    "npf_package = flopy.mf6.ModflowGwfnpf(model, pname='npf', save_flows=True, icelltype=1, k=k_template)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying MFArray Data\n",
    "\n",
    "MFArray data can also be specified as a numpy array, a list of values, or a single value.  Below strt (starting heads) are defined as a single value, 100.0, which is interpreted as an internal constant storage type of value 100.0.  Strt could also be defined as a list defining a value for every model cell:\n",
    "\n",
    "strt=[100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, \n",
    "      90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0]\n",
    "      \n",
    "Or as a list defining a value or values for each model layer:\n",
    "\n",
    "strt=[100.0, 90.0]\n",
    "\n",
    "or:\n",
    "\n",
    "strt=[[100.0], [90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0]]\n",
    "\n",
    "MFArray data can also be stored in an external file by using a dictionary using the keys 'filename' to specify the file name relative to the model folder and 'data' to specific the data.  The optional 'factor', 'iprn', and 'binary' keys may also be used.\n",
    "\n",
    "strt={'filename': 'strt.txt', 'factor':1.0, 'data':[100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, \n",
    "      90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0], 'binary': 'True'}\n",
    "      \n",
    "If the 'data' key is omitted from the dictionary flopy will try to read the data from an existing file 'filename'.  Any relative paths for loading data from a file should specified relative to the MF6 simulation folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strt={'filename': 'strt.txt', 'factor':1.0, 'data':[100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, \n",
    "      90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0, 90.0], 'binary': 'True'}\n",
    "ic_package = flopy.mf6.ModflowGwfic(model, pname='ic', strt=strt,\n",
    "                                    filename='{}.ic'.format(model_name))\n",
    "# move external file data into model folder\n",
    "icv_data_path = os.path.join('..', 'data', 'mf6', 'notebooks', 'iconvert.txt')\n",
    "copyfile(icv_data_path, os.path.join(sim_path, 'iconvert.txt'))\n",
    "# create storage package\n",
    "sto_package = flopy.mf6.ModflowGwfsto(model, pname='sto', save_flows=True, iconvert={'filename':'iconvert.txt'},\n",
    "                                      ss=[0.000001, 0.000002], \n",
    "                                      sy=[0.15, 0.14, 0.13, 0.12, 0.11, 0.11, 0.12, 0.13, 0.14, 0.15,\n",
    "                                         0.15, 0.14, 0.13, 0.12, 0.11, 0.11, 0.12, 0.13, 0.14, 0.15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MFList Templates\n",
    "\n",
    "Flopy supports specifying record and recarray \"MFList\" data in a number of ways.  Templates can be created that define the shape of the data.  The empty method for \"MFList\" data templates take up to 7 parameters.\n",
    "\n",
    "* model - The model object that the data is a part of.  A valid model object with a discretization package is required in order to build the proper array dimensions.  This parameter is required.\n",
    "\n",
    "* maxbound - The number of rows in the recarray.  If not specified one row is returned.\n",
    "\n",
    "* aux_vars - List of auxiliary variable names.  If not specified auxiliary variables are not used.\n",
    "\n",
    "* boundnames - True/False if boundnames is to be used.\n",
    "\n",
    "* nseg - Number of segments (only relevant for a few data types)\n",
    "\n",
    "* timeseries - True/False indicates that time series data will be used.\n",
    "\n",
    "* stress_periods - List of integer stress periods to be used (transient MFList data only).  If not specified for transient data, template will only be defined for stress period 1. \n",
    "\n",
    "MFList transient data templates are numpy recarrays stored in a dictionary with the dictionary key an integer zero based stress period value (stress period - 1).\n",
    "\n",
    "In the code below the well package is set up using a transient MFList template to help build the well's stress_periods.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxbound = 2\n",
    "# build a stress_period_data template with 2 wells over stress periods 1 and 2 with boundnames \n",
    "# and three aux variables\n",
    "wel_periodrec = flopy.mf6.ModflowGwfwel.stress_period_data.empty(model, maxbound=maxbound, boundnames=True, \n",
    "                                                             aux_vars=['var1', 'var2', 'var3'],\n",
    "                                                             stress_periods=[0,1])\n",
    "# define the two wells for stress period one\n",
    "wel_periodrec[0][0] = ((0,1,2), -50.0, -1, -2, -3, 'First Well')\n",
    "wel_periodrec[0][1] = ((1,1,4), -25.0, 2, 3, 4, 'Second Well')\n",
    "# define the two wells for stress period two\n",
    "wel_periodrec[1][0] = ((0,1,2), -200.0, -1, -2, -3, 'First Well')\n",
    "wel_periodrec[1][1] = ((1,1,4), -4000.0, 2, 3, 4, 'Second Well')\n",
    "# build the well package\n",
    "wel_package = flopy.mf6.ModflowGwfwel(model, pname='wel', print_input=True, print_flows=True,\n",
    "                                      auxiliary=['var1', 'var2', 'var3'], maxbound=maxbound,\n",
    "                                      stress_period_data=wel_periodrec, boundnames=True, save_flows=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cell IDs\n",
    "\n",
    "Cell IDs always appear as tuples in an MFList. For a structured grid cell IDs appear as:\n",
    "\n",
    "(&lt;layer&gt;, &lt;row&gt;, &lt;column&gt;)\n",
    "\n",
    "For vertice based grid cells IDs appear as:\n",
    "\n",
    "(&lt;layer&gt;, &lt;intralayer_cell_id&gt;)\n",
    "\n",
    "Unstructured grid cell IDs appear as:\n",
    "\n",
    "(&lt;cell_id&gt;)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying MFList Data\n",
    "\n",
    "MFList data can also be defined as a list of tuples, with each tuple being a row of the recarray.  For transient data the list of tuples can be stored in a dictionary with the dictionary key an integer zero based stress period value.  If only a list of tuples is specified for transient data, the data is assumed to apply to stress period 1.  Additional stress periods can be added with the add_transient_key method.  The code below defines saverecord and printrecord as a list of tuples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# printrecord data as a list of tuples.  since no stress\n",
    "# period is specified it will default to stress period 1\n",
    "printrec_tuple_list = [('HEAD', 'ALL'), ('BUDGET', 'ALL')]\n",
    "# saverecord data as a dictionary of lists of tuples for \n",
    "# stress periods 1 and 2.  \n",
    "saverec_dict = {0:[('HEAD', 'ALL'), ('BUDGET', 'ALL')],1:[('HEAD', 'ALL'), ('BUDGET', 'ALL')]}\n",
    "# create oc package\n",
    "oc_package = flopy.mf6.ModflowGwfoc(model, pname='oc', \n",
    "                                    budget_filerecord=[('{}.cbc'.format(model_name),)],\n",
    "                                    head_filerecord=[('{}.hds'.format(model_name),)],\n",
    "                                    saverecord=saverec_dict,\n",
    "                                    printrecord=printrec_tuple_list)\n",
    "# add stress period two to the print record\n",
    "oc_package.printrecord.add_transient_key(1)\n",
    "# set the data for stress period two in the print record\n",
    "oc_package.printrecord.set_data([('HEAD', 'ALL'), ('BUDGET', 'ALL')], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying MFList Data in an External File \n",
    "\n",
    "MFList data can be specified in an external file using a dictionary with the 'filename' key.  If the 'data' key is also included in the dictionary and is not None, flopy will create the file with the data contained in the 'data' key.  The code below creates a chd package which creates and references an external file containing data for stress period 1 and stores the data internally in the chd package file for stress period 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stress_period_data = {0: {'filename': 'chd_sp1.dat', 'data': [[(0, 0, 0), 70.]]},\n",
    "                      1: [[(0, 0, 0), 60.]]}\n",
    "chd = flopy.mf6.ModflowGwfchd(model, maxbound=1, stress_period_data=stress_period_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages that Support both List-based and Array-based Data\n",
    "\n",
    "The recharge and evapotranspiration packages can be specified using list-based or array-based input.  The array packages have an \"a\" on the end of their name:\n",
    "\n",
    "ModflowGwfrch - list based recharge package\n",
    "ModflowGwfrcha - array based recharge package\n",
    "ModflowGwfevt - list based evapotranspiration package\n",
    "ModflowGwfevta - array based evapotranspiration package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rch_recarray = {0:[((0,0,0), 'rch_1'), ((1,1,1), 'rch_2')],\n",
    "                1:[((0,0,0), 'rch_1'), ((1,1,1), 'rch_2')]}\n",
    "rch_package = flopy.mf6.ModflowGwfrch(model, pname='rch', fixed_cell=True, print_input=True, \n",
    "                                      maxbound=2, stress_period_data=rch_recarray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility Files (TS, TAS, OBS, TAB)\n",
    "\n",
    "Utility files, MF6 formatted files that reference by packages, include time series, time array series, observation, and tab files.  The file names for utility files are specified using the package that references them. The utility files can be created in several ways. A simple case is demonstrated below. More detail is given in the flopy3_mf6_obs_ts_tas notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build a time series array for the recharge package\n",
    "ts_data = [(0.0, 0.015, 0.0017), (1.0, 0.016, 0.0019), (2.0, 0.012, 0.0015),\n",
    "           (3.0, 0.020, 0.0014), (4.0, 0.015, 0.0021), (5.0, 0.013, 0.0012),\n",
    "           (6.0, 0.022, 0.0012), (7.0, 0.016, 0.0014), (8.0, 0.013, 0.0011),\n",
    "           (9.0, 0.021, 0.0011), (10.0, 0.017, 0.0016), (11.0, 0.012, 0.0015)]\n",
    "rch_package.ts.initialize(time_series_namerecord=['rch_1', 'rch_2'],                                     \n",
    "                          timeseries=ts_data, filename='recharge_rates.ts',\n",
    "                          interpolation_methodrecord=['stepwise', 'stepwise'])\n",
    "\n",
    "# build an recharge observation package that outputs the western recharge to a binary file and the eastern\n",
    "# recharge to a text file\n",
    "obs_data = {('rch_west.csv', 'binary'): [('rch_1_1_1', 'RCH', (0, 0, 0)),\n",
    "                                    ('rch_1_2_1', 'RCH', (0, 1, 0))],\n",
    "             'rch_east.csv': [('rch_1_1_5', 'RCH', (0, 0, 4)),\n",
    "                              ('rch_1_2_5', 'RCH', (0, 1, 4))]}\n",
    "rch_package.obs.initialize(filename='example_model.rch.obs', digits=10, \n",
    "                           print_input=True, continuous=obs_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving and Running a MF6 Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Saving and running a simulation are done with the MFSimulation class's write_simulation and run_simulation methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write simulation to new location\n",
    "sim.write_simulation()\n",
    "\n",
    "# run simulation\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exporting a MF6 Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exporting a MF6 model to a shapefile or netcdf is the same as exporting a MF2005 model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make directory\n",
    "pth = os.path.join('data', 'netCDF_export')\n",
    "if not os.path.exists(pth):\n",
    "    os.makedirs(pth)\n",
    "    \n",
    "# export the dis package to a netcdf file\n",
    "model.dis.export(os.path.join(pth, 'dis.nc'))\n",
    "\n",
    "# export the botm array to a shapefile\n",
    "model.dis.botm.export(os.path.join(pth, 'botm.shp'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading an Existing MF6 Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Loading a simulation can be done with the flopy.mf6.MFSimulation.load static method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the simulation\n",
    "loaded_sim = flopy.mf6.MFSimulation.load(sim_name, 'mf6', 'mf6', sim_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Retrieving Data and Modifying an Existing MF6 Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data can be easily retrieved from a simulation.  Data can be retrieved using two methods.  One method is to retrieve the data object from a master simulation dictionary that keeps track of all the data.  The master simulation dictionary is accessed by accessing a simulation's \"simulation_data\" property and then the \"mfdata\" property:\n",
    "\n",
    "sim.simulation_data.mfdata[<data path>]\n",
    "\n",
    "The data path is the path to the data stored as a tuple containing the model name, package name, block name, and data name.\n",
    "\n",
    "The second method is to get the data from the package object.  If you do not already have the package object, you can work your way down the simulation structure, from the simulation to the correct model, to the correct package, and finally to the data object.\n",
    "\n",
    "These methods are demonstrated in the code below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get hydraulic conductivity data object from the data dictionary\n",
    "hk = sim.simulation_data.mfdata[(model_name, 'npf', 'griddata', 'k')]\n",
    "\n",
    "# get specific yield data object from the storage package\n",
    "sy = sto_package.sy\n",
    "\n",
    "# get the model object from the simulation object using the get_model method, \n",
    "# which takes a string with the model's name and returns the model object\n",
    "mdl = sim.get_model(model_name)\n",
    "# get the package object from the model mobject using the get_package method,\n",
    "# which takes a string with the package's name or type\n",
    "ic = mdl.get_package('ic')\n",
    "# get the data object from the initial condition package object\n",
    "strt = ic.strt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have the appropriate data object there are a number methods to retrieve data from that object.  Data  retrieved can either be the data as it appears in the model file or the data with any factor specified in the model file applied to it.  To get the raw data without applying a factor use the get_data method.  To get the data with the factor already applied use .array.\n",
    "\n",
    "Note that MFArray data is always a copy of the data stored by flopy.  Modifying the copy of the flopy data will have no affect on the data stored in flopy.  Non-constant internal MFList data is returned as a reference to a numpy recarray.  Modifying this recarray will modify the data stored in flopy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the data without applying any factor\n",
    "hk_data_no_factor = hk.get_data()\n",
    "print('Data without factor:\\n{}\\n'.format(hk_data_no_factor))\n",
    "\n",
    "# get data with factor applied\n",
    "hk_data_factor = hk.array\n",
    "print('Data with factor:\\n{}\\n'.format(hk_data_factor))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data can also be retrieved from the data object using []. For unlayered data the [] can be used to slice the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# slice layer one row two\n",
    "print('SY slice of layer on row two\\n{}\\n'.format(sy[0,:,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For layered data specify the layer number within the brackets.  This will return a \"LayerStorage\" object which let's you change attributes of an individual layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get layer one LayerStorage object\n",
    "hk_layer_one = hk[0]\n",
    "# change the print code and factor for layer one\n",
    "hk_layer_one.iprn = '2'\n",
    "hk_layer_one.factor = 1.1\n",
    "print('Layer one data without factor:\\n{}\\n'.format(hk_layer_one.get_data()))\n",
    "print('Data with new factor:\\n{}\\n'.format(hk.array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying Data\n",
    "\n",
    "Data can be modified in several ways.  One way is to set data for a given layer within a LayerStorage object, like the one accessed in the code above.  Another way is to set the data attribute to the new data.  Yet another way is to call the data object's set_data method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set data within a LayerStorage object\n",
    "hk_layer_one.set_data([120.0, 100.0, 80.0, 70.0, 60.0, 50.0, 40.0, 30.0, 25.0, 20.0])\n",
    "print('New HK data no factor:\\n{}\\n'.format(hk.get_data()))\n",
    "# set data attribute to new data\n",
    "ic_package.strt = 150.0\n",
    "print('New strt values:\\n{}\\n'.format(ic_package.strt.array))\n",
    "# call set_data\n",
    "sto_package.ss.set_data([0.000003, 0.000004])\n",
    "print('New ss values:\\n{}\\n'.format(sto_package.ss.array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying the Simulation Path\n",
    "\n",
    "The simulation path folder can be changed by using the set_sim_path method in the MFFileMgmt object. The MFFileMgmt object can be obtained from the simulation object through properties:\n",
    "\n",
    "sim.simulation_data.mfpath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# create new path\n",
    "save_folder = os.path.join(sim_path, 'sim_modified')\n",
    "# change simulation path\n",
    "sim.simulation_data.mfpath.set_sim_path(save_folder)\n",
    "# create folder\n",
    "if not os.path.isdir(save_folder):\n",
    "    os.makedirs(save_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a Model Relative Path\n",
    "\n",
    "A model relative path lets you put all of the files associated with a model in a folder relative to the simulation folder.  Warning, this will override all of your file paths to model package files and will also override any relative file paths to external model data files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change path of model files relative to the simulation folder\n",
    "model.set_model_relative_path('model_folder')\n",
    "\n",
    "# create folder\n",
    "if not os.path.isdir(save_folder):\n",
    "    os.makedirs(os.path.join(save_folder,'model_folder'))\n",
    "\n",
    "# write simulation to new folder\n",
    "sim.write_simulation()\n",
    "\n",
    "# run simulation from new folder\n",
    "sim.run_simulation()   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-Processing the Results\n",
    "\n",
    "Results can be retrieved from the master simulation dictionary.  Results are retrieved from the master simulation dictionary with using a tuple key that identifies the data to be retrieved.  For head data use the key\n",
    "\n",
    "('&lt;model name&gt;', 'HDS', 'HEAD')\n",
    "\n",
    "where &lt;model name&gt; is the name of your model. For cell by cell budget data use the key\n",
    "\n",
    "('&lt;model name&gt;', 'CBC', '&lt;flow data name&gt;')\n",
    "\n",
    "where &lt;flow data name&gt; is the name of the flow data to be retrieved (ex. 'FLOW-JA-FACE').  All available output keys can be retrieved using the output_keys method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keys = sim.simulation_data.mfdata.output_keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entries in the list above are keys for data in the head file \"HDS\" and data in cell by cell flow file \"CBC\". Keys in this list are not gaurenteed to be in any particular order. The code below uses the head file key to retrieve head data and then plots head data using matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# get all head data\n",
    "head = sim.simulation_data.mfdata['example_model', 'HDS', 'HEAD']\n",
    "# get the head data from the end of the model run\n",
    "head_end = head[-1]\n",
    "# plot the head data from the end of the model run\n",
    "levels = np.arange(160,162,1)\n",
    "extent = (0.0, 1000.0, 2500.0, 0.0)\n",
    "plt.contour(head_end[0, :, :],extent=extent)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results can also be retrieved using the existing binaryfile method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get head data using old flopy method\n",
    "hds_path = os.path.join(sim_path, model_name + '.hds')\n",
    "hds = flopy.utils.HeadFile(hds_path)\n",
    "# get heads after 1.0 days\n",
    "head = hds.get_data(totim=1.0)\n",
    "# plot head data\n",
    "plt.contour(head[0, :, :],extent=extent)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}