diff --git a/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb new file mode 100644 index 00000000..c6002ff0 --- /dev/null +++ b/fvmbook/chap3/acousimple/acoustics_simple_waves.ipynb @@ -0,0 +1,109 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:78991e4acc93c32f483391d7a48b3dfe1559a0d2ee23b89d604c56147ff47ace" + }, + "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", + "from clawpack.visclaw import ianimate\n", + "\n", + "import matplotlib\n", + "matplotlib.rcParams.update({'font.size': 22})\n", + "import matplotlib.pyplot as plt" + ], + "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": [ + "ianimate.ianimate(claw,varname='Pressure')" + ], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file 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