{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Model Creation\nIn this example, we'll create an iMODFLOW model\nfrom scratch with complex boundary conditions\nand horizontal barriers.\n\nThere are two surface water systems:\nthe outer two edges of the grid feature ditches\nwith a rising stage, whereas the central ditch\nhas a periodic boundary conditions with a\nsummer and winter stage.\n\nThe model will be written as a projectfile\nwith a set of IDFs containing all the grid information\nand a .tim file containing the time discretization.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll start with the usual imports, supplied\nwith geopandas and shapely to specify vector data\nfor the `hfb` package.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import geopandas as gpd\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\nfrom shapely.geometry import LineString\n\nimport imod\nimport imod.flow as flow"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Discretization\n\nWe'll start off by creating a model discretization.\nThe model consists of a 3 by 9 by 9 three-dimensional grid.\n\nWe'll specify the grid first.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "shape = nlay, nrow, ncol = 3, 9, 9\n\ndx = 100.0\ndy = -100.0\ndz = np.array([5, 30, 100])\nxmin = 0.0\nxmax = dx * ncol\nymin = 0.0\nymax = abs(dy) * nrow\ndims = (\"layer\", \"y\", \"x\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we'll create the coordinates\nwhich set the grid dimensions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "layer = np.arange(1, nlay + 1)\ny = np.arange(ymax, ymin, dy) + 0.5 * dy\nx = np.arange(xmin, xmax, dx) + 0.5 * dx\ncoords = {\"layer\": layer, \"y\": y, \"x\": x}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The vertical grid discretization (tops and bottoms) are set with a 1D DataArray.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "surface = 0.0\ninterfaces = np.insert((surface - np.cumsum(dz)), 0, surface)\n\nbottom = xr.DataArray(interfaces[1:], coords={\"layer\": layer}, dims=\"layer\")\ntop = xr.DataArray(interfaces[:-1], coords={\"layer\": layer}, dims=\"layer\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll have to create a time discretization as well.\nCreate 1 year of monthly timesteps\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = pd.date_range(start=\"1/1/2018\", end=\"12/1/2018\", freq=\"MS\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll create our first 3 dimensional grid here,\nthe `ibound` grid, which sets where active cells are `(ibound = 1.0)`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ibound = xr.DataArray(np.ones(shape), coords=coords, dims=dims)\nibound"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Hydrogeology\n\nWe'll create a very simple hydrogeology,\nby specifying kh, kva, and sto as constants\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kh = 10.0\nkva = 1.0\nsto = 0.001"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initial conditions\n\nWe do not put much effort in the creation of the initial conditions\nin this example, instead we copy the ibounds.\nThis is a 3D grid filled with the value 1, and we can use it as a\ninital condition as well.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "starting_head = ibound.copy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions\n\nWe will put some more effort in creating some complex\nboundary conditions. We'll create both two outer ditches\nwith a rising stage, as well as a central ditch with periodic\n(summer-winter) stage.\n\n### Rising outer ditches\n\nFirst, we'll create rising trend,\nby creating a 1D array ones with the same size as the time dimension,\nand computing the cumulative sum over it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "trend = np.ones(times[:-1].shape)\ntrend = np.cumsum(trend)\ntrend_da = xr.DataArray(trend, coords={\"time\": times[:-1]}, dims=[\"time\"])\n\ntrend_da"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we assign values only to edges of model x domain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "is_x_edge = starting_head.x.isin([x[0], x[-1]])\nhead_edge = starting_head.where(is_x_edge)\nhead_edge"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's multiply our\n1D DataArray with dimension ``time``, with the static 3D grid\nwith dimension ``layer, y, x``,\nwhich xarray automatically broadcasts to a 4D array,\nwith dimensions ``time, layer, y, x``\nThis finishes our\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "head_edge_rising = trend_da * head_edge\nhead_edge_rising"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Periodic central ditch\n\nWe'll take only the central column of the grid with\n(``where``), the rest will be set to ``np.nan``,\nand from this we'll select only the upper layer,\nas the ditch will be located only in the upper layer.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "is_x_central = starting_head.x == x[4]\nhead_central = starting_head.where(is_x_central).sel(layer=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create period times, we let these times start before\nthe model starts.\nThis is necessary because `iMODFLOW` only forward fills periods.\nOtherwise, in this case there wouldn't be a\nperiodic boundary condition until april.\n\nWe will do this by selecting the months april and october,\nand then subtracting a year\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "period_times = times[[3, 9]] - np.timedelta64(365, \"D\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We are creating a summer and winter level.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "periods_da = xr.DataArray([4, 10], coords={\"time\": period_times}, dims=[\"time\"])\nhead_periodic = periods_da * head_central\n\nhead_periodic"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create dictionary to tell `iMOD`\nwhich period name corresponds to which date.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "timemap = {\n    period_times[0]: \"summer\",\n    period_times[1]: \"winter\",\n}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Wells\n\nWells are specified as a pandas dataframe.\nWe create a diagonal line of wells through the domain.\n\nBecause we can.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wel_df = pd.DataFrame()\nwel_df[\"id_name\"] = np.arange(len(x)).astype(str)\nwel_df[\"x\"] = x\nwel_df[\"y\"] = y\nwel_df[\"rate\"] = dx * dy * -1 * 0.5\nwel_df[\"time\"] = np.tile(times[:-1], 2)[: len(x)]\nwel_df[\"layer\"] = 2\n\nwel_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Horizontal Flow Barrier\n\nCreate barriers between ditches in layer 1 and 2 (but not 3).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "line1 = LineString([(x[2], ymin), (x[2], ymax)])\nline2 = LineString([(x[7], ymin), (x[7], ymax)])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll have to repeat each line for each layer\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lines = np.array([line1, line2, line1, line2], dtype=\"object\")\nhfb_layers = np.array([1, 1, 2, 2])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can specify names for our own bookkeeping\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "id_name = [\"left_upper\", \"right_upper\", \"left_lower\", \"right_lower\"]\n\n# The hfb has to specified as a geopandas `GeoDataFrame\n# <https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.html>`_\n\nhfb_gdf = gpd.GeoDataFrame(\n    geometry=lines, data=dict(id_name=id_name, layer=hfb_layers, resistance=100.0)\n)\n\nhfb_gdf"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build\n\nFinally, we build the model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = flow.ImodflowModel(\"my_first_imodflow_model\")\nm[\"pcg\"] = flow.PreconditionedConjugateGradientSolver()\n\nm[\"bnd\"] = flow.Boundary(ibound)\nm[\"top\"] = flow.Top(top)\nm[\"bottom\"] = flow.Bottom(bottom)\n\nm[\"khv\"] = flow.HorizontalHydraulicConductivity(kh)\nm[\"kva\"] = flow.VerticalAnisotropy(kva)\nm[\"sto\"] = flow.StorageCoefficient(sto)\n\nm[\"shd\"] = flow.StartingHead(starting_head)\n\nm[\"chd\"] = flow.ConstantHead(head=head_edge_rising)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create periodic boundary condition\nand specify it as a periodic stress package.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m[\"ghb\"] = flow.GeneralHeadBoundary(head=head_periodic, conductance=10.0)\nm[\"ghb\"].periodic_stress(timemap)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can specify a second stress package as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m[\"ghb2\"] = flow.GeneralHeadBoundary(head=head_periodic + 10.0, conductance=1.0)\n# You also need to specify periodic stresses for second system.\nm[\"ghb2\"].periodic_stress(timemap)\n\nm[\"wel\"] = flow.Well(**wel_df)\n\nm[\"hfb\"] = flow.HorizontalFlowBarrier(**hfb_gdf)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Specify output control; -1 specifies to save the output of interest for all layers.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m[\"oc\"] = flow.OutputControl(save_head=-1, save_flux=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "imod-python wants you to provide the model endtime to your time_discretization!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.create_time_discretization(additional_times=times[-1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we write the model\nWrites both .IDFs as well as projectfile, an inifile,\nand a .tim file that contains the time discretization.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "modeldir = imod.util.temporary_directory()\nm.write(directory=modeldir)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run\n\nYou can run the model using the comand prompt and the iMOD executables.\nThis is part of the iMOD v5 release, which can be downloaded here:\nhttps://oss.deltares.nl/web/imod/download-imod5 .\niMOD only works on Windows.\n\nTo run your model, open up a command prompt\nand run the following commands:\n\n```batch\ncd c:\\path\\to\\modeldir\nc:\\path\\to\\imod\\folder\\iMOD_v5_3_X64R.EXE config_run.ini\nc:\\path\\to\\imod\\folder\\iMODFLOW_V5_3_METASWAP_SVN1977_X64R.EXE my_first_imodflow_model.nam\n```\n<div class=\"alert alert-info\"><h4>Note</h4><p>iMODFLOW requires you to change directory into the model directory\n   before calling its executable. Otherwise the program will throw an error.</p></div>\n\n\n"
      ]
    }
  ],
  "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.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}