CEE314Problemsets4Q4Teaching.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/AlexCHEN-Engineer/CEE314HW/blob/main/CEE314Problemsets4Q4Teaching.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mediterranean-century"
      },
      "outputs": [],
      "source": [
        "%time\n",
        "# Do not revise here. \n",
        "try:\n",
        "    import google.colab\n",
        "except ImportError:\n",
        "    import ufl\n",
        "    import dolfin\n",
        "    import mshr\n",
        "else:\n",
        "    try:\n",
        "        import ufl\n",
        "        import dolfin\n",
        "        import mshr\n",
        "    except ImportError:\n",
        "        !wget \"https://fem-on-colab.github.io/releases/fenics-install.sh\" -O \"/tmp/fenics-install.sh\" && bash \"/tmp/fenics-install.sh\"\n",
        "        import ufl\n",
        "        import dolfin\n",
        "        import mshr"
      ],
      "id": "mediterranean-century"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Variational formulation\n",
        "\n",
        "In this part, we will solve a coupled poromechanics problems.\n",
        "\n",
        "For a domain denoted as $\\Omega$, the boundary is $\\Gamma = \\Gamma_u\\cup \\Gamma_t = \\Gamma_p\\cup\\Gamma_q$, $\\Gamma_u\\cap \\Gamma_t=\\phi$ and $\\Gamma_p\\cup\\Gamma_q=\\phi$. \n",
        "\n",
        "The strong form can be expressed as:\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "&\\nabla(\\boldsymbol{\\sigma}^\\prime-p\\boldsymbol{1})+\\rho\\boldsymbol{g}=\\boldsymbol{0},   \\\\\n",
        "&\\nabla⋅\\dot{\\boldsymbol{u}}+\\nabla⋅{\\boldsymbol{q}}=0, \\\\\n",
        "& \\boldsymbol{u} = \\overline{\\boldsymbol{u}}\\quad on \\quad \\Gamma_u, \\\\\n",
        "& \\boldsymbol{\\sigma}^\\prime\\cdot\\boldsymbol{n} = \\overline{\\boldsymbol{t}}\\quad on\\quad \\Gamma_t, \\\\\n",
        "& p = \\overline{p}\\quad on \\quad \\Gamma_p, \\\\\n",
        "& \\boldsymbol{q}\\cdot \\boldsymbol{n} = \\overline{q}\\quad on \\quad \\Gamma_q.\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "\n",
        "As for the constitutive equations, we need to refer to\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& \\boldsymbol{\\sigma}^\\prime = \\mathbb{C}^e:\\boldsymbol{\\epsilon} ,\\\\\n",
        "& \\boldsymbol{q} = -\\frac{\\boldsymbol{\\kappa}}{\\mu_f}\\cdot \\nabla p.\n",
        "\\end{aligned}\n",
        "\\end{equation}\n"
      ],
      "metadata": {
        "id": "P2bXKmcEqOJN"
      },
      "id": "P2bXKmcEqOJN"
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Derive the weak form\n",
        "\n",
        "The trial space is defined as \n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& S_u = \\{\\boldsymbol{u}:\\Omega→\\mathbb{R}^{ndim}|\\boldsymbol{u}→\\boldsymbol{H}^1,\\ \\boldsymbol{u}=\\overline{\\boldsymbol{u}}\\quad on \\quad \\Gamma_u \\} ,\\\\\n",
        "&S_p = \\{p:\\Omega→\\mathbb{R}|p→H^1,\\ p=\\overline{p}\\quad on \\quad \\Gamma_p \\}\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "\n",
        "For the weighting space,\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& V_u = \\{\\boldsymbol{\\eta}:\\Omega→\\mathbb{R}^{ndim}|\\boldsymbol{\\eta}→\\boldsymbol{H}^1,\\ \\boldsymbol{\\eta}={\\boldsymbol{0}}\\quad on \\quad \\Gamma_u \\} ,\\\\\n",
        "&V_\\psi = \\{\\psi:\\Omega→\\mathbb{R}|\\psi→H^1,\\ \\psi=0\\quad on \\quad \\Gamma_\\psi \\}\n",
        "\\end{aligned}\n",
        "\\end{equation}\n"
      ],
      "metadata": {
        "id": "MLHKWGvgv_Cg"
      },
      "id": "MLHKWGvgv_Cg"
    },
    {
      "cell_type": "markdown",
      "source": [
        "The weak form becomes:\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& \\mathcal{G}_1 = \\int_{\\Omega}(\\nabla^s\\boldsymbol{\\eta}:\\boldsymbol{\\sigma}^\\prime-p\\nabla\\cdot\\boldsymbol{\\eta}-\\boldsymbol{\\eta}\\cdot{\\rho}\\boldsymbol{g})dV - \\int_{\\Gamma_t}\\boldsymbol{\\eta}\\cdot\\overline{\\boldsymbol{t}}dA = \\boldsymbol{0} \\\\\n",
        "& \\mathcal{G}_2 = \\int_{\\Omega}(\\psi\\nabla\\cdot\\dot{\\boldsymbol{u}}-\\nabla \\psi\\cdot \\boldsymbol{q})dV + \\int_{\\Gamma_q} \\psi\\bar{q}dA\\\\\n",
        "\\end{aligned}\n",
        "\\end{equation}"
      ],
      "metadata": {
        "id": "62rA3t9lxyo6"
      },
      "id": "62rA3t9lxyo6"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yWGBMCZ5jAX2"
      },
      "source": [
        "# Import necessary library"
      ],
      "id": "yWGBMCZ5jAX2"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "infectious-train"
      },
      "outputs": [],
      "source": [
        "from fenics import *\n",
        "from mshr import *\n",
        "from ufl import nabla_div"
      ],
      "id": "infectious-train"
    },
    {
      "cell_type": "code",
      "source": [
        "# Create mesh and define expression\n",
        "# mesh = RectangleMesh.create([Point(0, 0), Point(1,1)],[2,2],CellType.Type.quadrilateral)\n",
        "# We can still create the mesh here.\n",
        "# But, we can't use plot to see the results."
      ],
      "metadata": {
        "id": "xXlk8W3E671m"
      },
      "id": "xXlk8W3E671m",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Parameters\n",
        "\n",
        "In this part, please use the following mechanical and hydraulic parameters. \n",
        "* $\\lambda = 0,\\quad \\mu = 500$\n",
        "* $k = 1.0\\times 10^{-11},\\quad \\mu_l = 1.0\\times 10^{-2}$\n",
        "All the parameters adopt the [international units](https://en.wikipedia.org/wiki/International_unit)."
      ],
      "metadata": {
        "id": "PsmZYVTj5Fd3"
      },
      "id": "PsmZYVTj5Fd3"
    },
    {
      "cell_type": "code",
      "source": [
        "#####################\n",
        "#### TODO ###########\n",
        "#####################\n",
        "########## Start your code here ###########\n",
        "T = None          # final time\n",
        "num_steps = None  # number of time steps\n",
        "dt = None         # time step size\n",
        "rho_l = None      # fluid density\n",
        "E = None          # Young's modulus\n",
        "nu = None         # poisson's ratio\n",
        "k = None   \t      # intrinsic permeability\n",
        "mu_l = None \t  # dynamic viscosity of IF\n",
        "lambda_= None     # Lame coefficient\n",
        "mu = None         # Lame coefficient\n",
        "########### End your code here #############"
      ],
      "metadata": {
        "id": "4Au2rvfqZMTc"
      },
      "id": "4Au2rvfqZMTc",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Create mesh and define expression\n",
        "mesh=RectangleMesh(Point(0.0,0.0), Point(10., 5.), 10, 5,'crossed')\n",
        "# mesh = RectangleMesh.create([Point(0, 0), Point(10, 5)], [20, 10], CellType.Type.quadrilateral)\n",
        "top_u =  CompiledSubDomain(\"near(x[1], side) && x[0] > -DOLFIN_EPS && x[0] < 2.+DOLFIN_EPS\", side = 5.)\n",
        "top_p =  CompiledSubDomain(\"near(x[1], side)\", side = 5.)\n",
        "bottom =  CompiledSubDomain(\"near(x[1], side)\", side = 0.0)\n",
        "left_right = CompiledSubDomain(\"(near(x[0], 0.0) ) || (near(x[0], 10.) )\")\n",
        "\n",
        "boundary = MeshFunction('size_t', mesh, mesh.topology().dim()-1)\n",
        "boundary.set_all(1)\n",
        "top_u.mark(boundary, 2)\n",
        "ds = Measure('ds', domain = mesh, subdomain_data = boundary)\n",
        "\n",
        "T  = Expression(('0.0', '-10'),degree=1)  # Load on the boundary\n",
        "\n",
        "#Define Mixed Space (R2,R) -> (u,p)\n",
        "V = VectorElement(\"CG\", mesh.ufl_cell(), 2)\n",
        "W = FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
        "L = FunctionSpace(mesh,W)\n",
        "MS = dolfin.FunctionSpace(mesh, MixedElement([V,W]))\n",
        "\n",
        "# Define boundary condition\n",
        "bcu1 = DirichletBC(MS.sub(0).sub(0), 0.0, left_right)  # slip condition\n",
        "bcu2 = DirichletBC(MS.sub(0).sub(1), 0.0, bottom)      # slip condition\n",
        "bcp = DirichletBC(MS.sub(1), 0.0, top_p)               # drained condition\n",
        "bc=[bcu1,bcu2,bcp]\n",
        "\n",
        "# Define strain.\n",
        "def epsilon(u):\n",
        "    return 0.5*(nabla_grad(u) + nabla_grad(u).T)\n",
        "# define stress.\n",
        "def sigma(u):\n",
        "   return lambda_*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)\n",
        "\n",
        "# Define variational problem and initial condition\n",
        "X0 = Function(MS)\n",
        "B = TestFunction(MS)\n",
        "\n",
        "e_u0 = Expression(('0.0', '0.0'), degree=1)\n",
        "e_p0 = Expression('0.0', degree=1)\n",
        "\n",
        "u0 = interpolate(e_u0, MS.sub(0).collapse())\n",
        "p0 = interpolate(e_p0, MS.sub(1).collapse())\n",
        "\n",
        "Xn = Function(MS)\n",
        "assign(Xn, [u0, p0])\n",
        "\n",
        "(u,p)=split(X0)\n",
        "(u_n,p_n)=split(Xn)\n",
        "(v,q)=split(B)\n",
        "\n",
        "\n",
        "#####################\n",
        "#### TODO ###########\n",
        "#####################\n",
        "############### Start your code here ###########\n",
        "F = None \n",
        "############### End your code here #############\n",
        "\n",
        "#solver tuning\n",
        "dX0 = TrialFunction(MS)\n",
        "J = derivative(F, X0, dX0)\n",
        "Problem = NonlinearVariationalProblem(F, X0, J = J, bcs = bc)\n",
        "Solver  = NonlinearVariationalSolver(Problem)\n",
        "Solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'\n",
        "Solver.parameters['newton_solver']['relative_tolerance'] = 1.e-11\n",
        "Solver.parameters['newton_solver']['absolute_tolerance'] = 5.e-10\n",
        "\n",
        "#if you want to store the results\n",
        "vtkfile_u = File('MC/u.pvd')\n",
        "vtkfile_p = File('MC/p.pvd')\n",
        "\n",
        "t = 0\n",
        "\n",
        "p_out = []\n",
        "\n",
        "from tqdm import tqdm\n",
        "for n in tqdm(range(1)):\n",
        "\tt += dt\n",
        "\tSolver.solve()    \n",
        "\tassign(Xn,X0)\n",
        "\n",
        "\t(_u,_p)=X0.split()\n",
        "\t\n",
        "\tvtkfile_u << (_u,t)\n",
        "\tvtkfile_p << (_p,t)\n",
        "\tp_out.append(_p(0., 3.))"
      ],
      "metadata": {
        "id": "pw-EBzjqQiha",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6b194bc0-5901-4112-ef7c-6169edbb1751"
      },
      "id": "pw-EBzjqQiha",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "<class 'ufl.indexed.Indexed'>\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 1/1 [00:00<00:00, 29.11it/s]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "u_result, p_result = X0.split()\n",
        "import matplotlib.pyplot as plt\n",
        "c=plot(p_result)\n",
        "plt.colorbar(c)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "G0dWBTs8TB3q"
      },
      "id": "G0dWBTs8TB3q",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Postprocess"
      ],
      "metadata": {
        "id": "bZFw87laa7Tp"
      },
      "id": "bZFw87laa7Tp"
    },
    {
      "cell_type": "code",
      "source": [
        "# Use the following method to zip the file.\n",
        "!zip -r 'MC.zip' '/content/MC' "
      ],
      "metadata": {
        "id": "joHFy9lSTN27"
      },
      "id": "joHFy9lSTN27",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 4(b)"
      ],
      "metadata": {
        "id": "rtxvdCroiuq_"
      },
      "id": "rtxvdCroiuq_"
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "time = np.linspace(0,1000,1000)\n",
        "pressure = np.array(p_out)\n",
        "plt.xlabel(r\"$t$\")\n",
        "plt.ylabel(r\"$p/w$\")\n",
        "plt.semilogx(100/9.8*1e-3*time/4, pressure/10, 'k')"
      ],
      "metadata": {
        "id": "BgAwNRBDfINY"
      },
      "id": "BgAwNRBDfINY",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "Zp1Rs8homNZ9"
      },
      "id": "Zp1Rs8homNZ9",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "machine_shape": "hm",
      "name": "CEE314Problemsets4Q4Teaching.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython"
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python"
    },
    "accelerator": "TPU"
  },
  "nbformat": 4,
  "nbformat_minor": 5
}