{ "metadata": { "name": "", "signature": "sha256:411f2e365ae2cc0d7b6cf94f2aee3a55dfc9154edd06de4a752f7c5cf53f3f43" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Elasticity equations and their Riemann solver in 2D and 3D\n", "by Mauricio J. Del Razo S. 2014\n", "\n", "Here we will give a brief introduction to the elasticity equations, and we will derive the normal Riemann solver for the two dimensional equations in a mapped grid. We will show how to recover the solver for a simple cartesian grid and derive the transverse solver. We will also introduce the normal Riemann solver for a three dimensional mapped grid (not yet implemented). \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. The elasticity equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Newton's second law \u2014the conservation of momentum\u2014 for a solid of mass $m$ is given by \n", "\n", "$$\n", " m\\frac{du}{dt} = F,\n", "$$\n", "\n", "where $u$ is the velocity and $F$ the force. We will combine this with Hooke's law for a spring, which states that the force exerted by the spring is proportional to the distance it has been displaced from its equilibrium position, i.e. $F=-kx$, so we obtain\n", "\n", "$$\n", " m\\frac{d^2x}{dt^2} = -kx,\n", "$$\n", "\n", "where $x$ is the displacement from the spring's resting state. This is the well known harmoni oscillator and the simplest equation of elasticity. \n", "\n", "Now we will consider a three dimensional elastic media with constant density $\\rho=m/V$, and velocity $\\bar{u}=(u,v,w)$ so Newton's law is \n", "\n", "$$\n", " \\rho\\frac{d\\bar{u}}{dt} = \\frac{F}{V},\n", "$$\n", "where $d/dt$ will now correspond to the material derivative. The force felt by the elastic media can be rewritten in terms of the pressure as: $\\frac{F}{V}=-\\nabla p$; however, the pressure is a scalar, so it can only model forces that act in the same direction than the deformation, so we need a more general concept than pressure to model elastic media. This is given by the stress tensor (same units as pressure: force per area),\n", "\n", "$$\n", "\\mathbf{\\sigma} = \\left( \\begin{array}{ccc}\n", "\\sigma_{11} & \\sigma_{12} & \\sigma_{13} \\\\\n", "\\sigma_{21} & \\sigma_{22} & \\sigma_{23} \\\\\n", "\\sigma_{31} & \\sigma_{32} & \\sigma_{33} \\end{array} \\right), \n", "\\hspace{20mm}\n", " \\frac{F}{V} = \\nabla \\cdot \\mathbf{\\sigma}. \n", "$$\n", "\n", "If $\\sigma_{ij}=-p\\delta_{ij}$, we obtain $\\nabla \\cdot \\mathbf{\\sigma} = \\nabla p$ as expected. Combining these two equation, we can actually write the elasticity equations as\n", "\n", "$$\n", " \\rho\\frac{d\\bar{u}}{dt} = \\nabla \\cdot \\mathbf{\\sigma}.\n", "$$\n", "\n", "For small displacements $\\bar{u} \\ll u$, the material derivative can be approximated by the partial derivative yielding the linear elasticity equations,\n", "\n", "$$\n", " \\rho\\frac{\\partial\\bar{u}}{\\partial t} = \\nabla \\cdot \\mathbf{\\sigma}.\n", "$$\n", "\n", "We have three equations and 12 unknown, so we still need some more equations. First of all, in most physical applications, the conservation of angular momentum imples the symmetry of the stress tensor, i.e. $\\sigma_{ij}=\\sigma_{ji}$, which eliminates three variables. The other six equations can be deduce from a generalized version of Hooke's law, the constitutive equation. It must also be linear in order to obtain linear equations. This relation relates the stresses to the strains (deformations) as the Hooke's law relates the forces to the elongation/compression of the spring, see [LeVeque 2002](http://depts.washington.edu/clawpack/book.html) for an intuituve deduction of the following linear constitutive equation. \n", "\n", "$$\n", " \\left[ \\begin{array}{c}\n", "\\sigma_{11} \\\\\n", "\\sigma_{22} \\\\\n", "\\sigma_{33} \\\\\n", "\\sigma_{12} \\\\\n", "\\sigma_{23} \\\\\n", "\\sigma_{13} \\end{array} \\right]\n", "= \\left[ \\begin{array}{cccccc}\n", "\\lambda + 2\\mu & \\lambda & \\lambda & 0 & 0 & 0 \\\\\n", "\\lambda & \\lambda + 2\\mu & \\lambda & 0 & 0 & 0 \\\\\n", "\\lambda & \\lambda & \\lambda + 2\\mu & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 2\\mu & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 2\\mu & 0\\\\\n", "0 & 0 & 0 & 0 & 0 & 2\\mu \\\\\n", "\\end{array} \\right]\n", "\\left[ \\begin{array}{c}\n", "\\epsilon_{11} \\\\\n", "\\epsilon_{22} \\\\\n", "\\epsilon_{33} \\\\\n", "\\epsilon_{12} \\\\\n", "\\epsilon_{23} \\\\\n", "\\epsilon_{13} \\end{array} \\right],\n", "$$\n", "where the strain tensor is expressed in term of the displacement vector $\\delta(x,y,z,t)$ as $\\mathbf{\\epsilon}=\\frac{1}{2}\\left[\\nabla \\bar{\\delta} + \\left(\\nabla \\bar{\\delta}\\right)^T\\right]$ and $\\lambda$ and $\\mu$ are the Lam\u00e9 parameters (analogous to the $x$ and $k$ in Hooke's law respectively). Differentiating by $t$ this stress-strain relation we obtain the other six equations to complete the system, \n", "\n", "$$\n", "\\left[ \\begin{array}{c}\n", "\\sigma_{11} \\\\\n", "\\sigma_{22} \\\\\n", "\\sigma_{33} \\\\\n", "\\sigma_{12} \\\\\n", "\\sigma_{23} \\\\\n", "\\sigma_{13} \\end{array} \\right]_t\n", "= \\left[ \\begin{array}{cccccc}\n", "\\lambda + 2\\mu & \\lambda & \\lambda & 0 & 0 & 0 \\\\\n", "\\lambda & \\lambda + 2\\mu & \\lambda & 0 & 0 & 0 \\\\\n", "\\lambda & \\lambda & \\lambda + 2\\mu & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 2\\mu & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 2\\mu & 0\\\\\n", "0 & 0 & 0 & 0 & 0 & 2\\mu \\\\\n", "\\end{array} \\right]\n", "\\left[ \\begin{array}{c}\n", "u_x \\\\\n", "v_y \\\\\n", "w_z \\\\\n", "(u_y + v_x)/2 \\\\\n", "(v_z + w_y)/2 \\\\\n", "(u_z + w_x)/2 \\end{array} \\right].\n", "$$\n", "\n", "We can rewrite the nine equations in the form \n", "\n", "$$\n", " \\bar{q}_t + \\mathbf{A}\\bar{q}_x + \\mathbf{B}\\bar{q}_x + \\mathbf{C}\\bar{q}_z = 0,\n", "$$\n", "\n", "with\n", "\n", "$$\n", "\\bar{q} = \\left[ \\begin{array}{c}\n", "\\sigma_{11} \\\\\n", "\\sigma_{22} \\\\\n", "\\sigma_{33} \\\\\n", "\\sigma_{12} \\\\\n", "\\sigma_{23} \\\\\n", "\\sigma_{13} \\\\\n", "u \\\\\n", "v \\\\\n", "w \\end{array} \\right]\n", "\\hspace{10mm}\n", "\\mathbf{A} = \n", "\\left[ \\begin{array}{ccccccccc}\n", "0 & 0 & 0 & 0 & 0 & 0 & -(\\lambda + 2\\mu) & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -\\lambda & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -\\lambda & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\mu & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\mu \\\\\n", "-1/\\rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 \\\\\n", "\\end{array} \\right]\n", "$$\n", "\n", "$$\n", "\\mathbf{B} = \n", "\\left[ \\begin{array}{ccccccccc}\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\lambda & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -(\\lambda + 2\\mu) & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\lambda & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -\\mu & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\mu \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & -1/\\rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 & 0 \\\\\n", "\\end{array} \\right]\n", "\\hspace{5mm}\n", "\\mathbf{C} = \n", "\\left[ \\begin{array}{ccccccccc}\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\lambda \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\lambda \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -(\\lambda + 2\\mu) \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\mu & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & 0 & -\\mu & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 \\\\\n", "0 & 0 & 0 & 0 & -1/\\rho & 0 & 0 & 0 & 0 \\\\\n", "0 & 0 & -1/\\rho & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n", "\\end{array} \\right]\n", "$$\n", "\n", "The three dimensional elasticity equations written in this form are ideal to solve the Riemann problem and implement numerically.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. The Riemann problem for 2D elasticity equations in a mapped grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equations obtained in the previous section are easily simplified to its two dimensional version. The stress tensor will only have 3 entries: $\\sigma_{11},\\sigma_{22}$ and $\\sigma_{12}$, and the velocity will have two: $u$ and $v$. The normal Riemann problem will solve $\\bar{q}_t + \\mathbf{A}\\bar{q}_x=0$ when sweeping in the $x$ direction and $ \\bar{q}_t + \\mathbf{B}\\bar{q}_y=0$ when sweeping in the $y$ direction. We can generalize this for a general quadrilateral mapped grid. The Riemann problem in the normal direction of a grid edge is $\\bar{q}_t + \\mathbf{A_n}\\bar{q}_n=0$, where $\\bar{q}_n$ is the derivative in the normal to the edge direction, $\\mathbf{A_n}=n_x \\mathbf{A} + n_y\\mathbf{B}$, and $\\hat{n}=(n_x,n_y)$ is the normal to the edge where the Riemann problem is beig solved. The jacobian matrix is,\n", "\n", "$$\n", "\\mathbf{A_n} = -\n", "\\left[ \\begin{array}{ccccc}\n", "0 & 0 & 0 & n_x(\\lambda + 2\\mu) & n_y\\lambda \\\\\n", "0 & 0 & 0 & n_x\\lambda & n_y(\\lambda + 2\\mu) \\\\\n", "0 & 0 & 0 & n_y\\mu & n_x\\mu \\\\\n", "n_x/\\rho & 0 & n_y/\\rho & 0 & 0 \\\\\n", "0 & n_y/\\rho & n_x/\\rho & 0 & 0 \\\\\n", "\\end{array} \\right].\n", "$$\n", "\n", "Using $\\hat{n}=(1,0)$ or $\\hat{n}=(0,1)$, we recover $\\mathbf{A}$ and $\\mathbf{B}$ respectively. In order to solve the Riemann problem, we need the eigenvalues and eigenvectors of this matrix. We can obtain them symblically:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Initialize sympy for symbolic computations\n", "import sympy as sy\n", "# Mathjax printing with simpy\n", "sy.init_printing(use_latex='mathjax')\n", "\n", "# Declare symbolic variables\n", "nx, ny, rho, lamda, mu = sy.symbols('n_x n_y rho lamda mu', real=True)\n", "R = sy.Matrix(5, 5, lambda i,j: 0)\n", "Rval = sy.Matrix(5, 1, lambda i,j: 0)\n", "\n", "# Define symbolic matrix\n", "An = -sy.Matrix([[0,0,0,nx*(lamda+2*mu),ny*lamda],\n", " [0,0,0,nx*lamda,ny*(lamda+2*mu)],\n", " [0,0,0,ny*mu,nx*mu],\n", " [nx/rho,0,ny/rho,0,0],\n", " [0,ny/rho,nx/rho,0,0]])\n", "\n", "# Calculate eigenvalues and eigenvectors of matrix An\n", "# eigvec[i][0] => ith eigenvalue\n", "# eigvec[i][2] => ith eigenvector\n", "eig_vec = An.eigenvects()\n", "\n", "\n", "# Simplify eigenvalues and eigenvectors and save in symbolic arrays\n", "for i in range(5):\n", " Rval[i] = sy.simplify(eig_vec[i][0].subs(nx**2,1-ny**2).cancel())\n", " for j in range(5):\n", " R[j,i] = sy.simplify(eig_vec[i][2][0][j].subs(nx**2,1-ny**2).cancel())\n", " R[j,i] = sy.simplify((ny*R[j,i]).subs(nx**2,1-ny**2).cancel())\n", " R[j,i] = sy.simplify(R[j,i].subs(nx**2+ny**2,1).cancel())" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "# Show eigenvalue and corresponding eigenvectors in a widget\n", "\n", "from IPython.html.widgets import interact, interactive\n", "from IPython.display import display, HTML, Latex\n", "\n", "# Function to print the eigenvalue and eigenvector\n", "def print_eigen(Eigenvalue):\n", " #print \"Eigenvalue, Eigenvector= \"\n", " #display((Rval[Eigenvalue], R[:,Eigenvalue]));\n", " display(Latex('Eigenvalue =' + '$' + sy.latex(Rval[Eigenvalue]) + '\\hspace{20mm} $' + \n", " 'Eigenvector =$' + sy.latex(R[:,Eigenvalue]) + '$'))\n", " \n", "# Interactive widget to display eigenvalue and eigenvector\n", "interact(print_eigen, Eigenvalue=(0,4));" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "Eigenvalue =$\\sqrt{\\frac{\\mu}{\\rho}}\\hspace{20mm} $Eigenvector =$\\left[\\begin{smallmatrix}2 n_{y}^{2} \\rho \\sqrt{\\frac{\\mu}{\\rho}}\\\\- 2 n_{y}^{2} \\rho \\sqrt{\\frac{\\mu}{\\rho}}\\\\\\frac{n_{x} n_{y} \\rho \\sqrt{\\frac{\\mu}{\\rho}} \\left(- 2 n_{y}^{2} + 1\\right)}{n_{y}^{2} -1}\\\\\\frac{n_{x} n_{y}^{2}}{n_{y}^{2} -1}\\\\n_{y}\\end{smallmatrix}\\right]$" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The eigenvalues and eigenvectors obtained are the same than those in [LeVeque 2002](http://depts.washington.edu/clawpack/book.html), except for sustitutions of $n_x^2 + n_y^2=1$ and constant factors in the eigenvectors. Note the wave with eigevaule zero doesn't causes any propagation, and it will not be relevant in solving the Riemann problem. The next two of this eigenvalues correspond to two shear waves (S-waves) with the shear speed as eigenvalues \n", "\n", "$$\n", " c_s = \\mp \\sqrt{\\mu/\\rho},\n", "$$\n", "\n", "and the next two correspond to two pressure wave (P-waves) with the usual sound speed as eigenvalues\n", "\n", "$$\n", " c_p = \\mp \\sqrt{\\frac{\\lambda + 2\\mu}{\\rho}}.\n", "$$\n", "\n", "The eigenvectors can be easily simplified multiplying each by some constant and rewriting them in terms of $c_s$ and $c_p$. We also rearrange them into the matrix of the five eigenvectors $\\mathbf{R}=[\\bar{r_1}, \\bar{r_2}, \\bar{r_3}, \\bar{r_4}, \\bar{r_5}]$,\n", "\n", "$$\n", "\\mathbf{R} = \n", "\\left[ \\begin{array}{ccccc}\n", "\\lambda_L + 2\\mu_L n_x^2 & \\lambda_R + 2\\mu_R n_x^2 & -2\\mu_L n_x n_y & -2\\mu_R n_x n_y & n_y^2 \\\\\n", "\\lambda_L + 2\\mu_L n_y^2 & \\lambda_R + 2\\mu_R n_y^2 & 2\\mu_L n_x n_y & 2\\mu_R n_x n_y & n_x^2 \\\\\n", "2\\mu_L n_x n_y & 2\\mu_R n_x n_y & \\mu_L [n_x^2 - n_y^2] & \\mu_R [n_x^2 - n_y^2] & -n_x n_y \\\\\n", "n_x c_{pL} & -n_x c_{pR} & -n_y c_{sL} & n_y c_{sR} & 0 \\\\\n", "n_y c_{pL} & -n_y c_{pR} & n_x c_{sL} & -n_x c_{sR} & 0 \\\\\n", "\\end{array} \\right],\n", "$$\n", "\n", "where each column correspond to the eigenvectors with the respective eigenvalues,\n", "\n", "$$\n", "s_1 = -c_{pL}, \\ \\ s_2=c_{pR}, \\ \\ s_3 = -c_{sL}, \\ \\ s_4 = c_{sR}, \\ \\ s_5 = 0,\n", "$$\n", "\n", "and the subindex for left and right, $L$ (and $R$), was added to the parameters of the eigenvector with negative (and positive) value respectively. This denote different values of the parameters on the left and right side of the discontinuity where the waves will be propagating correspondingly. If the parameters $\\lambda$ and $\\mu$ are different across the discontinuity of the Riemann problem; we could think there is an in an interface between two materials. This will also allow the Riemann solver to solve variable coefficients problems for inhomogeneous materials too. \n", "\n", "In order to solve the Riemann problem, we need to solve for $\\bar{\\alpha}$,\n", "\n", "$$\n", "\\mathbf{R} \\bar{\\alpha} = \\Delta \\bar{q},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\delta \\bar{Q} = [\\Delta \\sigma_{11}; \\Delta \\sigma_{22}; \\Delta \\sigma_{12}; \\Delta u; \\Delta v]$ is the jump of $\\bar{q}$ across the discontinuity. We could invert $\\mathbf{R}$ numerically to obtain $\\bar{\\alpha}$, or we could try to use sympy again. In the case of linear elasticity equation, it can be done analitically, and it yields\n", "\n", "$$\n", "\\alpha_1 = \\frac{c_{pR}\\left[\n", " \\Delta \\sigma_{11}n_x^2 + \n", " \\Delta \\sigma_{22}n_y^2 +\n", " 2n_x n_y \\Delta \\sigma_{12}\\right] + \n", " (\\lambda_R +2 \\mu_R)(n_x \\Delta u +\n", " n_y \\Delta v)]}{\\Delta_P}, \\\\\n", "\\alpha_2 = \\frac{c_{pL}\\left[\n", " \\Delta \\sigma_{11}n_x^2 + \n", " \\Delta \\sigma_{22}n_y^2 +\n", " 2n_x n_y \\Delta \\sigma_{12}\\right] - \n", " (\\lambda_L +2 \\mu_L)(n_x \\Delta u +\n", " n_y \\Delta v)}{\\Delta_P}, \\\\\n", "\\alpha_3 = \\frac{c_{sR}\\left[\n", " \\Delta \\sigma_{12}(n_x^2 - n_y^2) + \n", " n_x n_y(\\Delta \\sigma_{22} - \\Delta \\sigma_{11})\\right] +\n", " \\mu_R(n_x \\Delta v - n_y \\Delta u))}{\\Delta_S}, \\\\\n", "\\alpha_4 = \\frac{c_{sL}\\left[\n", " \\Delta \\sigma_{12}(n_x^2 - n_y^2) + \n", " n_x n_y(\\Delta \\sigma_{22} - \\Delta \\sigma_{11})\\right] -\n", " \\mu_L(n_x \\Delta v - n_y \\Delta u))}{\\Delta_S}, \\\\\n", "$$\n", "with\n", "$$\n", "\\Delta_P = c_{pR}(\\lambda_L + 2\\mu_L) + c_{pL}(\\lambda_R + 2\\mu_R), \\\\\n", "\\Delta_S = \\mu_L c_{sR} + \\mu_R c_{sL}.\n", "$$\n", "and corresponding speeds (eigenvalues):\n", "$$\n", "s_1 = -c_{pL} \\ \\ s_2 = c_{pR} \\ \\ s_3 = -c_{sL} \\ \\ s_4 = c_{sR} \\ \\ s_5 = 0.\n", "$$\n", "Note $\\alpha_5$ is not required, since its corresponding zero eigenvalue doesn't allow wave propagation. \n", "This information is enough to calculate the solution, since the vector $\\bar{\\alpha_i}$ represents the jump across the wave $i$ that moves along the eigenvector $\\bar{r_i}$ with speed $s_i$. \n", "Note when $\\hat{n} = (1,0)$ or $\\hat{n} = (0,1)$, we recover the normal Riemann solver for 2D Elasticity in a Cartesian grid." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Clawpack implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We solved the Riemann problem for a mapped grid, so we want it to solve a more complicated system than just a Riemann problem. Clawpack has all the framework required to solve hyperbolic equations by solving Riemann problems on each grid cell. If we can provide the Riemann solver to Clawpack, it can do all the hard work for us. In order to implement the solution into Clawpack, we need to rewrite it in its wave formulation, where we update the value at cell $\\bar{Q}_{i,j}^n$, with $i,j$ the indexes in space and $n$ in time following,\n", "\n", "$$\n", "\\bar{Q}_{i,j}^{n+1} = \\bar{Q}_{i,j}^{n} - \\frac{\\Delta t}{\\Delta x}\n", "\\left[A^+\\Delta Q_{i-1/2,j} + A^-\\Delta Q_{i+1/2,j}\\right] + \\text{High order corrections & limiters},\n", "$$ \n", "\n", "when sweeping in the $i$ direction of the computational domain. Note this should be equivalent to sweeping in the normal direction of the physical domain when using a mapped grid, that's why it's important solve the Riemann problem normal to the edge between the grid cells. The same process is analogous in the $j$ direction. Also note this equation requires output of two Riemann problems: the one in edge $i-1/2$ and the one in $i+1/2$. The output of the Riemann solver routine for Clawpack should be given by the waves speed and wave fluctuactions. The wave speed we already know from the eigenvalues in the previous section. In order to calculate the wave fluctuations $A^{\\pm}\\Delta Q_{i-1/2,j}$ from the left edge, we first need the waves given by\n", "\n", "$$\n", "W_i = \\alpha_i \\bar{r_i}\n", "$$\n", "\n", "with $i=1,2,3,4$ and $\\bar{r_i}$ the corresponding eigenvectors calculated previously. Solving the Riemann problem in the right edge is completely analogous. The positive and negative wave fluctuations $A^\\pm \\Delta Q_{i- 1/2,j}$ will be given in terms of the speeds(eigenvalues) and the waves,\n", "\n", "$$\n", "A^+ \\Delta Q_{i-1/2,j} = s_2W_2 + s_4W_4 \\\\\n", "A^- \\Delta Q_{i-1/2,j} = s_1W_1 + s_3W_3,\n", "$$\n", "When employing the mapped grid is important to remember to scale the wave speeds $s$ in the wave fluctuations by the factor $\\gamma_x = dy_{phy}/dy_{com}$ if sweeping in the $x$ direction or $\\gamma_y = dx_{phy}/dx_{com}$ when sweeping in the $y$ direction, where $dx_{phy}$ denotes mesh grid spacing in the physical domain and $dx_{com}$ in the computational one. The fluctuation in the mapped grid should then be,\n", "\n", "$$\n", "A^+ \\Delta Q_{i-1/2,j}^{map} = \\gamma c_{pR}W_2 + \\gamma c_{sR}W_4 \\\\\n", "A^- \\Delta Q_{i-1/2,j}^{map} = -\\gamma c_{pL}W_1 -\\gamma c_{sL}W_3.\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Implementing a sample solution into python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can implement a sample solution of the Riemann problem we just deduced into python. We will give an outline of the algorithm followed by the code:\n", "
    \n", "
  1. Define piecewise constant initial condition with the left and right states\n", "$$\n", "q_L = [\\sigma_{11L}, \\sigma_{22L}, \\sigma_{12L}, u_L, v_L ] \\ \\ \\text{with} \\ \\ \\mu_L \\ \\ \\text{and} \\ \\ \\lambda_L, \\\\\n", "q_R = [\\sigma_{11R}, \\sigma_{22R}, \\sigma_{12R}, u_R, v_R ] \\ \\ \\text{with} \\ \\ \\mu_R \\ \\ \\text{and} \\ \\ \\lambda_R. \n", "$$\n", "
  2. \n", "
  3. Set unitary normal direction, use $n_x=1$ and $n_y=0$ to obatin the normal solver in a cartesian grid
  4. \n", "
  5. Calculate the wave speeds given by the eigenvalues and the corresponding eigenvectors $\\mathbf{A_n}$ for the selected parametrs.\n", "
  6. Compute \\bar{\\alpha} that are the weights of the jumps in the direction of each of the eigenvectors
  7. \n", "
  8. Compute the middle states $q_{mL}$ and $q_{sL}$ as shown in the figure, note that $-c_p \\le -c_s$, so the jump from $q_L$ to $q_{mL}$ will be in the 1-wave (speed $-c_{pL})$, and the jump from $q_{mL}$ to $q_{sL}$ will be across the 3-wave (speed $-c_{sL}$)
  9. \n", "
  10. As we didn't calculate the jump across the 5-wave (speed 0), we can obtain $q_{mR}$ and $q_{sR}$ starting from the the right state $q_R$. The jump from $q_{mR}$ to $q_{R}$ is across the 2-wave (speed $c_{pR}$), and the jump from $q_{sR}$ to $q_{mR}$ is across the 4-wave (speed $c_{sR}$). These yields the following formulas for the middle states,\n", "\n", "$$\n", "q_{mL} = q_{L} + \\alpha_1 \\bar{r_1} \\\\\n", "q_{sL} = q_{mL} + \\alpha_3 \\bar{r_3} \\\\\n", "q_{mR} = q_{R} - \\alpha_2 \\bar{r_2} \\\\\n", "q_{sR} = q_{mR} - \\alpha_4\\bar{r_4}.\n", "$$\n", "
  11. \n", "\n", "Here is the figure showing the six different states of the Riemann problem and the corresponding waves:\n", "\n", "\n", "\n", "Now we implement this solution into a python script and plot the solutions interactively\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#import numpy\n", "import numpy as np\n", "\n", "# Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v])\n", "qL = np.array([0.0, 0.0, 0.0, 4.0, 0.0])\n", "qR = np.array([0.0, 0.0, 0.0, 0.0, 0.0])\n", "dq = np.array([qR[0]-qL[0], qR[1]-qL[1], qR[2]-qL[2], \n", " qR[3]-qL[3], qR[4]-qL[4]])\n", "\n", "# Set Lame parameters (\\mu, \\lambda) and density (left and right)\n", "muL = 1; muR = 4\n", "lamL = 1; lamR = 4\n", "rhoL = 1; rhoR = 2\n", "\n", "# Set normal (unitary) \n", "nx = 1.0/np.sqrt(2); \n", "ny = 1.0/np.sqrt(2)\n", "nx2 = nx*nx\n", "ny2 = ny*ny\n", "nxy = nx*ny\n", "\n", "# Define P and S speeds (left and right)\n", "bulkL = lamL + 2*muL\n", "bulkR = lamR + 2*muR\n", "cpL = np.sqrt(bulkL/rhoL)\n", "cpR = np.sqrt(bulkR/rhoR)\n", "csL = np.sqrt(muL/rhoL)\n", "csR = np.sqrt(muR/rhoR)\n", "\n", "# Define Eigenvalues \n", "s = np.array([-cpL, cpR, -csL, csR, 0])\n", "\n", "# Define the 4 eigenvectors (from columns of Matrix R)\n", "r1 = np.array([lamL + 2*muL*nx2, lamL + 2*muL*ny2,\n", " 2*muL*nxy, nx*cpL, ny*cpL])\n", "r2 = np.array([lamR + 2*muR*nx2, lamR + 2*muR*ny2, \n", " 2*muR*nxy, -nx*cpR, -ny*cpR])\n", "r3 = np.array([-2*muL*nxy, 2*muL*nxy, muL*(nx2 - ny2),\n", " -ny*csL, nx*csL])\n", "r4 = np.array([-2*muR*nxy, 2*muR*nxy, muR*(nx2 - ny2),\n", " ny*csR, -nx*csR])\n", "r5 = np.array([ny2, nx2, -nxy, 0, 0])\n", "\n", "# Compute the 4 alphas\n", "detP = cpR*bulkL + cpL*bulkR\n", "detS = csR*muL + csL*muR\n", "\n", "al1 = (cpR*(dq[0]*nx2 + dq[1]*ny2 + 2*nxy*dq[2]) +\n", " bulkR*(nx*dq[3] + ny*dq[4]))/detP\n", "al2 = (cpL*(dq[0]*nx2 + dq[1]*ny2 + 2*nxy*dq[2]) -\n", " bulkL*(nx*dq[3] + ny*dq[4]))/detP\n", "al3 = (csR*(dq[2]*(nx2 - ny2) + nxy*(dq[1] - dq[0])) + \n", " muR*(nx*dq[4] - ny*dq[3]))/detS\n", "al4 = (csL*(dq[2]*(nx2 - ny2) + nxy*(dq[1] - dq[0])) - \n", " muL*(nx*dq[4] - ny*dq[3]))/detS\n", "\n", "# Compute four middle states qmL, qsL, qsR, qmR --note order of characteristics is s3,s1,s5,s4,s2\n", "# Starting from left side\n", "qmL = qL + al1*r1 \n", "qsL = qmL + al3*r3\n", "# Starting from right side\n", "qmR = qR - al2*r2\n", "qsR = qmR - al4*r4\n", "\n", "## Calculate alpha 5 to test correctness\n", "## Note alpha 5 (al[4]) is not neccesary to calculate when using clawpack, since it doesn't propagate any wave\n", "#p1 = cpR*(lamL*(nx2-ny2)*(dq[0]-dq[1]) - 2*muL*(dq[0]*ny2 + dq[1]*nx2) + 4*nxy*dq[2]*(lamL + muL))\n", "#p2 = cpL*(lamR*(nx2-ny2)*(dq[0]-dq[1]) - 2*muR*(dq[0]*ny2 + dq[1]*nx2) + 4*nxy*dq[2]*(lamR + muR))\n", "#p3 = (nx*dq[3] + ny*dq[4])*(2*lamL*muR - 2*lamR*muL)\n", "#al5 = -(p1 + p2 + p3)/detP\n", "#qsR_test = qsL + al5*r5 # qsR_test must be equal to qsR\n", "#print(qsR_test-qsR) # Must be zero or almost zero in all entries\n", " \n", "# Compute waves characteristics for plotting\n", "x = np.linspace(-5,5,50)\n", "Wc = np.zeros((len(s),len(x)))\n", "Wc[0][:] = x/s[0]\n", "Wc[1][:] = x/s[1]\n", "Wc[2][:] = x/s[2]\n", "Wc[3][:] = x/s[3]\n", "# To plot vertical dashed line\n", "Wc[4][:] = x/(s[4]+0.000000000001)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 62 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the solution for different times and for each of the variables in an interactive widget:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Required for plotting\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pylab import *\n", "from IPython.html.widgets import interact\n", "from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget\n", "\n", "# Plot Riemann solution\n", "def plot_Riemann(time,variable):\n", " \n", " vardict = {'sigma 11':0, 'sigma 22':1, 'sigma 12':2, 'normal velocity u':3, 'transverse velocity v':4}\n", " qvar = vardict.get(variable)\n", " \n", " fig = plt.figure(figsize=(10, 4))\n", " \n", " # Create subplot for (x,t) plane\n", " tmax = 5\n", " ax = fig.add_subplot(1,2,1)\n", " ax.set_xlabel('x')\n", " ax.set_ylabel('time')\n", " ax.axis([min(x),max(x),0,tmax])\n", " \n", " # Plot characteristic lines\n", " # P-waves in red\n", " ax.plot(x,Wc[0][:], '-r', linewidth=2)\n", " ax.plot(x,Wc[1][:], '-r', linewidth=2)\n", " # S-waves in blue\n", " ax.plot(x,Wc[2][:], '-b', linewidth=2)\n", " ax.plot(x,Wc[3][:], '-b', linewidth=2)\n", " # Non-propagating wave in dashed\n", " ax.plot(x,Wc[4][:], '--k', linewidth=2)\n", " # Plot time-line in (x,t) plane\n", " ax.plot(x, 0*x+time, 'k', linewidth=3)\n", "\n", " # Create subplot for Riemann solution\n", " ax2 = fig.add_subplot(1,2,2)\n", " ax2.set_xlabel('x')\n", " ax2.set_ylabel('q [' + str(qvar+1) + ']')\n", " ax2.axis([min(x),max(x),-max(x),max(x)])\n", " \n", " # Create Riemann solution vector and plot\n", " sol = np.zeros(len(x))\n", " for i in range(len(x)):\n", " if x[i] < -cpL*time:\n", " sol[i] = qL[qvar]\n", " elif x[i] < -csL*time:\n", " sol[i] = qmL[qvar]\n", " elif x[i] < 0:\n", " sol[i] = qsL[qvar]\n", " elif x[i] < csR*time:\n", " sol[i] = qsR[qvar]\n", " elif x[i] < cpR*time:\n", " sol[i] = qmR[qvar]\n", " else:\n", " sol[i] = qR[qvar]\n", " ax2.plot(x,sol, 'k', linewidth = 3)\n", " ax2.fill_between(x,-20, sol, facecolor='blue', alpha=0.2)\n", " return fig\n", " \n", "# Create interactive widget to visualize solution \n", "#interact(plot_Riemann, time=(0,5,0.1), qvar={'sigma11':0, 'sigma22':1, 'sigma12':2, \n", " #'normal velocity u':3, 'transverse velocity v':4});\n", "StaticInteract(plot_Riemann, time=RangeWidget(0,5,0.25),\n", " variable=DropDownWidget(['sigma 11','sigma 22','sigma 12',\n", " 'normal velocity u','transverse velocity v'])) " ], "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", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", "
    \n", " \n", " time: \n", "
    \n", "variable: \n", "
    \n", " " ], "metadata": {}, "output_type": "pyout", "prompt_number": 63, "text": [ "" ] } ], "prompt_number": 63 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The characteristics of the P-Waves and S-waves are colored red and blue respectively. It is possible to observe the propagation of the P-Waves and S-Waves by moving the time slider for the different variables. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. The transverse Riemann problem for 2D elasticity equations in a mapped grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to obtain second order accuracy, we will need to solve a transverse Riemann solver in addition to the normal one. This will require to split the normal wave fluctuactions $A^\\pm \\Delta Q_{i- 1/2,j}$ at edge $i-1/2$ into transverse wave fluctuactions $B^\\pm A^+\\Delta Q_{i- 1/2,j}$ and $B^\\pm A^-\\Delta Q_{i- 1/2,j}$; the former descomposition is shown in the following figure and explained in detail on [LeVeque 2002](http://depts.washington.edu/clawpack/book.html).\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will only focuse on $B^\\pm A^+\\Delta Q_{i- 1/2,j}$, since the other case is analogous. We begin by decomposing the normal wave fluctuations into a linear combination of transverse waves, which will follow the same equations and jacobian as before, so we obtain\n", "\n", "$$\n", "A^+\\Delta Q_{i- 1/2,j} = \\mathbf{R} \\bar{\\beta} = \\beta_1 \\bar{r_1} + \\beta_2 \\bar{r_2} + \\beta_3 \\bar{r_3} + \\beta_4 \\bar{r_4}.\n", "$$\n", "\n", "In the context of mapped grids, note the eigenvectors $\\bar{r_i}$ depend on the normal $\\hat{n}$=($n_x,n_y$), so if the left edge normal was used for the normal Riemann solver, we need to use the bottom edge normal to calculate $\\bar{r_i}$. This normal might also be different for the up-going and down-going fluctuation. If we are on a cartesian grid, both up-going and down-going normal will be equal and if $\\hat{n}$=($1,0$) was used for the normal solver, the transverse solver will use the bottom normal $\\hat{n}$=$($0,1$) or viceversa. \n", "\n", "The eigenvectors will also depend on the parameters at each cell, and now we have three cells involved: the bottom one, the middle one and the upper one. We will denote the wave speeds and parameters with the subindexes $B,M$ and $U$ to denote bottom, middle and upper respectively. \n", "\n", "We solve the system in a very similar manner to the system for $\\bar{\\alpha}$. When we add the right input parameters for the eigenvectors, we obtain $\\bar{\\beta}$ is given by,\n", "\n", "$$\n", "\\beta_1 = \\frac{c_{pM}\\left[\n", " A^+_1 n_{xM}^2 + \n", " A^+_2 n_{yM}^2 +\n", " 2n_{xM} n_{yM} A^+_3\\right] + \n", " (\\lambda_M +2 \\mu_M)(n_{xM} A^+_4 +\n", " n_{yM} A^+_5)]}{\\Delta_P^-}, \\\\\n", "\\beta_2 = \\frac{c_{pM}\\left[\n", " A^+_1 n_{xU}^2 + \n", " A^+_2 n_{yU}^2 +\n", " 2n_{xU} n_{yU} A^+_3 \\right] - \n", " (\\lambda_M +2 \\mu_M)(n_{xU} A^+_4 +\n", " n_{yU} A^+_5)}{\\Delta_P^+}, \\\\\n", "\\beta_3 = \\frac{c_{sM}\\left[\n", " A^+_3 (n_{xB}^2 - n_{yB}^2) + \n", " n_{xB} n_{yB}(A^+_2 - A^+_1) \\right] +\n", " \\mu_M(n_{xB} A^+_5 - n_{yB} A^+_4))}{\\Delta_S^-}, \\\\\n", "\\beta_4 = \\frac{c_{sM}\\left[\n", " A^+_3 (n_{xU}^2 - n_{yU}^2) + \n", " n_{xU} n_{yU}(A^+_2 - A^+_1) \\right] -\n", " \\mu_M(n_{xU} A^+_5 - n_{yU} A^+_4 ))}{\\Delta_S^+}, \\\\\n", "$$\n", "with\n", "$$\n", "\\Delta_P^- = c_{pM}(\\lambda_B + 2\\mu_B) + c_{pB}(\\lambda_M + 2\\mu_M), \\\\\n", "\\Delta_P^+ = c_{pU}(\\lambda_M + 2\\mu_M) + c_{pM}(\\lambda_U + 2\\mu_U), \\\\\n", "\\Delta_S^- = \\mu_B c_{sM} + \\mu_M c_{sB}, \\\\\n", "\\Delta_S^+ = \\mu_M c_{sU} + \\mu_U c_{sM},\n", "$$\n", "\n", "where $A^+\\Delta Q_{i- 1/2,j}=[A^+_1, A^+_2,A^+_3,A^+_4,A^+_5]$, ($n_{xM},n_{yM}$) correspond to the bottom normal of the middle cell, ($n_{xU},n_{yU}$) correspond to the bottom normal of the upper cell and the speeds were calculated with the same formulas as before with the $\\lambda,\\mu$ and $\\rho$ parameters that correspond to each cell.\n", "\n", "The speeds (eigenvalues) are the same as before, but in the corresponding cell,\n", "$$\n", "s_1^T = -c_{pB} \\ \\ s_2^T = c_{pU} \\ \\ s_3^T = -c_{sB} \\ \\ s_4^T = c_{sU} \\ \\ s_5^T = 0.\n", "$$\n", "\n", "The transverse waves are given by $W_i^T = \\beta_i \\bar{r_i}$, so the fluctuations are given by,\n", "\n", "$$\n", "B^+A^+ \\Delta Q_{i-1/2,j} = s_2^T W_2^T + s_4^T W_4^T \\\\\n", "B^-A^+ \\Delta Q_{i-1/2,j} = s_1^T W_1^T + s_3^T W_3^T,\n", "$$\n", "\n", "Once again, when employing the mapped grid is important to remember to scale the wave speeds $s$ in the wave fluctuations by the factor $\\gamma$\n", "\n", "$$\n", "B^+A^+ \\Delta Q_{i-1/2,j}^{map} = \\gamma_U c_{pU} W_2^T + \\gamma_U c_{sU}W_4^T \\\\\n", "B^-A^+ \\Delta Q_{i-1/2,j}^{map} = -\\gamma_M c_{pB}W_1^T -\\gamma_M c_{sB}W_3^T,\n", "$$\n", "where $\\gamma_U$ correspond to the scaling ratio of the bottom edge of the upper cell and $\\gamma_M$ to the scaling ratio of the bottom edge of the middle cell. The same process can be repeated for the left normal fluctuation $A^-\\Delta Q_{i-1/2,j}$.\n" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }