{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kelvin-Helmholtz Instability\n", "\n", "\n", "\n", "(image: Chuck Doswell)\n", "\n", "We will simulate the incompressible Kelvin-Helmholtz problem. We non-dimensionalize the problem by taking the box height to be one and the jump in velocity to be one. Then the Reynolds number is given by\n", "\n", "$ {\\rm Re} = \\frac{U H}{\\nu} = \\frac{1}{\\nu}. $\n", "\n", "We use no slip boundary conditions, and a box with aspect ratio $L/H=2$. The initial velocity profile is given by a hyperbolic tangent, and only a single mode is initially excited. We will also track a passive scalar which will help us visualize the instability.\n", "\n", "First, we import the necessary modules." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import h5py\n", "import time" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from dedalus import public as de\n", "from dedalus.extras import flow_tools" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we set logging to `INFO` level. Currently, by default, Dedalus sets its logging output to `DEBUG`, which produces more info than we need here." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import logging\n", "root = logging.root\n", "for h in root.handlers:\n", " h.setLevel(\"INFO\")\n", " \n", "logger = logging.getLogger(__name__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Domain\n", "\n", "First, we will specify the domain. Domains are built by taking the direct product of bases. Here we are running a 2D simulation, so we will define $x$ and $y$ bases. From these, we build the domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#Aspect ratio 2\n", "Lx, Ly = (2., 1.)\n", "nx, ny = (128, 64)\n", "\n", "# Create bases and domain\n", "x_basis = de.Fourier('x', nx, interval=(0, Lx), dealias=3/2)\n", "y_basis = de.Chebyshev('y',ny, interval=(-Ly/2, Ly/2), dealias=3/2)\n", "domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Equations\n", "\n", "Next we will define the equations that will be solved on this domain. The equations are\n", "\n", "$$ \\partial_t u + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} u + \\partial_x p = \\frac{1}{{\\rm Re}} \\nabla^2 u $$\n", "$$ \\partial_t v + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} v + \\partial_y p = \\frac{1}{{\\rm Re}} \\nabla^2 v $$\n", "$$ \\boldsymbol{\\nabla}\\boldsymbol{\\cdot}\\boldsymbol{u} = 0 $$\n", "$$ \\partial_t s + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} s = \\frac{1}{{\\rm Re}{\\rm Sc}} \\nabla^2 s $$\n", "\n", "The equations are written such that the left-hand side (LHS) is treated implicitly, and the right-hand side (RHS) is treated explicitly. The LHS is limited to only linear terms, though linear terms can also be placed on the RHS. Since $y$ is our special direction in this example, we also restrict the LHS to be at most first order in derivatives with respect to $y$.\n", "\n", "We also set the parameters, the Reynolds and Schmidt numbers." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "problem = de.IVP(domain, variables=['p','u','v','uy','vy','s','sy'])\n", "\n", "Reynolds = 1e4\n", "Schmidt = 1.\n", "\n", "problem.parameters['Re'] = Reynolds\n", "problem.parameters['Sc'] = Schmidt\n", "\n", "problem.add_equation(\"dt(u) + dx(p) - 1/Re*(dx(dx(u)) + dy(uy)) = - u*dx(u) - v*uy\")\n", "problem.add_equation(\"dt(v) + dy(p) - 1/Re*(dx(dx(v)) + dy(vy)) = - u*dx(v) - v*vy\")\n", "problem.add_equation(\"dx(u) + vy = 0\")\n", "\n", "problem.add_equation(\"dt(s) - 1/(Re*Sc)*(dx(dx(s)) + dy(sy)) = - u*dx(s) - v*sy\")\n", "\n", "problem.add_equation(\"uy - dy(u) = 0\")\n", "problem.add_equation(\"vy - dy(v) = 0\")\n", "problem.add_equation(\"sy - dy(s) = 0\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "problem.add_bc(\"left(u) = 0.5\")\n", "problem.add_bc(\"right(u) = -0.5\")\n", "problem.add_bc(\"left(v) = 0\")\n", "problem.add_bc(\"right(v) = 0\", condition=\"(nx != 0)\")\n", "problem.add_bc(\"integ(p,'y') = 0\", condition=\"(nx == 0)\")\n", "problem.add_bc(\"left(s) = 0\")\n", "problem.add_bc(\"right(s) = 1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have a special boundary condition for the $k_x=0$ mode (singled out by `condition=\"(dx==0)\"`). This is because the continuity equation implies $\\partial_y v=0$ if $k_x=0$; thus, $v=0$ on the top and bottom are redundant boundary conditions. We replace one of these with a gauge choice for the pressure.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timestepping\n", "\n", "We have implemented about twenty implicit-explicit timesteppers in Dedalus. This contains both multi-stage and multi-step methods. For this problem, we will use a four-stage, fourth order Runge-Kutta integrator. Changing the timestepping algorithm is as easy as changing one line of code." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ts = de.timesteppers.RK443" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial Value Problem\n", "\n", "We now have the three ingredients necessary to set up our IVP:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "solver = problem.build_solver(ts)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = domain.grid(0)\n", "y = domain.grid(1)\n", "u = solver.state['u']\n", "uy = solver.state['uy']\n", "v = solver.state['v']\n", "vy = solver.state['vy']\n", "p = solver.state['p']\n", "s = solver.state['s']\n", "sy = solver.state['sy']\n", "\n", "a = 0.05\n", "sigma = 0.2\n", "flow = -0.5\n", "u['g'] = flow*np.tanh(y/a)\n", "s['g'] = 0.5*(1+np.tanh(y/a))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "amp = -0.2\n", "sigma = 0.2\n", "v['g'] = amp*np.exp(-y**2/sigma**2)*np.sin(2*np.pi*x/Lx)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.close()\n", "plt.plot(v['g'][int(nx/4),:],y[0])\n", "plt.ylabel('y')\n", "plt.xlabel('vertical velocity')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set integration parameters and the CFL." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "solver.stop_sim_time = 3.01\n", "solver.stop_wall_time = np.inf\n", "solver.stop_iteration = np.inf\n", "\n", "initial_dt = 0.2*Lx/nx\n", "cfl = flow_tools.CFL(solver,initial_dt,safety=0.8,threshold=0.05)\n", "cfl.add_velocities(('u','v'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "We have a sophisticated analysis framework in which the user specifies analysis tasks as strings. Users can output full data cubes, slices, volume averages, and more. Here we will only output a few 2D slices, and a 1D profile of the horizontally averaged concentration field. Data is output in the hdf5 file format." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "analysis = solver.evaluator.add_file_handler('analysis_tasks', sim_dt=0.1, max_writes=50)\n", "analysis.add_task('s')\n", "analysis.add_task('u')\n", "analysis.add_task('0.5*(u**2+v**2)',name='KE',scales=(3/2,3/2))\n", "analysis.add_task('0.5*(dx(v)-uy)**2',name='enstrophy')\n", "solver.evaluator.vars['Lx'] = Lx\n", "analysis.add_task(\"integ(s,'x')/Lx\", name='s profile')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Main Loop\n", "\n", "We now have everything set up for our simulation. In Dedalus, the user writes their own main loop." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "# Make plot of scalar field\n", "x = domain.grid(0,scales=domain.dealias)\n", "y = domain.grid(1,scales=domain.dealias)\n", "xm, ym = np.meshgrid(x,y)\n", "fig, axis = plt.subplots(figsize=(8,5))\n", "s.set_scales(domain.dealias)\n", "p = axis.pcolormesh(xm, ym, s['g'].T, cmap='RdBu_r');\n", "axis.set_title('t = %f' %solver.sim_time)\n", "axis.set_xlim([0,2.])\n", "axis.set_ylim([-0.5,0.5])\n", "\n", "logger.info('Starting loop')\n", "start_time = time.time()\n", "while solver.ok:\n", " dt = cfl.compute_dt()\n", " solver.step(dt)\n", " if solver.iteration % 10 == 0:\n", " # Update plot of scalar field\n", " p.set_array(np.ravel(s['g'][:-1,:-1].T))\n", " axis.set_title('t = %f' %solver.sim_time)\n", " fig.canvas.draw()\n", "\n", "end_time = time.time()\n", "\n", "# Print statistics\n", "logger.info('Run time: %f' %(end_time-start_time))\n", "logger.info('Iterations: %i' %solver.iteration)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "As an example of doing some analysis, we will load in the horizontally averaged profiles of the scalar field $s$ and plot them." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Read in the data\n", "f = h5py.File('analysis_tasks/analysis_tasks_s1/analysis_tasks_s1_p0.h5',flag='r')\n", "y = f['/scales/y/1.0'][:]\n", "t = f['scales']['sim_time'][:]\n", "s_ave = f['tasks']['s profile'][:]\n", "f.close()\n", "\n", "s_ave = s_ave[:,0,:] # remove length-one x dimension" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axis = plt.subplots(figsize=(8,6))\n", "for i in range(0,31,6):\n", " axis.plot(y,s_ave[i,:],label='t=%4.2f' %t[i])\n", "\n", "plt.xlim([-0.5,0.5])\n", "plt.ylim([0,1])\n", "plt.ylabel(r'$\\frac{\\int \\ s dx}{L_x}$',fontsize=24)\n", "plt.xlabel(r'$y$',fontsize=24)\n", "plt.legend(loc='lower right').draw_frame(False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }