From f1a25e9da83aac2a9aa8748b20c4c339d8779cf6 Mon Sep 17 00:00:00 2001 From: David Ketcheson Date: Mon, 3 Nov 2014 17:19:57 +0300 Subject: [PATCH 1/3] Add acoustics notebook for fvmbook Chapter 3. --- .../acousimple/acoustics_simple_waves.ipynb | 120 ++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb diff --git a/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb new file mode 100644 index 00000000..8815c49b --- /dev/null +++ b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb @@ -0,0 +1,120 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:fe4ac1911f7fa07461f66cbe9ed9f96aea142f18cac76b6289795e4be5a26f38" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "heading", + "level": 1, + "metadata": {}, + "source": [ + "Simple waves in the acoustics system" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The code in this notebook runs the simulation that produces the output shown in Figures 3.1 and 3.8." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%matplotlib inline\n", + "from clawpack import pyclaw\n", + "from clawpack import riemann\n", + "import numpy as np\n", + "\n", + "import matplotlib\n", + "matplotlib.rcParams.update({'font.size': 22})\n", + "import matplotlib.pyplot as plt\n", + "from clawpack.visclaw.JSAnimation import IPython_display" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Define the domain, equations, and boundary conditions\n", + "domain = pyclaw.Domain([-1.],[1.],[800])\n", + "solver = pyclaw.ClawSolver1D(riemann.acoustics_1D)\n", + "solver.bc_lower[0] = pyclaw.BC.wall\n", + "solver.bc_upper[0] = pyclaw.BC.extrap\n", + "\n", + "num_eqn = 2\n", + "solution = pyclaw.Solution(num_eqn,domain)\n", + "solver.limiters = pyclaw.limiters.tvd.MC\n", + "\n", + "# Define the material parameters\n", + "rho = 1.0 # Density\n", + "bulk = 0.25 # Bulk modulus\n", + "solution.problem_data['rho']=rho\n", + "solution.problem_data['bulk']=bulk\n", + "solution.problem_data['zz']=np.sqrt(rho*bulk) # Impedance\n", + "solution.problem_data['cc']=np.sqrt(bulk/rho) # Sound speed\n", + "\n", + "# Define the initial condition\n", + "xc=solution.p_centers[0] # Cell centers\n", + "solution.q[0,:] = 0.5*np.exp(-80 * xc**2) + 0.5*(np.abs(xc+0.2)<0.1)\n", + "solution.q[1,:] = 0.\n", + "\n", + "# Set up the simulation and output\n", + "claw = pyclaw.Controller()\n", + "claw.solution = solution\n", + "claw.solver = solver\n", + "claw.output_format = None\n", + "claw.keep_copy = True\n", + "claw.tfinal = 3.0\n", + "claw.num_output_times = 30\n", + "claw.verbosity = 1" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "status = claw.run()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig = plt.figure(figsize=[8,4])\n", + "ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 1.0))\n", + "line, = ax.plot([], [], lw=2)\n", + "plt.xlabel('x'); plt.ylabel('p')\n", + "\n", + "def fplot(i):\n", + " frame = claw.frames[i]\n", + " q = frame.q[0,:]\n", + " line.set_data(xc, q)\n", + " return line,\n", + "\n", + "matplotlib.animation.FuncAnimation(fig, fplot, frames=len(claw.frames))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file From d6faf1c8bd94c75ed1bd24062799f50044977307 Mon Sep 17 00:00:00 2001 From: David Ketcheson Date: Mon, 3 Nov 2014 21:26:54 +0300 Subject: [PATCH 2/3] Use new ianimate function from visclaw branch. See https://github.com/clawpack/visclaw/pull/100. --- .../acousimple/acoustics_simple_waves.ipynb | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb index 8815c49b..c6002ff0 100644 --- a/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb +++ b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb @@ -1,7 +1,7 @@ { "metadata": { "name": "", - "signature": "sha256:fe4ac1911f7fa07461f66cbe9ed9f96aea142f18cac76b6289795e4be5a26f38" + "signature": "sha256:78991e4acc93c32f483391d7a48b3dfe1559a0d2ee23b89d604c56147ff47ace" }, "nbformat": 3, "nbformat_minor": 0, @@ -31,11 +31,11 @@ "from clawpack import pyclaw\n", "from clawpack import riemann\n", "import numpy as np\n", + "from clawpack.visclaw import ianimate\n", "\n", "import matplotlib\n", "matplotlib.rcParams.update({'font.size': 22})\n", - "import matplotlib.pyplot as plt\n", - "from clawpack.visclaw.JSAnimation import IPython_display" + "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, @@ -96,18 +96,7 @@ "cell_type": "code", "collapsed": false, "input": [ - "fig = plt.figure(figsize=[8,4])\n", - "ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 1.0))\n", - "line, = ax.plot([], [], lw=2)\n", - "plt.xlabel('x'); plt.ylabel('p')\n", - "\n", - "def fplot(i):\n", - " frame = claw.frames[i]\n", - " q = frame.q[0,:]\n", - " line.set_data(xc, q)\n", - " return line,\n", - "\n", - "matplotlib.animation.FuncAnimation(fig, fplot, frames=len(claw.frames))" + "ianimate.ianimate(claw,varname='Pressure')" ], "language": "python", "metadata": {}, From 10093cb1a8087dd5021bd2c2717cbea54f9f0ca3 Mon Sep 17 00:00:00 2001 From: David Ketcheson Date: Mon, 3 Nov 2014 21:40:59 +0300 Subject: [PATCH 3/3] Add notebook for fvmbook chapter 6 example (limiters). --- fvmbook/chap6/limiters.ipynb | 9445 ++++++++++++++++++++++++++++++++++ 1 file changed, 9445 insertions(+) create mode 100644 fvmbook/chap6/limiters.ipynb diff --git a/fvmbook/chap6/limiters.ipynb b/fvmbook/chap6/limiters.ipynb new file mode 100644 index 00000000..936d5b9d --- /dev/null +++ b/fvmbook/chap6/limiters.ipynb @@ -0,0 +1,9445 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:56ad361e9cd3dcd84982a6446b4bfceb0481bc3c92f266e40c51b1586028740f" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demonstration of different limiters for advection\n", + "This notebook shows the effect of limiters on advection for various initial conditions. It includes the results for first and second order methods from Figures 6.1-6.3 of [Finite Volume Methods for Hyperbolic Problems](http://depts.washington.edu/clawpack/book.html), as well as results for WENO methods." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from clawpack import pyclaw\n", + "from clawpack import riemann\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import animation\n", + "from clawpack.visclaw.JSAnimation import IPython_display" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The function below sets up an advection simulation with the specified scheme, on the unit interval with periodic boundary conditions. Since output is written at integer times, the solution at each output time has been advected precisely back to its original location. Plotting each frame shows us how the solution is modified by numerical effects as it moves once through the domain." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def setup(scheme='minmod',cfl_max=0.9,IC='gauss_square',mx=100):\n", + " if 'weno' in scheme:\n", + " solver = pyclaw.SharpClawSolver1D(riemann.advection_1D)\n", + " else:\n", + " solver = pyclaw.ClawSolver1D(riemann.advection_1D)\n", + "\n", + " solver.bc_lower[0] = pyclaw.BC.periodic\n", + " solver.bc_upper[0] = pyclaw.BC.periodic\n", + " \n", + " if scheme in ('minmod','superbee','MC','vanleer'):\n", + " solver.limiters = getattr(pyclaw.limiters.tvd,scheme)\n", + " #elif scheme == 'CT':\n", + " #solver.limiters = pyclaw.limiters.tvd.cada_torrilhon_limiter\n", + " elif scheme == 'Lax-Wendroff':\n", + " solver.limiters = 0\n", + " elif scheme == 'first-order':\n", + " solver.order = 1\n", + " elif 'weno' in scheme:\n", + " solver.weno_order = int(scheme[4:]) #weno5, weno7, ...\n", + " else:\n", + " raise Exception('Unrecognized limiter')\n", + "\n", + " solver.cfl_max = cfl_max\n", + " solver.cfl_desired = cfl_max*0.9\n", + "\n", + " x = pyclaw.Dimension('x',0.0,1.0,mx)\n", + " domain = pyclaw.Domain(x)\n", + " num_eqn = 1\n", + " state = pyclaw.State(domain,num_eqn)\n", + " state.problem_data['u']=1.\n", + "\n", + " grid = state.grid\n", + " xc = grid.x.centers\n", + " if IC=='gauss_square':\n", + " beta=200.; x0=0.3\n", + " state.q[0,:] = np.exp(-beta * (xc-x0)**2) + (xc>0.6)*(xc<0.8)\n", + " elif IC=='wavepacket':\n", + " beta=100.; x0=0.5\n", + " state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.sin(80.*xc)\n", + " else:\n", + " raise Exception('Unrecognized initial condition.')\n", + "\n", + " claw = pyclaw.Controller()\n", + " claw.solution = pyclaw.Solution(state,domain)\n", + " claw.solver = solver\n", + " claw.keep_copy = True\n", + " claw.output_format = None\n", + " claw.verbosity = 1\n", + " claw.tfinal =10.0\n", + " return claw" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Gaussian and square wave" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "results = []\n", + "schemes = ('first-order','Lax-Wendroff','minmod','superbee','MC','vanleer','weno5','weno7','weno9')\n", + "for scheme in schemes:\n", + " claw = setup(scheme=scheme)\n", + " claw.run()\n", + " results.append(claw.frames)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def animate(results,ymin=-0.1):\n", + " fig = plt.figure(figsize=(12,8))\n", + "\n", + " N = len(results)\n", + " n = int(np.ceil(np.sqrt(N)))\n", + " axes = []\n", + " gs1 = matplotlib.gridspec.GridSpec(n, n)\n", + " gs1.update(wspace=0.,hspace=0.)\n", + "\n", + " for i in range(n):\n", + " for j in range(n):\n", + " k = n*i + j\n", + " if k0:\n", + " axes[-1].yaxis.set_ticklabels(())\n", + "\n", + " lines = [0]*len(schemes)\n", + " for i in range(len(lines)):\n", + " lines[i], = axes[i].plot([], [], lw=2)\n", + "\n", + " for i,ax in enumerate(axes):\n", + " ax.set_xlim(0,1); ax.set_ylim(ymin,1.3)\n", + " ax.set_title(schemes[i], x = 0.5, y=0.8, fontsize=15 )\n", + " \n", + " xc = results[0][0].p_centers[0]\n", + "\n", + " def fplot(frame_number):\n", + " fig.suptitle('Time t = '+str(results[0][frame_number].t),fontsize=20)\n", + " for i, line in enumerate(lines):\n", + " line.set_data(xc,results[i][frame_number].q[0,:])\n", + " return lines,\n", + "\n", + " return matplotlib.animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)\n", + "\n", + "animate(results)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It's easy to see the differences in dissipation and dispersion of the various schemes. Of course, these properties also depend on the CFL number. Try changing the CFL number in the call to `setup()` above and see how the results change." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Wavepacket\n", + "To see how the limiters affect intermediate frequencies, let's advect a wavepacket." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "results = []\n", + "\n", + "for scheme in schemes:\n", + " claw = setup(scheme=scheme,IC='wavepacket',mx=300)\n", + " claw.run()\n", + " results.append(claw.frames)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 13 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "animate(results,ymin=-1.2)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + "\n", + " \n", + " \n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "\n", + "
\n", + "\n", + "\n", + "\n" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "" + ] + } + ], + "prompt_number": 14 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file