new exercise notebooks

.gitignore

@@ -19,3 +19,5 @@ usr
 .env
 .iga-python
 *_autosummary*
+
+.DS_Store

_toc.yml

@@ -66,6 +66,10 @@ parts:
       - file: chapter3/navier-stokes-steady
         sections:
           - file: chapter3/navier-stokes-steady-streamfunction-velocity
+  - caption: Exercises
+    chapters:
+      - file: exercises/harmonic-fem
+      - file: exercises/harmonic-feec
   - caption: Advanced subjects 
     chapters:
       - file: chapter4/subdomains

exercises/electro_static_2d.ipynb

@@ -0,0 +1,285 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Note: This notebook must be changed by the student. As it is now, it does not approximate the solution to the Poisson problem!**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 2D Solver for Electro-Static Problem on the Unit Square\n",
+    "\n",
+    "In this exercise we write a solver for the 2D electro-static problem on the unit square, which can be shown to be equivalent to an $\\mathcal{L}^1$ Hodge-Laplace problem.\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    &\\text{Find }\\boldsymbol{E} \\in D(\\mathcal{L}^1) \\\\\n",
+    "    &\\mathcal{L}^1\\boldsymbol{E} = (\\boldsymbol{curl}\\ curl - \\boldsymbol{grad}\\ div)\\boldsymbol{E} = -\\boldsymbol{grad}\\ \\rho \\quad \\text{ in }\\Omega=]0,1[^2,\n",
+    "\\end{align*}\n",
+    "for the specific choice\n",
+    "\\begin{align*}\n",
+    "    \\rho(x, y) = 2\\sin(\\pi x)\\sin(\\pi y)\n",
+    "\\end{align*}\n",
+    "\n",
+    "## The Discrete Variational Formulation\n",
+    "\n",
+    "The corresponding discrete variational formulation reads\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    &\\text{Find }\\boldsymbol{E} \\in V_h^1 \\subset V^1 = H_0(curl;\\Omega)\\text{ such that } \\\\\n",
+    "    &a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V_h^1\n",
+    "\\end{align*}\n",
+    "\n",
+    "where \n",
+    "- $a(\\boldsymbol{E},\\boldsymbol{v}) \\coloneqq \\int_{\\Omega} (curl\\ \\boldsymbol{E})(curl\\ \\boldsymbol{v}) + (\\widetilde{div}_h\\boldsymbol{E})(\\widetilde{div}_h\\boldsymbol{v}) ~ d\\Omega$ ,\n",
+    "- $l(\\boldsymbol{v}) := -\\int_{\\Omega} \\boldsymbol{grad}\\ \\rho\\cdot\\boldsymbol{v} ~ d\\Omega$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Discrete Model using de Rham objects"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympde.calculus      import dot\n",
+    "from sympde.topology      import elements_of, Square, Derham\n",
+    "from sympde.expr.expr     import BilinearForm, integral\n",
+    "\n",
+    "domain  = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n",
+    "derham  = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n",
+    "\n",
+    "V0      = derham.V0\n",
+    "V1      = derham.V1\n",
+    "V2      = derham.V2\n",
+    "\n",
+    "u0, v0  = elements_of(V0, names='u0, v0')\n",
+    "u1, v1  = elements_of(V1, names='u1, v1')\n",
+    "u2, v2  = elements_of(V2, names='u2, v2')\n",
+    "\n",
+    "# bilinear forms corresponding to V0, V1 and V2 mass matrices\n",
+    "m0      = BilinearForm((u0, v0), integral(domain,      u0 * v0))\n",
+    "m1      = BilinearForm((u1, v1), integral(domain,      dot(u1, v1)))\n",
+    "m2      = BilinearForm((u2, v2), integral(domain,      u2 * v2))\n",
+    "\n",
+    "# callable source term \\rho\n",
+    "from sympy  import pi, sin, cos, lambdify, Tuple, Matrix\n",
+    "x,y             = domain.coordinates\n",
+    "rho             = 2*sin(pi*x)*sin(pi*y)\n",
+    "rho_lambdified  = lambdify((x, y), rho)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from psydac.api.discretization  import discretize\n",
+    "from psydac.api.settings        import PSYDAC_BACKEND_GPYCCEL\n",
+    "\n",
+    "backend = PSYDAC_BACKEND_GPYCCEL"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ncells  = [32, 32]  # Bspline cells\n",
+    "degree  = [4, 4]    # Bspline degree"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# discretize domain and derham\n",
+    "# <-\n",
+    "# <-\n",
+    "\n",
+    "# define FEM spaces V0_h, V1_h and V2_h\n",
+    "#V0_h        = derham_h.V0\n",
+    "#V1_h        = derham_h.V1\n",
+    "#V2_h        = derham_h.V2\n",
+    "\n",
+    "# Commuting projection operators\n",
+    "#P0, P1, P2  = derham_h.projectors()\n",
+    "\n",
+    "# Exterior Derivative operators (grad and curl)\n",
+    "# <-\n",
+    "\n",
+    "# Mass matrices\n",
+    "# <-\n",
+    "# <-\n",
+    "# <-\n",
+    "\n",
+    "# <-\n",
+    "# <-\n",
+    "# <-\n",
+    "\n",
+    "# Boundary Projectors\n",
+    "from utils                  import H1BoundaryProjector2D, HcurlBoundaryProjector2D\n",
+    "from psydac.linalg.basic    import IdentityOperator\n",
+    "\n",
+    "#P_H1        = H1BoundaryProjector2D(V0, V0_h.coeff_space)\n",
+    "#P_Hcurl     = HcurlBoundaryProjector2D(V1, V1_h.coeff_space)\n",
+    "\n",
+    "#I0          = IdentityOperator(V0_h.coeff_space)\n",
+    "#I1          = IdentityOperator(V1_h.coeff_space)\n",
+    "\n",
+    "#P_H1_Gamma      = I0 - P_H1\n",
+    "#P_Hcurl_Gamma   = I1 - P_Hcurl\n",
+    "\n",
+    "# The inner product < div v, div E > poses a difficulty for the projection method.\n",
+    "# In order for us to obtain the correct result, we must modify the transposed gradient\n",
+    "# as well as the inverse M0 mass matrix\n",
+    "from psydac.linalg.solvers  import inverse\n",
+    "\n",
+    "#Gt_0 = P_H1 @ G.T\n",
+    "\n",
+    "#M0_0 = P_H1 @ M0 @ P_H1 + P_H1_Gamma\n",
+    "#M0_inv = inverse(M0_0, 'cg', maxiter=1000, tol=1e-11)\n",
+    "\n",
+    "# System Matrix A and A_bc\n",
+    "#A           = # <-\n",
+    "#A_bc        = P_Hcurl @ A @ P_Hcurl + P_Hcurl_Gamma\n",
+    "\n",
+    "# rhs vector f and f_bc\n",
+    "#rho_coeffs  = P0(rho_lambdified).coeffs\n",
+    "#f           = # <-\n",
+    "#f_bc        = P_Hcurl @ f"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import  time\n",
+    "\n",
+    "tol     = 1e-10\n",
+    "maxiter = 1000\n",
+    "\n",
+    "# We use a preconditioner for acceleration\n",
+    "from utils import get_M1_block_kron_solver_2D\n",
+    "#pc          = get_M1_block_kron_solver_2D(V1_h.coeff_space, ncells, degree, periodic=[False, False])\n",
+    "#M1_0        = P_Hcurl @ M1 @ P_Hcurl + P_Hcurl_Gamma\n",
+    "#M1_0_inv    = inverse(M1_0, 'pcg', pc=pc, tol=1e-11, maxiter=1000)\n",
+    "#A_bc_inv    = inverse(A_bc, 'pcg', pc=M1_0_inv, tol=tol, maxiter=maxiter)\n",
+    "\n",
+    "t0      = time.time()\n",
+    "#E_h     = A_bc_inv @ f_bc\n",
+    "t1      = time.time()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Computing the error norm\n",
+    "\n",
+    "As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n",
+    "In this example, the analytical solution is given by\n",
+    "\n",
+    "$$\n",
+    "\\boldsymbol{E}_{ex}(x, y) = -\\frac{1}{\\pi}\\left(\\begin{matrix} \n",
+    "                \\cos(\\pi x) \\sin(\\pi y) \\\\\n",
+    "                \\cos(\\pi y) \\sin(\\pi x)\n",
+    "                \\end{matrix}\\right)\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from psydac.fem.basic   import FemField\n",
+    "\n",
+    "from sympde.expr        import Norm\n",
+    "\n",
+    "E_ex            = Tuple( (-1/pi)*cos(pi*x) * sin(pi*y),\n",
+    "                         (-1/pi)*cos(pi*y) * sin(pi*x) )\n",
+    "#E_h_FemField    = FemField(V1_h, E_h)\n",
+    "\n",
+    "error           = Matrix([u1[0] - E_ex[0], u1[1] - E_ex[1]])\n",
+    "\n",
+    "# create the formal Norm object\n",
+    "l2norm          = Norm(error, domain, kind='l2')\n",
+    "\n",
+    "# discretize the norm\n",
+    "#l2norm_h        = discretize(l2norm, domain_h, V1_h, backend=backend)\n",
+    "\n",
+    "# assemble the norm\n",
+    "#l2_error        = l2norm_h.assemble(u1=E_h_FemField)\n",
+    "\n",
+    "# print the result\n",
+    "#print( '> Grid          :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
+    "#print( '> Degree        :: [{p1},{p2}]'  .format( p1=degree[0], p2=degree[1] ) )\n",
+    "#print( '> CG info       :: ',A_bc_inv.get_info() )\n",
+    "#print( '> L2 error      :: {:.2e}'.format( l2_error ) )\n",
+    "#print( '' )\n",
+    "#print( '> Solution time :: {:.3g}'.format( t1-t0 ) )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualization\n",
+    "\n",
+    "We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from utils import plot\n",
+    "\n",
+    "E_ex_x   = lambdify((x, y), E_ex[0])\n",
+    "E_ex_y   = lambdify((x, y), E_ex[1])\n",
+    "#E_h_x    = E_h_FemField[0]\n",
+    "#E_h_y    = E_h_FemField[1]\n",
+    "#error_x  = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n",
+    "#error_y  = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n",
+    "\n",
+    "#plot(gridsize_x     = 100, \n",
+    "#    gridsize_y     = 100, \n",
+    "#    title          = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n",
+    "#    funs           = [E_ex_x, E_h_x, error_x], \n",
+    "#    titles         = [r'$(\\boldsymbol{E}_{ex})_1(x,y)$', r'$(\\boldsymbol{E}_h)_1(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_1(x,y)|$'],\n",
+    "#    surface_plot   = True\n",
+    "#)\n",
+    "\n",
+    "#plot(gridsize_x     = 100,\n",
+    "#    gridsize_y     = 100,\n",
+    "#    title          = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n",
+    "#    funs           = [E_ex_y, E_h_y, error_y],\n",
+    "#    titles         = [r'$(\\boldsymbol{E}_{ex})_2(x,y)$', r'$(\\boldsymbol{E}_h)_2(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_2(x,y)|$'],\n",
+    "#    surface_plot   = True\n",
+    "#)"
+   ]
+  }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 4
+}