Skip to content

Commit

Permalink
examples: more explaination and addressed comments in curvilinear not…
Browse files Browse the repository at this point in the history 
…ebook
  • Loading branch information
@EdCaunt
EdCaunt committed Jun 21, 2023
1 parent 9d70670 commit 64e1d7d
Showing 1 changed file with 129 additions and 66 deletions.
195 changes: 129 additions & 66 deletions examples/seismic/tutorials/16_variable_z_topography.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,36 @@
"\n",
"In this notebook, we will create a simple implementation of a cosine hill topography based on coordinate transforms. This approach to irregular topography is commonly used in industry and academia, and is based upon creating some coordinate transform which maps the curved free-surface to a horizontal, grid-coincident plane in the iteration space. The most straightforward way of achieving this is by stretching the grid in the z axis, often referred to as \"variable-z\". This yields a transformed wave equation which can be solved with an explict FD solver.\n",
"\n",
"The isotropic acoustic wave equation in 2D is given as\n",
"\n",
"$\\frac{\\partial^2p}{\\partial t^2} = c^2\\left(\\frac{\\partial^2}{\\partial x^2}+\\frac{\\partial^2}{\\partial z^2}\\right)$,\n",
"\n",
"where $p$ is pressure and $c$ is wavespeed. To transform this equation from Cartesian dimensions $x$ and $z$ into a domain based on curvilinear dimensions $x'$ and $\\zeta$, spatial derivatives in this equation must be replaced with their transformed equivalents. Given that only the vertical dimension is stretched, these transformed derivatives are of the form\n",
"\n",
"$\\frac{\\partial^2p}{\\partial x^2} = \\frac{\\partial^2p}{\\partial x'^2} + \\frac{\\partial^2p}{\\partial x'\\partial \\zeta}\\frac{\\partial\\zeta}{\\partial x} + \\left(\\frac{\\partial^2 p}{\\partial x' \\partial\\zeta}+\\frac{\\partial^2 p}{\\partial\\zeta^2}\\frac{\\partial\\zeta}{\\partial x}\\right)\\frac{\\partial\\zeta}{\\partial x}+\\frac{\\partial p}{\\partial\\zeta}\\frac{\\partial^2\\zeta}{\\partial x^2}$,\n",
"\n",
"and\n",
"\n",
"$\\frac{\\partial^2p}{\\partial z^2} = \\frac{\\partial^2 p}{\\partial\\zeta^2}\\left(\\frac{\\partial\\zeta}{\\partial z}\\right)^2$.\n",
"\n",
"Note that we now have several derivatives pertaining to the transformation, $\\frac{\\partial\\zeta}{\\partial x}$ for example, which will require evaluation.\n",
"\n",
"## The topography\n",
"The surface topography for this particular example corresponds to\n",
"\n",
"$z_M = z_b + \\frac{h\\left(1 - cos\\left(\\frac{2\\pi(x-x_m)}{w}\\right)\\right)}{2}$,\n",
"\n",
"where $z_b$ is the baseline from which the topography rises, $h$ is the height of the topography, $x$ is the horizontal dimension with $x_m$ being its lower bound, and $w$ is the width of the domain. This surface represents the maximum extent of $z$, which will be mapped to the maximum extent of the vertical iteration dimension $\\zeta$. The minimum $\\zeta$, $\\zeta_m$ will be mapped to $z_m$.\n",
"\n",
"By constructing a linear function $\\zeta(x, z)$ using these two points, derivatives of $\\zeta$ with respect to these variables can be calculated. Note that the horizontal component $x'(x, z) = x$. In this case, since $z_m$ and $z_M$ are defined as above, these derivatives can be calculated analytically. In practice, more complex topographies may require calculation via finite-difference approximations.\n",
"\n",
"Note that that this transformation is not a conformal mapping, simplifying this process as some of these derivatives go to zero or unity.\n",
"\n",
"## The transformation\n",
"\n",
"To implement our curvilinear solver, we need to know the transformation required to map the surface onto a horizontal plane, in this case the centre of the iteration space. For simplicity, we apply a linear warping of the z dimension."
"To implement our curvilinear solver, we need to know the transformation required to map the surface onto a horizontal plane, in this case the maximum of the iteration space. For simplicity, we apply a linear warping of the z dimension.\n",
"\n",
"Firstly, we define a function to create a symbolic linear fit between two points:"
]
},
{
Expand All @@ -19,11 +46,34 @@
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import sympy as sp\n",
"\n",
"def linear_fit(A, b, z):\n",
" \"\"\"Fit a linear function symbolically\"\"\"\n",
" sol = A**-1*b\n",
" \n",
" # Form the function\n",
" poly = sol[0] + sol[1]*z\n",
" return poly"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will use this to create a function to generate symbolic expressions for $\\zeta(x, z)$ and derivatives, its inverse $z(x', \\zeta)$, and the equation of the surface $z_M(x)$, and return them as Python functions."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def zeta_derivs(lambdify=False):\n",
"def zeta_derivs():\n",
" \"\"\"\n",
" Generates an expression for Z position of\n",
" cosine hill of height h, with a baseline\n",
Expand Down Expand Up @@ -67,53 +117,38 @@
" b = sp.Matrix([[zetam],\n",
" [zetaM]])\n",
"\n",
" sol = A**-1*b\n",
" \n",
" # Form the function\n",
" poly = sol[0] + sol[1]*z\n",
" poly = linear_fit(A, b, z)\n",
" \n",
" # Make a list of functions to return\n",
" funcs = []\n",
"\n",
" if lambdify: # Return as Python functions\n",
" params = [x, z, xM, xm, zb, zm, zetaM, zetam, h]\n",
" funcs.append(sp.lambdify(params, poly)) # No derivs taken\n",
" for i in range(1, 3): # Derivs\n",
" funcs.append(sp.lambdify(params, sp.diff(poly, x, i)))\n",
" funcs.append(sp.lambdify(params, sp.diff(poly, z, i)))\n",
" else: # Return as SymPy expressions\n",
" funcs.append(poly) # No derivs taken\n",
" for i in range(1, 3): # Derivs\n",
" funcs.append(sp.diff(poly, x, i))\n",
" funcs.append(sp.diff(poly, z, i))\n",
"\n",
" # Inverse of the linear warp (used for plotting)\n",
" C = sp.Matrix([[1, zetam],\n",
" [1, zetaM]])\n",
" d = sp.Matrix([[zm],\n",
" [zs]])\n",
"\n",
" polyinv = linear_fit(C, d, zeta)\n",
" \n",
" solinv = C**-1*d\n",
" # Inverse of the polynomial\n",
" polyinv = solinv[0] + solinv[1]*zeta\n",
" # Make a list of functions to return\n",
" funcs = []\n",
"\n",
" # zeta(x, z) and its derivatives\n",
" params = [x, z, xM, xm, zb, zm, zetaM, zetam, h]\n",
" funcs.append(sp.lambdify(params, poly)) # No derivs taken\n",
" for i in range(1, 3): # Symbolic derivatives\n",
" funcs.append(sp.lambdify(params, sp.diff(poly, x, i)))\n",
" funcs.append(sp.lambdify(params, sp.diff(poly, z, i)))\n",
" \n",
" if lambdify: # Return as Python function\n",
" paramsinv = [x, zeta, xM, xm, zb, zm, zetaM, zetam, h]\n",
" funcs.append(sp.lambdify(paramsinv, polyinv))\n",
" else: # Return as SymPy expression\n",
" funcs.append(polyinv)\n",
" # z(x, zeta) to be used for plotting purposes\n",
" paramsinv = [x, zeta, xM, xm, zb, zm, zetaM, zetam, h]\n",
" funcs.append(sp.lambdify(paramsinv, polyinv))\n",
" \n",
" # Function of the surface\n",
" if lambdify: # Return as Python function\n",
" paramssurf = [x, xM, xm, zb, h]\n",
" funcs.append(sp.lambdify(paramssurf, zs))\n",
" else: # Return as SymPy expression\n",
" funcs.append(zs)\n",
" paramssurf = [x, xM, xm, zb, h]\n",
" funcs.append(sp.lambdify(paramssurf, zs))\n",
"\n",
" funcs = tuple(funcs)\n",
" return funcs\n",
"\n",
"zeta_xz, zeta_dx1, zeta_dz1, zeta_dx2, zeta_dz2, z_xzeta, z_s = zeta_derivs(lambdify=True)"
"zeta_xz, zeta_dx1, zeta_dz1, zeta_dx2, zeta_dz2, z_xzeta, z_s = zeta_derivs()"
]
},
{
Expand All @@ -122,44 +157,44 @@
"source": [
"## Set up subdomains for the free surface\n",
"\n",
"For an MPI-safe free surface implementation, it is necessary to set up subdomains."
"For an MPI-safe free surface implementation, it is necessary to set up subdomains to avoid the need for hardcoded indices. See `userapi/03_subdomains.ipynb` for more detail on this functionality."
]
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from devito import SubDomain\n",
"\n",
"s_o = 8 # Space order\n",
"so = 8 # Space order\n",
"\n",
"\n",
"class MainDomain(SubDomain): # Main model domain\n",
" name = 'main'\n",
" def __init__(self, s_o):\n",
" def __init__(self, so):\n",
" super().__init__()\n",
" self.s_o = s_o\n",
" self.so = so\n",
" \n",
" def define(self, dimensions):\n",
" x, zeta = dimensions\n",
" return {x: x, zeta: ('middle', 0, 1+self.s_o//2)}\n",
" return {x: x, zeta: ('middle', 0, 1+self.so//2)}\n",
"\n",
"\n",
"class FreeSurface(SubDomain): # Free surface region\n",
" name = 'freesurface'\n",
" def __init__(self, s_o):\n",
" def __init__(self, so):\n",
" super().__init__()\n",
" self.s_o = s_o\n",
" self.so = so\n",
" \n",
" def define(self, dimensions):\n",
" x, zeta = dimensions\n",
" return {x: x, zeta: ('right', 1+self.s_o//2)}\n",
" return {x: x, zeta: ('right', 1+self.so//2)}\n",
" \n",
"\n",
"main_domain = MainDomain(s_o)\n",
"freesurface = FreeSurface(s_o)"
"main_domain = MainDomain(so)\n",
"freesurface = FreeSurface(so)"
]
},
{
Expand All @@ -177,7 +212,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 4,
"metadata": {
"scrolled": false
},
Expand Down Expand Up @@ -210,12 +245,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating derivatives of zeta with respect to z for the transformation:"
"The functions defined earlier can now be used to evaluate the requisite derivatives pertaining to the transformation. These values will be discretized as `Function`s onto the finite-difference grid."
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -243,6 +278,7 @@
"\n",
"# dzeta/dx\n",
"zeta_dx = Function(name='zeta_dx', grid=grid)\n",
"\n",
"zeta_dx.data[:] = zeta_dx1(x_msh, z_msh, x_M, x_m, z_b, z_m, zeta_M, zeta_m, h)\n",
"\n",
"# dzeta2/dx2\n",
Expand All @@ -259,7 +295,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 6,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -296,12 +332,14 @@
"source": [
"## The transformed equation\n",
"\n",
"Now we can write out our transformed equation. Both derivatives in the Laplacian will have modified expressions reflecting the new coordinate system."
"Now we can write out our transformed equation. Both derivatives in the Laplacian will have modified expressions reflecting the new coordinate system.\n",
"\n",
"These are substituted into the acoustic wave equation."
]
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -310,13 +348,13 @@
"# Wavespeed\n",
"c = 2. # km/s\n",
"\n",
"# Set up the function and equations\n",
"p = TimeFunction(name=\"p\", grid=grid, time_order=2, space_order=s_o)\n",
"# Set up the TimeFunction and equations\n",
"p = TimeFunction(name=\"p\", grid=grid, time_order=2, space_order=so)\n",
"# Transformed derivative expressions (note these are not the shortcuts)\n",
"pdx2 = p.dx2+p.dx.dzeta*zeta_dx+(p.dx.dzeta+p.dzeta2*zeta_dx)*zeta_dx+p.dzeta*zeta_dxx\n",
"pdz2 = p.dzeta2*zeta_dz**2\n",
"\n",
"# Looks just like the normal equation\n",
"# Superficially resembles the untransformed equation\n",
"eq = Eq(p.forward, 2*p - p.backward + c**2*(pdx2 + pdz2),\n",
" subdomain=grid.subdomains['main'])"
]
Expand All @@ -332,7 +370,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -372,9 +410,13 @@
" zetaind = f.indices[-1]\n",
" if (zetaind - zeta).as_coeff_Mul()[0] > 0:\n",
" # Flip the sign on folded-back points\n",
" s = sign(zetafs.symbolic_max - zetaind.subs({zeta: zetafs, zeta.spacing: 1}))\n",
" s = sign(zetafs.symbolic_max\n",
" - zetaind.subs({zeta: zetafs, zeta.spacing: 1}))\n",
" # Antisymmetric mirror\n",
" mapper.update({f: s*f.subs({zetaind: INT(zeta.symbolic_max*zeta.spacing - abs(zeta.symbolic_max*zeta.spacing - zetaind))})})\n",
" mapper.update({f:\n",
" s*f.subs({zetaind: INT(zeta.symbolic_max*zeta.spacing\n",
" - abs(zeta.symbolic_max*zeta.spacing\n",
" - zetaind))})})\n",
" return Eq(lhs, rhs.subs(mapper), subdomain=subdomain)\n",
" \n",
"fs = free_surface_top(eq, grid.subdomains['freesurface'], p)"
Expand All @@ -389,7 +431,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -414,7 +456,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 10,
"metadata": {},
"outputs": [
{
Expand All @@ -430,12 +472,12 @@
"PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
" PerfEntry(time=1e-06, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section1', rank=None),\n",
" PerfEntry(time=0.016621000000000028, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" PerfEntry(time=0.016083, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
" (PerfKey(name='section2', rank=None),\n",
" PerfEntry(time=1.7000000000000003e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
" PerfEntry(time=1.3000000000000003e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
]
},
"execution_count": 9,
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -458,7 +500,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 11,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -493,7 +535,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 12,
"metadata": {
"scrolled": false
},
Expand Down Expand Up @@ -523,9 +565,30 @@
"plt.ylabel(\"y (m)\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"154.39851"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.linalg.norm(p.data)"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
Expand Down

0 comments on commit 64e1d7d

Please sign in to comment.