phasefield-gravity-wave.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import time\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "from pystencils.session import *\n",
    "from lbmpy.session import *\n",
    "\n",
    "from pystencils.simp import sympy_cse\n",
    "from pystencils.boundaries import BoundaryHandling\n",
    "\n",
    "from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters\n",
    "from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle\n",
    "from lbmpy.phasefield_allen_cahn.kernel_equations import *\n",
    "\n",
    "from lbmpy.advanced_streaming import LBMPeriodicityHandling\n",
    "from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `pycuda` is installed the simulation automatically runs on GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No pycuda installed\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    import pycuda\n",
    "except ImportError:\n",
    "    pycuda = None\n",
    "    gpu = False\n",
    "    target = ps.Target.CPU\n",
    "    print('No pycuda installed')\n",
    "\n",
    "if pycuda:\n",
    "    gpu = True\n",
    "    target = ps.Target.GPU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gravity wave simulated with a phase-field model for immiscible fluids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry Setup\n",
    "\n",
    "First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stencil_phase = LBStencil(Stencil.D2Q9)\n",
    "stencil_hydro = LBStencil(Stencil.D2Q9)\n",
    "assert(stencil_hydro.D == stencil_phase.D)\n",
    "\n",
    "dimensions = stencil_phase.D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "timesteps = 238800\n",
      "Re = 10\n",
      "Pe = 0.41876046901171776\n",
      "Cn = 0.1\n",
      "domain_width = 50\n",
      "fluid_depth = 25.0\n",
      "amplitude = 0.5\n",
      "relaxation_rate_heavy = 1.99\n",
      "mobility = 0.02\n",
      "interface_width = 5\n",
      "density_heavy = 1.0\n",
      "density_light = 0.001\n",
      "density_ratio = 1000\n",
      "kinematic_viscosity_heavy = 0.0008375209380234356\n",
      "kinematic_viscosity_light = 0.0008375209380234356\n",
      "kinematic_viscosity_ratio = 1\n",
      "dynamic_viscosity_heavy = 0.0008375209380234356\n",
      "dynamic_viscosity_light = 8.375209380234357e-07\n",
      "dynamic_viscosity_ratio = 1000.0\n",
      "wavelength = 50\n",
      "wavenumber = 0.12566370614359174\n",
      "wave_frequency = 0.00033500837520937423\n",
      "surface_tension = 0\n",
      "gravitational_acceleration = -8.964447065459421e-07\n",
      "data_extract_frequency = 298\n",
      "vtk_output_frequency = 2985\n"
     ]
    }
   ],
   "source": [
    "# user defined input\n",
    "\n",
    "Re = 10 # Reynolds number\n",
    "domain_width = 50\n",
    "fluid_depth = 0.5 * domain_width\n",
    "amplitude = 0.01 * domain_width\n",
    "relaxation_rate_heavy = 1.99\n",
    "mobility = 0.02 # phase field mobility\n",
    "interface_width = 5 # phase field interface width\n",
    "density_heavy = 1.0 # density of heavy phase\n",
    "density_ratio = 1000\n",
    "density_light = density_heavy / density_ratio # density of light phase\n",
    "kinematic_viscosity_ratio = 1\n",
    "\n",
    "kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)\n",
    "kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio\n",
    "wavelength = domain_width\n",
    "wavenumber = 2.0 * np.pi / domain_width\n",
    "wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency\n",
    "surface_tension = 0\n",
    "gravitational_acceleration = - wave_frequency**2 / wavenumber / np.tanh(wavenumber * fluid_depth)\n",
    "Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number\n",
    "Cn = interface_width / domain_width # Cahn number\n",
    "dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy\n",
    "relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy\n",
    "kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio\n",
    "dynamic_viscosity_light = kinematic_viscosity_light * density_light\n",
    "relaxation_time_light = 3.0 * kinematic_viscosity_light\n",
    "\n",
    "timesteps = int(80 / wave_frequency)\n",
    "\n",
    "data_extract_frequency = int(0.1 / wave_frequency)\n",
    "vtk_output_frequency = int(1 / wave_frequency)\n",
    "vtk_output_path = \"vtk_out/gravity-wave\"\n",
    "vtk_base_directory = vtk_output_path.split(\"/\")[0] # create directory for vtk-output if it does not yet exist\n",
    "if not os.path.exists(vtk_base_directory):\n",
    "    os.mkdir(os.getcwd() + \"/\" + vtk_base_directory)\n",
    "\n",
    "domain_size = (domain_width, domain_width)\n",
    "filename = \"pf-re-\" + str(Re) + \"-resolution-\" + str(domain_width) + \".txt\"\n",
    "\n",
    "print(\"timesteps =\", timesteps)\n",
    "print(\"Re =\", Re)\n",
    "print(\"Pe =\", Pe)\n",
    "print(\"Cn =\", Cn)\n",
    "print(\"domain_width =\", domain_width)\n",
    "print(\"fluid_depth =\", fluid_depth)\n",
    "print(\"amplitude =\", amplitude)\n",
    "print(\"relaxation_rate_heavy =\", relaxation_rate_heavy)\n",
    "print(\"mobility =\", mobility)\n",
    "print(\"interface_width =\", interface_width)\n",
    "print(\"density_heavy =\", density_heavy)\n",
    "print(\"density_light =\", density_light)\n",
    "print(\"density_ratio =\", density_ratio)\n",
    "print(\"kinematic_viscosity_heavy =\", kinematic_viscosity_heavy)\n",
    "print(\"kinematic_viscosity_light =\", kinematic_viscosity_light)\n",
    "print(\"kinematic_viscosity_ratio =\", kinematic_viscosity_ratio)\n",
    "print(\"dynamic_viscosity_heavy =\", dynamic_viscosity_heavy)\n",
    "print(\"dynamic_viscosity_light =\", dynamic_viscosity_light)\n",
    "print(\"dynamic_viscosity_ratio =\", dynamic_viscosity_heavy/dynamic_viscosity_light)\n",
    "print(\"wavelength =\", wavelength)\n",
    "print(\"wavenumber =\", wavenumber)\n",
    "print(\"wave_frequency =\", wave_frequency)\n",
    "print(\"surface_tension =\", surface_tension)\n",
    "print(\"gravitational_acceleration =\", gravitational_acceleration)\n",
    "print(\"data_extract_frequency =\", data_extract_frequency)\n",
    "print(\"vtk_output_frequency =\", vtk_output_frequency)"
   ]