{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4abcbdae",
   "metadata": {},
   "source": [
    "# Introduction to pyISSM\n",
    "\n",
    "This notebook provides a first introduction to **pyISSM**, the Python interface to the Ice-sheet and Sea-level System Model (ISSM).\n",
    "\n",
    "This notebook focuses on helping you understand how pyISSM is structured, how its syntax works, and how to interact with ISSM models in a Python-native way.\n",
    "\n",
    "By the end of this notebook, you should be comfortable navigating pyISSM, inspecting model objects, and understanding where different pieces of functionality live.\n",
    "\n",
    "## Scope of this notebook\n",
    "\n",
    "This notebook focuses on **structure and syntax**, not on glaciology or ice sheet model theory, or model development. If you are new to pyISSM, we strongly recommend working through this notebook before moving on."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc396bf9",
   "metadata": {},
   "source": [
    "```{warning}\n",
    "**MATLAB -> Python Differences** \n",
    "\n",
    "For those users who are familiar with the existing MATLAB API, pyISSM will feel familiar, but not identical. pyISSM preserves ISSMs core concepts (e.g. the model object and general workflows), while adopting Pythonic syntax and conventions.\n",
    "\n",
    "Important differences include:\n",
    "- Adoption of keyword arguments, rather than positional argument pairs.\n",
    "- Function names use snake_case (underscores)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0334a24a",
   "metadata": {},
   "source": [
    "## Import pyISSM\n",
    "\n",
    "We'll start by importing `pyISSM` and getting a feel for it's structure. Instructions for installing `pyISSM` are provided [here](https://github.com/ACCESS-NRI/pyISSM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "04ae730a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyISSM contains the following submodules: ['analysis', 'data', 'model', 'plot', 'tools']\n"
     ]
    }
   ],
   "source": [
    "# Import pyISSM\n",
    "import pyissm\n",
    "\n",
    "# Check available submodules\n",
    "import pkgutil\n",
    "print(f'pyISSM contains the following submodules: {[m.name for m in pkgutil.iter_modules(pyissm.__path__)]}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cf87739",
   "metadata": {},
   "source": [
    "We can see that pyISSM contains five sub-modules.\n",
    "- `analysis` contains tailored post-processing tools\n",
    "- `data` contains tailored tools for interacting with external datasets, including interpolation routines. \n",
    "- `model` contains the central functionality, including the core `Model()` class.\n",
    "- `plot` conatains tailored tools for visualising ISSM models.\n",
    "- `tools` contains diverse tools for interacting with the ISSM Model Object.\n",
    "\n",
    "This notebook focuses on the structure of the `model` submodule. Other modules are primarily used for parameterising, visualising, and analysing ISSM models and will be explored in other notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225b550c",
   "metadata": {},
   "source": [
    "## The ISSM Model Object (`md`)\n",
    "\n",
    "`pyISSM` is developed around a central `Model()` class. Each model contains a series of sub-classes that define model components and parameters (e.g. `mesh`, `geometry`, etc.). Each model and sub-class contains detailed information that can be interrogated/explored interactively.\n",
    "\n",
    "Across the ISSM ecosystem (within the core C++ code, in MATLAB, and here in pyISSM), ISSM models are commonly named `md`.\n",
    "\n",
    "Let's create an empty ISSM model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dc9dbbe3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   ISSM Model Class                         \n",
       "                                            \n",
       "               mesh:  mesh properties         \n",
       "               mask:  defines grounded and floating elements \n",
       "           geometry:  surface elevation, bedrock topography, ice thickness, ... \n",
       "          constants:  physical constants      \n",
       "                smb:  surface mass balance    \n",
       "      basalforcings:  bed forcings            \n",
       "          materials:  material properties     \n",
       "             damage:  damage propagation laws \n",
       "           friction:  basal friction / drag properties \n",
       "       flowequation:  flow equations          \n",
       "       timestepping:  timestepping for transient models \n",
       "     initialization:  initial guess / state   \n",
       "              rifts:  rifts properties        \n",
       "         solidearth:  solidearth inputs and settings \n",
       "                dsl:  dynamic sea level       \n",
       "              debug:  debugging tools (valgrind, gprof \n",
       "            verbose:  verbosity level in solve \n",
       "           settings:  settings properties     \n",
       "           toolkits:  PETSc options for each solution \n",
       "            cluster:  cluster parameters (number of CPUs...) \n",
       "   balancethickness:  parameters for balancethickness solution \n",
       "      stressbalance:  parameters for stressbalance solution \n",
       "      groundingline:  parameters for groundingline solution \n",
       "          hydrology:  parameters for hydrology solution \n",
       "      masstransport:  parameters for masstransport solution \n",
       "            thermal:  parameters for thermal solution \n",
       "        steadystate:  parameters for steadystate solution \n",
       "          transient:  parameters for transient solution \n",
       "           levelset:  parameters for moving boundaries (level-set method) \n",
       "            calving:  parameters for calving  \n",
       "    frontalforcings:  parameters for frontalforcings \n",
       "                esa:  parameters for elastic adjustment solution \n",
       "           sampling:  parameters for stochastic sampler \n",
       "               love:  parameters for love solution \n",
       "           autodiff:  automatic differentiation parameters \n",
       "          inversion:  parameters for inverse methods \n",
       "                qmu:  Dakota properties       \n",
       "                amr:  adaptive mesh refinement properties \n",
       "   outputdefinition:  output definition       \n",
       "            results:  model results           \n",
       "       radaroverlay:  radar image for plot overlay \n",
       "      miscellaneous:  miscellaneous fields    \n",
       "  stochasticforcing:  stochasticity applied to model forcings "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create an empty model\n",
    "md = pyissm.model.Model()\n",
    "\n",
    "# View all model sub-classes by viewing the __repr__ of the model\n",
    "md"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7587545f",
   "metadata": {},
   "source": [
    "We can access any model sub-class using standard Python dot notation. By default, interactively interrogating `md` will print a summary of the model, or specific model sub-class. When interrogating a model object, we see all available sub-classes (like the above example). when interrigating a model sub-class, we see all available fields/parameters within that sub-class. If a field/parameter is empty, `N/A` is displayed.\n",
    "\n",
    "Let's have a look at pre-defined parameters in the `md.constants` sub-class, and empty fields in the `md.geometry` sub-class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0eb4c7c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   constants parameters:\n",
       "         g                      : 9.81            -- gravitational acceleration [m/s^2]\n",
       "         omega                  : 7.292e-05       -- angular velocity of Earth [rad/s]\n",
       "         yts                    : 31536000.0      -- number of seconds in a year [s/yr]\n",
       "         referencetemperature   : 223.15          -- reference temperature used in the enthalpy model [K]\n",
       "         gravitational_constant : 6.67259e-11     -- Newtonian constant of gravitation [m^3/kg/s^2]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md.constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f91b1d2a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   geometry parameters:\n",
       "         surface                : N/A             -- ice upper surface elevation [m]\n",
       "         thickness              : N/A             -- ice thickness [m]\n",
       "         base                   : N/A             -- ice base elevation [m]\n",
       "         bed                    : N/A             -- bed elevation [m]\n",
       "         hydrostatic_ratio      : N/A             -- hydrostatic ratio for floating ice"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "md.geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4edb27",
   "metadata": {},
   "source": [
    "In the next cell, all available sub-classes are listed and commented out. Spend a few minutes interrogating each of them and familiarizing yourself with the structure of a pyISSM model. You can simply uncomment one line at a time and re-run the cell to view the fields/parameters included in each sub-class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "115a7ec4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# md.amr\n",
    "# md.autodiff\n",
    "# md.balancethickness\n",
    "# md.basalforcings\n",
    "# md.calving\n",
    "# md.cluster\n",
    "# md.constants\n",
    "# md.damage\n",
    "# md.debris\n",
    "# md.debug\n",
    "# md.dsl\n",
    "# md.esa\n",
    "# md.flowequation\n",
    "# md.friction\n",
    "# md.frontalforcings\n",
    "# md.geometry\n",
    "# md.groundingline\n",
    "# md.hydrology\n",
    "# md.initialization\n",
    "# md.inversion\n",
    "# md.levelset\n",
    "# md.love\n",
    "# md.mask\n",
    "# md.masstransport\n",
    "# md.materials\n",
    "# md.mesh\n",
    "# md.miscellaneous\n",
    "# md.outputdefinition\n",
    "# md.private\n",
    "# md.qmu\n",
    "# md.radaroverlay\n",
    "# md.results\n",
    "# md.rifts,\n",
    "# md.sampling\n",
    "# md.settings\n",
    "# md.smb\n",
    "# md.solidearth\n",
    "# md.steadystate\n",
    "# md.stochasticforcing\n",
    "# md.stressbalance\n",
    "# md.thermal\n",
    "# md.timestepping\n",
    "# md.toolkits\n",
    "# md.transient\n",
    "# md.verbose"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a84cab0a",
   "metadata": {},
   "source": [
    "## Alternative parameterisations\n",
    "\n",
    "Within ISSM, there are multiple ways to represent various processes, depending on the use case. For example, modellers may wish to represent basal friction using the Budd friction law, or the Schoof friction law. Or, modellers may wish to use the Positive Degree Day (PDD) surface mass balance (SMB) approach, rather than a prescribed SMB. All of this flexibility is embedded in pyISSM and accessible via different model classes.\n",
    "\n",
    "Let's take a look at a few different friction classes, accessible via `pyissm.model.classes.friction....`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d553b337",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n",
       "(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n",
       "         coefficient            : N/A             -- friction coefficient [SI]\n",
       "         p                      : N/A             -- p exponent\n",
       "         q                      : N/A             -- q exponent\n",
       "         coupling               : 0               -- Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)\n",
       "         linearize              : 0               -- 0: not linearized, 1: interpolated linearly, 2: constant per element (default is 0)\n",
       "         effective_pressure     : N/A             -- Effective Pressure for the forcing if not coupled [Pa]\n",
       "         effective_pressure_l...: 0               -- Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# By default, md.friction uses the default friction law, which is a Budd friction law\n",
    "md.friction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "99c844a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Weertman sliding law parameters: Sigma_b = C^(- 1 / m) * |u_b|^(1 / m - 1) * u_b\n",
       "         C                      : N/A             -- friction coefficient [SI]\n",
       "         m                      : N/A             -- m exponent\n",
       "         linearize              : 0               -- 0: not linearized, 1: interpolated linearly, 2: constant per element (default is 0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To change the friction law to a Weertman friction law, we can do:\n",
    "md.friction = pyissm.model.classes.friction.weertman()\n",
    "md.friction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "453cb12c",
   "metadata": {},
   "source": [
    "The same appraoch allows users to adjust the SMB parameterisation (as well as other model sub-classes). Let's look at a few different SMB options, accessible via `pyissm.model.classes.smb....` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "543debdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   surface forcings parameters:\n",
       "         isdelta18o             : 0               -- is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)\n",
       "         ismungsm               : 0               -- is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)\n",
       "         issetpddfac            : 0               -- is user passing in defined pdd factors (0 or 1, default is 0)\n",
       "         desfac                 : 0.5             -- desertification elevation factor (between 0 and 1, default is 0.5) [m]\n",
       "         s0p                    : N/A             -- should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]\n",
       "         s0t                    : N/A             -- should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]\n",
       "         rlaps                  : 6.5             -- present day lapse rate [degree/km]\n",
       "         rlapslgm               : 6.5             -- LGM lapse rate [degree/km]\n",
       "         monthlytemperatures    : N/A             -- monthly surface temperatures [K], required if pdd is activated and delta18o not activated\n",
       "         precipitation          : N/A             -- monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated\n",
       "         delta18o               : 0               -- delta18o [per mil], required if pdd is activated and delta18o activated\n",
       "         delta18o_surface       : N/A             -- surface elevation of the delta18o site, required if pdd is activated and delta18o activated\n",
       "         temperatures_presentday: N/A             -- monthly present day surface temperatures [K], required if delta18o/mungsm is activated\n",
       "         temperatures_lgm       : N/A             -- monthly LGM surface temperatures [K], required if delta18o or mungsm is activated\n",
       "         precipitations_prese...: N/A             -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated\n",
       "         precipitations_lgm     : N/A             -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated\n",
       "         Tdiff                  : N/A             -- time interpolation parameter for temperature, 1D(year), required if mungsm is activated\n",
       "         sealev                 : N/A             -- sea level [m], 1D(year), required if mungsm is activated\n",
       "         temperatures_presentday: N/A             -- monthly present day surface temperatures [K], required if delta18o/mungsm is activated\n",
       "         temperatures_lgm       : N/A             -- monthly LGM surface temperatures [K], required if delta18o or mungsm is activated\n",
       "         precipitations_prese...: N/A             -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated\n",
       "         precipitations_lgm     : N/A             -- monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated\n",
       "         Pfac                   : N/A             -- time interpolation parameter for precipitation, 1D(year), required if mungsm is activated\n",
       "         Tdiff                  : N/A             -- time interpolation parameter for temperature, 1D(year), required if mungsm is activated\n",
       "         sealev                 : N/A             -- sea level [m], 1D(year), required if mungsm is activated\n",
       "         steps_per_step         : 1               -- number of smb steps per time step\n",
       "         averaging              : 0               -- averaging methods from short to long steps\n",
       "\t\t0: Arithmetic (default)\n",
       "\t\t1: Geometric\n",
       "\t\t2: Harmonic\n",
       "         requested_outputs      : ['default',]    -- additional outputs requested"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# By default, md.smb uses the default surface mass balance model which uses a prescribed smb field\n",
    "md.smb\n",
    "\n",
    "# To change the smb model to a positive degree day smb model, we can do:\n",
    "md.smb = pyissm.model.classes.smb.pdd()\n",
    "md.smb"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-26.01] *",
   "language": "python",
   "name": "conda-env-analysis3-26.01-py"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}