08_coupled_problems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 8 - Coupled Problems, monolithic and staggered schemes\n",
    "\n",
    "Consider two coupled ODEs:\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{y} &= k_1 (y_\\text{eq}-y)(z_\\text{eq}-c_1z)\n",
    "    \\\\\n",
    "    \\dot{z} &= k_2 (z_\\text{eq}-z)^3 + c_2\\dot{y}\n",
    "\\end{align}\n",
    "\n",
    "First, we import some libraries for numerical basics and for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Some plot settings\n",
    "plt.style.use('seaborn-deep')\n",
    "plt.rcParams['lines.linewidth']= 2.0\n",
    "plt.rcParams['lines.color']= 'black'\n",
    "plt.rcParams['legend.frameon']=True\n",
    "plt.rcParams['font.family'] = 'serif'\n",
    "plt.rcParams['legend.fontsize']=14\n",
    "plt.rcParams['font.size'] = 14\n",
    "plt.rcParams['axes.spines.right'] = False\n",
    "plt.rcParams['axes.spines.top'] = False\n",
    "plt.rcParams['axes.spines.left'] = True\n",
    "plt.rcParams['axes.spines.bottom'] = True\n",
    "plt.rcParams['axes.axisbelow'] = True\n",
    "plt.rcParams['figure.figsize'] = (8, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, let's consider a Backward Euler scheme, so that\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} &= k_1 (y_\\text{eq} - y_{n+1})(z_\\text{eq} - c_1 z_{n+1})\n",
    "    \\\\\n",
    "    \\frac{z_{n+1}- z_n}{t_{n+1} - t_n} &= k_2 (z_\\text{eq} - z_{n+1})^3 + c_2\\frac{y_{n+1}- y_n}{t_{n+1} - t_n}\n",
    "\\end{align}\n",
    "\n",
    "This can be re-cast in a residual form, i.e based on an approximation error:\n",
    "\n",
    "\\begin{align}\n",
    "    r_y &= \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} - k_1 (y_\\text{eq} - y_{n+1})(z_\\text{eq} - c_1 z_{n+1})\n",
    "    \\\\\n",
    "    r_z &= \\frac{z_{n+1}- z_n}{t_{n+1} - t_n} - k_2 (z_\\text{eq} - z_{n+1})^3 - c_2 \\frac{y_{n+1}- y_n}{t_{n+1} - t_n}\n",
    "\\end{align}\n",
    "\n",
    "where we would like to find the solutions $y_{n+1}$ and $z_{n+1}$ such that $r_y(y_{n+1},z_{n+1}) = 0$ and $r_z(y_{n+1},z_{n+1}) = 0$ *simultaneously*. To achieve that, we have several options.\n",
    "\n",
    "## Monolithic scheme\n",
    "\n",
    "The monolithic scheme iteratively solves both equations simultaneously, i.e. in a monolithic entity. Choosing a Newton-Raphson approach to resolve the non-linearities results in \n",
    "\n",
    "\\begin{align}\n",
    "    r_y(y_{n+1},z_{n+1}) &\\approx r_y(y_n,z_n) + \\left. \\frac{\\partial r_y}{\\partial y} \\right|_n \\underbrace{(y_{n+1} - y_n)}_{\\displaystyle \\Delta y_{n+1}} + \\left. \\frac{\\partial r_y}{\\partial z} \\right|_n \\underbrace{(z_{n+1} - z_n)}_{\\displaystyle \\Delta z_{n+1}} \\overset{!}{=} 0\n",
    "    \\\\\n",
    "    r_z(y_{n+1},z_{n+1}) &\\approx r_z(y_n,z_n) + \\left. \\frac{\\partial r_z}{\\partial y} \\right|_n \\underbrace{(y_{n+1} - y_n)}_{\\displaystyle \\Delta y_{n+1}} + \\left. \\frac{\\partial r_z}{\\partial z} \\right|_n \\underbrace{(z_{n+1} - z_n)}_{\\displaystyle \\Delta z_{n+1}} \\overset{!}{=} 0\n",
    "\\end{align}\n",
    "\n",
    "To make the monolithic structure more clearly visible, let's re-arrange this to $-r_i = J_{ij} \\Delta u_j$ form and cast it into a vector-matrix format (making the time step index $n+1$ implicit and introducing the iteration counter $k$):\n",
    "\n",
    "$$\n",
    "-\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "    r_y(y_k,z_k)\n",
    "    \\\\\n",
    "    r_z(y_k,z_k)\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left.\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "    \\displaystyle \\frac{\\partial r_y}{\\partial y} & \\displaystyle \\frac{\\partial r_y}{\\partial z}\n",
    "    \\\\\n",
    "    \\displaystyle \\frac{\\partial r_z}{\\partial y} & \\displaystyle \\frac{\\partial r_z}{\\partial z}\n",
    "\\end{array}\n",
    "\\right)\n",
    " \\right|_k\n",
    " \\left(\n",
    "\\begin{array}{c}\n",
    "    \\Delta y_{k+1}\n",
    "    \\\\\n",
    "    \\Delta z_{k+1}\n",
    "\\end{array}\n",
    "\\right)\n",
    "= \n",
    "\\left.\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "    J_{yy} & J_{yz}\n",
    "    \\\\\n",
    "    J_{zy} & J_{zz}\n",
    "\\end{array}\n",
    "\\right)\n",
    " \\right|_k\n",
    " \\left(\n",
    "\\begin{array}{c}\n",
    "    \\Delta y_{k+1}\n",
    "    \\\\\n",
    "    \\Delta z_{k+1}\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_1 = 2.\n",
    "y_eq = 1.\n",
    "k_2 = 2\n",
    "z_eq = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residual_y(y_i,y_n,z_i,dt,c_1):\n",
    "    return (y_i - y_n)/dt - k_1 * (y_eq - y_i)*(z_eq - c_1 * z_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residual_z(y_i,y_n,z_i,z_n,dt,c_2):\n",
    "    return (z_i - z_n)/dt - k_2 * (z_eq - z_i)**3 - c_2 * (y_i - y_n)/dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_yy(z_i,dt, c_1):\n",
    "    return 1./dt + k_1*(z_eq - c_1 * z_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_yz(y_i, c_1):\n",
    "    return k_1 * (y_eq - y_i) * c_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_zz(z_i,dt):\n",
    "    return 1./dt + 3*k_2*(z_eq - z_i)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_zy(dt, c_2):\n",
    "    return -c_2/dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for every time step we advance the system by applying the generalized midpoint Newton-Raphson update rule. Each new time step value then serves as the initial condition for the following linearization step.\n",
    "\n",
    "We will *recursively* apply the Newton update until the residual is below a specified tolerance:\n",
    "\n",
    "$\\left| r_{n+1}^{i+1}\\right| < \\epsilon_\\text{abs}$\n",
    "\n",
    "with a suitably chosen absolute tolerance set here to $10^{-6}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_monolithic(dt,c1,c2):\n",
    "    #Startwerte\n",
    "    t_end = 5.\n",
    "    absolute_tolerance = 1.e-6\n",
    "    max_iter = 1000\n",
    "    iteration_counter = np.array([0])\n",
    "    u = [np.array([0., 0.])]#initial values for y and z\n",
    "    times = np.array([0.])\n",
    "    #\n",
    "    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point\n",
    "        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt\n",
    "        u_old = u[-1] #Starting value for recursive update\n",
    "        i = 0\n",
    "        #\n",
    "        while True:\n",
    "            #evaluate residual\n",
    "            res = np.array([residual_y(u_old[0],u[-1][0],u_old[1],dt,c1),\n",
    "                            residual_z(u_old[0],u[-1][0],u_old[1],u[-1][1],dt,c2)])\n",
    "            #if residual is below tolerance, above maximum iterations, stop iterations\n",
    "            if (np.linalg.norm(res) < absolute_tolerance or i > max_iter): \n",
    "                break\n",
    "            #evaluate Jacobian\n",
    "            Jac_yy = J_yy(u_old[1],dt,c1)\n",
    "            Jac_yz = J_yz(u_old[0],c1)\n",
    "            Jac_zy = J_zy(dt,c2)\n",
    "            Jac_zz = J_zz(u_old[1],dt)\n",
    "            Jac = np.array([[Jac_yy,Jac_yz],[Jac_zy,Jac_zz]])\n",
    "            #perform linear step\n",
    "            u_new = u_old + np.linalg.solve(Jac,-res)\n",
    "            #update counter\n",
    "            i += 1\n",
    "            u_old = u_new #preparation of next recursion\n",
    "\n",
    "        u.append(u_new) #append the new found solution to the solution vector\n",
    "        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge\n",
    "        \n",
    "    return times, u, iteration_counter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a29e9894a8804fedbf8d9e0303362b5f"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(dt=widgets.BoundedFloatText(value=0.2,min=1e-3,max=0.8,description='$\\Delta t$ / s'))\n",
    "\n",
    "def plot_monolithic(dt=0.01):\n",
    "    time_n, u_n, iters11 = integrate_monolithic(dt,1.,1.)\n",
    "    u11 = np.array(u_n)\n",
    "    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))\n",
    "    #solution\n",
    "    ax[0][0].plot(time_n,u11.T[0],marker='d',label='11')\n",
    "    ax[0][0].legend()\n",
    "    ax[0][0].set_xlabel('$t$ / s')\n",
    "    ax[0][0].set_ylabel('$y$')\n",
    "    #error\n",
    "    ax[0][1].plot(time_n,u11.T[1],marker='d',label='11')\n",
    "    ax[0][1].legend()\n",
    "    ax[0][1].set_xlabel('$t$ / s')\n",
    "    ax[0][1].set_ylabel('$z$')\n",
    "    #iterations\n",
    "    ax[1][0].plot(iters11,marker='d',label='11')\n",
    "    ax[1][0].set_xlabel('time step')\n",
    "    ax[1][0].set_ylabel('number of iterations')\n",
    "    ax[1][0].legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tasks\n",
    "* How does the number of iterations change if the system if you a) uncouple the system or b) choose a uni-lateral coupling?\n",
    "* How does the system evolve for other parameter settings? Try to understand why. Examples: $(c_1,c_2) = (0,2),(1,-5)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Staggered scheme\n",
    "\n",
    "A staggered scheme works by solving one equation for one unknown while keeping the other unknown(s) fixed. With the new solution, one proceeds to the other equation(s) to solve for the next unknown, while keeping the other unknown(s) fixed. In this way, one iterates between the equations until the solution is stationary.\n",
    "\n",
    "Choosing a Newton-Raphson approach to resolve the non-linearities results in \n",
    "\n",
    "\\begin{align}\n",
    "    r_y(y_{n+1},z_n) &\\approx r_y(y_n,z_n) + \\left. \\frac{\\partial r_y}{\\partial y} \\right|_n \\underbrace{(y_{n+1} - y_n)}_{\\displaystyle \\Delta y_{n+1}} \\overset{!}{=} 0\n",
    "    \\\\\n",
    "    r_z(y_n,z_{n+1}) &\\approx r_z(y_n,z_n) + \\left. \\frac{\\partial r_z}{\\partial z} \\right|_n \\underbrace{(z_{n+1} - z_n)}_{\\displaystyle \\Delta z_{n+1}} \\overset{!}{=} 0\n",
    "\\end{align}\n",
    "\n",
    "To make the staggered scheme more clearly visible, let's re-arrange this to $-r_i = J_{ij} \\Delta u_j$ form, cast it into a vector-matrix format, and introduce the Newton iteration counters $k$ and $l$ as well as the coupling iteration counter $o$:\n",
    "\n",
    "\\begin{align}\n",
    "            - r_y({}^oy^k_{n+1},{}^o z_{n+1}) &= {}^o J_{yy}|_{n+1}^k \\Delta {}^oy_{n+1}^{k+1}\n",
    "            \\\\\n",
    "            - r_z({}^oy_{n+1},{}^o z_{n+1}^k) &= {}^o J_{zz}|_{n+1}^k \\Delta {}^oz_{n+1}^{k+1}\n",
    "\\end{align}\n",
    "\n",
    "The coupling iterations $o$ are the outer loop, while the inner loop consists of the Newton iterations. The overall solution is found when the outer solutions of the individual unknowns $z$ and $y$ no longer change between coupling iterations:\n",
    "\n",
    "$$\n",
    "    |{}^oy^k_{n+1} - {}^{o-1}y^k_{n+1}| + |{}^oz^k_{n+1} - {}^{o-1}z^k_{n+1}| < \\epsilon_\\text{abs}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_staggered(dt,c1,c2):\n",
    "    #Startwerte\n",
    "    t_end = 5.\n",
    "    absolute_tolerance = 1.e-6\n",
    "    max_iter = 100\n",
    "    iteration_counter = np.array([0])\n",
    "    y = np.array([0.])#initial values for y\n",
    "    z = np.array([0.])#initial values for z\n",
    "    times = np.array([0.])\n",
    "    #\n",
    "    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point\n",
    "        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt\n",
    "        y_old = y[-1] #Starting value for recursive update\n",
    "        z_old = z[-1] #Starting value for recursive update\n",
    "        y_coup = y_old\n",
    "        z_coup = z_old\n",
    "        i = 0\n",
    "        #\n",
    "        while True:\n",
    "            while True:\n",
    "                #evaluate residual\n",
    "                res = residual_y(y_old,y[-1],z_old,dt,c1)\n",
    "                #if residual is below tolerance, above maximum iterations, stop iterations\n",
    "                if (np.abs(res) < absolute_tolerance or i > max_iter): \n",
    "                    break\n",
    "                #perform linear step\n",
    "                y_new = y_old - res/J_yy(z_old,dt,c1)\n",
    "                #update counter\n",
    "                i += 1\n",
    "                y_old = y_new #preparation of next recursion\n",
    "\n",
    "            while True:\n",
    "                #evaluate residual\n",
    "                res = residual_z(y_old,y[-1],z_old,z[-1],dt,c2)\n",
    "                #if residual is below tolerance, above maximum iterations, stop iterations\n",
    "                if (np.abs(res) < absolute_tolerance or i > max_iter): \n",
    "                    break\n",
    "                #perform linear step\n",
    "                z_new = z_old - res/J_zz(z_old,dt)\n",
    "                #update counter\n",
    "                i += 1\n",
    "                z_old = z_new #preparation of next recursion\n",
    "\n",
    "            if ((np.abs(y_coup - y_new) + np.abs(z_coup - z_new)) < absolute_tolerance or i > max_iter):\n",
    "                break\n",
    "            y_coup = y_new\n",
    "            z_coup = z_new\n",
    "\n",
    "        y = np.append(y,y_new) #append the new found solution to the solution vector\n",
    "        z = np.append(z,z_new) #append the new found solution to the solution vector\n",
    "        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge\n",
    "        \n",
    "    return times, y, z, iteration_counter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "96e36ba86a494c09aaf9a46f02ad197a"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(dt=widgets.BoundedFloatText(value=0.2,min=1e-3,max=0.8,description='$\\Delta t$ / s'))\n",
    "\n",
    "def plot_staggered(dt=0.01):\n",
    "    time_n, y, z, iters11 = integrate_staggered(dt,1.,1.)\n",
    "    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))\n",
    "    #solution\n",
    "    ax[0][0].plot(time_n,y,marker='d',label='11')\n",
    "    ax[0][0].legend()\n",
    "    ax[0][0].set_xlabel('$t$ / s')\n",
    "    ax[0][0].set_ylabel('$y$')\n",
    "    #error\n",
    "    ax[0][1].plot(time_n,z,marker='d',label='11')\n",
    "    ax[0][1].legend()\n",
    "    ax[0][1].set_xlabel('$t$ / s')\n",
    "    ax[0][1].set_ylabel('$z$')\n",
    "    #iterations\n",
    "    ax[1][0].plot(iters11,marker='d',label='11')\n",
    "    ax[1][0].set_xlabel('time step')\n",
    "    ax[1][0].set_ylabel('number of iterations')\n",
    "    ax[1][0].legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1efc83c7f58c4a5e87074cea569e408b"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(dt=widgets.BoundedFloatText(value=0.2,min=1e-3,max=0.8,description='$\\Delta t$ / s'))\n",
    "\n",
    "def plot_comparison(dt=0.01):\n",
    "    \n",
    "    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))\n",
    "    #set values for coupling parameters\n",
    "    c1 = [0.,1.]\n",
    "    c2 = [0.,1.,-5.]\n",
    "    # run all combinations\n",
    "    for i in range(len(c1)):\n",
    "        for j in range(len(c2)):\n",
    "            #run both schemes\n",
    "            time_n, u_n, iters_m = integrate_monolithic(dt,c1[i],c2[j])\n",
    "            u_m = np.array(u_n)\n",
    "            time_n, y_s, z_s, iters_s = integrate_staggered(dt,c1[i],c2[j])\n",
    "            #plot color\n",
    "            color = next(ax[0][0]._get_lines.prop_cycler)['color']\n",
    "            #solution y\n",
    "            ax[0][0].plot(time_n,u_m.T[0],ls='--',color=color)\n",
    "            ax[0][0].plot(time_n,y_s,color=color)\n",
    "            #\n",
    "\n",
    "            #solution z\n",
    "            ax[0][1].plot(time_n,u_m.T[1],ls='--',color=color)\n",
    "            ax[0][1].plot(time_n,z_s,color=color)\n",
    "\n",
    "            #iterations\n",
    "            ax[1][0].plot(iters_m,ls='--',color=color)\n",
    "            ax[1][0].plot(iters_s,color=color)\n",
    "            \n",
    "            #dummy\n",
    "            ax[1][1].plot(0,0,color=color,label=\"$c_1 = %.1f$, $c_2 = %.1f$\" %(c1[i],c2[j]))\n",
    "    \n",
    "    ax[0][0].set_xlabel('$t$ / s')\n",
    "    ax[0][0].set_ylabel('$y$')\n",
    "    ax[0][1].set_xlabel('$t$ / s')\n",
    "    ax[0][1].set_ylabel('$z$')\n",
    "    ax[1][0].set_xlabel('time step')\n",
    "    ax[1][0].set_ylabel('number of iterations')\n",
    "    ax[1][1].legend()\n",
    "    ax[1][1].set_title('monolithic --; staggered __')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tasks\n",
    "\n",
    "* What main difference do you observe between both schemes?\n",
    "* Run the combination $c_1=-3$, $c_2 = -5$. What do you observe? What happens if you a) increase and b) decrease the time steps?"
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}