rbspacetest.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from nutils import mesh, function as fn, plot, matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_elements = 10\n",
    "degree = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create the original mesh and basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pts = np.linspace(0, 1, n_elements + 1)\n",
    "domain, geom = mesh.rectilinear([pts, pts])\n",
    "basis = domain.basis('spline', degree=(degree, degree))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's form the stiffness matrix in the original basis. I call `.toarray()` to get a Numpy matrix object. (The matrix is sparse but small so I'm doing it this way for convenience.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "integrand = fn.outer(basis.grad(geom)).sum([-1])\n",
    "orig_mx = domain.integrate(integrand, geometry=geom, ischeme='gauss9').toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we proceed we need the indices of the boundary and the interior DoFs. To do this we project a boundary function (zero in this case) onto the basis. All the interior DoFs should be NaN in the resulting vector, while the boundary DoFs should have actual numbers. I use 1-point quadrature because the function is constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "project > constrained 48/169 dofs, error 0.00e+00/area\n"
     ]
    }
   ],
   "source": [
    "cons = domain.boundary.project(0, onto=basis, geometry=geom, ischeme='gauss1')\n",
    "interior = np.isnan(cons)             # A vector of booleans\n",
    "boundary = np.logical_not(interior)   # ... same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do the eigenvalue decomposition of the interior part of the matrix. Numpy's `eigh` function returns eigenvalues in ascending order, whereas we would like them to be descending, so let's also reverse that order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "interior_mx = orig_mx[np.ix_(interior, interior)]\n",
    "eigvals, Vred = np.linalg.eigh(interior_mx)\n",
    "eigvals = eigvals[::-1]\n",
    "Vred = Vred[:,::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, `Vred` is square, being $n \\times n$, where $n$ is the number of interior DoFs. We would like to expand the columns so that it's $N \\times n$, where $N$ is the number of total DoFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "V = np.zeros((len(basis), Vred.shape[1]))\n",
    "V[interior,:] = Vred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to compute the stiffness matrix in the \"reduced\" space. It is ${\\boldsymbol A}_\\text{red} = {\\boldsymbol V}^\\intercal {\\boldsymbol A}_\\text{orig} {\\boldsymbol V}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "red_mx = V.T.dot(orig_mx).dot(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should be equal to the diagonal eigenvalue matrix. Why? Well, if we reorder the basis to represent $\\boldsymbol A_\\text{orig}$ as\n",
    "\n",
    "$$ {\\boldsymbol A}_\\text{orig} = \\begin{pmatrix} \n",
    "{\\boldsymbol A}_\\text{ii} & {\\boldsymbol A}_\\text{ib} \\\\\n",
    "{\\boldsymbol A}_\\text{bi} & {\\boldsymbol A}_\\text{bb}\n",
    "\\end{pmatrix} $$\n",
    "\n",
    "where subscript i and b refers to *interior* and *boundary*. We know that ${\\boldsymbol A}_\\text{ii} = {\\boldsymbol V}_\\text{red} {\\boldsymbol \\Lambda} {\\boldsymbol V}^\\intercal_\\text{red}$ and \n",
    "\n",
    "$${\\boldsymbol V} = \\begin{pmatrix} {\\boldsymbol V}_\\text{red} \\\\ {\\boldsymbol 0} \\end{pmatrix}$$\n",
    "\n",
    "so we get\n",
    "\n",
    "$$ {\\boldsymbol V}^\\intercal {\\boldsymbol A}_\\text{orig} {\\boldsymbol V}\n",
    "= \\begin{pmatrix} {\\boldsymbol V}^\\intercal_\\text{red} & {\\boldsymbol 0} \\end{pmatrix}\n",
    "\\begin{pmatrix} \n",
    "{\\boldsymbol V}_\\text{red} {\\boldsymbol \\Lambda} {\\boldsymbol V}^\\intercal_\\text{red} & {\\boldsymbol A}_\\text{ib} \\\\\n",
    "{\\boldsymbol A}_\\text{bi} & {\\boldsymbol A}_\\text{bb}\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix} {\\boldsymbol V}_\\text{red} \\\\ {\\boldsymbol 0} \\end{pmatrix}\n",
    "= \\begin{pmatrix} {\\boldsymbol V}^\\intercal_\\text{red} & {\\boldsymbol 0} \\end{pmatrix}\n",
    "\\begin{pmatrix} {\\boldsymbol V}_\\text{red} {\\boldsymbol \\Lambda} \\\\ {\\boldsymbol A}_\\text{bi} {\\boldsymbol V}_\\text{red} \\end{pmatrix}\n",
    "= {\\boldsymbol \\Lambda}\n",
    "$$\n",
    "\n",
    "(Remember, $\\boldsymbol{V}_\\text{red}$ is an orthogonal matrix.) Let's check that this holds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(red_mx, np.diag(eigvals))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great. Now let's create a random right hand side for solving the system. We'll be using homogeneous Dirichlet boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "orig_rhs = np.random.rand(len(basis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, by using Nutils' own method with the `constrain` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_lhs = matrix.NumpyMatrix(orig_mx).solve(orig_rhs, constrain=cons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, by doing the constraining manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "manual_lhs = np.zeros((len(basis),))\n",
    "manual_lhs[interior] = np.linalg.solve(interior_mx, orig_rhs[interior])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solutions should be equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(orig_lhs, manual_lhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great. Now let's compute the reduced basis representation of this solution. (It's the reduced solution that we expect to get.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "expected_red_lhs = np.linalg.solve(Vred, orig_lhs[interior])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's compute the *actual* reduced solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "red_rhs = V.T.dot(orig_rhs)\n",
    "actual_red_lhs = np.linalg.solve(red_mx, red_rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are they equal, as hoped?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(expected_red_lhs, actual_red_lhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good. So when moved back into the original basis the reduced and original solution should be identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(orig_lhs, V.do)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}