{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4 - Hybdrid Absorbing Boundary Condition (HABC)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.1 - Introduction\n",
"\n",
"In this notebook we describe absorbing boundary conditions and their use combined with the *Hybdrid Absorbing Boundary Condition* (*HABC*). The common points to the previous notebooks Introduction to Acoustic Problem, Damping and PML will be used here, with brief descriptions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.2 - Absorbing Boundary Conditions\n",
"\n",
" We initially describe absorbing boundary conditions, the so called A1 and A2 Clayton's conditions and\n",
" the scheme from Higdon. These methods can be used as pure boundary conditions, designed to reduce reflections,\n",
" or as part of the Hybrid Absorbing Boundary Condition, in which they are combined with an absorption layer in a manner to be described ahead. \n",
" \n",
" In the presentation of these boundary conditions we initially consider the wave equation to be solved on\n",
" the spatial domain $\\Omega=\\left[x_{I},x_{F}\\right] \\times\\left[z_{I},z_{F}\\right]$ as show in the figure bellow. More details about the equation and domain definition can be found in the Introduction to Acoustic Problem notebook. \n",
"\n",
"\n",
" \n",
"## 4.2.1 - Clayton's A1 Boundary Condition\n",
"\n",
"Clayton's A1 boundary condition is based on a one way wave equation (OWWE). This simple condition\n",
" is such that outgoing waves normal to the border would leave without reflection. At the $\\partial \\Omega_1$ part of the boundary\n",
" we have, \n",
"\n",
"- $\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial t}-c(x,z)\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial x}=0.$\n",
"\n",
"while at $\\partial \\Omega_3$ the condition is\n",
" \n",
"- $\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial t}+c(x,z)\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial x}=0.$\n",
"\n",
"and at $\\partial \\Omega_2$\n",
"\n",
"- $\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial t}+c(x,z)\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial z}=0.$\n",
"\n",
"## 4.2.2 - Clayton's A2 Boundary Condition \n",
"\n",
"The A2 boundary condition aims to impose a boundary condition that would make outgoing waves leave the domain without being reflected. This condition is approximated (using a Padé approximation in the wave dispersion relation) by the following equation to be imposed on the boundary part $\\partial \\Omega_1$\n",
"\n",
"- $\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial t^{2}}+c(x,z)\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial x \\partial t}+\\frac{c^2(x,z)}{2}\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial z^{2}}=0.$\n",
"\n",
"At $\\partial \\Omega_3$ we have\n",
"\n",
"- $\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial t^{2}}-c(x,z)\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial z \\partial t}+\\frac{c^2(x,z)}{2}\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial x^{2}}=0.$\n",
"\n",
"while at $\\partial \\Omega_2$ the condition is\n",
"\n",
"- $\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial t^{2}}-c(x,z)\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial x \\partial t}+\\frac{c^2(x,z)}{2}\\displaystyle\\frac{\\partial^{2} u(x,z,t)}{\\partial z^{2}}=0.$\n",
"\n",
"At the corner points the condition is \n",
"\n",
"- $\\displaystyle\\frac{\\sqrt{2}\\partial u(x,z,t)}{\\partial t}+c(x,z)\\left(\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial x}+\\displaystyle\\frac{\\partial u(x,z,t)}{\\partial z}\\right)=0.$\n",
"\n",
"## 4.2.3 - Higdon Boundary Condition\n",
"\n",
"The Higdon Boundary condition of order p is given at $\\partial \\Omega_1$ and $\\partial \\Omega_3$n by:\n",
"\n",
"- $\\Pi_{j=1}^p(\\cos(\\alpha_j)\\left(\\displaystyle\\frac{\\partial }{\\partial t}-c(x,z)\\displaystyle\\frac{\\partial }{\\partial x}\\right)u(x,z,t)=0.$\n",
"\n",
"and at $\\partial \\Omega_2$\n",
"\n",
"- $\\Pi_{j=1}^p(\\cos(\\alpha_j)\\left(\\displaystyle\\frac{\\partial}{\\partial t}-c(x,z)\\displaystyle\\frac{\\partial}{\\partial z}\\right)u(x,z,t)=0.$\n",
"\n",
" This method would make that outgoing waves with angle of incidence at the boundary equal to $\\alpha_j$ would\n",
" present no reflection. The method we use in this notebook employs order 2 ($p=2$) and angles $0$ and $\\pi/4$.\n",
" \n",
" Observation: There are similarities between Clayton's A2 and the Higdon condition. If one chooses $p=2$ and\n",
" both angles equal to zero in Higdon's method, this leads to the condition:\n",
" $ u_{tt}-2cu_{xt}+c^2u_{xx}=0$. But, using the wave equation, we have that $c^2u_{xx}=u_{tt}-c^2u_{zz}$. Replacing this relation in the previous equation, we get: $2u_{tt}-2cu_{xt}-c^2u_{zz}=0$ which is Clayton's A2\n",
" boundary condition. In this sence, Higdon's method would generalize Clayton's scheme. But the discretization of\n",
" both methods are quite different, since in Higdon's scheme the boundary operators are unidirectional, while\n",
" in Clayton's A2 not."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.3 - Acoustic Problem with HABC\n",
"\n",
"In the hybrid absorption boundary condition (HABC) scheme we will also extend the spatial domain as $\\Omega=\\left[x_{I}-L,x_{F}+L\\right] \\times\\left[z_{I},z_{F}+L\\right]$.\n",
"\n",
"We added to the target domain $\\Omega_{0}=\\left[x_{I},x_{F}\\right]\\times\\left[z_{I},z_{F}\\right]$ an extension zone, of length $L$ in both ends of the direction $x$ and at the end of the domain in the direction $z$, as represented in the figure bellow.\n",
"\n",
"\n",
"\n",
" The difference with respect to previous schemes, is that this extended region will now be considered as the union of several gradual extensions. As represented in the next figure, we define a region $A_M=\\Omega_{0}$. The regions $A_k, k=M-1,\\cdots,1$ will be defined as the previous region $A_{k+1}$ to which we add one extra grid line to the left,\n",
" right and bottom sides of it, such that the final region $A_1=\\Omega$ (we thus have $M=L+1$).\n",
" \n",
" \n",
" \n",
" We now consider the temporal evolution\n",
" of the solution of the HABC method. Suppose that $u(x,z,t-1)$ is the solution at a given instant $t-1$ in all the \n",
" extended $\\Omega$ domain. We update it to instant $t$, using one of the absorbing boundary conditions described in the previous section (A1, A2 or Higdon) producing a preliminar new function $u(x,z,t)$. Now, call $u_{1}(x,z,t)$ the solution at instant $t$ constructed in the extended region, by applying the same absorbing boundary condition at the border of each of the domains $A_k,k=1,..,M$. The HABC solution will be constructed as a convex combination of $u(x,z,t)$ and $u_{1}(x,z,t)$:\n",
" \n",
"- $u(x,z,t) = (1-\\omega)u(x,z,t)+\\omega u_{1}(x,z,t)$.\n",
"\n",
"The function $u_{1}(x,z,t)$ is defined (and used) only in the extension of the domain. The function $w$ is a \n",
"weight function growing from zero at the boundary $\\partial\\Omega_{0}$ to one at $\\partial\\Omega$. The particular weight function to be used could vary linearly, as when the scheme was first proposed by Liu and Sen. But HABC produces better results with a non-linear weight function to be described ahead.\n",
"\n",
"The wave equation employed here will be the same as in the previous notebooks, with same velocity model, source term and initial conditions.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4.3.1 The weight function $\\omega$\n",
"\n",
"One can choose a *linear* weight function as \n",
"\n",
"\\begin{equation}\n",
"\\omega_{k} = \\displaystyle\\frac{M-k}{M};\n",
"\\end{equation}\n",
"\n",
"or preferably a *non linear*\n",
"\n",
"\\begin{equation}\n",
"\\omega_{k}=\\left\\{ \\begin{array}{ll}\n",
"1, & \\textrm{if $1\\leq k \\leq P+1$,} \\\\ \\left(\\displaystyle\\frac{M-k}{M-P}\\right)^{\\alpha} , & \\textrm{if $P+2 \\leq k \\leq M-1.$} \\\\ 0 , & \\textrm{if $k=M$.}\\end{array}\\right.\n",
"\\label{eq:elo8}\n",
"\\end{equation} \n",
"\n",
"In general we take $P=2$ and we choose $\\alpha$ as follows:\n",
"\n",
"- $\\alpha = 1.5 + 0.07(npt-P)$, in the case of A1 and A2;\n",
"- $\\alpha = 1.0 + 0.15(npt-P)$, in the case of Higdon.\n",
"\n",
"The value *npt* designates the number of discrete points that define the length of the blue band in the direction $x$ and/or $z$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.4 - Finite Difference Operators and Discretization of Spatial and Temporal Domains\n",
"\n",
" We employ the same methods as in the previous notebooks. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.5 - Standard Problem\n",
"\n",
"Redeeming the Standard Problem definitions discussed on the notebook Introduction to Acoustic Problem we have that:\n",
"\n",
"- $x_{I}$ = 0.0 Km;\n",
"- $x_{F}$ = 1.0 Km = 1000 m;\n",
"- $z_{I}$ = 0.0 Km;\n",
"- $z_{F}$ = 1.0 Km = 1000 m;\n",
"\n",
"The spatial discretization parameters are given by:\n",
"- $\\Delta x$ = 0.01 km = 10m;\n",
"- $\\Delta z$ = 0.01 km = 10m;\n",
"\n",
"Let's consider a $I$ the time domain with the following limitations:\n",
"\n",
"- $t_{I}$ = 0 s = 0 ms;\n",
"- $t_{F}$ = 1 s = 1000 ms;\n",
"\n",
"The temporal discretization parameters are given by:\n",
"\n",
"- $\\Delta t$ $\\approx$ 0.0016 s = 1.6 ms;\n",
"- $NT$ = 626.\n",
"\n",
"The source term, velocity model and positioning of receivers will be as in the previous notebooks."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4.6 - Numerical Simulations\n",
"\n",
"For the numerical simulations of this notebook we use several of the notebook codes presented in Damping e PML. The new features will be described in more detail.\n",
"\n",
"So, we import the following Python and Devito packages:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# NBVAL_IGNORE_OUTPUT\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plot\n",
"import math as mt\n",
"import matplotlib.ticker as mticker \n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"from matplotlib import cm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From Devito's library of examples we import the following structures:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# NBVAL_IGNORE_OUTPUT\n",
"\n",
"%matplotlib inline\n",
"from examples.seismic import TimeAxis\n",
"from examples.seismic import RickerSource\n",
"from examples.seismic import Receiver\n",
"from examples.seismic import plot_velocity\n",
"from devito import SubDomain, Grid, NODE, TimeFunction, Function, Eq, solve, Operator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mesh parameters that we choose define the domain $\\Omega_{0}$ plus the absorption region. For this, we use the following data:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"nptx = 101\n",
"nptz = 101\n",
"x0 = 0.\n",
"x1 = 1000. \n",
"compx = x1-x0\n",
"z0 = 0.\n",
"z1 = 1000.\n",
"compz = z1-z0;\n",
"hxv = (x1-x0)/(nptx-1)\n",
"hzv = (z1-z0)/(nptz-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we saw previously, HABC has three approach possibilities (A1, A2 and Higdon) and two types of weights (linear and non-linear). So, we insert two control variables. The variable called *habctype* chooses the type of HABC approach and is such that:\n",
"\n",
"- *habctype=1* is equivalent to choosing A1;\n",
"- *habctype=2* is equivalent to choosing A2;\n",
"- *habctype=3* is equivalent to choosing Higdon;\n",
"\n",
"Regarding the weights, we will introduce the variable *habcw* that chooses the type of weight and is such that:\n",
"\n",
"- *habcw=1* is equivalent to linear weight;\n",
"- *habcw=2* is equivalent to non-linear weights;\n",
"\n",
"In this way, we make the following choices:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"habctype = 3\n",
"habcw = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The number of points of the absorption layer in the directions $x$ and $z$ are given, respectively, by:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"npmlx = 20\n",
"npmlz = 20"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The lengths $L_{x}$ and $L_{z}$ are given, respectively, by:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"lx = npmlx*hxv\n",
"lz = npmlz*hzv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the construction of the *grid* we have:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"nptx = nptx + 2*npmlx\n",
"nptz = nptz + 1*npmlz\n",
"x0 = x0 - hxv*npmlx\n",
"x1 = x1 + hxv*npmlx\n",
"compx = x1-x0\n",
"z0 = z0\n",
"z1 = z1 + hzv*npmlz\n",
"compz = z1-z0\n",
"origin = (x0,z0)\n",
"extent = (compx,compz)\n",
"shape = (nptx,nptz)\n",
"spacing = (hxv,hzv)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the case of the acoustic equation with Damping and in the acoustic equation with PML, we can define specific regions in our domain, since the solution $u_{1}(x,z,t)$ is only calculated in the blue region. We will soon follow a similar scheme for creating *subdomains* as was done on notebooks Damping and PML.\n",
"\n",
"First, we define a region corresponding to the entire domain, naming this region as *d0*. In the language of *subdomains* *d0* it is written as:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"class d0domain(SubDomain):\n",
" name = 'd0'\n",
" def define(self, dimensions):\n",
" x, z = dimensions\n",
" return {x: x, z: z}\n",
"d0_domain = d0domain()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The blue region will be built with 3 divisions:\n",
"\n",
"- *d1* represents the left range in the direction *x*, where the pairs $(x,z)$ satisfy: $x\\in\\{0,npmlx\\}$ and $z\\in\\{0,nptz\\}$;\n",
"- *d2* represents the rigth range in the direction *x*, where the pairs $(x,z)$ satisfy: $x\\in\\{nptx-npmlx,nptx\\}$ and $z\\in\\{0,nptz\\}$;\n",
"- *d3* represents the left range in the direction *y*, where the pairs $(x,z)$ satisfy: $x\\in\\{npmlx,nptx-npmlx\\}$ and $z\\in\\{nptz-npmlz,nptz\\}$;\n",
"\n",
"Thus, the regions *d1*, *d2* and *d3* aare described as follows in the language of *subdomains*:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"class d1domain(SubDomain):\n",
" name = 'd1'\n",
" def define(self, dimensions):\n",
" x, z = dimensions\n",
" return {x: ('left',npmlx), z: z}\n",
"d1_domain = d1domain()\n",
"\n",
"class d2domain(SubDomain):\n",
" name = 'd2'\n",
" def define(self, dimensions):\n",
" x, z = dimensions\n",
" return {x: ('right',npmlx), z: z}\n",
"d2_domain = d2domain()\n",
"\n",
"class d3domain(SubDomain):\n",
" name = 'd3'\n",
" def define(self, dimensions):\n",
" x, z = dimensions\n",
" if((habctype==3)&(habcw==1)):\n",
" return {x: x, z: ('right',npmlz)}\n",
" else:\n",
" return {x: ('middle', npmlx, npmlx), z: ('right',npmlz)}\n",
"d3_domain = d3domain()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The figure below represents the division of domains that we did previously:\n",
"\n",
"\n",
"\n",
"After we defining the spatial parameters and constructing the *subdomains*, we then generate the *spatial grid* and set the velocity field:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"grid = Grid(origin=origin, extent=extent, shape=shape, subdomains=(d0_domain,d1_domain,d2_domain,d3_domain), dtype=np.float64)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"v0 = np.zeros((nptx,nptz)) \n",
"X0 = np.linspace(x0,x1,nptx)\n",
"Z0 = np.linspace(z0,z1,nptz)\n",
" \n",
"x10 = x0+lx\n",
"x11 = x1-lx\n",
" \n",
"z10 = z0\n",
"z11 = z1 - lz\n",
"\n",
"xm = 0.5*(x10+x11)\n",
"zm = 0.5*(z10+z11)\n",
" \n",
"pxm = 0\n",
"pzm = 0\n",
" \n",
"for i in range(0,nptx):\n",
" if(X0[i]==xm): pxm = i\n",
" \n",
"for j in range(0,nptz):\n",
" if(Z0[j]==zm): pzm = j\n",
" \n",
"p0 = 0 \n",
"p1 = pzm\n",
"p2 = nptz\n",
" \n",
"v0[0:nptx,p0:p1] = 1.5\n",
"v0[0:nptx,p1:p2] = 2.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Previously we introduce the local variables *x10,x11,z10,z11,xm,zm,pxm* and *pzm* that help us to create a specific velocity field, where we consider the whole domain (including the absorpion region). Below we include a routine to plot the velocity field."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def graph2dvel(vel):\n",
" plot.figure()\n",
" plot.figure(figsize=(16,8))\n",
" fscale = 1/10**(3)\n",
" scale = np.amax(vel[npmlx:-npmlx,0:-npmlz])\n",
" extent = [fscale*(x0+lx),fscale*(x1-lx), fscale*(z1-lz), fscale*(z0)]\n",
" fig = plot.imshow(np.transpose(vel[npmlx:-npmlx,0:-npmlz]), vmin=0.,vmax=scale, cmap=cm.seismic, extent=extent)\n",
" plot.gca().xaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f km'))\n",
" plot.gca().yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f km'))\n",
" plot.title('Velocity Profile')\n",
" plot.grid()\n",
" ax = plot.gca()\n",
" divider = make_axes_locatable(ax)\n",
" cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
" cbar = plot.colorbar(fig, cax=cax, format='%.2e')\n",
" cbar.set_label('Velocity [km/s]')\n",
" plot.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below we include the plot of velocity field."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"