examples/notebooks/Helmholtz_non_periodic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving Helmholtz's equation on non-periodic topological domains\n",
    "In this notebook, we will show how to solve Helmholtz's equation on single patch and multipatch domains. As we are running this in a notebook, everything will be run in serial and hence we are limiting ourselves to a fairly coarse discretization to avoid taking too much time. However, Psydac allows for hybrid MPI + OpenMP parallelization with barely any changed to the code. The lines that are impacted by that will be preceeded by their MPI equivalent\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 : Building the domain\n",
    "Psydac uses the powerful topological tools of SymPDE to build a large variety of domains. Here we show how to create a Square."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sympde.topology import Square\n",
    "Omega=Square()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Defining the Abstract PDE model using SymPDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympde.calculus import grad, dot\n",
    "from sympde.calculus import minus, plus\n",
    "\n",
    "from sympde.expr.expr          import LinearForm, BilinearForm\n",
    "from sympde.expr.expr          import integral              \n",
    "from sympde.expr.expr          import Norm                       \n",
    "from sympde.expr               import find, EssentialBC\n",
    "\n",
    "from sympde.topology import ScalarFunctionSpace\n",
    "from sympde.topology import elements_of\n",
    "from sympde.topology import NormalVector\n",
    "from sympde.topology import Union\n",
    "from sympy import pi, cos, sin, symbols, conjugate, exp, sqrt\n",
    "\n",
    "# Define the abstract model to solve Helmholtz's equation using the manufactured solution method\n",
    "\n",
    "x, y = Omega.coordinates\n",
    "kappa = 2*pi\n",
    "solution = exp(1j * kappa * x) * sin(kappa * y)\n",
    "e_w_0 = sin(kappa * y) # value of incoming wave at x=0, forall y\n",
    "dx_e_w_0 = 1j*kappa*sin(kappa * y) # derivative wrt. x of incoming wave at x=0, forall y\n",
    "\n",
    "V   = ScalarFunctionSpace('V', Omega, kind=None)\n",
    "V.codomain_type='complex'\n",
    "\n",
    "u, v = elements_of(V, names='u, v')\n",
    "\n",
    "error  = u - solution\n",
    "\n",
    "expr   = dot(grad(u),grad(v)) - 2 * kappa ** 2 * u * v\n",
    "\n",
    "#we impose an incoming wave from the left and absorbing boundary conditions at the right.\n",
    "boundary_expr = - 1j * kappa * u * v\n",
    "x_boundary = Union(Omega.get_boundary(axis=0, ext=-1), Omega.get_boundary(axis=0, ext=1))\n",
    "boundary_source_expr = - dx_e_w_0 * v - 1j * kappa * e_w_0 * v    \n",
    "\n",
    "a = BilinearForm((u, v), integral(Omega, expr) + integral(x_boundary, boundary_expr))\n",
    "l = LinearForm(v, integral(Omega.get_boundary(axis=0, ext=-1), boundary_source_expr))\n",
    "\n",
    "equation = find(u, forall=v, lhs=a(u,v), rhs=l(v))\n",
    "\n",
    "l2norm = Norm(error, Omega, kind='l2')\n",
    "h1norm = Norm(error, Omega, kind='h1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Discretizing the domain, spaces and equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Discretize the geometry and equation\n",
    "from psydac.api.discretization import discretize\n",
    "from psydac.api.settings import PSYDAC_BACKENDS\n",
    "\n",
    "# The backends are as follows and can all be run with OpenMP by setting backend['omp'] = True.\n",
    "# For all but the Python backend, compilation flags can be accessed and changed using backend['flags'].\n",
    "# - Python backend: Generate and runs python files. Accessed via PSYDAC_BACKENDS['python']\n",
    "# - Pyccel GCC backend: Accessed via PSYDAC_BACKENDS['pyccel-gcc']\n",
    "# - Pyccel Intel backend: Accessed via PSYDAC_BACKENDS['pyccel-intel']\n",
    "# - Pyccel PGI backend: Accessed via PSYDAC_BACKENDS['pyccel-pgi']\n",
    "# - Pyccel Numba backend: Accessed via PSYDAC_BACKENDS['numba']\n",
    "\n",
    "backend = PSYDAC_BACKENDS['python']\n",
    "\n",
    "# Uncomment to use OpenMp\n",
    "# import os\n",
    "# os.environ['OMP_NUM_THREADS'] = \"4\"\n",
    "# backend['omp'] = True\n",
    "\n",
    "ncells = [30, 30]\n",
    "degree = [3, 3]\n",
    "periodic = [False, True]\n",
    "\n",
    "nquads = [p + 1 for p in degree]\n",
    "\n",
    "# MPI version\n",
    "# from mpi4py import MPI\n",
    "# comm = MPI.COMM_WORLD\n",
    "# Omega_h = discretize(Omega, ncells=ncells, comm=comm)\n",
    "Omega_h = discretize(Omega, ncells=ncells, periodic=periodic)\n",
    "\n",
    "Vh         = discretize(V, Omega_h, degree=degree)\n",
    "equation_h = discretize(equation, Omega_h, [Vh, Vh], nquads=nquads, backend=backend)\n",
    "l2norm_h   = discretize(l2norm, Omega_h, Vh, nquads=nquads, backend=backend)\n",
    "h1norm_h   = discretize(h1norm, Omega_h, Vh, nquads=nquads, backend=backend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Solving the equation and computing the error norms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> CG info       ::  {'niter': 46, 'success': True, 'res_norm': 5.330478803717119e-15}\n",
      "> L2 error      :: 2.07e-07\n",
      "> H1 error      :: 3.80e-05\n",
      "> Solution time :: 3.76e+00s\n",
      "> Evaluat. time :: 4.00e-01s \n"
     ]
    }
   ],
   "source": [
    "# Set the solver parameters\n",
    "# 'cbig' -> Biconjugate gradient method\n",
    "equation_h.set_solver('bicg', info=True, tol=1e-14)\n",
    "\n",
    "import time\n",
    "\n",
    "t0_s     = time.time()\n",
    "uh, info = equation_h.solve()\n",
    "t1_s     = time.time()\n",
    "\n",
    "\n",
    "t0_d     = time.time()\n",
    "l2_error = l2norm_h.assemble(u=uh)\n",
    "h1_error = h1norm_h.assemble(u=uh)\n",
    "t1_d     = time.time()\n",
    "\n",
    "print( '> CG info       :: ',info )\n",
    "print( '> L2 error      :: {:.2e}'.format(l2_error))\n",
    "print( '> H1 error      :: {:.2e}'.format(h1_error))\n",
    "print( '> Solution time :: {:.2e}s'.format(t1_s - t0_s))\n",
    "print( '> Evaluat. time :: {:.2e}s '.format(t1_d - t0_d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Saving the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the results using OutputManager\n",
    "from psydac.api.postprocessing import OutputManager\n",
    "import os\n",
    "\n",
    "os.makedirs('results_Helmholtz', exist_ok=True)\n",
    "\n",
    "Om = OutputManager(\n",
    "    f'results_Helmholtz/space_info_{Omega.name}',\n",
    "    f'results_Helmholtz/field_info_{Omega.name}',\n",
    "# MPI version\n",
    "#   comm=comm,    \n",
    ")\n",
    "\n",
    "Om.add_spaces(V=Vh)\n",
    "Om.export_space_info()\n",
    "\n",
    "Om.set_static()\n",
    "Om.export_fields(u=uh)\n",
    "\n",
    "Om.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Exporting the results to VTK and visualizing them with Paraview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the results to VTK using PostProcessManager\n",
    "from psydac.api.postprocessing import PostProcessManager\n",
    "from sympy import lambdify\n",
    "\n",
    "Pm = PostProcessManager(\n",
    "    domain=Omega,\n",
    "    space_file=f'results_Helmholtz/space_info_{Omega.name}.yaml',\n",
    "    fields_file=f'results_Helmholtz/field_info_{Omega.name}.h5',\n",
    "# MPI version\n",
    "#   comm=comm,    \n",
    ")\n",
    "\n",
    "# The complex fields will be exported in two field: a real and an imaginary field\n",
    "Pm.export_to_vtk(\n",
    "    f'results_Helmholtz/visu_{Omega.name}',\n",
    "    grid=None,\n",
    "    npts_per_cell=3,\n",
    "    fields='u',\n",
    "    additional_physical_functions={'exact_solution': lambdify(Omega.coordinates, solution, modules='numpy')}\n",
    ")\n",
    "\n",
    "Pm.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Results on the square domain\n",
    "\n",
    "<p float=\"left\">\n",
    "<table><tr>\n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Real_Computed.png\" width=\"100%\" /> \n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Imag_Computed.png\" width=\"100%\" />\n",
    "</tr></table>\n",
    "</p>\n",
    "\n",
    "<p float=\"left\">\n",
    "<table><tr>\n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Imag_exact.png\" width=\"90%\" /> \n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Real_exact.png\" width=\"100%\" /> \n",
    "</tr></table> \n",
    "</p>\n",
    "      \n",
    "<p float=\"left\">\n",
    "<table><tr>\n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Real_error.png\" width=\"100%\" />\n",
    "  <td><img src=\"./paraview_images/Helmholtz_Square_Imag_error.png\" width=\"100%\" />\n",
    "</tr></table> \n",
    "</p>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "v_psydac",
   "language": "python",
   "name": "v_psydac"
  },
  "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.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}