sww2d_central_finite_media.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11a20daf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator\n",
    "from matplotlib import ticker"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6731f4c",
   "metadata": {},
   "source": [
    "## 01. Parameters\n",
    "### 01-01. Physical parameters\n",
    "#### - Computational domain: 2D square region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80efbba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# g?\n",
    "g = 9.81\n",
    "\n",
    "# rho0\n",
    "rho0 = 1e-2\n",
    "\n",
    "# Temporal\n",
    "t_sta = 0.\n",
    "t_end = 150.\n",
    "\n",
    "# Spatial\n",
    "x_min = -10.\n",
    "x_max = 10.\n",
    "\n",
    "y_min = copy.copy(x_min) # shallow copy\n",
    "y_max = copy.copy(x_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf9b4aa",
   "metadata": {},
   "source": [
    "### 01-02. Computational parameters\n",
    "#### - Sizes of timestep and grids\n",
    "#### - Sampling frequency "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78be28ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling\n",
    "sample_freq = 10 # sampling data per 10 timesteps\n",
    "\n",
    "# Temporal\n",
    "dt = 5e-2\n",
    "nt = int((t_end - t_sta) / dt) + 1\n",
    "\n",
    "# Spatial\n",
    "nx = 128\n",
    "ny = copy.copy(nx)\n",
    "\n",
    "dx = np.abs(x_max - x_min) / nx\n",
    "dy = np.abs(y_max - y_min) / ny"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6acc34fc",
   "metadata": {},
   "source": [
    "### 01-03. Allocation of memory on each variable (n, u1, u2) and mesh (X, Y) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4d01907",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variables\n",
    "#eta  = np.zeros(shape=(2, nx, ny)) #(time(old(0), new(1)), x, y)\n",
    "#velx = copy.deepcopy(eta) # deep copy\n",
    "#vely = copy.deepcopy(eta)\n",
    "eta  = np.zeros(shape=(3, nx, ny)) #(time(old(0), new(1), half(2)), x, y) for Leapfrog integration\n",
    "velx = copy.deepcopy(eta) # deep copy\n",
    "vely = copy.deepcopy(eta)\n",
    "\n",
    "# Grids\n",
    "x = np.linspace(x_min, x_max, nx)\n",
    "y = copy.copy(x) # shallow copy\n",
    "X, Y = np.meshgrid(x, y) # For visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e60dd408",
   "metadata": {},
   "source": [
    "\n",
    "### 01-04. Setup immersed boundary (IB): IB method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "706634be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wall: circle\n",
    "wall_radius = 0.8 * x_max\n",
    "\n",
    "# Distance (from wall) Field; (-): internal, (+): external\n",
    "distance = wall_radius - np.sqrt(X**2 + Y**2)\n",
    "\n",
    "# Epsilon Field\n",
    "sml = 1e-10\n",
    "thickness = 0.1\n",
    "eps = 0.5 * (1.0 - sml) * np.tanh(distance / thickness) + 0.5 * (1.0 + sml)\n",
    "eps[:, 0] = 0.0\n",
    "eps[:, ny-1] = 0.0\n",
    "eps[0, :] = 0.0\n",
    "eps[nx-1, :] = 0.0\n",
    "#eps = np.ones(shape=(nx, ny))\n",
    "plt.imshow(eps)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5441f19e",
   "metadata": {},
   "source": [
    "## 02. Time integration\n",
    "### 02-01. Initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f18baa48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Continuous I.C.\n",
    "eta[0,:,:] = rho0 + 0.001 * np.exp(-0.5 * (X**2 + Y**2))\n",
    "\n",
    "## Discontinuous I.C.\n",
    "#for i in range(0, nx):\n",
    "#    for j in range(0, ny):\n",
    "#        eta[0, i, j] = np.where((X[i, j]**2 + Y[i, j]**2) < (0.2*x_max)**2, rho0*1.1, rho0)\n",
    "\n",
    "velx[0,:,:] = 0.\n",
    "vely[0,:,:] = 0.\n",
    "\n",
    "fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"})\n",
    "surf = ax.plot_surface(X, Y, eta[0, :, :], cmap=cm.coolwarm,\n",
    "                       vmin=(1.0-0.08)*rho0, vmax=(1.0+0.08)*rho0, \n",
    "                       linewidth=0, antialiased=False)\n",
    "\n",
    "# ax.view_init(elev=10,azim=10)\n",
    "\n",
    "ax.set_xlim(x_min, x_max)\n",
    "ax.set_ylim(y_min, y_max)\n",
    "ax.set_zlim((1.0-0.08)*rho0, (1.0+0.08)*rho0)\n",
    "\n",
    "formatter = ticker.ScalarFormatter(useMathText=True)\n",
    "# Formatting of scientific representation\n",
    "formatter.set_scientific(True) \n",
    "# Formatting of floating point representation\n",
    "formatter.set_powerlimits((-2,2)) \n",
    "formatter.set_scientific\n",
    "ax.zaxis.set_major_formatter(formatter)\n",
    "\n",
    "# Make the panes transparent(1)\n",
    "# ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "# ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "# ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# # Make the panes transparent(2)\n",
    "ax.xaxis.pane.fill = False\n",
    "ax.yaxis.pane.fill = False\n",
    "ax.zaxis.pane.fill = False\n",
    "# ## Set color to white\n",
    "# ax.xaxis.pane.set_edgecolor('w')\n",
    "# ax.yaxis.pane.set_edgecolor('w')\n",
    "# ax.zaxis.pane.set_edgecolor('w')\n",
    "\n",
    "# Make the grid lines transparent(1)\n",
    "# ax.xaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "# ax.yaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "# ax.zaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "\n",
    "# Make the grid lines transparent(2): get rid of grid\n",
    "# ax.grid(False)\n",
    "\n",
    "# Add a color bar which maps values to colors.\n",
    "fig.colorbar(surf, shrink=0.5, aspect=15, pad=0.15, extend='both')\n",
    "\n",
    "ax.set_title('Elapsed Time: ' + str(0) + '[s]', loc='left')\n",
    "ax.set_xlabel('x [m]')\n",
    "ax.set_ylabel('y [m]')\n",
    "\n",
    "# Save figure\n",
    "id_snapshot = 0\n",
    "fname = 'central2d/img_finite_media/leapfrog/ib/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#fname = 'central2d/img_finite_media/leapfrog/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#fname = 'central2d/img_finite_media/leapfrog/discontinuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "##fname = 'central2d/img_finite_media/euler/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#fname = 'central2d/img_finite_media/euler/discontinuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "plt.savefig(fname, format='tiff', bbox_inches='tight', dpi=200)\n",
    "id_snapshot += 1\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc822aa",
   "metadata": {},
   "source": [
    "### 02-02. Time integration\n",
    "#### - Time: Explicit Euler\n",
    "#### - Space: Central difference "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc9c9dbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#for l in range(1, nt-1):\n",
    "#    for i in range(1, nx-1):\n",
    "#        for j in range(1, ny-1):\n",
    "#            #01. n\n",
    "#            eta[1, i, j] = eta[0, i, j] - dt * (\n",
    "#                                #01-01. rho * dui / dxi: Central difference\n",
    "#                                rho0 * (\n",
    "#                                    #01-01-01. du1 / dx\n",
    "#                                    0.5 * (velx[0, i+1, j] - velx[0, i-1, j]) / dx\n",
    "#\n",
    "#                                    #01-01-02. du2 / dx2\n",
    "#                                    + 0.5 * (vely[0, i, j+1] - vely[0, i, j-1]) / dy\n",
    "#                                )\n",
    "#\n",
    "#                                #01-02. ui * dn / dxi: Central difference \n",
    "#                                + (\n",
    "#                                    #01-02-01. u1 * dn/dx1\n",
    "#                                    velx[0, i, j] * (0.5 * (eta[0, i+1, j] - eta[0, i-1, j]) / dx)\n",
    "#\n",
    "#                                    #01-02-02. u2 * dn/dx2\n",
    "#                                    + vely[0, i, j] * (0.5 * (eta[0, i, j+1] - eta[0, i, j-1]) / dy)\n",
    "#                                ) \n",
    "#\n",
    "##                                #01-01. dui / dxi: Upwinding\n",
    "##                                rho0 * (\n",
    "##                                    #01-01-01. du1 / dx1\n",
    "##                                    np.where(velx[0, i, j] > 0, 1., 0.) * (velx[0, i, j] - velx[0, i-1, j]) / dx\n",
    "##                                    + np.where(velx[0, i, j] < 0, 1., 0.) * (velx[0, i+1, j] - velx[0, i, j]) / dx\n",
    "##    \n",
    "##                                    #01-01-02. du2 / dx2\n",
    "##                                    + np.where(vely[0, i, j] > 0, 1., 0.) * (vely[0, i, j] - vely[0, i, j-1]) / dy\n",
    "##                                    + np.where(vely[0, i, j] < 0, 1., 0.) * (vely[0, i, j+1] - vely[0, i, j]) / dy\n",
    "##                                )\n",
    "##\n",
    "##                                #01-02. ui * dn / dxi: Upwinding \n",
    "##                                + (\n",
    "##                                    #01-02-01. u1 * dn/dx1\n",
    "##                                    velx[0, i, j] * (\n",
    "##                                        np.where(velx[0, i, j] > 0, 1., 0.) * (eta[0, i, j] - eta[0, i-1, j]) / dx\n",
    "##                                        + np.where(velx[0, i, j] < 0, 1., 0.) * (eta[0, i+1, j] - eta[0, i, j]) / dx\n",
    "##                                    )\n",
    "##\n",
    "##                                    #01-02-02. u2 * dn/dx2\n",
    "##                                    + vely[0, i, j] * (\n",
    "##                                        np.where(vely[0, i, j] > 0, 1., 0.) * (eta[0, i, j] - eta[0, i, j-1]) / dy\n",
    "##                                        + np.where(vely[0, i, j] < 0, 1., 0.) * (eta[0, i, j+1] - eta[0, i, j]) / dy\n",
    "##                                    )\n",
    "##                                )\n",
    "#                            )\n",
    "#            #02. u1    \n",
    "#            velx[1, i, j] = velx[0, i, j] - dt * (\n",
    "#                                #02-01. uj * du1 / dxj; Central difference\n",
    "#                                #02-01-01. u1 * du1 / dx1\n",
    "#                                velx[0, i, j] * (\n",
    "#                                    0.5 * (velx[0, i+1, j] - velx[0, i-1, j]) / dx\n",
    "#                                )\n",
    "#                                \n",
    "#                                #02-01-02. u2 * du1 / dx2\n",
    "#                                + vely[0, i, j] * (\n",
    "#                                    0.5 * (velx[0, i, j+1] - velx[0, i, j-1]) / dy\n",
    "#                                )\n",
    "#                                \n",
    "#                                #02-02. g * dn / dx1; Central difference\n",
    "#                                + g * (\n",
    "#                                    0.5 * (eta[0, i+1, j] - eta[0, i-1, j]) / dx\n",
    "#                                )\n",
    "#\n",
    "##                                #02-01. uj * du1 / dxj; Upwinding \n",
    "##                                #02-01-01. u1 * du1 / dx1\n",
    "##                                velx[0, i, j] * (\n",
    "##                                    np.where(velx[0, i, j] > 0, 1., 0.) * (velx[0, i, j] - velx[0, i-1, j]) / dx\n",
    "##                                    + np.where(velx[0, i, j] < 0, 1., 0.) * (velx[0, i+1, j] - velx[0, i, j]) / dx\n",
    "##                                )\n",
    "##                                \n",
    "##                                #02-01-02. u2 * du1 / dx2\n",
    "##                                + vely[0, i, j] * (\n",
    "##                                    np.where(vely[0, i, j] > 0, 1., 0.) * (velx[0, i, j] - velx[0, i, j-1]) / dy\n",
    "##                                    + np.where(vely[0, i, j] < 0, 1., 0.) * (velx[0, i, j+1] - velx[0, i, j]) / dy\n",
    "##                                )\n",
    "##                                \n",
    "##                                #02-02. g * dn / dx1; Upwinding\n",
    "##                                + g * (\n",
    "##                                    np.where(velx[0, i, j] > 0, 1., 0.) * (eta[0, i, j] - eta[0, i-1, j]) / dx\n",
    "##                                    + np.where(velx[0, i, j] < 0, 1., 0.) * (eta[0, i+1, j] - eta[0, i, j]) / dx                                   \n",
    "##                                )\n",
    "#                            )\n",
    "#            #03. u2    \n",
    "#            vely[1, i, j] = vely[0, i, j] - dt * (\n",
    "#                                #03-01. uj * du2 / dxj; Central difference\n",
    "#                                #03-01-01. u1 * du2 / dx1\n",
    "#                                velx[0, i, j] * 0.5 * (vely[0, i+1, j] - vely[0, i-1, j]) / dx\n",
    "#                               \n",
    "#                                #03-01-02. u2 * du2 / dx2\n",
    "#                                + vely[0, i, j] * 0.5 * (vely[0, i, j+1] - vely[0, i, j-1]) / dy\n",
    "#                               \n",
    "#                                #03-02. g * dn / dx2; Central difference\n",
    "#                                + g * 0.5 * (eta[0, i, j+1] - eta[0, i, j-1]) / dy\n",
    "#\n",
    "##                                #03-01. uj * du2 / dxj; Upwinding\n",
    "##                                #03-01-01. u1 * du2 / dx1\n",
    "##                                velx[0, i, j] * (\n",
    "##                                    np.where(velx[0, i, j] > 0, 1., 0.) * (vely[0, i, j] - vely[0, i-1, j]) / dx\n",
    "##                                    + np.where(velx[0, i, j] < 0, 1., 0.) * (vely[0, i+1, j] - vely[0, i, j]) / dx\n",
    "##                                )\n",
    "##                               \n",
    "##                               #03-01-02. u2 * du2 / dx2\n",
    "##                                + vely[0, i, j] * (\n",
    "##                                    np.where(velx[0, i, j] > 0, 1., 0.) * (vely[0, i, j] - vely[0, i, j-1]) / dx\n",
    "##                                    + np.where(vely[0, i, j] < 0, 1., 0.) * (vely[0, i, j+1] - vely[0, i, j]) / dy\n",
    "##                                )\n",
    "##                               \n",
    "##                               #03-02. g * dn / dx2; Upwinding\n",
    "##                                + g * (\n",
    "##                                    np.where(vely[0, i, j] > 0, 1., 0.) * (eta[0, i, j] - eta[0, i-1, j]) / dy\n",
    "##                                    + np.where(vely[0, i, j] < 0, 1., 0.) * (eta[0, i+1, j] - eta[0, i, j]) / dy                                   \n",
    "##                                )\n",
    "#                            )\n",
    "#\n",
    "#    # Apply bouncary condition (B.C.) \n",
    "#    # Neumann B.C. for n         \n",
    "#    eta[1, 0, :] = eta[1, 1, :]\n",
    "#    eta[1, nx-1, :] = eta[1, nx-2, :]\n",
    "#    eta[1, :, 0] = eta[1, :, 1]\n",
    "#    eta[1, :, ny-1] = eta[1, :, ny-2]\n",
    "#   \n",
    "#    # Dirichlet B.C. for u, v        \n",
    "#    velx[1, 0, :] = 0.\n",
    "#    velx[1, nx-1, :] = 0.\n",
    "#    velx[1, :, 0] = 0.\n",
    "#    velx[1, :, ny-1] = 0.\n",
    "#   \n",
    "#    vely[1, 0, :] = 0.\n",
    "#    vely[1, nx-1, :] = 0.\n",
    "#    vely[1, :, 0] = 0.\n",
    "#    vely[1, :, ny-1] = 0.\n",
    "#    \n",
    "#    # Update (n)th timestep data with (n+1) timestep data \n",
    "#    # 0: (n)th timestep (Present)\n",
    "#    # 1: (n+1)th timestep (Future)\n",
    "#    eta[0, :, :] = copy.deepcopy(eta[1, :, :])\n",
    "#    velx[0, :, :] = copy.deepcopy(velx[1, :, :])\n",
    "#    vely[0, :, :] = copy.deepcopy(vely[1, :, :])\n",
    "#\n",
    "#    # Sampling/Visualizing the data\n",
    "#    if(l % sample_freq) == 0:\n",
    "#        fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"})\n",
    "#        surf = ax.plot_surface(X, Y, eta[0, :, :], cmap=cm.coolwarm,\n",
    "#                                vmin=(1.0-0.08)*rho0, vmax=(1.0+0.08)*rho0, \n",
    "#                                linewidth=0, antialiased=False)\n",
    "#\n",
    "##        ax.view_init(elev=10,azim=10)\n",
    "#\n",
    "#        ax.set_xlim(x_min, x_max)\n",
    "#        ax.set_ylim(y_min, y_max)\n",
    "#        ax.set_zlim((1.0-0.08)*rho0, (1.0+0.08)*rho0)\n",
    "#\n",
    "#        formatter = ticker.ScalarFormatter(useMathText=True)\n",
    "#\n",
    "#        # Formatting of scientific representation\n",
    "#        formatter.set_scientific(True) \n",
    "#\n",
    "#        # Formatting of floating point representation\n",
    "#        formatter.set_powerlimits((-2,2)) \n",
    "#        formatter.set_scientific\n",
    "#        ax.zaxis.set_major_formatter(formatter)\n",
    "#\n",
    "##        # Make the panes transparent(1)\n",
    "##        ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "##        ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "##        ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "#\n",
    "#        # Make the panes transparent(2)\n",
    "#        ax.xaxis.pane.fill = False\n",
    "#        ax.yaxis.pane.fill = False\n",
    "#        ax.zaxis.pane.fill = False\n",
    "#        # Set color to white\n",
    "#        ax.xaxis.pane.set_edgecolor('w')\n",
    "#        ax.yaxis.pane.set_edgecolor('w')\n",
    "#        ax.zaxis.pane.set_edgecolor('w')\n",
    "#\n",
    "#        # Make the grid lines transparent(1)\n",
    "##        ax.xaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "##        ax.yaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "##        ax.zaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "#\n",
    "#        # Make the grid lines transparent(2): get rid of grid\n",
    "##        ax.grid(False)\n",
    "#\n",
    "#        # Add a color bar which maps values to colors.\n",
    "#        fig.colorbar(surf, shrink=0.5, aspect=15, pad=0.15, extend='both')\n",
    "#\n",
    "#        # Set title\n",
    "#        present_time = t_sta + l * dt \n",
    "#        ax.set_title('Elapsed Time: ' + str(present_time) + '[s]', loc='left')\n",
    "#        ax.set_xlabel('x [m]')\n",
    "#        ax.set_ylabel('y [m]')\n",
    "#\n",
    "#        # Save figure\n",
    "#        #fname = 'central2d/img_finite_media/euler/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#        fname = 'central2d/img_finite_media/euler/discontinuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#        plt.savefig(fname, format='tiff', bbox_inches='tight', dpi=200)\n",
    "#        id_snapshot += 1\n",
    "#\n",
    "#        plt.show()\n",
    "#            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41b36bb1",
   "metadata": {},
   "source": [
    "\n",
    "### 02-02. Time integration\n",
    "#### - Time: Leapfrog integration**\n",
    "#### - Space: Central difference \n",
    "\n",
    "###### **Unlike Euler integration, it is stable for oscillatory motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c5191de",
   "metadata": {},
   "outputs": [],
   "source": [
    "for l in range(1, nt-1):\n",
    "\n",
    "    # Predict velocity at (l+1/2) timestep: u1^(l+1/2) and u2^(l+1/2)\n",
    "    for i in range(1, nx-1):\n",
    "        for j in range(1, ny-1):\n",
    "\n",
    "            #01. u1^(l+1/2)    \n",
    "            velx[2, i, j] = velx[0, i, j] - 0.5 * dt * (\n",
    "                                #01-01. uj * du1 / dxj; Central difference\n",
    "                                #01-01-01. u1 * du1 / dx1\n",
    "                                velx[0, i, j] * 0.5 * (velx[0, i+1, j] - velx[0, i-1, j]) / dx\n",
    "                                \n",
    "                                #01-01-02. u2 * du1 / dx2\n",
    "                                + vely[0, i, j] * 0.5 * (velx[0, i, j+1] - velx[0, i, j-1]) / dy\n",
    "                                \n",
    "                                #01-02. g * dn / dx1; Central difference\n",
    "#                                + g * 0.5 * (eta[0, i+1, j] - eta[0, i-1, j]) / dx\n",
    "                                + eps[i, j] * g * 0.5 * (eta[0, i+1, j] - eta[0, i-1, j]) / dx\n",
    "                            )\n",
    "            #02. u2^(l+1/2)    \n",
    "            vely[2, i, j] = vely[0, i, j] - 0.5 * dt * (\n",
    "                                #02-01. uj * du2 / dxj; Central difference\n",
    "                                #02-01-01. u1 * du2 / dx1\n",
    "                                velx[0, i, j] * 0.5 * (vely[0, i+1, j] - vely[0, i-1, j]) / dx\n",
    "                               \n",
    "                                #02-01-02. u2 * du2 / dx2\n",
    "                                + vely[0, i, j] * 0.5 * (vely[0, i, j+1] - vely[0, i, j-1]) / dy\n",
    "                               \n",
    "                                #02-02. g * dn / dx2; Central difference\n",
    "#                                + g * 0.5 * (eta[0, i, j+1] - eta[0, i, j-1]) / dy\n",
    "                                + eps[i, j] * g * 0.5 * (eta[0, i, j+1] - eta[0, i, j-1]) / dy\n",
    "                            )\n",
    "\n",
    "    # Apply bouncary condition (B.C.) \n",
    "    # Dirichlet B.C. for u, v        \n",
    "    velx[2, 0, :] = 0.\n",
    "    velx[2, nx-1, :] = 0.\n",
    "    velx[2, :, 0] = 0.\n",
    "    velx[2, :, ny-1] = 0.\n",
    "   \n",
    "    vely[2, 0, :] = 0.\n",
    "    vely[2, nx-1, :] = 0.\n",
    "    vely[2, :, 0] = 0.\n",
    "    vely[2, :, ny-1] = 0.\n",
    "\n",
    "    velx[2, :, :] = eps[:, :] * velx[2, :, :]\n",
    "    vely[2, :, :] = eps[:, :] * vely[2, :, :]\n",
    "\n",
    "    # Obtain n^(l+1) using u1^(l+1/2) and u2^(l+1/2)\n",
    "    for i in range(1, nx-1):\n",
    "        for j in range(1, ny-1):\n",
    "\n",
    "            #03. n^(l+1)\n",
    "            eta[1, i, j] = eta[0, i, j] - dt * (\n",
    "                                #03-01. rho * dui / dxi: Central difference\n",
    "                                rho0 * (\n",
    "                                    #03-01-01. du1 / dx\n",
    "                                    0.5 * (velx[2, i+1, j] - velx[2, i-1, j]) / dx\n",
    "\n",
    "                                    #03-01-02. du2 / dx2\n",
    "                                    + 0.5 * (vely[2, i, j+1] - vely[2, i, j-1]) / dy\n",
    "                                )\n",
    "\n",
    "                                #03-02. ui * dn / dxi: Central difference \n",
    "                                + eps[i, j] * (\n",
    "#                                + (\n",
    "                                    #03-02-01. u1 * dn/dx1\n",
    "                                    velx[2, i, j] * (0.5 * (eta[0, i+1, j] - eta[0, i-1, j]) / dx)\n",
    "\n",
    "                                    #03-02-02. u2 * dn/dx2\n",
    "                                    + vely[2, i, j] * (0.5 * (eta[0, i, j+1] - eta[0, i, j-1]) / dy)\n",
    "                                ) \n",
    "                            )\n",
    "\n",
    "    # Apply bouncary condition (B.C.) \n",
    "    # Neumann B.C. for n         \n",
    "    eta[1, 0, :] = eta[1, 1, :]\n",
    "    eta[1, nx-1, :] = eta[1, nx-2, :]\n",
    "    eta[1, :, 0] = eta[1, :, 1]\n",
    "    eta[1, :, ny-1] = eta[1, :, ny-2]\n",
    "\n",
    "    # Obtain velocity (u1, u2) using n^(l+1) and predited value of u1^(l+1/2) and u2^(l+1/2)\n",
    "    for i in range(1, nx-1):\n",
    "        for j in range(1, ny-1):\n",
    "            #02. u1    \n",
    "            velx[1, i, j] = velx[2, i, j] - 0.5 * dt * (\n",
    "                                #02-01. uj * du1 / dxj; Central difference\n",
    "                                #02-01-01. u1 * du1 / dx1\n",
    "                                velx[2, i, j] * 0.5 * (velx[2, i+1, j] - velx[2, i-1, j]) / dx\n",
    "                                \n",
    "                                #02-01-02. u2 * du1 / dx2\n",
    "                                + vely[2, i, j] * 0.5 * (velx[2, i, j+1] - velx[2, i, j-1]) / dy\n",
    "                                \n",
    "                                #02-02. g * dn / dx1; Central difference\n",
    "#                                + g * 0.5 * (eta[1, i+1, j] - eta[1, i-1, j]) / dx\n",
    "                                + eps[i, j] * g * 0.5 * (eta[1, i+1, j] - eta[1, i-1, j]) / dx\n",
    "                            )\n",
    "            #03. u2    \n",
    "            vely[1, i, j] = vely[2, i, j] - 0.5 * dt * (\n",
    "                                #03-01. uj * du2 / dxj; Central difference\n",
    "                                #03-01-01. u1 * du2 / dx1\n",
    "                                velx[2, i, j] * 0.5 * (vely[2, i+1, j] - vely[2, i-1, j]) / dx\n",
    "                               \n",
    "                                #03-01-02. u2 * du2 / dx2\n",
    "                                + vely[2, i, j] * 0.5 * (vely[2, i, j+1] - vely[2, i, j-1]) / dy\n",
    "                               \n",
    "                                #03-02. g * dn / dx2; Central difference\n",
    "#                                + g * 0.5 * (eta[1, i, j+1] - eta[1, i, j-1]) / dy\n",
    "                                + eps[i, j] * g * 0.5 * (eta[1, i, j+1] - eta[1, i, j-1]) / dy\n",
    "                            )\n",
    "\n",
    "    # Apply bouncary condition (B.C.) \n",
    "    # Dirichlet B.C. for u, v        \n",
    "    velx[1, 0, :] = 0.\n",
    "    velx[1, nx-1, :] = 0.\n",
    "    velx[1, :, 0] = 0.\n",
    "    velx[1, :, ny-1] = 0.\n",
    "   \n",
    "    vely[1, 0, :] = 0.\n",
    "    vely[1, nx-1, :] = 0.\n",
    "    vely[1, :, 0] = 0.\n",
    "    vely[1, :, ny-1] = 0.\n",
    "\n",
    "    velx[1, :, :] = eps[:, :] * velx[1, :, :]\n",
    "    vely[1, :, :] = eps[:, :] * vely[1, :, :]\n",
    "\n",
    "    # Update (l)th timestep data with (l+1) timestep data \n",
    "    # 0: (l)th timestep (Present)\n",
    "    # 1: (l+1)th timestep (Future)\n",
    "    eta[0, :, :] = copy.deepcopy(eta[1, :, :])\n",
    "    velx[0, :, :] = copy.deepcopy(velx[1, :, :])\n",
    "    vely[0, :, :] = copy.deepcopy(vely[1, :, :])\n",
    "\n",
    "    # Sampling/Visualizing the data\n",
    "    if(l % sample_freq) == 0:\n",
    "        fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"})\n",
    "        surf = ax.plot_surface(X, Y, eta[0, :, :], cmap=cm.coolwarm,\n",
    "                                vmin=(1.0-0.08)*rho0, vmax=(1.0+0.08)*rho0, \n",
    "                                linewidth=0, antialiased=False)\n",
    "\n",
    "#        ax.view_init(elev=10,azim=10)\n",
    "\n",
    "        ax.set_xlim(x_min, x_max)\n",
    "        ax.set_ylim(y_min, y_max)\n",
    "        ax.set_zlim((1.0-0.08)*rho0, (1.0+0.08)*rho0)\n",
    "\n",
    "        formatter = ticker.ScalarFormatter(useMathText=True)\n",
    "\n",
    "        # Formatting of scientific representation\n",
    "        formatter.set_scientific(True) \n",
    "\n",
    "        # Formatting of floating point representation\n",
    "        formatter.set_powerlimits((-2,2)) \n",
    "        formatter.set_scientific\n",
    "        ax.zaxis.set_major_formatter(formatter)\n",
    "\n",
    "#        # Make the panes transparent(1)\n",
    "#        ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "#        ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "#        ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "        # Make the panes transparent(2)\n",
    "        ax.xaxis.pane.fill = False\n",
    "        ax.yaxis.pane.fill = False\n",
    "        ax.zaxis.pane.fill = False\n",
    "        # Set color to white\n",
    "        ax.xaxis.pane.set_edgecolor('w')\n",
    "        ax.yaxis.pane.set_edgecolor('w')\n",
    "        ax.zaxis.pane.set_edgecolor('w')\n",
    "\n",
    "        # Make the grid lines transparent(1)\n",
    "#        ax.xaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "#        ax.yaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "#        ax.zaxis._axinfo['grid']['color'] = (1.0, 1.0, 1.0, 0.0)\n",
    "\n",
    "        # Make the grid lines transparent(2): get rid of grid\n",
    "#        ax.grid(False)\n",
    "\n",
    "        # Add a color bar which maps values to colors.\n",
    "        fig.colorbar(surf, shrink=0.5, aspect=15, pad=0.15, extend='both')\n",
    "\n",
    "        # Set title\n",
    "        present_time = t_sta + l * dt \n",
    "        ax.set_title('Elapsed Time: ' + str(present_time) + '[s]', loc='left')\n",
    "        ax.set_xlabel('x [m]')\n",
    "        ax.set_ylabel('y [m]')\n",
    "\n",
    "        # Save figure\n",
    "        fname = 'central2d/img_finite_media/leapfrog/ib/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#        fname = 'central2d/img_finite_media/leapfrog/continuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "#        fname = 'central2d/img_finite_media/leapfrog/discontinuous_ic/ww2d_height_' + str(id_snapshot).zfill(5) + '.tiff'\n",
    "        plt.savefig(fname, format='tiff', bbox_inches='tight', dpi=200)\n",
    "        id_snapshot += 1\n",
    "\n",
    "        plt.show()\n",
    "        print('Total variation = ' + str(dx * dx * np.sum((eta[1, :, :] - rho0) * eps[:, :])))\n",
    "            "
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}