Workshop_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e099f5b2-4b8c-4406-814d-2cd85d46dbe8",
   "metadata": {},
   "source": [
    "# Workshop 1\n",
    "Hi and welcome to the first of a series of workshops on using MODFLOW-6 with Flopy.\n",
    "To complete this workshop you will need the following:\n",
    "1. The USGS Modflow-6 suite - get it [here](https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6.4.2.zip)\n",
    "2. A Python environment with flopy, numpy, pandas, os, shutil, sys and matplotlib\n",
    "\n",
    "Ideally the Modflow6 execuatble (located in the bin folder of the downloaded suite) should be avaialble system wide otherwise it needs to be placed in the folder where this notebook will running.\n",
    "\n",
    "Open the mf6io.pdf document that comes with the MF6 suite in the documentation folder. This is the reference document that you should be using to check your model input files and it explains all the different options available to you. This point needs to be emphasized. As the modeller, ***it is your responsibility to check your input files***. Do not presume that because Flopy executed without error that your input files are fine. Flopy is a tool to help you create those files – it can still produce input files that will cause your model to fail.\n",
    "\n",
    "Only the base flopy installation is used - not the optional dependencies.\n",
    "\n",
    "This notebook was developed with flopy version 3.4.2\n",
    "\n",
    "We'll get straight into it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e8b1d94-3571-48a2-97f2-2d450384d32e",
   "metadata": {},
   "source": [
    "# Imports\n",
    "These are some of the most common libraries that I use when working on building a model with Flopy. You will notice that some libraries are imported with an alias. Hopefully you are all familiar with the reasons for this. If not don't worry it will become clearer later. I find it useful to always check which version of the library I am using in case I get any issues with errors that may be specific to that version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d7c2d5f-9ec4-4bf7-bd09-5881293bbca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import flopy\n",
    "\n",
    "\n",
    "print(f\"Pandas version = {pd.__version__}\")\n",
    "print(f\"Numpy version = {np.__version__}\")\n",
    "print(f\"Flopy version = {flopy.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae086923-9dce-46c2-8c0a-8bc951cede45",
   "metadata": {},
   "source": [
    "# Setting up some folders to dump things in later\n",
    "We will be creating files and perhaps using files from other utilities, so it is nice to have things organised. You can do this however you want but as a simple example you can follow what I do below. Let’s first look at where we are on the computer by using the os library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a37b4b2-d7d9-42ce-8f37-e1e1a64f55b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a24be431-1c9f-46d0-aaf8-605d193d5357",
   "metadata": {},
   "source": [
    "The path above is where we are working right now. Lets make a new folder specific to this workshop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbae53c1-d936-4e6b-9cdc-52570e829bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws1 = os.path.join('workshop_1') # here we are making a path not creating the folder\n",
    "if os.path.exists(ws1): # here we are asking if the path exists on the computer. \n",
    "    shutil.rmtree(ws1)# if it does exist, delete it and all the files in it\n",
    "    os.mkdir(ws1) # then remake it\n",
    "else:\n",
    "    os.mkdir(ws1) # if it doesn't exist then make the folder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a520089-3218-4ae7-835b-6e16122c5365",
   "metadata": {},
   "source": [
    "Let's make a few extra paths to folders that we will use later. Notice how the workshop_1 folder path (which is ws1) is now the first entry in the join chain. This means that the paths we are making will be sub-directories of the workshop_1 folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f8a654d-3be9-4dda-8bb5-6545b4ef5770",
   "metadata": {},
   "outputs": [],
   "source": [
    "gis_f = os.path.join(ws1,'GIS') # creating a sub-directory path for our gis input/output\n",
    "model_f = os.path.join(ws1,'model') # creating a sub-directory path for our model input/output\n",
    "plots_f = os.path.join(ws1,'plots') # creating a sub-directory path for our plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75caf1a0-8e34-4c4c-a7c0-4ebdda6d1ddd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at the paths we made\n",
    "print(gis_f)\n",
    "print(model_f)\n",
    "print(plots_f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1fee469-22c5-4a7c-803d-aac16bc22d51",
   "metadata": {},
   "source": [
    "Now we will make a few sub-directories inside our workshop folder by using a for loop from a list of strings called subdirectory_list. Once you execute the code below check you folder for the subdirectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4dc1031-859e-4c00-b6f3-82bef183db4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "subdirectory_pathlist = [gis_f,model_f,plots_f] # this is our list of sub-directory paths, each one is a path string\n",
    "for path in subdirectory_pathlist: # for each path in the list above\n",
    "    if os.path.exists(path): # check if the folder already exists\n",
    "        pass # if it does exist then do nothing\n",
    "    else:\n",
    "        os.mkdir(path) # otherwise make the folder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2c2458a-4295-415f-8e8d-73a63fa43466",
   "metadata": {},
   "source": [
    "# Setting up a MF6 simulation\n",
    "The simulation namefile is required for all MF6 modelling work. What you are creating next is the mfsim.nam file. This is the file at the highest level when it comes to all modelling work in you specific project using MF6. The mfsim.nam file is what the MF6 executable will look for as soon as you try to run a model in you project folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdfd06ea-3fc4-4a6d-a513-77952bc3b97b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", # this is a path to your exe for mf6, if it is available system wide then all you need is \"mf6\"\n",
    "                             verbosity_level=1, # this is optional but can be handy\n",
    "                             sim_ws=model_f) # this one is mandatory if you want to make sure that the model files are created in the folder we want"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f4319a5-b470-43c5-bff1-f54f9c54e01d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "344ef337-c26d-47f6-82bd-4a4d47c64d24",
   "metadata": {},
   "source": [
    "Let's see if we can write the mfsim.nam file using Flopy's package write function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1602822b-e141-4d1f-857d-b7e83471d37a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19951c2b-1f52-4990-b117-4ec7054e4808",
   "metadata": {},
   "source": [
    "Hmm, looks like it failed trying to find the tdis package object but it did succeed in writing the mfsim.nam file before failing. The write_simualtion() method expects a timing package object to be present or it will fail. We'll discuss this more later but for now let's take a quick peek at the mfsim.nam file. You can execute the code below to view the contents of the file or alternatively open the file in its folder with your favourite text editor. Do you know where to look for the file?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1d59e5c-7760-44fb-aa15-a559c22ee8b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"mfsim.nam\"))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fab33e9-c476-4956-a66d-34a069437f60",
   "metadata": {},
   "source": [
    "Open the mf6io.pdf document to see what else could be in this file if we did things a bit differently.\n",
    "We will discuss what the other sections are in the document and whether you add them explicitly to the sim object of if Flopy will handle that for you."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16df3d23-3b54-4f90-a9fa-97956f404980",
   "metadata": {},
   "source": [
    "Things to think about.\n",
    "\n",
    "Can I remake the object? Can a simulation have more than one simulation object? Are there any conflicts with options? What happens if I enter something incorrectly? Try to remake the sim object with all options turned on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "289ea2cf-f338-43db-bbe4-c5aa0c62e15a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", # this is a path to your exe for mf6, if it is available system wide then all you need is \"mf6\"\n",
    "                             verbosity_level=1, # this is optional but can be handy\n",
    "                             sim_ws=model_f, # this one is mandatory if you want to make sure that the model files are created in the folder we want\n",
    "                             continue_=True,\n",
    "                             nocheck=True,\n",
    "                             lazy_io=True,\n",
    "                             memory_print_option=\"ALL\",\n",
    "                             write_headers=False) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b2421c-f931-46ab-adf1-289ecf965d04",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdebd91d-269a-44a6-9235-1182ac74f83a",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"mfsim.nam\"))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cca8bf3-a60e-460d-9018-7aa5b9e44132",
   "metadata": {},
   "source": [
    "Wait a minute what about maxerrors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6cae287-9e86-4cab-8781-b13ef45381bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", # this is a path to your exe for mf6, if it is available system wide then all you need is \"mf6\"\n",
    "                             verbosity_level=1, # this is optional but can be handy\n",
    "                             sim_ws=model_f, # this one is mandatory if you want to make sure that the model files are created in the folder we want\n",
    "                             continue_=False,\n",
    "                             nocheck=False,\n",
    "                             lazy_io=False,\n",
    "                             memory_print_option=\"ALL\",\n",
    "                             write_headers=False,\n",
    "                             maxerrors=20) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a4d47b7-a727-4c36-9e8c-ce67faff5727",
   "metadata": {},
   "source": [
    "Out of options. Time to go to GitHub or dig into the source code of the library installed on your computer. What you will find is that the installed version does not have the maxerrors option available but the online develop branch version does. This is likely to change in the future. This does not mean you can't use the option because it already exists in MF6. It does mean that you can't use it through Flopy because it is in development. For the purposes of this workshop we don't need it but this brief diversion is handy to demonstrate that the software you are using is being developed on a continuous basis. This includes both MF6 and Flopy. New features added to MF6 may take a while to become available in Flopy. So you need to know how to edit the input files manually, which is where the mf6io.pdf becomes your primary resource. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d82b83-3f33-4bfa-b188-1e5671f158f7",
   "metadata": {},
   "source": [
    "What if I want to change something in the sim object after it has already been created. There are methods that allow you to access settings directly but for these workshops we're going to keep it as simple as possible so just recreate the object with new settings. but be sure to keep the object names the same. Otherwise, you may not be replacing the previous object but attempting to add another one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da4ebb79-6c19-43db-886d-c7235e0b55cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", # this is a path to your exe for mf6, if it is available system wide then all you need is \"mf6\"\n",
    "                             verbosity_level=1, # this is optional but can be handy\n",
    "                             sim_ws=model_f) # this one is mandatory if you want to make sure that the model files are created in the folder we want"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6461b491-2fd9-4e1e-9e09-9b013e855d28",
   "metadata": {},
   "source": [
    "Thats basically all you need to know for the sim object in order to create it. But what if you have an existing model and you want to load it instead. The Flopy GitHub introduction model is available as in import for use with the sim package that you just created. Execute the following block of code to import the function and then run it by passing in your sim object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe8ac413-3373-4851-84dc-ebbd5acd0ad4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers import ws1_mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b230e23-f715-4e88-ba74-0edddfb2e83c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ws1_mod(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb19ab70-a17d-4c1d-8971-9ecbb0b16cba",
   "metadata": {},
   "source": [
    "So you've just built a complete model. Take a quick look at the sim object now. Then open up the simulation workspace folder. Note that you haven't written your model files yet. Let's take a look at what the sim object looks like now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "122f3a36-0162-4547-a7cd-13865fcf3cd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a02de8b-20cb-40c6-9eba-7b116cbfd690",
   "metadata": {},
   "source": [
    "Now to write your simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e99e370-dcf4-41ab-bcba-55d0c0718ffa",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "863c7823-74b7-491d-ae81-eee29b7bfe3d",
   "metadata": {},
   "source": [
    "Have a look at your model folder now. The model input files should be there - but no output files just yet.\n",
    "Now lets run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27f426c3-261a-4b00-aad1-6ef2615d27e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ec4a2cd-a41f-42d0-808a-55aa68ce4ef0",
   "metadata": {},
   "source": [
    "Have a look at the folder again. If your model completed successfully. You should now see output files in addition to your input files. Great now we can look at how to load an existing model using the simulation object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a418e6b9-fbf4-44e8-b5ed-e5e4fcd75085",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = flopy.mf6.MFSimulation.load(sim_ws=model_f) \n",
    "# note that for MF6 you generally only have to provide the path to the model\n",
    "# It will automatically search for the mfsim.nam file in that location \n",
    "# then load the model from that."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2385f868-f858-478c-999b-7d92554585e5",
   "metadata": {},
   "source": [
    "If the model was not created by Flopy but some other GUI then you may have trouble loading all packages.\n",
    "Again this is where you need the mf6io.pdf to help you out because you may be able to load the package with a minor tweaks to the problematic input file.\n",
    "If it requires major modifications then perhaps consuder a load without that package using the load only options as demonstrated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6abb9bd-a293-4b01-84b8-9b1e78871fb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "load_list = [\"npf\"]\n",
    "sim = flopy.mf6.MFSimulation.load(sim_ws=model_f, load_only=load_list) \n",
    "#note that there are certain packages that it will always load irrespective of what you declare as load only."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2d655de-d486-42a6-99eb-33e1c7683e61",
   "metadata": {},
   "source": [
    "Once you load the model and whatever packages you want, then you can \"get\" the simulation level package objects from the simulation object. Note only simulation level package objects. Accessing loaded model level package objects requires that you first \"get\" the loaded model object from the loaded simulation. We will do more with the model object later but for completeness we will access both simulation and model level package objects here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18aca910-6576-477b-b1e2-8172c5d7aa6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis = sim.get_package('tdis') # a simulation level package object\n",
    "gwf = sim.get_model(\"MyModel\") # get the model object in order to access model level package\n",
    "npf = gwf.get_package(\"NPF\") # get the loaded npf package object and assign it to a variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea64d8fd-cca9-47e4-b54a-588e14724570",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that the variable type mtaches the loaded package\n",
    "type(npf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac27f76-f2e4-4524-8402-ec98b3be0c41",
   "metadata": {},
   "source": [
    "Now we are going to clear all that and start again but this time we will make our own time discretization package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a199a84-ec09-4b9a-8834-de2be61769f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in [sim,gwf,npf,tdis]: # loop through a list of objects we want to remove\n",
    "    del x # use the del statement to remove the objects, this clears them from memory\n",
    "flist = [x for x in os.listdir(model_f)] # gets a list of all the files in the model folder\n",
    "_ = [os.remove(os.path.join(model_f,f)) for f in flist] # removes all the files in the model folder, faster as list comprehension? Maybe.\n",
    "# now check that it worked\n",
    "flist = [x for x in os.listdir(model_f)] # get the files again\n",
    "flist # should see an empty list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "538354df-67ca-4412-9b3d-3cdb1433cc21",
   "metadata": {},
   "source": [
    "# The timing file or TDIS\n",
    "\n",
    "The important thing to note with the creation of the tdis object is that you are passing in the simulation object. Only the simulation level package objects require that you pass in a sim object.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aafe9b15-763a-4a6b-a071-63978ccdf6b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", # this is a path to your exe for mf6, if it is available system wide then all you need is \"mf6\"\n",
    "                             verbosity_level=1, # this is optional but can be handy\n",
    "                             sim_ws=model_f) # this one is mandatory if you want to make sure that the model files are created in the folder we want\n",
    "tdis = flopy.mf6.ModflowTdis(sim) # will give you a basic tdis object set for single steady-state stress period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fab7f218-5632-44eb-a50a-8b1658772583",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e15321-e71b-4217-9cda-26b3a77f1b57",
   "metadata": {},
   "source": [
    "That may be all you need but let’s try to make a slightly better one specifying time units and our model start datetime as ISO 8601 format. Note also that it will warn you that it detects an existing timing package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d71b2c7f-9481-44f9-9261-a4f573bef033",
   "metadata": {},
   "outputs": [],
   "source": [
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\")\n",
    "tdis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25950b5c-84f6-4e4b-8e1d-2470d35fa838",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380c19dc-5fbb-4a17-8490-477171d9b204",
   "metadata": {},
   "source": [
    "What if we have a transient model? Lest start with a transient single stress period model. Its going to run for a year and have approximately weekly time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77d9d518-9cd3-4a9f-85ba-d5171090007f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# start by creating a new tuple for the first stress period\n",
    "# the tuple needs to be in a sequence of stress period length, number of steps and mutiplier.\n",
    "# you can change it directly or rebuild the object\n",
    "tdis.perioddata.array[0] = (365.25,52,1)\n",
    "tdis\n",
    "sim.write_simulation()\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"MySim.tdis\"))]\n",
    "#tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d67ab45-e552-4700-a27a-5b8fb5c475a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# alternativley just rebuild it\n",
    "del tdis\n",
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\",perioddata=[(365.25,52,1)])\n",
    "tdis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78a12380-1eda-4336-8cd9-2eae80af463b",
   "metadata": {},
   "source": [
    "Okay now we will try something a bit tricky with time discretization. We are going to build our stress period data using actual dates. We want an initial steady-state period followed by 20 years of monthly stress periods with monthly time steps followed by a single stress period of 100 years with approximately monthly time steps. Our model will start on the 31/12/2023 with the steady-state period of 1 day. This means we want our transient period to start on the 01/01/2024."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ac32d0e-f5ee-4fc1-9906-fbfba1ab72b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we start with the 20 years of monthly stress periods, note the starting date used in the range.\n",
    "dates = pd.date_range('2024-01-01','2043-01-01', freq='MS').tolist() # use pandas to create a list of timestamps in monthly frquency\n",
    "dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2709f0fe-d63b-473a-b88f-8e5d79a72564",
   "metadata": {},
   "outputs": [],
   "source": [
    "perlens = [(dates[x]-dates[x-1]).days for x in range(1,len(dates))] # gives me the monthly period lengths in days. note use of .days\n",
    "perlens # this will help to check if we have the correct starting values.\n",
    "# Question: why do I start my range at 1 and not zero? Why not for x in dates\n",
    "# Recall that we need to provide a list of tuples to the tdis object where the first item in the tuple is the stress period length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a03d598-4539-414f-a089-484881c135db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the number of entries in period lengths and assgn it to a variable because we will need it later\n",
    "numper = len(perlens) # checks the length of the list. It should be 228 because we have 20 * 12 stress periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b368bafa-ea15-431a-a19c-a3c0b34c466b",
   "metadata": {},
   "outputs": [],
   "source": [
    "stp = 1 # we'll see how the model stability goes with 1 time step per period to start off with \n",
    "# build the transient stress period data tuples for the tdis [ackage\n",
    "_ = [(x,stp,1) for x in perlens] # note the final 1 in the tuple is the multiplier. We won't use it to begin with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "648b0c31-5968-49da-944a-74c80d2cb38c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in the steady state stress period at the start as (1,1,1) and 100 year period at the end\n",
    "# note my 100 year period has 1200 steps so sticking with approximatley monthly time steps\n",
    "pdata = [(1,1,1), *_, (36525, 1200, 1)] # note the use of unpacking the _ tuple list\n",
    "# now we have a single recovery period that is 100 years in length with approximatley mothly time steps\n",
    "numper = len(pdata) # we are updating the value of numper here to include our steady-state and 100 year periods\n",
    "print(numper) # print the value for viewing\n",
    "pdata # view all the period data. Note that this is what we will provide to the "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a522a2bb-43ca-4107-b78a-5eab58e3b0c0",
   "metadata": {},
   "source": [
    "As a reference we are going to create a dataframe for use later. This will provide a check for any transient behaviour for boundary conditions. For example, if we know pumping is scheduled to commence in February 2025 then we might need to know which stress period this is. I say might here because there are alternative ways to control transient boundary conditions in MF6 that were not available in previous versions of MODFLOW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f3608e2-d60a-4882-9fe0-bcf4fb69416d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ss = pd.Timestamp(\"2023-12-31\") # set the date for the SS period\n",
    "dates = [ss,*dates] # add the date\n",
    "# Hmm, could we have done this originally and included the steady-state stress period when we got the tuple lengths? Try it yourself.\n",
    "len(dates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7541f3c-f83f-484b-9861-4a42324feb18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now lets build the dataframe\n",
    "df = pd.DataFrame() # start with en empty dataframe\n",
    "df['Date'] = dates # create a column with dates in it\n",
    "df['SP'] = range(1,len(dates)+1) # create a column with the actual MF6 stress period number in it\n",
    "df['Flopy_SP'] = range(len(dates)) # create a column with the Flopy zero-based stress period numbers in it\n",
    "df.to_csv('model_timing.csv',index=None) # export the csv to the project folder\n",
    "# if you want to keep this dataframe available for use later in the script then make a copy of it using a descriptive name\n",
    "modtime_df = df.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89067074-669c-46cb-b981-1e0b8c3534c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we create the new tdis object\n",
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\",\n",
    "                             perioddata=pdata,nper=numper)\n",
    "tdis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4184c58e-22f5-433a-8e35-c10ef5e03e70",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets write our simulation\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "294a12b9-25e7-4efe-96bb-d8402c8d61ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets look at the tdis file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"MySim.tdis\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e331d35-96ca-49a2-b772-032287c56903",
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers import ws1_mod_trans \n",
    "# now we are importing the same model as before but with one extra package being built. \n",
    "# What is the new package? Why is it needed now?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "effbc288-8a31-4f97-82f9-b2490a3dac56",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now lets run the transient version of the same model \n",
    "ws1_mod_trans(sim,tdis) # remember this is only putting the objects in memory you havent created the files yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bb5b613-3f83-4bad-ac86-ade68c6a98f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation() # write the simulation files\n",
    "sim.run_simulation() # run the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf94985-aa72-43db-ae66-fcd64db2e814",
   "metadata": {},
   "source": [
    "Ideally you should be building models that are very stable but in some instances you may have to use an unstable model or you have advanced stress packages that are causing you some strife so you may opt to use adaptive time stepping. This will automatically change the time step length in an attempt to solve the model for you. Note this is one strategy to get to a solution for your model there are other options too. You may also choose to set this up immediately from the get go but there are some caveats to doing this with regards to model output that you should be aware of. Take some time to read the ATS section in the mf6io.pdf to make sure that you understand how it all works. I'll provide some brief descriptions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a743a3c2-20cb-4a6a-89aa-df84d5958eb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need to provide the tdis object with some data for the ATS package in order for it to be created for us.\n",
    "# in this example we are only going to want ATS active during the monthly stress periods\n",
    "# first we define which stress perios we want the ATS active for\n",
    "iperats = [x for x in range(1,numper-1)] # what am i doing here exactly?\n",
    "dt0 = 0 # Now we say what we want the initial time step for the stress period to be. An entry of zero is special. Read up why.\n",
    "dtmin = 0.01 # next we need to define the minimum time step length we will allow. Recall our time units are days.\n",
    "# This next one requires the most thought and depending on how much your model is struggling to converge may require adjusting.\n",
    "# We will use a recommended value for now\n",
    "dtadj = 3.0\n",
    "# The last variable is also worth reading up about, it controls how the time step length is reduced when your model fails to converge\n",
    "# for brevity we will use the recommended value\n",
    "dtfailadj = 5.0\n",
    "# now we create our input data list of tuples. Check the mf6io.pdf for the entries to make sure we have them ordered correctly.\n",
    "ats_pdata = [[x,dt0,dtmin,y,dtadj,dtfailadj] for x,y in zip(iperats,perlens)] # What are x and y here?\n",
    "# now we need to pass this on to the tdis package and re-write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdd745a8-69d5-4f28-a44d-f1d91dcca962",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we initialize the ats in the existing tdis object\n",
    "tdis.ats.initialize(maxats=len(ats_pdata),perioddata=ats_pdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e2b9dd9-6d7d-4f88-8cb7-0f072f87775d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.write_simulation() # look for the ats "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f832f9d8-d129-48b8-a472-834db733ed86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets look at our tdis input file now\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"MySim.tdis\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6912803-0b00-47b7-af25-8b6a9e8188fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets look at our ATS input file now\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"MySim.tdis.ats\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea6d3ee6-1538-4051-a02e-2fd5d7aef163",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets run it\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbdcdfa4-8ca5-4d07-958d-90000b3f45bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets clear everything again\n",
    "for x in [sim,tdis]: # loop through a list of objects we created that we want to remove\n",
    "    del x # use the del statement to remove the objects, this clears them from memory\n",
    "flist = [x for x in os.listdir(model_f)] # gets a list of all the files in the model folder\n",
    "_ = [os.remove(os.path.join(model_f,f)) for f in flist] # removes all the files in the model folder, faster as list comprehension? Maybe.\n",
    "# now check that it worked\n",
    "flist = [x for x in os.listdir(model_f)] # get the files again\n",
    "flist # should see an empty list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43c79f3-505c-4c7c-9118-f9c98e22026a",
   "metadata": {},
   "source": [
    "# The solver or IMS package\n",
    "\n",
    "This package is mandatory. But the input file is reasonably simple to construct with Flopy. There are lots of options for you to use if you really want to customise the solver but the defaults are generally fine to ensure convergence. The most common settings that you will change are your convergence criteria and you number of iterations. There may be occasions where you want to optimise your solver for speed then the other settings could be of benefit but if you stick with the \"SIMPLE\", \"MODERATE\" and \"COMPLEX\" you should be fine. If not, consider investigating model design issues first before resorting to exotic solver configurations just to converge. The mf6io.pdf document has some information on the benefit of each solver setting plus recommendations for values. These are worth a read at least once. You have to have a sim object in memory before you create this object because it is needed as part of the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d69e32eb-1e43-4107-8614-5fee4e03af20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets create our sim and tdis objects again we'll drop ATS for now. But stick with the transient model configuration.\n",
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", \n",
    "                             verbosity_level=1,\n",
    "                             sim_ws=model_f)\n",
    "dates = pd.date_range('2024-01-01','2043-01-01', freq='MS').tolist()\n",
    "ss_date = pd.Timestamp(\"2023-12-31\")\n",
    "dates = [ss_date,*dates] # using all the dates\n",
    "perlens = [(dates[x]-dates[x-1]).days for x in range(1,len(dates))]\n",
    "stp = 1\n",
    "_ = [(x,stp,1) for x in perlens] # now they all have a step of 1 but we want 1200 steps for the last one\n",
    "pdata = [*_,(36525,1200,1)] # adding the final stress period as a tuple\n",
    "numper = len(pdata)\n",
    "# Could we have added the final date to our dates and then just changed the number of time steps?\n",
    "# Look up mutable and imutable objects in Pyhton to find out why.\n",
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\",\n",
    "                             perioddata=pdata,nper=numper)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43500709-9285-4413-9497-fdc3ec60910f",
   "metadata": {},
   "source": [
    "The ims package object below is an example. As you can see the settings are very straight forward. Most of the input variables are either a sinlge string, single integer or float variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb86fdc8-c12d-4eba-b37e-a7d38421a29f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ims = flopy.mf6.ModflowIms(sim, complexity='MODERATE', # using the moderate solver settings option Note passing in sim object\n",
    "                           csv_inner_output_filerecord='inner.csv', # asking to output a csv of the inner iterations\n",
    "                           csv_outer_output_filerecord='outer.csv', # asking for an output of the outer itertions\n",
    "                           outer_maximum=500, # Overiding the default MODERATE value (50) setting the maximum number of outer iterations\n",
    "                           inner_maximum=500, #  Overiding the default MODERATE value (100) setting the maximum number of inner iterations\n",
    "                           outer_dvclose=0.01, # This is the default MODERATE value setting the convergence criteria for the outer iteration\n",
    "                           inner_dvclose=0.001) #  Overiding the default MODERATE value (0.01) setting the convergence criteria for the inner iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88888a72-c870-4913-af85-abe69288adad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's write the simulation\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e437093-82db-49ff-8248-fa9766a32ffe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's look at the file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"MySim.ims\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73b9e1bb-c507-44ba-b687-dc813e82a66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import a new model that needs the sim, tdis and ims objects\n",
    "from helpers import ws1_mod_trans2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbfbbc75-a91a-4a70-a2d1-43454cfce98a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now run the imported function and pass\n",
    "ws1_mod_trans2(sim,tdis,ims) # recall we are only building package objects here not writing files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3899dda-f978-405e-96bc-a98a796f4897",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now write the files then check you model folder\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4294b312-70d5-48c4-8769-681f60640da7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the simulation\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c67753f9-1724-4175-b437-e97e9bf52d85",
   "metadata": {},
   "source": [
    "Now we should have some information on the model convergence behaviour in two csv files named \"inner.csv\" and \"outer.csv\". Take a look at them and see if you can decipher how your chosen convergence criteria influence the solver behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8db1eb55-d3c0-4cc0-a521-a102a9834b1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets clear everything again\n",
    "for x in [sim,tdis,ims]: # loop through a list of objects we created that we want to remove\n",
    "    del x # use the del statement to remove the objects, this clears them from memory\n",
    "flist = [x for x in os.listdir(model_f)] # gets a list of all the files in the model folder\n",
    "_ = [os.remove(os.path.join(model_f,f)) for f in flist] # removes all the files in the model folder, faster as list comprehension? Maybe.\n",
    "# now check that it worked\n",
    "flist = [x for x in os.listdir(model_f)] # get the files again\n",
    "flist # should see an empty list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "553fcc60-851a-4953-8ba5-8b9dc1f9ed94",
   "metadata": {},
   "source": [
    "# The model name-file\n",
    "This is the last simulation level package object that you need to create and is also generally very straight forward to build. It is at this point that you need to provide a model-name. Recall that we set a simulation name previously and a simulation folder to write our model files in. With the model object you can change where the model level files are written using a relative path (optional). The example below is generally all you need. The use of newtonoptions is probably the only thing that requires some consideration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1afe4cb5-548b-430a-b87e-460c76a0277b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we build our simulation from scratch again (I'm channeling Zed Shaw here)\n",
    "# Lets create our sim and tdis objects again we'll drop ATS for now. But stick with the transient model configuration.\n",
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", \n",
    "                             verbosity_level=1,\n",
    "                             sim_ws=model_f)\n",
    "dates = pd.date_range('2024-01-01','2043-01-01', freq='MS').tolist()\n",
    "ss_date = pd.Timestamp(\"2023-12-31\")\n",
    "dates = [ss_date,*dates] # using all the dates\n",
    "perlens = [(dates[x]-dates[x-1]).days for x in range(1,len(dates))]\n",
    "stp = 1\n",
    "_ = [(x,stp,1) for x in perlens] # now they all have a step of 1 but we want 1200 steps for the last one\n",
    "pdata = [*_,(36525,1200,1)] # adding the final stress period as a tuple\n",
    "numper = len(pdata)\n",
    "# Could we have added the final date to our dates and then just changed the number of time steps?\n",
    "# Look up mutable and imutable objects in Pyhton to find out why.\n",
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\",\n",
    "                             perioddata=pdata,nper=numper)\n",
    "ims = flopy.mf6.ModflowIms(sim, complexity='MODERATE', # using the moderate solver settings option Note passing in sim object\n",
    "                           csv_inner_output_filerecord='inner.csv', # asking to output a csv of the inner iterations\n",
    "                           csv_outer_output_filerecord='outer.csv', # asking for an output of the outer itertions\n",
    "                           outer_maximum=500, # Overiding the default MODERATE value (50) setting the maximum number of outer iterations\n",
    "                           inner_maximum=500, #  Overiding the default MODERATE value (100) setting the maximum number of inner iterations\n",
    "                           outer_dvclose=0.01, # This is the default MODERATE value setting the convergence criteria for the outer iteration\n",
    "                           inner_dvclose=0.001) #  Overiding the default MODERATE value (0.01) setting the convergence criteria for the inner iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "812dae62-ebfe-430b-8824-c2850c198169",
   "metadata": {},
   "source": [
    "We will discuss some of the other options here. They tend to be provided for flexibility and are more use case specific. The example should still provide sufficient guidance on how to use the different options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb0e70f0-d212-41a0-9861-857849387b8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now first set the model name\n",
    "model_name = 'flow' \n",
    "# build the object without a specified relative path for the model files\n",
    "gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, save_flows=True, newtonoptions=\"under_relaxation\")\n",
    "# The save flows option just indicates that we want to save a budget file\n",
    "# What is the budget file? Do you need it? \n",
    "# This is worth thinking about because it can make a big different to your model speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c423c60-b8f7-410c-9cd5-b19cbfef0130",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now look at the sim object and note the inclusion of the model in the sim at the bottom of the output.\n",
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bfc2a4b-5e5d-4e5e-a8b0-acff585621be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now write the simulation.\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "949a3174-e793-4056-958f-f26b692ce6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's look at the model name-file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"flow.nam\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb451e96-860d-45d7-bbf8-d8a83dfe64ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the simulation name-file\n",
    "_ = [print(line.rstrip()) for line in open(os.path.join(model_f,\"mfsim.nam\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "422afabe-86c5-49e9-99a5-aeb05646c11a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import a new model that needs the sim, tdis, ims and gwf objects plus the model name\n",
    "from helpers import ws1_mod_trans3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0f31620-e4e9-4ecb-b285-75bfd534389c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now run the imported function and pass in the required objects and variables\n",
    "ws1_mod_trans3(sim,tdis,ims,gwf,model_name) # recall we are only building package objects here not writing files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0549d11b-a97d-4ad1-a652-544133e5c743",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the model files\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c8cd14-8a86-476b-83e3-03425145707f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model and check the model directory for the output files\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a5bc1c-1f4c-4050-958e-348582a00ed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's clear all the files from the model folder\n",
    "shutil.rmtree(model_f) # remove the model folder (this will also remove subdirectories in the folder)\n",
    "os.mkdir(model_f) # remake a new model folder\n",
    "# This is a different way to get rid of the files. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eef827d-9633-46db-8e18-9fcc33e9fe65",
   "metadata": {},
   "source": [
    "# Now repeat with a new option\n",
    "\n",
    "Here we use a relative path for the model which will place the model specific input files in a subdirectory. Keep an eye on your files and folders when you write the simualtion to get a feel for what this does. Also have a look at the mfsim.nam file to see how the path to the modle name fiel is altered. There is also a slight change to how we go about creating the tdis object's period data. Can you see where?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "093e4146-c5ab-4de8-a195-421ecc2249c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we build our simulation from scratch again (I'm channeling Zed Shaw here)\n",
    "# Lets create our sim and tdis objects again we'll drop ATS for now. But stick with the transient model configuration.\n",
    "sim_name = \"MySim\" # give your simulation a name, usually I use a project name but this is entirley up to you\n",
    "sim = flopy.mf6.MFSimulation(sim_name=sim_name, \n",
    "                             exe_name=\"mf6\", \n",
    "                             verbosity_level=1,\n",
    "                             sim_ws=model_f)\n",
    "dates = pd.date_range('2024-01-01','2043-01-01', freq='MS').tolist()\n",
    "ss_date = pd.Timestamp(\"2023-12-31\")\n",
    "dates = [ss_date,*dates] # using all the dates\n",
    "perlens = [(dates[x]-dates[x-1]).days for x in range(1,len(dates))]\n",
    "stp = 1\n",
    "_ = [(x,stp,1) for x in perlens] # now they all have a step of 1 but we want 1200 steps for the last one\n",
    "pdata = [*_,(36525,1200,1)] # adding the final stress period as a tuple\n",
    "numper = len(pdata)\n",
    "# Could we have added the final date to our dates and then just changed the number of time steps?\n",
    "# Look up mutable and imutable objects in Pyhton to find out why.\n",
    "tdis = flopy.mf6.ModflowTdis(sim,time_units='days',start_date_time=\"2023-12-31T00:00:00\",\n",
    "                             perioddata=pdata,nper=numper)\n",
    "ims = flopy.mf6.ModflowIms(sim, complexity='MODERATE', # using the moderate solver settings option Note passing in sim object\n",
    "                           csv_inner_output_filerecord='inner.csv', # asking to output a csv of the inner iterations\n",
    "                           csv_outer_output_filerecord='outer.csv', # asking for an output of the outer itertions\n",
    "                           outer_maximum=500, # Overiding the default MODERATE value (50) setting the maximum number of outer iterations\n",
    "                           inner_maximum=500, #  Overiding the default MODERATE value (100) setting the maximum number of inner iterations\n",
    "                           outer_dvclose=0.01, # This is the default MODERATE value setting the convergence criteria for the outer iteration\n",
    "                           inner_dvclose=0.001) #  Overiding the default MODERATE value (0.01) setting the convergence criteria for the inner iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0b15e65-f609-42af-b5e3-87f4c3be0a29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now first set the model name\n",
    "model_name = 'flow' \n",
    "# build the object again this time with a specified relative path for the model files\n",
    "gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, \n",
    "                           save_flows=True, \n",
    "                           newtonoptions=\"under_relaxation\",\n",
    "                           model_rel_path='.\\\\flow')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56729d92-f3ca-402c-a3bb-7e1c0081e648",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets take a look at the sim object again\n",
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0609b90e-5d64-47fc-bd67-216be8d55ef3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import a new model that needs the sim, tdis, ims and gwf objects plus the model name\n",
    "from helpers import ws1_mod_trans3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed4f8577-0b95-4816-bdce-56f26d1a2333",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now run the imported function and pass in the required objects and variables\n",
    "ws1_mod_trans3(sim,tdis,ims,gwf,model_name) # recall we are only building package objects here not writing files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad52d9e5-7211-460d-8915-dbeab949320b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the model files\n",
    "sim.write_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38e7ce79-3183-4fb1-9052-303db46b52aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model and check the model directory for the output files\n",
    "sim.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84d0c27-dcb9-4f3f-ae22-4a6ec21163a9",
   "metadata": {},
   "source": [
    "# Thats all folks\n",
    "\n",
    "Today we covered the main simualtion level packages. Next session we will focus purley on spatial discretisation. Hope it was worthwhile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "794c9bc2-f487-4245-a181-be2c3919af31",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}