diff --git a/docs/source/Advanced/boundary/boundary.rst b/docs/source/Advanced/boundary/boundary.rst new file mode 100644 index 0000000..c20e97f --- /dev/null +++ b/docs/source/Advanced/boundary/boundary.rst @@ -0,0 +1,29 @@ +:author: 朱邓达 +:date: 2025-01-08 + +顶底界面的边界条件 +===================== + +默认情况下，一维层状模型的顶界面为 **自由边界** ，即过 z 平面的牵引力为 0； +底界面为 **半无限空间** ，其对应模型文件的最后一行。 + +自版本 0.15.0 起， **PyGRT** 支持设置顶底界面的边界条件，可设置为 **自由边界，刚性边界和半无限空间** ， 公式推导详见 :doc:`/Formula/boundary` ，示例详见 :doc:`/Gallery/ex16/ex16` 。 +为了不改变默认的模型输入，顶、底界面对边界条件的处理会有些区别，详见以下示意图。 + +.. figure:: plot_boundary.svg + :align: center + +通过以下可选参数，可设置顶底界面的边界条件。 + +.. tabs:: + + .. group-tab:: C + + 详见 :doc:`/Module/greenfn` 和 :doc:`/Module/static_greenfn` 模块的 **-B** 选项。 + + .. group-tab:: Python + + :func:`PyModel1D() ` 创建模型时可设置边界条件： + + + ``topbound:Literal['free', 'rigid', 'halfspace']`` 顶层边界条件 + + ``botbound:Literal['free', 'rigid', 'halfspace']`` 底层边界条件 diff --git a/docs/source/Advanced/boundary/plot_boundary.svg b/docs/source/Advanced/boundary/plot_boundary.svg new file mode 100644 index 0000000..4024927 --- /dev/null +++ b/docs/source/Advanced/boundary/plot_boundary.svg @@ -0,0 +1,767 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/Formula/boundary.rst b/docs/source/Formula/boundary.rst new file mode 100644 index 0000000..acce38e --- /dev/null +++ b/docs/source/Formula/boundary.rst @@ -0,0 +1,90 @@ +:author: 朱邓达 +:date: 2026-01-07 + +自由边界和刚性边界的反射系数 +============================================ + +本质上与 |yao2026| 中推导顶层的自由界面反射系数的过程基本一致，只是替换了相应的矩阵表达式。 + +具体而言， 以动态 P-SV 为例， + ++ 自由边界满足过 z 平面的牵引力为 0， 即 + +.. math:: + + \begin{bmatrix} + q_m \\ + w_m \\ + 0 \\ + 0 \\ + \end{bmatrix} = + \begin{bmatrix} + k & b & k & -b \\ + a & k & -a & k \\ + 2\mu\Omega & 2k\mu b & 2\mu\Omega & -2k\mu b \\ + 2k\mu a & 2\mu\Omega & -2k\mu a & 2\mu\Omega \\ + \end{bmatrix} + \begin{bmatrix} + \phi_m^- \\ + \psi_m^- \\ + \phi_m^+ \\ + \psi_m^+ \\ + \end{bmatrix} + ++ 刚性边界满足位移为 0 ，即 + +.. math:: + + \begin{bmatrix} + 0 \\ + 0 \\ + \sigma_{Rm} \\ + \tau_{Rm} \\ + \end{bmatrix} = + \begin{bmatrix} + k & b & k & -b \\ + a & k & -a & k \\ + 2\mu\Omega & 2k\mu b & 2\mu\Omega & -2k\mu b \\ + 2k\mu a & 2\mu\Omega & -2k\mu a & 2\mu\Omega \\ + \end{bmatrix} + \begin{bmatrix} + \phi_m^- \\ + \psi_m^- \\ + \phi_m^+ \\ + \psi_m^+ \\ + \end{bmatrix} + +根据以上结果为 0 的表达式，可以构建出以下关系： + ++ 对于顶层界面的 z+ 侧 + +.. math:: + + \begin{bmatrix} + \phi_m^+ \\ + \psi_m^+ \\ + \end{bmatrix} = + \mathbf{R}_U + \begin{bmatrix} + \phi_m^- \\ + \psi_m^- \\ + \end{bmatrix} + ++ 对于底层界面的 z- 侧 + +.. math:: + + \begin{bmatrix} + \phi_m^- \\ + \psi_m^- \\ + \end{bmatrix} = + \mathbf{R}_D + \begin{bmatrix} + \phi_m^+ \\ + \psi_m^+ \\ + \end{bmatrix} + + +具体公式我使用 `SymPy `_ 做了推导，相应的 ``.ipynb`` 文件如下，以供读者参考。 + +:download:`boundary_condition.ipynb` （ :doc:`预览 ` ） diff --git a/docs/source/Formula/boundary_condition.ipynb b/docs/source/Formula/boundary_condition.ipynb new file mode 100644 index 0000000..3cf37db --- /dev/null +++ b/docs/source/Formula/boundary_condition.ipynb @@ -0,0 +1,393 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6ca6ddd0", + "metadata": {}, + "source": [ + "# 自由边界和刚性边界的反射系数\n", + "\n", + "以 P-SV 为例\n", + "\n", + "+ Author: Zhu Dengda \n", + "+ Email: zhudengda@mail.iggcas.ac.cn" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "b5a7a686", + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sp" + ] + }, + { + "cell_type": "markdown", + "id": "5811a951", + "metadata": {}, + "source": [ + "## 动态解" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "ec9378bd", + "metadata": {}, + "outputs": [], + "source": [ + "# 定义基本变量\n", + "k, b, d, a, mu, Omg, kb, ka, w, Delta = sp.symbols(r'k b d, a mu \\Omega k_{b} k_{a} \\omega \\Delta')\n", + "lam, rho = sp.symbols(r'\\lambda \\rho')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "789ac906", + "metadata": {}, + "outputs": [], + "source": [ + "D = sp.Matrix([\n", + " [k, b, k, -b],\n", + " [a, k, -a, k],\n", + " [2*mu*Omg, 2*k*mu*b, 2*mu*Omg, -2*k*mu*b],\n", + " [2*k*mu*a, 2*mu*Omg, -2*k*mu*a, 2*mu*Omg]\n", + "])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "239e2dd2", + "metadata": {}, + "outputs": [], + "source": [ + "def get_R(D, at_bottom:bool):\n", + " if not at_bottom:\n", + " return - (D[:, 2:])**(-1) * (D[:, :2])\n", + " else:\n", + " return - (D[:, :2])**(-1) * (D[:, 2:])" + ] + }, + { + "cell_type": "markdown", + "id": "8feffbc4", + "metadata": {}, + "source": [ + "### 自由界面（z平面牵引力为0）" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "eafaff07", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}} & \\frac{2 \\Omega b k}{- \\Omega^{2} + a b k^{2}}\\\\\\frac{2 \\Omega a k}{- \\Omega^{2} + a b k^{2}} & \\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[(\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2), 2*\\Omega*b*k/(-\\Omega**2 + a*b*k**2)],\n", + "[ 2*\\Omega*a*k/(-\\Omega**2 + a*b*k**2), (\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2)]])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z+ 侧， 例如顶层顶界面\n", + "R = get_R(D[2:, :], False)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "002464b0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}\\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}} & \\frac{2 \\Omega b k}{\\Omega^{2} - a b k^{2}}\\\\\\frac{2 \\Omega a k}{\\Omega^{2} - a b k^{2}} & \\frac{\\Omega^{2} + a b k^{2}}{- \\Omega^{2} + a b k^{2}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[(\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2), 2*\\Omega*b*k/(\\Omega**2 - a*b*k**2)],\n", + "[ 2*\\Omega*a*k/(\\Omega**2 - a*b*k**2), (\\Omega**2 + a*b*k**2)/(-\\Omega**2 + a*b*k**2)]])" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z- 侧， 例如底层底界面\n", + "R = get_R(D[2:, :], True)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "markdown", + "id": "b171d3f9", + "metadata": {}, + "source": [ + "### 刚性界面（位移为0）" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "1f167bde", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}\\frac{a b + k^{2}}{a b - k^{2}} & \\frac{2 b k}{a b - k^{2}}\\\\\\frac{2 a k}{a b - k^{2}} & \\frac{a b + k^{2}}{a b - k^{2}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[(a*b + k**2)/(a*b - k**2), 2*b*k/(a*b - k**2)],\n", + "[ 2*a*k/(a*b - k**2), (a*b + k**2)/(a*b - k**2)]])" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z+ 侧\n", + "R = get_R(D[:2, :], False)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8dd46e32", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}\\frac{a b + k^{2}}{a b - k^{2}} & - \\frac{2 b k}{a b - k^{2}}\\\\- \\frac{2 a k}{a b - k^{2}} & \\frac{a b + k^{2}}{a b - k^{2}}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[(a*b + k**2)/(a*b - k**2), -2*b*k/(a*b - k**2)],\n", + "[ -2*a*k/(a*b - k**2), (a*b + k**2)/(a*b - k**2)]])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z- 侧\n", + "R = get_R(D[:2, :], True)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "markdown", + "id": "76a678aa", + "metadata": {}, + "source": [ + "## 静态解" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "2c6d0581", + "metadata": {}, + "outputs": [], + "source": [ + "# 在静态解中，矩阵 D 和相对深度 z-zj 有关, 因此 z+ 侧和 z- 侧需使用两种 D 矩阵\n", + "D1 = sp.Matrix([\n", + " [1, -(1 + 2*k*Delta*d), 1, -(1 - 2*k*Delta*d)],\n", + " [1, (1 - 2*k*Delta*d), -1, -(1 + 2*k*Delta*d)],\n", + " [2*mu*k, 2*mu*Delta*k*(1 - 2*k*d), 2*mu*k, 2*mu*Delta*k*(1 + 2*k*d)],\n", + " [2*mu*k, -2*mu*Delta*k*(1 + 2*k*d), -2*mu*k, 2*mu*Delta*k*(1 - 2*k*d)],\n", + "])\n", + "\n", + "z = 0\n", + "D2 = sp.Matrix([\n", + " [1, -(1 + 2*k*Delta*z), 1, -(1 - 2*k*Delta*z)],\n", + " [1, (1 - 2*k*Delta*z), -1, -(1 + 2*k*Delta*z)],\n", + " [2*mu*k, 2*mu*Delta*k*(1 - 2*k*z), 2*mu*k, 2*mu*Delta*k*(1 + 2*k*z)],\n", + " [2*mu*k, -2*mu*Delta*k*(1 + 2*k*z), -2*mu*k, 2*mu*Delta*k*(1 - 2*k*z)],\n", + "])" + ] + }, + { + "cell_type": "markdown", + "id": "15d392c1", + "metadata": {}, + "source": [ + "### 自由界面（z平面牵引力为0）" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "8afb9768", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & - \\Delta\\\\- \\frac{1}{\\Delta} & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 0, -\\Delta],\n", + "[-1/\\Delta, 0]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z+ 侧\n", + "R = get_R(D2[2:, :], False)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "0faa5c1f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}- 2 d k & \\Delta \\left(- 4 d^{2} k^{2} - 1\\right)\\\\- \\frac{1}{\\Delta} & - 2 d k\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ -2*d*k, \\Delta*(-4*d**2*k**2 - 1)],\n", + "[-1/\\Delta, -2*d*k]])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z- 侧\n", + "R = get_R(D1[2:, :], True)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "markdown", + "id": "415c011b", + "metadata": {}, + "source": [ + "### 刚性界面（位移为0）" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "849e8723", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}0 & 1\\\\1 & 0\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[0, 1],\n", + "[1, 0]])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z+ 侧\n", + "R = get_R(D2[:2, :], False)\n", + "R.simplify()\n", + "R" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "ef82bd64", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}2 \\Delta d k & 4 \\Delta^{2} d^{2} k^{2} + 1\\\\1 & 2 \\Delta d k\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[2*\\Delta*d*k, 4*\\Delta**2*d**2*k**2 + 1],\n", + "[ 1, 2*\\Delta*d*k]])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# z- 侧\n", + "R = get_R(D1[:2, :], True)\n", + "R.simplify()\n", + "R" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py310", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/Formula/index.rst b/docs/source/Formula/index.rst index 65da52c..d1667b1 100644 --- a/docs/source/Formula/index.rst +++ b/docs/source/Formula/index.rst @@ -22,6 +22,7 @@ - :doc:`RT_liquid` - :doc:`uiz` - :doc:`DS_zero` + - :doc:`boundary` **静力学震源参数** @@ -42,6 +43,7 @@ RT_liquid uiz DS_zero + boundary static_uniform static_source \ No newline at end of file diff --git a/docs/source/Gallery/ex16/ex16.rst b/docs/source/Gallery/ex16/ex16.rst new file mode 100644 index 0000000..e027efe --- /dev/null +++ b/docs/source/Gallery/ex16/ex16.rst @@ -0,0 +1,50 @@ +:author: 朱邓达 +:date: 2025-01-08 + +(16) 不同边界条件的测试 +------------------------------------- + +下载示例： :download:`ex16.tar.gz` + +在 :doc:`/Advanced/boundary/boundary` 部分和 :doc:`/Formula/boundary` 部分对边界条件和公式进行了介绍， +尽管顶、底界面对边界条件的处理有些区别，但经过适当配置源点场点深度以及模型参数，可形成两个 **上下翻转** 的系统， +从而进行对比验证。 + +.. figure:: flip1.SVG + :align: center + + 红点为源点，蓝点为场点 + + +如上图所示，两个系统呈现上下翻转的关系，这样二者计算的结果之间最多只有正负号的符号差异， +这为顶底界面的边界条件提供了一个互相验证的测试。 +以下进行一些对比测试，可见对分量符号进行校正后，结果完美吻合。 + +示例一 +~~~~~~~~~~~~~~~ + +.. figure:: free_halfspace.svg + :align: center + + 动态解已卷积阶跃函数，下同 + +示例二 +~~~~~~~~~~~~~~~ + +.. figure:: rigid_halfspace.svg + :align: center + +示例三 +~~~~~~~~~~~~~~~ + +以上实例都是在层状介质中进行。若仅设置单层介质并把顶底界面都设置为半空间， +则相当于求解 **无限均匀空间中的格林函数** ，这可用解析解验证。 + +.. figure:: halfspace_halfspace.svg + :align: center + +以下为解析解，其中动态解公式参考《Quantitative Seismology》， +静态解公式参考 :doc:`/Formula/static_uniform` 。与上图比较可见结果吻合。 + +.. figure:: theoretical.svg + :align: center \ No newline at end of file diff --git a/docs/source/Gallery/ex16/flip1.SVG b/docs/source/Gallery/ex16/flip1.SVG new file mode 100644 index 0000000..b588b3c --- /dev/null +++ b/docs/source/Gallery/ex16/flip1.SVG @@ -0,0 +1,523 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/Gallery/ex16/infinite_space.py b/docs/source/Gallery/ex16/infinite_space.py new file mode 100644 index 0000000..b123532 --- /dev/null +++ b/docs/source/Gallery/ex16/infinite_space.py @@ -0,0 +1,120 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.figure import Figure + +nt = 801 # 总点数，不要求2的幂次 +dt = 0.01 # 采样时间间隔(s) +t = np.arange(0, nt)*dt +rs = np.arange(0.01, 10.01, 0.01) + +Vp, Vs, Rho = 5.0, 3.0, 2.2 +mu = Rho * Vs**2 +lam = Rho * Vp**2 - 2*mu +Delta = (lam + mu) / (lam + 3.0*mu) + +depsrc = 3. # 震源深度 km +deprcv = 1. # 台站深度 km + +def static_theoretical(ij:str): + D = np.hypot(depsrc - deprcv, rs) + theta = np.arctan((depsrc - deprcv) / rs) + gamz = np.sin(theta) + gamr = np.cos(theta) + + delta_ij = 1.0 if ij[0]==ij[1] else 0.0 + gam_ij = np.ones_like(D) + for s in list(ij): + if s == 'z': + gam_ij *= gamz + elif s == 'r': + gam_ij *= gamr + elif s == 't': + gam_ij *= 0.0 + + return 1/(4.0*np.pi*mu) * Delta/(1 + Delta) / D * (delta_ij / Delta + gam_ij) + +def dynamic_theoretical(ij:str): + r = 10 + R = np.hypot(depsrc - deprcv, r) + theta = np.arctan((depsrc - deprcv) / r) + gamz = np.sin(theta) + gamr = np.cos(theta) + + delta_ij = 1.0 if ij[0]==ij[1] else 0.0 + gam_ij = 1.0 + for s in list(ij): + if s == 'z': + gam_ij *= gamz + elif s == 'r': + gam_ij *= gamr + elif s == 't': + gam_ij *= 0.0 + + def heaviside(t0): + arr = np.zeros_like(t) + arr[t >= t0] = 1 + return arr + + # 积分后的结果 + near = 1/(4.0*np.pi*Rho) * (3.0*gam_ij - delta_ij) / R**3 * 0.5 * heaviside(R/Vp)*(np.fmin(t, R/Vs)**2 - (R/Vp)**2) + Pfar = 1/(4.0*np.pi*Rho*Vp**2) * gam_ij / R * heaviside(R/Vp) + Sfar = - 1/(4.0*np.pi*Rho*Vs**2) * (gam_ij - delta_ij) / R * heaviside(R/Vs) + + return near + Pfar + Sfar + +# ============================================================= +chLst = ['VFZ', 'VFR', 'HFZ', 'HFR', 'HFT'] + +staticLst = [] +staticLst.append(static_theoretical('zz') * (-1)) +staticLst.append(static_theoretical('zr') * (-1)) +staticLst.append(static_theoretical('zr')) +staticLst.append(static_theoretical('rr')) +staticLst.append(static_theoretical('tt')) + +dynamicLst = [] +dynamicLst.append(dynamic_theoretical('zz') * (-1)) +dynamicLst.append(dynamic_theoretical('zr') * (-1)) +dynamicLst.append(dynamic_theoretical('zr')) +dynamicLst.append(dynamic_theoretical('rr')) +dynamicLst.append(dynamic_theoretical('tt')) + +fig = plt.figure(figsize=(8, 6)) +fig1, fig2 = fig.subfigures(1, 2, wspace=0.1) + +axs1 = fig1.subplots(5, 1, sharex=True, gridspec_kw=dict(hspace=0)) +axs2 = fig2.subplots(5, 1, sharex=True, gridspec_kw=dict(hspace=0)) +for i in range(5): + prop1 = dict(ls='-', c='b', lw=2) + + ax = axs1[i] + ax.plot(t, dynamicLst[i], **prop1) + ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, ha='right', va='top', bbox=dict(fc='w')) + ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + ax.grid() + ax.set_xmargin(0) + ax.set_ymargin(0.3) + + ax = axs2[i] + ax.plot(rs, staticLst[i], **prop1) + ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, ha='right', va='top', bbox=dict(fc='w')) + ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + ax.grid() + ax.set_xmargin(0) + ax.set_ymargin(0.3) + +ax = axs1[0] +ax.set_title("Dynamic displacements") +ax = axs1[-1] +ax.set_xlabel("Time (s)") + +ax = axs2[0] +ax.set_title("Static displacements") +ax = axs2[-1] +ax.set_xlabel("Distance (km)") + +fig.text(0.0, 1.0, f"theoretical solution in an infinite space", + ha='left', va='bottom', fontsize=12, bbox=dict(fc='w')) + +fig.tight_layout() +fig.savefig(f"theoretical.svg", bbox_inches='tight') diff --git a/docs/source/Gallery/ex16/plot.py b/docs/source/Gallery/ex16/plot.py new file mode 100644 index 0000000..422efff --- /dev/null +++ b/docs/source/Gallery/ex16/plot.py @@ -0,0 +1,117 @@ +import pygrt +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.figure import Figure +import sys + +yarr = np.arange(0.01, 10.01, 0.01) +rs = np.array([10]) # 震中距数组，km +nt = 801 # 总点数，不要求2的幂次 +dt = 0.01 # 采样时间间隔(s) +Vp, Vs, Rho = 5.0, 3.0, 2.2 +t = np.arange(0, nt)*dt + +bound1 = sys.argv[1] +bound2 = sys.argv[2] +isuniform = sys.argv[3] == 'y' + +# ============================================================= +depsrc = 3. # 震源深度 km +deprcv = 1. # 台站深度 km +if isuniform: + modarr = np.array([ + [2.0, Vp, Vs, Rho], + ]) +else: + modarr = np.array([ + [2.0, Vp, Vs, Rho], + [1.6, Vp+0.1, Vs+0.1, Rho+0.1], + [2.8, Vp+0.4, Vs+0.4, Rho+0.2], + [0, Vp+0.6, Vs+0.6, Rho+0.3], + ]) +pymod1 = pygrt.PyModel1D(modarr, depsrc, deprcv, topbound=bound1, botbound=bound2) +st1 = pymod1.compute_grn(distarr=rs, nt=nt, dt=dt, keepAllFreq=True)[0] +pygrt.utils.stream_integral(st1) +static1 = pymod1.compute_static_grn(xarr=[0.0], yarr=yarr) + +# ============================================================= +# 设置上下翻转模型 +if bound2 == 'halfspace': + modarr2 = np.flipud(modarr) + H = np.sum(modarr[:-1, 0]) + modarr2[0,0] = max(max((depsrc, deprcv)) - H, 0) + 1.0 # 顶部半空间加 1km 厚度作为参考偏移 +else: + modarr2 = np.flipud(modarr[:-1, :]) + +H = np.sum(modarr2[:, 0]) +if bound1 != 'halfspace': + modarr2 = np.vstack((modarr2, modarr[0,:])) + +depsrc2 = H - depsrc # 震源深度 km +deprcv2 = H - deprcv # 台站深度 km +print(modarr2) +print(modarr2.shape, depsrc2, deprcv2) + +pymod2 = pygrt.PyModel1D(modarr2, depsrc2, deprcv2, topbound=bound2, botbound=bound1) # 整理好的模型对象 +st2 = pymod2.compute_grn(distarr=rs, nt=nt, dt=dt, keepAllFreq=True)[0] +static2 = pymod2.compute_static_grn(xarr=[0.0], yarr=yarr) +pygrt.utils.stream_integral(st2) + + +# ============================================================= +chLst = ['VFZ', 'VFR', 'HFZ', 'HFR', 'HFT'] + +fig = plt.figure(figsize=(8, 6)) +fig1, fig2 = fig.subfigures(1, 2, wspace=0.1) + +axs1 = fig1.subplots(5, 1, sharex=True, gridspec_kw=dict(hspace=0)) +axs2 = fig2.subplots(5, 1, sharex=True, gridspec_kw=dict(hspace=0)) +for i in range(5): + ch = chLst[i] + sgn = -1 if ch in ['VFR', 'HFZ'] else 1 + + prop1 = dict(ls='-', c='0.7', lw=3.5) + prop2 = dict(ls='--', c='k', lw=1.2) + + ax = axs1[i] + data = st1.select(channel=chLst[i])[0].data + ax.plot(t, data, **prop1) + data = st2.select(channel=chLst[i])[0].data * sgn + ax.plot(t, data, **prop2) + ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, + ha='right', va='top', bbox=dict(fc='w')) + + ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + ax.grid() + ax.set_xmargin(0) + ax.set_ymargin(0.3) + + ax = axs2[i] + ax.plot(yarr, static1[chLst[i]][0], **prop1) + ax.plot(yarr, static2[chLst[i]][0] * sgn, **prop2) + ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, + ha='right', va='top', bbox=dict(fc='w')) + ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + ax.grid() + ax.set_xmargin(0) + ax.set_ymargin(0.3) + +ax = axs1[0] +ax.set_title("Dynamic displacements") +ax = axs1[-1] +ax.set_xlabel("Time (s)") + +ax = axs2[0] +ax.set_title("Static displacements") +ax = axs2[-1] +ax.set_xlabel("Distance (km)") + +fig.text(0.0, 1.0, f"gray-solid: {bound1} + {bound2}\n" + f"black-dash: {bound2} + {bound1}", + ha='left', va='bottom', fontsize=12, bbox=dict(fc='w')) + +fig.tight_layout() +# fig.savefig(f"test_flip1.svg", bbox_inches='tight') +fig.savefig(f"{bound1}_{bound2}.svg", bbox_inches='tight') + + diff --git a/docs/source/Gallery/ex16/run.sh b/docs/source/Gallery/ex16/run.sh new file mode 100644 index 0000000..ba2affa --- /dev/null +++ b/docs/source/Gallery/ex16/run.sh @@ -0,0 +1,17 @@ +#!/bin/bash + +set -euo pipefail + +rm -rf *.svg *.tar.gz + +python plot.py free halfspace n +python plot.py rigid halfspace n +python plot.py free free n +python plot.py halfspace halfspace y +python infinite_space.py + +cp free_halfspace.svg cover.svg + +ex=$(basename $(pwd)) +cd .. && tar -czvf ${ex}.tar.gz ${ex} && mv ${ex}.tar.gz ${ex} && cd - + diff --git a/docs/source/Gallery/gallery.rst b/docs/source/Gallery/gallery.rst index 793e1d6..9432971 100644 --- a/docs/source/Gallery/gallery.rst +++ b/docs/source/Gallery/gallery.rst @@ -14,7 +14,7 @@ .. jinja:: - {% for i in range(1, 16) %} + {% for i in range(1, 17) %} {% set i = '%02d' % i %} .. grid-item-card:: :doc:`ex{{i}}/ex{{i}}` :padding: 1 diff --git a/docs/source/Module/explain_-Bbound.rst_ b/docs/source/Module/explain_-Bbound.rst_ new file mode 100644 index 0000000..950bee8 --- /dev/null +++ b/docs/source/Module/explain_-Bbound.rst_ @@ -0,0 +1,7 @@ +.. _-B: + +**-B**\ **f|F|r|R|h|H** + 设置边界条件，小写和大写字母分别控制顶界面和底界面的边界条件， + 分别为自由边界(**f|F**)，刚性边界(**r|R**)和半空间(**h|H**)， + 默认顶界面自由边界，底界面半空间。 + 用法类似 **-BrH** , **-BhF** 等。 diff --git a/docs/source/Module/greenfn.rst b/docs/source/Module/greenfn.rst index 263622d..bff7b0b 100644 --- a/docs/source/Module/greenfn.rst +++ b/docs/source/Module/greenfn.rst @@ -18,6 +18,7 @@ greenfn |-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**] |-R|\ *file*\|\ *r1,r2,...* |-O|\ *outdir* +[ |-B|\ **f|F|r|R|h|H** ] [ |-H|\ *f1/f2* ] [ |-L|\ *length*\ [**+l**\ *Flength*][**+a**\ *Ftol*][**+o**\ *offset*] ] [ |-C|\ **d|p|n** ] @@ -81,6 +82,8 @@ greenfn 可选选项 -------- +.. include:: explain_-Bbound.rst_ + .. _-H: **-H**\ *f1/f2* diff --git a/docs/source/Module/static_greenfn.rst b/docs/source/Module/static_greenfn.rst index 92848e0..ac9f1b6 100644 --- a/docs/source/Module/static_greenfn.rst +++ b/docs/source/Module/static_greenfn.rst @@ -18,6 +18,7 @@ static_greenfn |-X|\ *x1/x2/dx* |-Y|\ *y1/y2/dy* |-O|\ *outdir* +[ |-B|\ **f|F|r|R|h|H** ] [ |-L|\ *length*\ [**+l**\ *Flength*][**+a**\ *Ftol*][**+o**\ *offset*] ] [ |-C|\ **d|p|n** ] [ |-K|\ [**+k**\ *k0*][**+e**\ *keps*] ] @@ -59,6 +60,8 @@ static_greenfn 可选选项 -------- +.. include:: explain_-Bbound.rst_ + .. include:: explain_-L.rst_ .. _-K: diff --git a/docs/source/index.rst b/docs/source/index.rst index ffd6253..8b884ce 100755 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -25,6 +25,7 @@ - :doc:`Advanced/k_integ/index` - :doc:`Advanced/integ_converg/index` - :doc:`Advanced/kernel/kernel` + - :doc:`Advanced/boundary/boundary` - :doc:`Formula/index` - :doc:`Module/index` @@ -65,6 +66,7 @@ Advanced/k_integ/index Advanced/integ_converg/index Advanced/kernel/kernel + Advanced/boundary/boundary Formula/index Module/index diff --git a/test/greenfn/test_greenfn.py b/test/greenfn/test_greenfn.py index d015e4e..0b3cdca 100644 --- a/test/greenfn/test_greenfn.py +++ b/test/greenfn/test_greenfn.py @@ -35,3 +35,12 @@ # multi distances stgrn = pymod.compute_grn([6, 8, 10], nt, dt)[0] + + +# boundary condition +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv, topbound='free', botbound='free') +stgrn = pymod.compute_grn(dist, nt, dt)[0] +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv, topbound='halfspace', botbound='free') +stgrn = pymod.compute_grn(dist, nt, dt)[0] +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv, topbound='rigid', botbound='rigid') +stgrn = pymod.compute_grn(dist, nt, dt)[0] diff --git a/test/greenfn/test_greenfn.sh b/test/greenfn/test_greenfn.sh index 97ecefb..b458ef5 100755 --- a/test/greenfn/test_greenfn.sh +++ b/test/greenfn/test_greenfn.sh @@ -22,6 +22,12 @@ grt greenfn -M../milrow -D0.1/0 -N600/0.02 -R10 -OGRN -Cn grt greenfn -M../milrow -D2/3 -N600/0.02 -L20 -R10 -S -OGRN grt greenfn -M../milrow -D2/3 -N600/0.02 -L20 -R10 -S1,10,20 -OGRN +# boundary +grt greenfn -M../milrow -D2/3 -N600/0.02 -R10 -BrF -OGRN +grt greenfn -M../milrow -D2/3 -N600/0.02 -R10 -BhR -OGRN +grt greenfn -M../milrow -D2/3 -N600/0.02 -R10 -BrH -OGRN + + # multi distance grt greenfn -M../milrow -D2/0 -N600/0.02 -R6,8,10 -OGRN cat > dists <