The coupled equation in the time domain.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": [],
      "authorship_tag": "ABX9TyN0cY6Dqpal4FkyyfXXQdFk",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/GEORMC/Nnumerical_Methods_Course/blob/main/The%20coupled%20equation%20in%20the%20time%20domain.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": "boju8AYKJmLI"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "\n",
        "# Define parameters\n",
        "theta = 0.5  # Time-stepping parameter\n",
        "delta_t = 0.01  # Time step size\n",
        "num_steps = 100  # Number of time steps\n",
        "\n",
        "# Example matrices (replace with your specific problem matrices)\n",
        "n_dof = 3  # Degrees of freedom (example size)\n",
        "K_e = np.array([[10, -2, 0], [-2, 10, -2], [0, -2, 10]])  # Example stiffness matrix\n",
        "Q = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])  # Coupling matrix\n",
        "S = np.eye(n_dof)  # Storage matrix\n",
        "H = 0.1 * np.eye(n_dof)  # Hydraulic conductivity matrix\n",
        "\n",
        "# Initial conditions\n",
        "u_n = np.zeros(n_dof)  # Initial displacement\n",
        "p_n = np.zeros(n_dof)  # Initial pore pressure\n",
        "f_u = np.zeros(n_dof)  # External force\n",
        "f_p = np.zeros(n_dof)  # Fluid source\n",
        "\n",
        "# Time-stepping loop\n",
        "u_history = [u_n]  # Displacement history\n",
        "p_history = [p_n]  # Pressure history\n",
        "\n",
        "# Assemble constant matrices\n",
        "LHS = np.block([\n",
        "    [theta * K_e, -theta * Q],\n",
        "    [Q.T, S + delta_t * theta * H]\n",
        "])\n",
        "\n",
        "for step in range(num_steps):\n",
        "    # Calculate RHS\n",
        "    RHS_u = (1 - theta) * K_e @ u_n + f_u\n",
        "    RHS_p = S @ p_n - (1 - theta) * delta_t * H @ p_n + delta_t * f_p\n",
        "    RHS = np.hstack([RHS_u, RHS_p])\n",
        "\n",
        "    # Solve for next step\n",
        "    solution = np.linalg.solve(LHS, RHS)\n",
        "    u_n1, p_n1 = solution[:n_dof], solution[n_dof:]\n",
        "\n",
        "    # Update history\n",
        "    u_history.append(u_n1)\n",
        "    p_history.append(p_n1)\n",
        "\n",
        "    # Update for next iteration\n",
        "    u_n, p_n = u_n1, p_n1\n",
        "\n",
        "# Post-process results\n",
        "u_history = np.array(u_history)\n",
        "p_history = np.array(p_history)\n",
        "\n",
        "# Example: plot results\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.figure(figsize=(12, 6))\n",
        "plt.plot(u_history, label=\"Displacement (u)\")\n",
        "plt.plot(p_history, label=\"Pore Pressure (p)\", linestyle=\"--\")\n",
        "plt.xlabel(\"Time Step\")\n",
        "plt.ylabel(\"Value\")\n",
        "plt.legend()\n",
        "plt.title(\"Coupled System Solution\")\n",
        "plt.show()\n"
      ]
    }
  ]
}