nodal_flux_tests_1D.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false
      },
      "source": [
        "<div style=\"display: flex; align-items: center;\">\n",
        "    <img src=\"https://github.com/nagelt/Teaching_Scripts/raw/9d9e29ecca4b04eaf7397938eacbf116d37ddc93/Images/TUBAF_Logo_blau.png\" width=\"500\" height=\"auto\" height=\"auto\" style=\"margin-right: 100px;\" />\n",
        "    <div>\n",
        "        <p><strong>Prof. Dr. Thomas Nagel</strong></p>\n",
        "        <p>Chair of Soil Mechanics and Foundation Engineering<br>Geotechnical Institute<br>Technische Universität Bergakademie Freiberg.</p>\n",
        "        <p><a href=\"https://tu-freiberg.de/en/soilmechanics\">https://tu-freiberg.de/en/soilmechanics</a></p>\n",
        "    </div>\n",
        "</div>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import scipy as sp\n",
        "from pypardiso import spsolve\n",
        "\n",
        "#Some plot settings\n",
        "%run plot_functions/plot_settings.py"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Quick tests on nodal fluxes\n",
        "\n",
        "## Divergence-free flow\n",
        "\n",
        "Consider a linear concentration gradient of $\\nabla c = -1$ on $x \\in [0,1]$ for $v = 1$ and $D = 1$. We calculate nodal contributions based on the following terms (no units for the moment):\n",
        "\n",
        "\\begin{align}\n",
        "    f_1 &= \\int \\limits_0^1 \\nabla N v c \\text{d}x\n",
        "    \\\n",
        "    f_2 &= \\int \\limits_0^1 N v \\nabla c \\text{d}x\n",
        "    \\\n",
        "    f_3 &= -\\int \\limits_0^1 \\nabla N v \\nabla c \\text{d}x\n",
        "\\end{align}\n",
        "\n",
        "$f_1$ is a treatment used for advection in total flux models, $f_2$ is motivated by the implementation of advection with substituted mass balance, $f_3$ for diffusion."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Finite elements in 1D\n",
        "\n",
        "We first create an element class. An element knows the number of nodes it has, their IDs in the global node vector, and the coordinates of its nodes. Linear elements have 2 nodes and 2 quadrature points, quadratic elements 3 nodes and 3 quadrature points. The natural coordinates of the element run from -1 to 1, and the quadrature points and weights are directly taken from Numpy."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [],
      "source": [
        "#element class\n",
        "class line_element():#local coordinates go from -1 to 1\n",
        "    #takes number of nodes, global nodal coordinates, global node ids\n",
        "    def __init__(self, nnodes=2, ncoords=[0.,1.], nids=[0,1]):\n",
        "        self.__nnodes = nnodes\n",
        "        if (len(ncoords) != self.__nnodes):\n",
        "            raise Exception(\"Number of coordinates does not match number \\\n",
        "                            of nodes of element (%i vs of %i)\" %(self.__nnodes,len(ncoords)))\n",
        "        else:\n",
        "            self.__coords = np.array(ncoords)\n",
        "        \n",
        "        self.__natural_coords = (self.__coords-self.__coords[0])/(self.__coords[-1]-self.__coords[0])*2. - 1.\n",
        "        \n",
        "        if (len(nids) != self.__nnodes):\n",
        "            raise Exception(\"Number of node IDs does not match number \\\n",
        "                            of nodes of element (%i vs of %i)\" %(self.__nnodes,len(nids)))\n",
        "        else:\n",
        "            self.__global_ids = np.array(nids)\n",
        "        self.__quad_degree = self.__nnodes\n",
        "        self.__quad_points, self.__quad_weights = np.polynomial.legendre.leggauss(self.__quad_degree)\n",
        "                "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [],
      "source": [
        "#N\n",
        "def shape_function(element_order,xi):\n",
        "    if (element_order == 1):\n",
        "            return np.array([1])\n",
        "    elif (element_order == 2): #-1,1\n",
        "            return np.array([(1.-xi)/2., (1.+xi)/2.])\n",
        "    elif (element_order == 3): #-1, 0, 1\n",
        "            return np.array([(xi - 1.)*xi/2., (1-xi)*(1+xi), (1+xi)*xi/2.])\n",
        "    elif (element_order == 4): #-1, -1/3, 1/3, 1\n",
        "            return np.array([9/16*(1-xi)*(xi**2 - 1/9),\n",
        "                            27/16*(xi**2-1)*(xi-1/3),\n",
        "                            27/16*(1-xi**2)*(xi+1/3),\n",
        "                            9/16*(xi+1)*(xi**2-1/9)])\n",
        "        \n",
        "#dN_dxi\n",
        "def dshape_function_dxi(element_order,xi):\n",
        "    if (element_order == 1):\n",
        "            return np.array([0])\n",
        "    elif (element_order == 2): #-1,1\n",
        "        return np.array([-0.5, 0.5])  #xi only later for plotting dimensions\n",
        "    elif (element_order == 3):#-1,0,1\n",
        "        return np.array([xi - 0.5,-2.*xi,xi + 0.5])\n",
        "    elif (element_order == 4): #-1, -1/3, 1/3, 1\n",
        "            return np.array([-27*xi**2/16 + 9*xi/8 + 1/16,\n",
        "                            81/16*xi**2 - 9/8 * xi - 27/16, \n",
        "                            -81/16*xi**2 - 9/8 * xi +27/16,\n",
        "                            27*xi**2/16 + 9*xi/8 - 1/16])\n",
        "\n",
        "#dz_dxi\n",
        "def element_jacobian(element,xi):\n",
        "    element_order = element._line_element__nnodes\n",
        "    Jacobian = 0.\n",
        "    Jacobian += dshape_function_dxi(element_order,xi).dot(element._line_element__coords)\n",
        "    return Jacobian\n",
        "\n",
        "#dN_dz\n",
        "def grad_shape_function(element,xi):\n",
        "    element_order = element._line_element__nnodes\n",
        "    Jac = element_jacobian(element,xi)\n",
        "    return dshape_function_dxi(element_order,xi)/Jac"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Local assember"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 99,
      "metadata": {},
      "outputs": [],
      "source": [
        "def test_assembler(elem, c_nodes, v_nodes):\n",
        "    element_order = elem._line_element__nnodes\n",
        "    b_1 = np.zeros(element_order)\n",
        "    b_2 = np.zeros(element_order)\n",
        "    b_3 = np.zeros(element_order)\n",
        "    b_4 = np.zeros(element_order)\n",
        "    #z_nodes = elem._line_element__coords\n",
        "    for i in range(elem._line_element__quad_degree):\n",
        "        #local integration point coordinate\n",
        "        xi = elem._line_element__quad_points[i]\n",
        "        #shape function\n",
        "        N = shape_function(element_order,xi)\n",
        "        #gradient of shape function\n",
        "        dN_dX = grad_shape_function(elem,xi)\n",
        "        #determinant of Jacobian\n",
        "        detJ = np.abs(element_jacobian(elem,xi))\n",
        "        #integration weight\n",
        "        w = elem._line_element__quad_weights[i]\n",
        "        \n",
        "        #global integration point coordinate (for spatially varying properties)\n",
        "        c_glob = np.dot(N,c_nodes)\n",
        "        v_glob = np.dot(N,v_nodes)\n",
        "        grad_c = np.dot(dN_dX,c_nodes)\n",
        "        grad_v = np.dot(dN_dX,v_nodes)\n",
        "        #evaluation of local material/structural properties\n",
        "        #assembly of local RHS\n",
        "        b_1 += dN_dX * c_glob * v_glob * w * detJ\n",
        "        b_2 += N * grad_c * v_glob * w * detJ\n",
        "        b_3 -= dN_dX * grad_c * w * detJ\n",
        "        b_4 += N * c_glob * grad_v * w * detJ\n",
        "    return b_1, b_2, b_3, b_4"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 103,
      "metadata": {
        "tags": []
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "element order:  1\n",
            "f_1:  [-0.5  0.5]\n",
            "f_2:  [-0.5 -0.5]\n",
            "f_1 + f_2:  [-1.  0.]\n",
            "f_3:  [-1.  1.]\n",
            "---\n",
            "\n",
            "element order:  2\n",
            "f_1:  [-0.83333333  0.66666667  0.16666667]\n",
            "f_2:  [-0.16666667 -0.66666667 -0.16666667]\n",
            "f_1 + f_2:  [-1.00000000e+00  2.22044605e-16 -1.11022302e-16]\n",
            "f_3:  [-1.  0.  1.]\n",
            "---\n",
            "\n",
            "element order:  3\n",
            "f_1:  [-0.875  0.375  0.375  0.125]\n",
            "f_2:  [-0.125 -0.375 -0.375 -0.125]\n",
            "f_1 + f_2:  [-1.00000000e+00 -3.88578059e-16  4.44089210e-16 -1.38777878e-17]\n",
            "f_3:  [-1.00000000e+00 -7.77156117e-16  8.88178420e-16  1.00000000e+00]\n",
            "---\n",
            "\n"
          ]
        }
      ],
      "source": [
        "#generate elements of order 1 to 5 and calculate fluxes\n",
        "for i in range(1,4):\n",
        "    nodes = np.linspace(0,i,i+1)\n",
        "    ncoords = np.linspace(0,1,i+1)\n",
        "    elem = line_element(i+1,ncoords,nodes)\n",
        "    nconc = np.flip(ncoords)\n",
        "    f1, f2, f3, f4 = test_assembler(elem,nconc,np.ones(len(nodes)))\n",
        "    print(\"element order: \", i)\n",
        "    print(\"f_1: \", f1)\n",
        "    print(\"f_2: \", f2)\n",
        "    print(\"f_1 + f_2: \", f1+f2)\n",
        "    print(\"f_3: \", f3)\n",
        "    print(\"---\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Diffusion is always correctly quantified by $f_3$: 1 influx, 1 outflux, no internal fluxes in the domain.\n",
        "\n",
        "Advection is approached (1 influx and 0 outflux) by $f_1$ as the element order increases, but some averaging remains also at higher order. Internal fluxes show the redistribution of mass from left to right.\n",
        "\n",
        "The sum $f_1 + f_2$ recovers the advective flux with expected values on the boundaries, without internal redistribution information."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Divergent flux\n",
        "\n",
        "Now we add a linear velocity profile with mean value of 1 and repeat the test."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 109,
      "metadata": {
        "tags": []
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1.25, 0.75])"
            ]
          },
          "execution_count": 109,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "ncoords = np.linspace(0,1,2)\n",
        "v = np.ones(len(ncoords))-(ncoords-0.5)*0.5\n",
        "v"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 105,
      "metadata": {
        "tags": []
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "element order:  1\n",
            "f_1:  [-0.54166667  0.54166667]\n",
            "f_2:  [-0.54166667 -0.45833333]\n",
            "f_1 + f_2:  [-1.08333333  0.08333333]\n",
            "f_1 + f_2 + f_4:  [-1.25000000e+00  1.38777878e-16]\n",
            "f_3:  [-1.  1.]\n",
            "---\n",
            "\n",
            "element order:  2\n",
            "f_1:  [-0.95833333  0.83333333  0.125     ]\n",
            "f_2:  [-0.20833333 -0.66666667 -0.125     ]\n",
            "f_1 + f_2:  [-1.16666667e+00  1.66666667e-01 -1.24900090e-16]\n",
            "f_1 + f_2 + f_4:  [-1.25000000e+00  2.77555756e-16 -1.24900090e-16]\n",
            "f_3:  [-1.  0.  1.]\n",
            "---\n",
            "\n",
            "element order:  3\n",
            "f_1:  [-1.04791667  0.58125     0.35625     0.11041667]\n",
            "f_2:  [-0.14791667 -0.43125    -0.31875    -0.10208333]\n",
            "f_1 + f_2:  [-1.19583333  0.15        0.0375      0.00833333]\n",
            "f_1 + f_2 + f_4:  [-1.25000000e+00 -1.94289029e-16  2.22044605e-16 -3.98986399e-17]\n",
            "f_3:  [-1.00000000e+00 -7.77156117e-16  8.88178420e-16  1.00000000e+00]\n",
            "---\n",
            "\n"
          ]
        }
      ],
      "source": [
        "#generate elements of order 1 to 5 and calculate fluxes\n",
        "for i in range(1,4):\n",
        "    nodes = np.linspace(0,i,i+1)\n",
        "    ncoords = np.linspace(0,1,i+1)\n",
        "    elem = line_element(i+1,ncoords,nodes)\n",
        "    nconc = np.flip(ncoords)\n",
        "    f1, f2, f3, f4 = test_assembler(elem,nconc,np.ones(len(nodes))-(ncoords-0.5)*0.5)\n",
        "    print(\"element order: \", i)\n",
        "    print(\"f_1: \", f1)\n",
        "    print(\"f_2: \", f2)\n",
        "    print(\"f_1 + f_2: \", f1+f2)\n",
        "    print(\"f_1 + f_2 + f_4: \", f1+f2+f4)\n",
        "    print(\"f_3: \", f3)\n",
        "    print(\"---\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Der outflux von 0 wird in Ordnung 2 erreicht, in 3 aber wieder verfehlt (dennoch eine Ordnung besser als im linearen Element. Im linearen Element hätten wir allerdings konstante Geschwindigkeit, wenn wir mit einem hydraulischen Potenzial an den Knoten arbeiten.). Der Influx von 1.25 wird zunehmend besser approximiert, aber Unterschiede bleiben. Wenn man den Term verfollständigt um\n",
        "\n",
        "$$\n",
        "    f_4 = \\int \\limits_{0}^1 N c \\nabla v \\text{d}x\n",
        "$$\n",
        "\n",
        "geht die Bilanz wieder auf. Im Standard FE allerdings nicht so easy."
      ]
    }
  ],
  "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.9.7"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}