{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Freshwater Lens\n\nThis 2D examples illustrates the growth of a fresh water lens in an initially\nfully saline domain.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll start with the usual imports\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport xarray as xr\n\nimport imod"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Discretization\n\nWe'll start off by creating a model discretization, since\nthis is a simple conceptual model.\nThe model is a 2D cross-section, hence ``nrow = 1``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nrow = 1  # number of rows\nncol = 40  # number of columns\nnlay = 15  # number of layers\n\ndz = 10\ndx = 250\ndy = -dx"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set up tops and bottoms\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "top1D = xr.DataArray(\n    np.arange(nlay * dz, 0.0, -dz), {\"layer\": np.arange(1, nlay + 1)}, (\"layer\")\n)\n\nbot = top1D - dz"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set up ibound, which sets where active cells are `(ibound = 1.0)`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bnd = xr.DataArray(\n    data=np.full((nlay, nrow, ncol), 1.0),\n    coords={\n        \"y\": [0.5],\n        \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n        \"layer\": np.arange(1, 1 + nlay),\n        \"dx\": dx,\n        \"dy\": dy,\n    },\n    dims=(\"layer\", \"y\", \"x\"),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary Conditions\n\nSet the constant heads by specifying a negative value in iboud,\nthat is: ``bnd[index] = -1```\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bnd[0, :, 0:12] = -1\nbnd[0, :, 28:40] = -1\n\nfig, ax = plt.subplots()\nbnd.plot(y=\"layer\", yincrease=False, ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the recharge rates\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_rate = xr.DataArray(\n    data=np.full((nrow, ncol), 0.0),\n    coords={\"y\": [0.5], \"x\": np.arange(0.5 * dx, dx * ncol, dx), \"dx\": dx, \"dy\": dy},\n    dims=(\"y\", \"x\"),\n)\nrch_rate[:, 13:27] = 0.001\n\nfig, ax = plt.subplots()\nrch_rate.plot(ax=ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The model is recharged with fresh water\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rch_conc = xr.full_like(rch_rate, fill_value=0.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initial Conditions\n\nDefining the starting concentrations\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sconc = xr.DataArray(\n    data=np.full((nlay, nrow, ncol), 35.0),\n    coords={\n        \"y\": [0.5],\n        \"x\": np.arange(0.5 * dx, dx * ncol, dx),\n        \"layer\": np.arange(1, nlay + 1),\n        \"dx\": dx,\n        \"dy\": dy,\n    },\n    dims=(\"layer\", \"y\", \"x\"),\n)\n\nsconc[:, 13:27, 0] = 0.0\n\nfig, ax = plt.subplots()\nsconc.plot(y=\"layer\", yincrease=False, ax=ax)"
      ]
    },
    {
      "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 = imod.wq.SeawatModel(\"FreshwaterLens\")\nm[\"bas\"] = imod.wq.BasicFlow(ibound=bnd, top=150.0, bottom=bot, starting_head=0.0)\nm[\"lpf\"] = imod.wq.LayerPropertyFlow(\n    k_horizontal=10.0, k_vertical=20.0, specific_storage=0.0\n)\nm[\"btn\"] = imod.wq.BasicTransport(\n    icbund=bnd, starting_concentration=sconc, porosity=0.35\n)\nm[\"adv\"] = imod.wq.AdvectionTVD(courant=1.0)\nm[\"dsp\"] = imod.wq.Dispersion(longitudinal=0.0, diffusion_coefficient=0.0)\nm[\"vdf\"] = imod.wq.VariableDensityFlow(density_concentration_slope=0.71)\nm[\"rch\"] = imod.wq.RechargeHighestActive(rate=rch_rate, concentration=0.0)\nm[\"pcg\"] = imod.wq.PreconditionedConjugateGradientSolver(\n    max_iter=150, inner_iter=30, hclose=0.0001, rclose=0.1, relax=0.98, damp=1.0\n)\nm[\"gcg\"] = imod.wq.GeneralizedConjugateGradientSolver(\n    max_iter=150,\n    inner_iter=30,\n    cclose=1.0e-6,\n    preconditioner=\"mic\",\n    lump_dispersion=True,\n)\nm[\"oc\"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)\nm.create_time_discretization(additional_times=[\"1900-01-01T00:00\", \"2000-01-01T00:00\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we write the model, including runfile:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "modeldir = imod.util.temporary_directory()\nm.write(modeldir, resultdir_is_workdir=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run\n\nYou can run the model using the comand prompt and the iMOD-WQ executable.\nThis is part of the iMOD v5 release, which can be downloaded here:\nhttps://oss.deltares.nl/web/imod/download-imod5 .\nThis only works on Windows.\n\nTo run your model, open up a command prompt\nand run the following commands:\n\n.. code-block:: batch\n\n   cd c:\\path\\to\\modeldir\n   c:\\path\\to\\imod\\folder\\iMOD-WQ_V5_3_SVN359_X64R.exe FreshwaterLens.run\n\nNote that the version name of your executable might differ.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualise results\n\nAfter succesfully running the model, you can\nplot results as follows:\n\n.. code:: python\n\n   head = imod.idf.open(modeldir / \"results/head/*.idf\")\n\n   fig, ax = plt.subplots()\n   head.plot(yincrease=False, ax=ax)\n\n   conc = imod.idf.open(modeldir / \"results/conc/*.idf\")\n\n   fig, ax = plt.subplots()\n   conc.plot(levels=range(0, 35, 5), yincrease=False, ax=ax)\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.10.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}