{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Offline Phase: SGREEDY algorithm for sensor placement\n",
    "\n",
    "**Aim of the tutorial:** this notebook shows how to use the SGREEDY algorithm for sensor placement starting from a reduced basis, e.g., the POD modes.\n",
    "\n",
    "----\n",
    "\n",
    "*To execute this notebook* it is necessary to have the POD modes stored in `Offline_results/BasisFunctions` folder, placed in this directory (otherwise modify `path_off` variable)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "from IPython.display import clear_output\n",
    "import pickle\n",
    "\n",
    "from dolfinx.fem import FunctionSpace\n",
    "\n",
    "from pyforce.tools.write_read import ImportH5, StoreFunctionsList\n",
    "from pyforce.tools.functions_list import FunctionsList\n",
    "\n",
    "path_off ='./Offline_results/'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geometry is imported from \"ANL11A2_octave.geo\", generated with GMSH. Then, the mesh is created with the gmsh module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neutronics import create_anl11a2_mesh\n",
    "\n",
    "domain, _, _ = create_anl11a2_mesh(use_msh=True, save_mesh=False)\n",
    "\n",
    "fuel1_marker    = 1\n",
    "fuel2_marker    = 2\n",
    "fuel_rod_marker = 3\n",
    "refl_marker     = 4\n",
    "\n",
    "void_marker     = 10\n",
    "sym_marker      = 20\n",
    "\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us import the POD modes, using the `ImportH5` function from the `pyforce` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the functional space\n",
    "V = FunctionSpace(domain, (\"Lagrange\", 1))\n",
    "\n",
    "# Define the variables to load\n",
    "var_names = [\n",
    "             'phi_1',\n",
    "             'phi_2'\n",
    "             ]\n",
    "\n",
    "tex_var_names = [\n",
    "                 r'\\phi_1',\n",
    "                 r'\\phi_2'\n",
    "                 ]\n",
    "\n",
    "bf = dict()\n",
    "\n",
    "for field_i in range(len(var_names)):\n",
    "    \n",
    "    bf[var_names[field_i]] = ImportH5(V, path_off+'/BasisFunctions/basisPOD_'+var_names[field_i], 'POD_'+var_names[field_i])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SGREEDY for POD basis\n",
    "The SGREEDY algorithm, [Maday et al. (2014)](https://doi.org/10.1002/nme.4747), aims at maximising the *inf-sup* constant $\\beta_{N,M}$: roughly speaking, this constant measures how the sensors can include unmodelled physics with respect to the one embedded in the reduced basis.\n",
    "\n",
    "The overall algorithm is summarised in the following figure from [Riva et al. (2024)](https://doi.org/10.1016/j.apm.2024.06.040)\n",
    "\n",
    "![SGREEDY](sgreedy-algo.png)\n",
    "\n",
    "As for the GEIM method, three different values of the point spread will be considered: 0.1, 1 and 2.5.\n",
    "The SGREEDY method is implemented in the `pyforce` package, and it is called by `SGREEDY`: the class must be initialised with the domain, the basis functions (the POD modes in this case), the functional space in which the snapshots live, the name of the variable and the value of the point spread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyforce.offline.sensors import SGREEDY\n",
    "\n",
    "s = 2.\n",
    "sgreedy_pod = dict()\n",
    "\n",
    "for field_i in range(len(var_names)):\n",
    "    sgreedy_pod[var_names[field_i]] = [SGREEDY(domain, bf[var_names[field_i]], V, var_names[field_i], s)]*2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the snapshots belong to $H^1\\subset L^2$, the Riesz representation of a functional is sought in this space, endowed with the inner product\n",
    "\\begin{equation*}\n",
    "\\left(u,v\\right)_{H^1}=\\int_\\Omega \\nabla u\\cdot \\nabla v\\,d\\Omega + \\int_\\Omega u\\cdot v\\,d\\Omega\n",
    "\\end{equation*}\n",
    "The results will be compared with the approximation properties of the case with the Riesz representation in $L^2$, whose inner product is \n",
    "\\begin{equation*}\n",
    "\\left(u,v\\right)_{L^2} = \\int_\\Omega u\\cdot v\\,d\\Omega\n",
    "\\end{equation*}\n",
    "\n",
    "The SGREEDY algorithm is called by the method `generate` from the `SGREEDY` class: it requires the dimension of the reduced space to use, the maximum number of sensors to place and a tolerance value for the stability-loop of the algorithm. Moreover, the method can have the optional parameters `verbose` and `sampleEvery`, which allow to print the progress and telling how many cells are to be used, respectively. In the end, there is a switch option for the use of the Riesz representation in $H^1$ or $L^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SGREEDY for phi_1 with s=2.00 and Riesz representation in L2\n",
      " \n",
      "SGREEDY for phi_1 with s=2.00 and Riesz representation in H1\n",
      " \n",
      "SGREEDY for phi_2 with s=2.00 and Riesz representation in L2\n",
      " \n",
      "SGREEDY for phi_2 with s=2.00 and Riesz representation in H1\n",
      " \n"
     ]
    }
   ],
   "source": [
    "Nmax = 10\n",
    "Mmax = 20\n",
    "\n",
    "sam_every = 4\n",
    "\n",
    "is_H1 = [False, True]\n",
    "fun_space_label = ['L2', 'H1']\n",
    "\n",
    "for field_i in range(len(var_names)):\n",
    "        for kk in range(len(is_H1)):\n",
    "            \n",
    "            print('SGREEDY for '+var_names[field_i]+' with s={:.2f}'.format(s)+' and Riesz representation in '+fun_space_label[kk])\n",
    "            \n",
    "            sgreedy_pod[var_names[field_i]][kk].generate(   Nmax, Mmax, tol = 0.2, \n",
    "                                                            sampleEvery = sam_every, \n",
    "                                                            verbose=False, is_H1 = is_H1[kk])\n",
    "            \n",
    "            print(' ')\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can store the sensors using the `StoreFunctionsList` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for field_i in range(len(var_names)):\n",
    "        for kk in range(len(is_H1)):\n",
    "                StoreFunctionsList( domain, sgreedy_pod[var_names[field_i]][kk].basis_sens, \n",
    "                                    'SGREEDYPOD_' +var_names[field_i]+'_s_{:.2e}'.format(s),\n",
    "                                    path_off+'/BasisSensors/sensorsSGREEDYPOD_' + var_names[field_i]+'_s_{:.2e}_'.format(s)+fun_space_label[kk])"
   ]
  }
 ],
 "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.14"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e5b3af111de828d84cfea01c8e4cc990d7262e947155f31e208c22ad6b78199b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}