From 164efb05895ae9c88484e51b720b53e5681f5563 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 8 Jan 2026 12:37:01 +0800 Subject: [PATCH 1/8] DOC: add formula of free and rigid boundary --- docs/source/Formula/bound.rst | 90 +++++ docs/source/Formula/boundary_condition.ipynb | 393 +++++++++++++++++++ docs/source/Formula/index.rst | 2 + 3 files changed, 485 insertions(+) create mode 100644 docs/source/Formula/bound.rst create mode 100644 docs/source/Formula/boundary_condition.ipynb diff --git a/docs/source/Formula/bound.rst b/docs/source/Formula/bound.rst new file mode 100644 index 0000000..acce38e --- /dev/null +++ b/docs/source/Formula/bound.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..ef9e673 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:`bound` **静力学震源参数** @@ -42,6 +43,7 @@ RT_liquid uiz DS_zero + bound static_uniform static_source \ No newline at end of file From fc9744af85b30fc439d534bfd95e94b959716b94 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 8 Jan 2026 12:29:24 +0800 Subject: [PATCH 2/8] support boundary condition setting --- .../include/grt/common/RT_matrix.h | 22 +--- pygrt/C_extension/include/grt/common/const.h | 6 + pygrt/C_extension/include/grt/common/model.h | 15 ++- pygrt/C_extension/include/grt/dynamic/layer.h | 15 ++- .../include/grt/static/static_layer.h | 13 ++- pygrt/C_extension/src/common/RT_matrix.c | 19 +++ pygrt/C_extension/src/common/model.c | 20 ++++ pygrt/C_extension/src/dynamic/grt_greenfn.c | 44 ++++++- pygrt/C_extension/src/dynamic/layer.c | 109 +++++++++++++++--- .../src/integral/kernel_template.c_ | 26 ++++- .../src/static/grt_static_greenfn.c | 43 ++++++- pygrt/C_extension/src/static/static_layer.c | 88 ++++++++++++-- pygrt/c_interfaces.py | 5 + pygrt/c_structures.py | 3 + pygrt/pymod.py | 26 ++++- 15 files changed, 398 insertions(+), 56 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h index d92b01c..ed6478e 100644 --- a/pygrt/C_extension/include/grt/common/RT_matrix.h +++ b/pygrt/C_extension/include/grt/common/RT_matrix.h @@ -33,22 +33,12 @@ typedef struct { } RT_MATRIX; -/** 初始化 R/T 矩阵 */ -#define grt_init_RT_matrix(M) \ - RT_MATRIX *M = &(RT_MATRIX){\ - .RD = GRT_INIT_ZERO_2x2_MATRIX,\ - .RU = GRT_INIT_ZERO_2x2_MATRIX,\ - .TD = GRT_INIT_IDENTITY_2x2_MATRIX,\ - .TU = GRT_INIT_IDENTITY_2x2_MATRIX,\ - .RDL = 0.0,\ - .RUL = 0.0,\ - .TDL = 1.0,\ - .TUL = 1.0,\ - .invT = GRT_INIT_ZERO_2x2_MATRIX,\ - .invTL = 0.0,\ - .stats = GRT_INVERSE_SUCCESS \ - };\ - +/** + * 重新设置 R/T 矩阵为默认值 + * + * @param[out] M RT_MATRIX 指针 + */ +void grt_reset_RT_matrix(RT_MATRIX *M); /** * 合并 recursion_RD_PSV(SH) ，仅计算RD/RDL diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index def2bd5..a645846 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -141,6 +141,12 @@ typedef size_t sizeIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM]; for(int im = 0; im < GRT_SRC_M_NUM; ++im) \ for(int v = 0; v < GRT_INTEG_NUM; ++v) +/** 边界条件 */ +typedef enum { + GRT_BOUND_FREE = 0, + GRT_BOUND_RIGID, + GRT_BOUND_HALFSPACE, +} GRT_BOUND_TYPE; /** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */ enum { diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 6dbfe67..416ad0a 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -58,8 +58,13 @@ typedef struct { RT_MATRIX M_FA; RT_MATRIX M_FB; - /* 自由表面的反射系数矩阵，仅 RU, RUL 有用 */ + + /* 顶界面的反射系数矩阵，仅 RU, RUL 有用 */ + GRT_BOUND_TYPE topbound; RT_MATRIX M_top; + /* 底界面的反射系数矩阵，仅 RD, RDL 有用 */ + GRT_BOUND_TYPE botbound; + RT_MATRIX M_bot; /* 接收点处的接收矩阵 (转为位移u和位移导数uiz的(B_m, C_m, P_m)系分量) */ cplx_t R_EV[2][2]; @@ -146,6 +151,14 @@ void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n); */ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); +/** + * 设置模型的边界条件，并对底界面做检查 + * + * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in] topbound 顶层边界条件 + * @param[in] botbound 底层边界条件 + */ +void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound); /** * 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数 diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index ca4e938..75af04b 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -18,13 +18,22 @@ #include "grt/common/RT_matrix.h" /** - * 计算自由表面的反射系数，公式(5.3.10-14) + * 计算不同边界条件下顶界面的反射系数RU，其中自由表面的公式见(5.3.10-14) * * - * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_topfree_RU(GRT_MODEL1D *mod1d); +void grt_topbound_RU(GRT_MODEL1D *mod1d); + +/** + * 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14) + * + * + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * + */ +void grt_botbound_RD(GRT_MODEL1D *mod1d); /** diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index 9115697..7be2e83 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -20,13 +20,20 @@ #include "grt/common/model.h" /** - * 计算自由表面的静态反射系数，公式(6.3.12) + * 计算不同边界条件下顶界面的静态反射系数RU，其中自由表面的公式见(6.3.12) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_static_topfree_RU(GRT_MODEL1D *mod1d); +void grt_static_topbound_RU(GRT_MODEL1D *mod1d); +/** + * 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12) + * + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * + */ +void grt_static_botbound_RD(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35) diff --git a/pygrt/C_extension/src/common/RT_matrix.c b/pygrt/C_extension/src/common/RT_matrix.c index 9be9896..d0f7d23 100644 --- a/pygrt/C_extension/src/common/RT_matrix.c +++ b/pygrt/C_extension/src/common/RT_matrix.c @@ -13,12 +13,31 @@ #include #include #include +#include #include "grt/common/const.h" #include "grt/common/matrix.h" #include "grt/common/RT_matrix.h" +void grt_reset_RT_matrix(RT_MATRIX *M) +{ + RT_MATRIX *M0 = &(RT_MATRIX){ + .RD = GRT_INIT_ZERO_2x2_MATRIX, + .RU = GRT_INIT_ZERO_2x2_MATRIX, + .TD = GRT_INIT_IDENTITY_2x2_MATRIX, + .TU = GRT_INIT_IDENTITY_2x2_MATRIX, + .RDL = 0.0, + .RUL = 0.0, + .TDL = 1.0, + .TUL = 1.0, + .invT = GRT_INIT_ZERO_2x2_MATRIX, + .invTL = 0.0, + .stats = GRT_INVERSE_SUCCESS + }; + memcpy(M, M0, sizeof(*M)); +} + void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { grt_recursion_RD_PSV(M1, M2, M); diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 853e0ca..c380502 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -421,10 +421,30 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea // 先指定负频率，仅填充弹性模量 grt_attenuate_mod1d(mod1d, -1); + // 设置一个默认边界条件 + mod1d->topbound = GRT_BOUND_FREE; + mod1d->botbound = GRT_BOUND_HALFSPACE; + return mod1d; } +void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound) +{ + mod1d->topbound = topbound; + mod1d->botbound = botbound; + + size_t n = mod1d->n; + bool atBottom = mod1d->isrc+1 == n || mod1d->ircv+1 == n; + if (atBottom && botbound != GRT_BOUND_HALFSPACE){ + GRTRaiseError( + "When source/recevier is at the bottom, the bottom layer must be a halfspace. " + "Please check your boundary condition setting." + ); + } +} + + void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]){ FILE *fp = GRTCheckOpenFile(modelpath, "r"); size_t len; diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index c20e23f..2ccce70 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -54,6 +54,12 @@ typedef struct { char *s_depsrc; char *s_deprcv; } D; + /** 顶层和底层的边界条件 */ + struct { + bool active; + GRT_BOUND_TYPE topbound; + GRT_BOUND_TYPE botbound; + } B; /** 波形时窗 */ struct { bool active; @@ -212,6 +218,9 @@ static void print_Ctrl(const GRT_MODULE_CTRL *Ctrl){ printf(format, "applyDCM", (Ctrl->C.applyDCM)? "true" : "false"); printf(format, "applyPTAM", (Ctrl->C.applyPTAM)? "true" : "false"); + printf(format_size, "topbound", (size_t)Ctrl->B.topbound); + printf(format_size, "botbound", (size_t)Ctrl->B.botbound); + printf(format_size, "nt", Ctrl->N.nt); printf(format_real, "dt", Ctrl->N.dt); printf(format_real, "winT", Ctrl->N.winT); @@ -289,7 +298,7 @@ printf("\n" " -R,[,...] -O [-H/] \n" " [-L] [-C[d|p|n]] [-E[/]] \n" " [-K[+k][+s][+e][+v]]\n" -" [-P] [-Ge|v|h|s] \n" +" [-P] [-Ge|v|h|s] [-Bf|F|r|R|h|H]\n" " [-S[,,...]] [-e] [-s]\n" "\n\n" "Options:\n" @@ -397,6 +406,13 @@ printf("\n" " : Horizontal Force (HF)\n" " : Shear (DC)\n" "\n" +" -Bf|F|r|R|h|H\n" +" Boundary condition of top layer (lowercase) and\n" +" bottom layer (uppercase).\n" +" f|F: Free boundary.\n" +" r|R: Rigid boundary.\n" +" h|H: Halfspace.\n" +"\n" " -S[,,...]\n" " Frequency (index) of statsfile in wavenumber\n" " integration to be output, require 0 <= i <= nf-1,\n" @@ -429,6 +445,9 @@ printf("\n" /** 从命令行中读取选项，处理后记录到全局变量中 */ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 先为个别参数设置非0初始值 + Ctrl->B.topbound = GRT_BOUND_FREE; + Ctrl->B.botbound = GRT_BOUND_HALFSPACE; + Ctrl->N.zeta = GRT_GREENFN_N_ZETA; Ctrl->N.upsample_n = GRT_GREENFN_N_UPSAMPLE; Ctrl->H.freq1 = GRT_GREENFN_H_FREQ1; @@ -441,7 +460,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ Ctrl->G.doDC = GRT_GREENFN_G_DC; int opt; - while ((opt = getopt(argc, argv, ":M:D:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -471,6 +490,24 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 顶层和底层的边界条件 -Bf|F|r|R|h|H + case 'B': + Ctrl->B.active = true; + if(strlen(optarg) == 0 || strlen(optarg) > 2) GRTBadOptionError(B, ""); + for(size_t i = 0; i < strlen(optarg); ++i) { + switch(optarg[i]) { + case 'f': Ctrl->B.topbound = GRT_BOUND_FREE; break; + case 'r': Ctrl->B.topbound = GRT_BOUND_RIGID; break; + case 'h': Ctrl->B.topbound = GRT_BOUND_HALFSPACE; break; + case 'F': Ctrl->B.botbound = GRT_BOUND_FREE; break; + case 'R': Ctrl->B.botbound = GRT_BOUND_RIGID; break; + case 'H': Ctrl->B.botbound = GRT_BOUND_HALFSPACE; break; + default: + GRTBadOptionError(B, "unsupported -B%s.", optarg); + } + } + break; + // 点数,采样间隔,虚频率 -Nnt/dt[+w][+n][+a] case 'N': Ctrl->N.active = true; @@ -870,6 +907,9 @@ int greenfn_main(int argc, char **argv) { } GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + // 边界条件 + grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + // 当震源位于液体层中时，仅允许计算爆炸源对应的格林函数 // 程序结束前会输出对应警告 if(mod1d->Vb[mod1d->isrc]==0.0){ diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index b34cf77..7e31703 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -29,37 +29,118 @@ T N##2 = mod1d->N[iy];\ -void grt_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) { - cplx_t cbcb0 = mod1d->cbcb[0]; if(cbcb0 != 0.0){ // 固体表面 // 公式(5.3.10-14) cplx_t Delta = 0.0; cplx_t cbcb02 = 0.25*cbcb0*cbcb0; - cplx_t xa0 = mod1d->xa[0]; - cplx_t xb0 = mod1d->xb[0]; - // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 Delta = -1.0 + xa0*xb0 + cbcb0 - cbcb02; if(Delta == 0.0){ - mod1d->M_top.stats = GRT_INVERSE_FAILURE; - return; + return GRT_INVERSE_FAILURE; } Delta = 1.0 / Delta; - mod1d->M_top.RU[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; - mod1d->M_top.RU[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][0]; - mod1d->M_top.RUL = 1.0; + M[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; + M[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = 1.0; } else { // 液体表面 - mod1d->M_top.RU[0][0] = -1.0; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][1] = mod1d->M_top.RU[1][0] = 0.0; + M[0][0] = -1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) +{ + if(cbcb0 != 0.0){ + // 固体表面 + cplx_t Delta = 0.0; + + Delta = xa0*xb0 - 1.0; + if(Delta == 0.0){ + return GRT_INVERSE_FAILURE; + } + Delta = 1.0 / Delta; + M[0][0] = (xa0*xb0 + 1.0) * Delta; + M[0][1] = 2.0 * xb0 * Delta * adsgn; + M[1][0] = 2.0 * xa0 * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = - 1.0; + } + else { + // 液体表面 + M[0][0] = 1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +void grt_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t cbcb = mod1d->cbcb[0]; + cplx_t xa = mod1d->xa[0]; + cplx_t xb = mod1d->xb[0]; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t cbcb = mod1d->cbcb[nlay-2]; + cplx_t xa = mod1d->xa[nlay-2]; + cplx_t xb = mod1d->xb[nlay-2]; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + real_t k = mod1d->k; + real_t thk = mod1d->Thk[nlay-2]; + + cplx_t exa, exb, ex2a, ex2b, exab; + exa = exp(- k*thk*xa); + exb = exp(- k*thk*xb); + ex2a = exa * exa; + ex2b = exb * exb; + exab = exa * exb; + + mod1d->M_bot.RD[0][0] *= ex2a; mod1d->M_bot.RD[0][1] *= exab; + mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b; + + mod1d->M_bot.RDL *= ex2b; } diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_ index caca5de..4438412 100644 --- a/pygrt/C_extension/src/integral/kernel_template.c_ +++ b/pygrt/C_extension/src/integral/kernel_template.c_ @@ -47,7 +47,8 @@ #define __grt_RT_matrix_SH CONCAT_PREFIX(RT_matrix_SH) #define __grt_delay_RT_matrix CONCAT_PREFIX(delay_RT_matrix) #define __grt_source_coef CONCAT_PREFIX(source_coef) -#define __grt_topfree_RU CONCAT_PREFIX(topfree_RU) +#define __grt_topbound_RU CONCAT_PREFIX(topbound_RU) +#define __grt_botbound_RD CONCAT_PREFIX(botbound_RD) #define __grt_wave2qwv_REV_PSV CONCAT_PREFIX(wave2qwv_REV_PSV) #define __grt_wave2qwv_REV_SH CONCAT_PREFIX(wave2qwv_REV_SH) #define __grt_wave2qwv_z_REV_PSV CONCAT_PREFIX(wave2qwv_z_REV_PSV) @@ -72,14 +73,22 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) RT_MATRIX *M_BL = &mod1d->M_BL; RT_MATRIX *M_RS = &mod1d->M_RS; RT_MATRIX *M_FA = &mod1d->M_FA; - // 自由表面 + // 顶底界面 RT_MATRIX *M_top = &mod1d->M_top; + RT_MATRIX *M_bot = &mod1d->M_bot; + + grt_reset_RT_matrix(M_BL); + grt_reset_RT_matrix(M_RS); + grt_reset_RT_matrix(M_FA); + grt_reset_RT_matrix(M_top); + grt_reset_RT_matrix(M_bot); // 从顶到底进行矩阵递推, 公式(5.5.3) - for(size_t iy=1; iyn; ++iy){ // 因为n>=3, 故一定会进入该循环 + for(size_t iy = 1; iy < mod1d->n - 1; ++iy){ // 因为n>=3, 故一定会进入该循环 // 定义物理层内的反射透射系数矩阵 - grt_init_RT_matrix(M); + RT_MATRIX *M = &(RT_MATRIX){}; + grt_reset_RT_matrix(M); // 只有动态解才可以使用这个 if 简化， // 对于静态解，即使是震源层、接收层这种虚拟层位也需要显式计算R/T矩阵 @@ -143,11 +152,18 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) //=================================================================================== // 递推RU_FA - __grt_topfree_RU(mod1d); + __grt_topbound_RU(mod1d); if(M_top->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + // 递推RD_BL, 在模型定义阶段已检查震源/台站位于底层时的底层边界条件 + __grt_botbound_RD(mod1d); + if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + grt_recursion_RD(M_BL, M_bot, M_BL); + if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + + BEFORE_RETURN: return; } diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 29a92f9..4b6270d 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -41,6 +41,12 @@ typedef struct { char *s_depsrc; char *s_deprcv; } D; + /** 顶层和底层的边界条件 */ + struct { + bool active; + GRT_BOUND_TYPE topbound; + GRT_BOUND_TYPE botbound; + } B; /** 波数积分间隔以及方法 */ struct { bool active; @@ -143,9 +149,9 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static greenfn -M -D/ -X// \n" +" grt static greenfn -M -D/ -X// \n" " -Y// -O [-L] [-C[d|p|n]] \n" -" [-K[+k][+e]] [-S] [-e]\n" +" [-Bf|F|r|R|h|H] [-K[+k][+e]] [-S] [-e]\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" @@ -198,6 +204,13 @@ printf("\n" " + n: None.\n" " Default use +cd when fabs(depsrc-deprcv) <= %.1f.\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( "\n" +" -Bf|F|r|R|h|H\n" +" Boundary condition of top layer (lowercase) and\n" +" bottom layer (uppercase).\n" +" f|F: Free boundary.\n" +" r|R: Rigid boundary.\n" +" h|H: Halfspace.\n" +"\n" " -K[+k][+e]\n" " Several parameters designed to define the\n" " behavior in wavenumber integration. The upper\n" @@ -236,11 +249,14 @@ printf("\n" /** 从命令行中读取选项，处理后记录到全局变量中 */ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 先为个别参数设置非0初始值 + Ctrl->B.topbound = GRT_BOUND_FREE; + Ctrl->B.botbound = GRT_BOUND_HALFSPACE; + Ctrl->K.k0 = GRT_GREENFN_K_K0; int opt; - while ((opt = getopt(argc, argv, ":M:D:L:C:K:X:Y:O:Seh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:O:Seh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -270,6 +286,24 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 顶层和底层的边界条件 -Bf|F|r|R|h|H + case 'B': + Ctrl->B.active = true; + if(strlen(optarg) == 0 || strlen(optarg) > 2) GRTBadOptionError(B, ""); + for(size_t i = 0; i < strlen(optarg); ++i) { + switch(optarg[i]) { + case 'f': Ctrl->B.topbound = GRT_BOUND_FREE; break; + case 'r': Ctrl->B.topbound = GRT_BOUND_RIGID; break; + case 'h': Ctrl->B.topbound = GRT_BOUND_HALFSPACE; break; + case 'F': Ctrl->B.botbound = GRT_BOUND_FREE; break; + case 'R': Ctrl->B.botbound = GRT_BOUND_RIGID; break; + case 'H': Ctrl->B.botbound = GRT_BOUND_HALFSPACE; break; + default: + GRTBadOptionError(B, "unsupported -B%s.", optarg); + } + } + break; + // 波数积分间隔 -L[][+l][+a][+o] case 'L': Ctrl->L.active = true; @@ -482,6 +516,9 @@ int static_greenfn_main(int argc, char **argv){ } GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + // 边界条件 + grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + // 判断是否要自动使用收敛方法 if( ! Ctrl->C.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) { Ctrl->C.applyDCM = true; diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 687f927..a614a62 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -14,6 +14,7 @@ #include #include #include +#include #include "grt/static/static_layer.h" #include "grt/common/model.h" @@ -25,14 +26,87 @@ T N##2 = mod1d->N[iy];\ -void grt_static_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) { - cplx_t delta1 = mod1d->delta[0]; - // 公式(6.3.12) - mod1d->M_top.RU[0][0] = mod1d->M_top.RU[1][1] = 0.0; - mod1d->M_top.RU[0][1] = -delta1; - mod1d->M_top.RU[1][0] = -1.0/delta1; - mod1d->M_top.RUL = 1.0; + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = -delta1; + M[1][0] = -1.0/delta1; + } else { + M[0][0] = M[1][1] = -2.0*d*k; + M[0][1] = delta1 * (-4.0*d*d*k*k - 1.0); + M[1][0] = -1.0/delta1; + } + *ML = 1.0; + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) +{ + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = M[1][0] = 1.0; + } else { + M[0][0] = M[1][1] = 2.0*delta1*d*k; + M[0][1] = (4.0*GRT_SQUARE(delta1*d*k) + 1.0); + M[1][0] = 1.0; + } + *ML = -1.0; + return GRT_INVERSE_SUCCESS; +} + +void grt_static_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t delta = mod1d->delta[0]; + real_t k = mod1d->k; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); + mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_static_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t delta = mod1d->delta[nlay-2]; + real_t thk = mod1d->Thk[nlay-2]; + real_t k = mod1d->k; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_static_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + cplx_t ex, ex2; + ex = exp(- k*thk); + ex2 = ex * ex; + + mod1d->M_bot.RD[0][0] *= ex2; mod1d->M_bot.RD[0][1] *= ex2; + mod1d->M_bot.RD[1][0] *= ex2; mod1d->M_bot.RD[1][1] *= ex2; + + mod1d->M_bot.RDL *= ex2; } void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 96ae845..62f1ad4 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -85,6 +85,11 @@ def set_num_threads(n): C_grt_read_mod1d_from_file.restype = POINTER(c_GRT_MODEL1D) C_grt_read_mod1d_from_file.argtypes = [c_char_p, c_double, c_double, c_bool] +C_grt_set_mod1d_boundary = libgrt.grt_set_mod1d_boundary +"""设置模型边界条件并检查底界面""" +C_grt_set_mod1d_boundary.restype = None +C_grt_set_mod1d_boundary.argtypes = [POINTER(c_GRT_MODEL1D), c_int, c_int] + C_grt_free_mod1d = libgrt.grt_free_mod1d """释放C程序中申请的 GRT_MODEL1D 结构体内存""" C_grt_free_mod1d.restype = None diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index ce17a51..f34aa83 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -123,7 +123,10 @@ class c_GRT_MODEL1D(Structure): ('M_FA', c_RT_MATRIX), ('M_FB', c_RT_MATRIX), + ('topbound', c_int), ('M_top', c_RT_MATRIX), + ('botbound', c_int), + ('M_bot', c_RT_MATRIX), ('R_EV', CPLX * 4), ('R_EVL', CPLX), diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 253008d..23228d2 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -15,7 +15,7 @@ from obspy import read, Stream, Trace, UTCDateTime from scipy.fft import irfft, ifft from obspy.core import AttribDict -from typing import List, Dict, Union +from typing import List, Dict, Union, Literal import tempfile from time import time @@ -33,7 +33,9 @@ class PyModel1D: - def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False): + def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False, + topbound:Literal['free', 'rigid', 'halfspace']='free', + botbound:Literal['free', 'rigid', 'halfspace']='halfspace'): ''' Create 1D model instance, and insert the imaginary layer of source and receiver. @@ -41,12 +43,27 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b :param depsrc: source depth (km) :param deprcv: receiver depth (km) :param allowLiquid: whether liquid layers are allowed + :param topbound: boundary condition of the top layer + :param botbound: boundary condition of the bottom layer ''' self.depsrc:float = depsrc self.deprcv:float = deprcv self.c_mod1d:c_GRT_MODEL1D self.hasLiquid:bool = allowLiquid # 传入的模型是否有液体层 + self.topbound:str = topbound + self.botbound:str = botbound + + boundDct = { + 'free': 0, + 'rigid': 1, + 'halfspace': 2 + } + + if topbound not in boundDct: + raise ValueError(f"Unsupported topbound={topbound}.") + if botbound not in boundDct: + raise ValueError(f"Unsupported botbound={botbound}.") # 将modarr写入临时数组 with tempfile.NamedTemporaryFile(mode='w', delete=False) as tmpfile: @@ -57,10 +74,15 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b try: c_mod1d_ptr = C_grt_read_mod1d_from_file(tmp_path.encode("utf-8"), depsrc, deprcv, allowLiquid) self.c_mod1d = c_mod1d_ptr.contents # 这部分内存在C中申请，需由C函数释放。占用不多，这里跳过 + C_grt_set_mod1d_boundary(self.c_mod1d, boundDct[topbound], boundDct[botbound]) finally: if os.path.exists(tmp_path): os.unlink(tmp_path) + # 设置边界条件 + self.c_mod1d.topbound = boundDct[topbound] + self.c_mod1d.botbound = boundDct[botbound] + self.isrc = self.c_mod1d.isrc self.ircv = self.c_mod1d.ircv From f3ac52ce1b7fbee25c454ed5c441849422030c9e Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 8 Jan 2026 12:29:24 +0800 Subject: [PATCH 3/8] support boundary condition setting --- .../include/grt/common/RT_matrix.h | 22 +--- pygrt/C_extension/include/grt/common/const.h | 6 + pygrt/C_extension/include/grt/common/model.h | 15 ++- pygrt/C_extension/include/grt/dynamic/layer.h | 15 ++- .../include/grt/static/static_layer.h | 13 ++- pygrt/C_extension/src/common/RT_matrix.c | 19 +++ pygrt/C_extension/src/common/model.c | 20 ++++ pygrt/C_extension/src/dynamic/grt_greenfn.c | 44 ++++++- pygrt/C_extension/src/dynamic/layer.c | 109 +++++++++++++++--- .../src/integral/kernel_template.c_ | 26 ++++- .../src/static/grt_static_greenfn.c | 43 ++++++- pygrt/C_extension/src/static/static_layer.c | 88 ++++++++++++-- pygrt/c_interfaces.py | 5 + pygrt/c_structures.py | 3 + pygrt/pymod.py | 26 ++++- 15 files changed, 398 insertions(+), 56 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h index d92b01c..ed6478e 100644 --- a/pygrt/C_extension/include/grt/common/RT_matrix.h +++ b/pygrt/C_extension/include/grt/common/RT_matrix.h @@ -33,22 +33,12 @@ typedef struct { } RT_MATRIX; -/** 初始化 R/T 矩阵 */ -#define grt_init_RT_matrix(M) \ - RT_MATRIX *M = &(RT_MATRIX){\ - .RD = GRT_INIT_ZERO_2x2_MATRIX,\ - .RU = GRT_INIT_ZERO_2x2_MATRIX,\ - .TD = GRT_INIT_IDENTITY_2x2_MATRIX,\ - .TU = GRT_INIT_IDENTITY_2x2_MATRIX,\ - .RDL = 0.0,\ - .RUL = 0.0,\ - .TDL = 1.0,\ - .TUL = 1.0,\ - .invT = GRT_INIT_ZERO_2x2_MATRIX,\ - .invTL = 0.0,\ - .stats = GRT_INVERSE_SUCCESS \ - };\ - +/** + * 重新设置 R/T 矩阵为默认值 + * + * @param[out] M RT_MATRIX 指针 + */ +void grt_reset_RT_matrix(RT_MATRIX *M); /** * 合并 recursion_RD_PSV(SH) ，仅计算RD/RDL diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index def2bd5..a645846 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -141,6 +141,12 @@ typedef size_t sizeIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM]; for(int im = 0; im < GRT_SRC_M_NUM; ++im) \ for(int v = 0; v < GRT_INTEG_NUM; ++v) +/** 边界条件 */ +typedef enum { + GRT_BOUND_FREE = 0, + GRT_BOUND_RIGID, + GRT_BOUND_HALFSPACE, +} GRT_BOUND_TYPE; /** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */ enum { diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 6dbfe67..416ad0a 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -58,8 +58,13 @@ typedef struct { RT_MATRIX M_FA; RT_MATRIX M_FB; - /* 自由表面的反射系数矩阵，仅 RU, RUL 有用 */ + + /* 顶界面的反射系数矩阵，仅 RU, RUL 有用 */ + GRT_BOUND_TYPE topbound; RT_MATRIX M_top; + /* 底界面的反射系数矩阵，仅 RD, RDL 有用 */ + GRT_BOUND_TYPE botbound; + RT_MATRIX M_bot; /* 接收点处的接收矩阵 (转为位移u和位移导数uiz的(B_m, C_m, P_m)系分量) */ cplx_t R_EV[2][2]; @@ -146,6 +151,14 @@ void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n); */ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); +/** + * 设置模型的边界条件，并对底界面做检查 + * + * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in] topbound 顶层边界条件 + * @param[in] botbound 底层边界条件 + */ +void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound); /** * 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数 diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index ca4e938..75af04b 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -18,13 +18,22 @@ #include "grt/common/RT_matrix.h" /** - * 计算自由表面的反射系数，公式(5.3.10-14) + * 计算不同边界条件下顶界面的反射系数RU，其中自由表面的公式见(5.3.10-14) * * - * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_topfree_RU(GRT_MODEL1D *mod1d); +void grt_topbound_RU(GRT_MODEL1D *mod1d); + +/** + * 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14) + * + * + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * + */ +void grt_botbound_RD(GRT_MODEL1D *mod1d); /** diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index 9115697..7be2e83 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -20,13 +20,20 @@ #include "grt/common/model.h" /** - * 计算自由表面的静态反射系数，公式(6.3.12) + * 计算不同边界条件下顶界面的静态反射系数RU，其中自由表面的公式见(6.3.12) * - * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_top * */ -void grt_static_topfree_RU(GRT_MODEL1D *mod1d); +void grt_static_topbound_RU(GRT_MODEL1D *mod1d); +/** + * 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12) + * + * @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot + * + */ +void grt_static_botbound_RD(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35) diff --git a/pygrt/C_extension/src/common/RT_matrix.c b/pygrt/C_extension/src/common/RT_matrix.c index 9be9896..d0f7d23 100644 --- a/pygrt/C_extension/src/common/RT_matrix.c +++ b/pygrt/C_extension/src/common/RT_matrix.c @@ -13,12 +13,31 @@ #include #include #include +#include #include "grt/common/const.h" #include "grt/common/matrix.h" #include "grt/common/RT_matrix.h" +void grt_reset_RT_matrix(RT_MATRIX *M) +{ + RT_MATRIX *M0 = &(RT_MATRIX){ + .RD = GRT_INIT_ZERO_2x2_MATRIX, + .RU = GRT_INIT_ZERO_2x2_MATRIX, + .TD = GRT_INIT_IDENTITY_2x2_MATRIX, + .TU = GRT_INIT_IDENTITY_2x2_MATRIX, + .RDL = 0.0, + .RUL = 0.0, + .TDL = 1.0, + .TUL = 1.0, + .invT = GRT_INIT_ZERO_2x2_MATRIX, + .invTL = 0.0, + .stats = GRT_INVERSE_SUCCESS + }; + memcpy(M, M0, sizeof(*M)); +} + void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M) { grt_recursion_RD_PSV(M1, M2, M); diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 853e0ca..c380502 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -421,10 +421,30 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea // 先指定负频率，仅填充弹性模量 grt_attenuate_mod1d(mod1d, -1); + // 设置一个默认边界条件 + mod1d->topbound = GRT_BOUND_FREE; + mod1d->botbound = GRT_BOUND_HALFSPACE; + return mod1d; } +void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound) +{ + mod1d->topbound = topbound; + mod1d->botbound = botbound; + + size_t n = mod1d->n; + bool atBottom = mod1d->isrc+1 == n || mod1d->ircv+1 == n; + if (atBottom && botbound != GRT_BOUND_HALFSPACE){ + GRTRaiseError( + "When source/recevier is at the bottom, the bottom layer must be a halfspace. " + "Please check your boundary condition setting." + ); + } +} + + void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]){ FILE *fp = GRTCheckOpenFile(modelpath, "r"); size_t len; diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index c20e23f..2ccce70 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -54,6 +54,12 @@ typedef struct { char *s_depsrc; char *s_deprcv; } D; + /** 顶层和底层的边界条件 */ + struct { + bool active; + GRT_BOUND_TYPE topbound; + GRT_BOUND_TYPE botbound; + } B; /** 波形时窗 */ struct { bool active; @@ -212,6 +218,9 @@ static void print_Ctrl(const GRT_MODULE_CTRL *Ctrl){ printf(format, "applyDCM", (Ctrl->C.applyDCM)? "true" : "false"); printf(format, "applyPTAM", (Ctrl->C.applyPTAM)? "true" : "false"); + printf(format_size, "topbound", (size_t)Ctrl->B.topbound); + printf(format_size, "botbound", (size_t)Ctrl->B.botbound); + printf(format_size, "nt", Ctrl->N.nt); printf(format_real, "dt", Ctrl->N.dt); printf(format_real, "winT", Ctrl->N.winT); @@ -289,7 +298,7 @@ printf("\n" " -R,[,...] -O [-H/] \n" " [-L] [-C[d|p|n]] [-E[/]] \n" " [-K[+k][+s][+e][+v]]\n" -" [-P] [-Ge|v|h|s] \n" +" [-P] [-Ge|v|h|s] [-Bf|F|r|R|h|H]\n" " [-S[,,...]] [-e] [-s]\n" "\n\n" "Options:\n" @@ -397,6 +406,13 @@ printf("\n" " : Horizontal Force (HF)\n" " : Shear (DC)\n" "\n" +" -Bf|F|r|R|h|H\n" +" Boundary condition of top layer (lowercase) and\n" +" bottom layer (uppercase).\n" +" f|F: Free boundary.\n" +" r|R: Rigid boundary.\n" +" h|H: Halfspace.\n" +"\n" " -S[,,...]\n" " Frequency (index) of statsfile in wavenumber\n" " integration to be output, require 0 <= i <= nf-1,\n" @@ -429,6 +445,9 @@ printf("\n" /** 从命令行中读取选项，处理后记录到全局变量中 */ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 先为个别参数设置非0初始值 + Ctrl->B.topbound = GRT_BOUND_FREE; + Ctrl->B.botbound = GRT_BOUND_HALFSPACE; + Ctrl->N.zeta = GRT_GREENFN_N_ZETA; Ctrl->N.upsample_n = GRT_GREENFN_N_UPSAMPLE; Ctrl->H.freq1 = GRT_GREENFN_H_FREQ1; @@ -441,7 +460,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ Ctrl->G.doDC = GRT_GREENFN_G_DC; int opt; - while ((opt = getopt(argc, argv, ":M:D:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -471,6 +490,24 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 顶层和底层的边界条件 -Bf|F|r|R|h|H + case 'B': + Ctrl->B.active = true; + if(strlen(optarg) == 0 || strlen(optarg) > 2) GRTBadOptionError(B, ""); + for(size_t i = 0; i < strlen(optarg); ++i) { + switch(optarg[i]) { + case 'f': Ctrl->B.topbound = GRT_BOUND_FREE; break; + case 'r': Ctrl->B.topbound = GRT_BOUND_RIGID; break; + case 'h': Ctrl->B.topbound = GRT_BOUND_HALFSPACE; break; + case 'F': Ctrl->B.botbound = GRT_BOUND_FREE; break; + case 'R': Ctrl->B.botbound = GRT_BOUND_RIGID; break; + case 'H': Ctrl->B.botbound = GRT_BOUND_HALFSPACE; break; + default: + GRTBadOptionError(B, "unsupported -B%s.", optarg); + } + } + break; + // 点数,采样间隔,虚频率 -Nnt/dt[+w][+n][+a] case 'N': Ctrl->N.active = true; @@ -870,6 +907,9 @@ int greenfn_main(int argc, char **argv) { } GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + // 边界条件 + grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + // 当震源位于液体层中时，仅允许计算爆炸源对应的格林函数 // 程序结束前会输出对应警告 if(mod1d->Vb[mod1d->isrc]==0.0){ diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index b34cf77..7e31703 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -29,37 +29,118 @@ T N##2 = mod1d->N[iy];\ -void grt_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) { - cplx_t cbcb0 = mod1d->cbcb[0]; if(cbcb0 != 0.0){ // 固体表面 // 公式(5.3.10-14) cplx_t Delta = 0.0; cplx_t cbcb02 = 0.25*cbcb0*cbcb0; - cplx_t xa0 = mod1d->xa[0]; - cplx_t xb0 = mod1d->xb[0]; - // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 Delta = -1.0 + xa0*xb0 + cbcb0 - cbcb02; if(Delta == 0.0){ - mod1d->M_top.stats = GRT_INVERSE_FAILURE; - return; + return GRT_INVERSE_FAILURE; } Delta = 1.0 / Delta; - mod1d->M_top.RU[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; - mod1d->M_top.RU[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][0]; - mod1d->M_top.RUL = 1.0; + M[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; + M[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = 1.0; } else { // 液体表面 - mod1d->M_top.RU[0][0] = -1.0; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][1] = mod1d->M_top.RU[1][0] = 0.0; + M[0][0] = -1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) +{ + if(cbcb0 != 0.0){ + // 固体表面 + cplx_t Delta = 0.0; + + Delta = xa0*xb0 - 1.0; + if(Delta == 0.0){ + return GRT_INVERSE_FAILURE; + } + Delta = 1.0 / Delta; + M[0][0] = (xa0*xb0 + 1.0) * Delta; + M[0][1] = 2.0 * xb0 * Delta * adsgn; + M[1][0] = 2.0 * xa0 * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = - 1.0; + } + else { + // 液体表面 + M[0][0] = 1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +void grt_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t cbcb = mod1d->cbcb[0]; + cplx_t xa = mod1d->xa[0]; + cplx_t xb = mod1d->xb[0]; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t cbcb = mod1d->cbcb[nlay-2]; + cplx_t xa = mod1d->xa[nlay-2]; + cplx_t xb = mod1d->xb[nlay-2]; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + real_t k = mod1d->k; + real_t thk = mod1d->Thk[nlay-2]; + + cplx_t exa, exb, ex2a, ex2b, exab; + exa = exp(- k*thk*xa); + exb = exp(- k*thk*xb); + ex2a = exa * exa; + ex2b = exb * exb; + exab = exa * exb; + + mod1d->M_bot.RD[0][0] *= ex2a; mod1d->M_bot.RD[0][1] *= exab; + mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b; + + mod1d->M_bot.RDL *= ex2b; } diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_ index caca5de..4438412 100644 --- a/pygrt/C_extension/src/integral/kernel_template.c_ +++ b/pygrt/C_extension/src/integral/kernel_template.c_ @@ -47,7 +47,8 @@ #define __grt_RT_matrix_SH CONCAT_PREFIX(RT_matrix_SH) #define __grt_delay_RT_matrix CONCAT_PREFIX(delay_RT_matrix) #define __grt_source_coef CONCAT_PREFIX(source_coef) -#define __grt_topfree_RU CONCAT_PREFIX(topfree_RU) +#define __grt_topbound_RU CONCAT_PREFIX(topbound_RU) +#define __grt_botbound_RD CONCAT_PREFIX(botbound_RD) #define __grt_wave2qwv_REV_PSV CONCAT_PREFIX(wave2qwv_REV_PSV) #define __grt_wave2qwv_REV_SH CONCAT_PREFIX(wave2qwv_REV_SH) #define __grt_wave2qwv_z_REV_PSV CONCAT_PREFIX(wave2qwv_z_REV_PSV) @@ -72,14 +73,22 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) RT_MATRIX *M_BL = &mod1d->M_BL; RT_MATRIX *M_RS = &mod1d->M_RS; RT_MATRIX *M_FA = &mod1d->M_FA; - // 自由表面 + // 顶底界面 RT_MATRIX *M_top = &mod1d->M_top; + RT_MATRIX *M_bot = &mod1d->M_bot; + + grt_reset_RT_matrix(M_BL); + grt_reset_RT_matrix(M_RS); + grt_reset_RT_matrix(M_FA); + grt_reset_RT_matrix(M_top); + grt_reset_RT_matrix(M_bot); // 从顶到底进行矩阵递推, 公式(5.5.3) - for(size_t iy=1; iyn; ++iy){ // 因为n>=3, 故一定会进入该循环 + for(size_t iy = 1; iy < mod1d->n - 1; ++iy){ // 因为n>=3, 故一定会进入该循环 // 定义物理层内的反射透射系数矩阵 - grt_init_RT_matrix(M); + RT_MATRIX *M = &(RT_MATRIX){}; + grt_reset_RT_matrix(M); // 只有动态解才可以使用这个 if 简化， // 对于静态解，即使是震源层、接收层这种虚拟层位也需要显式计算R/T矩阵 @@ -143,11 +152,18 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) //=================================================================================== // 递推RU_FA - __grt_topfree_RU(mod1d); + __grt_topbound_RU(mod1d); if(M_top->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + // 递推RD_BL, 在模型定义阶段已检查震源/台站位于底层时的底层边界条件 + __grt_botbound_RD(mod1d); + if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + grt_recursion_RD(M_BL, M_bot, M_BL); + if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + + BEFORE_RETURN: return; } diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 29a92f9..4b6270d 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -41,6 +41,12 @@ typedef struct { char *s_depsrc; char *s_deprcv; } D; + /** 顶层和底层的边界条件 */ + struct { + bool active; + GRT_BOUND_TYPE topbound; + GRT_BOUND_TYPE botbound; + } B; /** 波数积分间隔以及方法 */ struct { bool active; @@ -143,9 +149,9 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static greenfn -M -D/ -X// \n" +" grt static greenfn -M -D/ -X// \n" " -Y// -O [-L] [-C[d|p|n]] \n" -" [-K[+k][+e]] [-S] [-e]\n" +" [-Bf|F|r|R|h|H] [-K[+k][+e]] [-S] [-e]\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" @@ -198,6 +204,13 @@ printf("\n" " + n: None.\n" " Default use +cd when fabs(depsrc-deprcv) <= %.1f.\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( "\n" +" -Bf|F|r|R|h|H\n" +" Boundary condition of top layer (lowercase) and\n" +" bottom layer (uppercase).\n" +" f|F: Free boundary.\n" +" r|R: Rigid boundary.\n" +" h|H: Halfspace.\n" +"\n" " -K[+k][+e]\n" " Several parameters designed to define the\n" " behavior in wavenumber integration. The upper\n" @@ -236,11 +249,14 @@ printf("\n" /** 从命令行中读取选项，处理后记录到全局变量中 */ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 先为个别参数设置非0初始值 + Ctrl->B.topbound = GRT_BOUND_FREE; + Ctrl->B.botbound = GRT_BOUND_HALFSPACE; + Ctrl->K.k0 = GRT_GREENFN_K_K0; int opt; - while ((opt = getopt(argc, argv, ":M:D:L:C:K:X:Y:O:Seh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:O:Seh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -270,6 +286,24 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 顶层和底层的边界条件 -Bf|F|r|R|h|H + case 'B': + Ctrl->B.active = true; + if(strlen(optarg) == 0 || strlen(optarg) > 2) GRTBadOptionError(B, ""); + for(size_t i = 0; i < strlen(optarg); ++i) { + switch(optarg[i]) { + case 'f': Ctrl->B.topbound = GRT_BOUND_FREE; break; + case 'r': Ctrl->B.topbound = GRT_BOUND_RIGID; break; + case 'h': Ctrl->B.topbound = GRT_BOUND_HALFSPACE; break; + case 'F': Ctrl->B.botbound = GRT_BOUND_FREE; break; + case 'R': Ctrl->B.botbound = GRT_BOUND_RIGID; break; + case 'H': Ctrl->B.botbound = GRT_BOUND_HALFSPACE; break; + default: + GRTBadOptionError(B, "unsupported -B%s.", optarg); + } + } + break; + // 波数积分间隔 -L[][+l][+a][+o] case 'L': Ctrl->L.active = true; @@ -482,6 +516,9 @@ int static_greenfn_main(int argc, char **argv){ } GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + // 边界条件 + grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + // 判断是否要自动使用收敛方法 if( ! Ctrl->C.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) { Ctrl->C.applyDCM = true; diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 687f927..a614a62 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -14,6 +14,7 @@ #include #include #include +#include #include "grt/static/static_layer.h" #include "grt/common/model.h" @@ -25,14 +26,87 @@ T N##2 = mod1d->N[iy];\ -void grt_static_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) { - cplx_t delta1 = mod1d->delta[0]; - // 公式(6.3.12) - mod1d->M_top.RU[0][0] = mod1d->M_top.RU[1][1] = 0.0; - mod1d->M_top.RU[0][1] = -delta1; - mod1d->M_top.RU[1][0] = -1.0/delta1; - mod1d->M_top.RUL = 1.0; + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = -delta1; + M[1][0] = -1.0/delta1; + } else { + M[0][0] = M[1][1] = -2.0*d*k; + M[0][1] = delta1 * (-4.0*d*d*k*k - 1.0); + M[1][0] = -1.0/delta1; + } + *ML = 1.0; + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) +{ + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = M[1][0] = 1.0; + } else { + M[0][0] = M[1][1] = 2.0*delta1*d*k; + M[0][1] = (4.0*GRT_SQUARE(delta1*d*k) + 1.0); + M[1][0] = 1.0; + } + *ML = -1.0; + return GRT_INVERSE_SUCCESS; +} + +void grt_static_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t delta = mod1d->delta[0]; + real_t k = mod1d->k; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); + mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_static_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t delta = mod1d->delta[nlay-2]; + real_t thk = mod1d->Thk[nlay-2]; + real_t k = mod1d->k; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_static_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + cplx_t ex, ex2; + ex = exp(- k*thk); + ex2 = ex * ex; + + mod1d->M_bot.RD[0][0] *= ex2; mod1d->M_bot.RD[0][1] *= ex2; + mod1d->M_bot.RD[1][0] *= ex2; mod1d->M_bot.RD[1][1] *= ex2; + + mod1d->M_bot.RDL *= ex2; } void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 96ae845..62f1ad4 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -85,6 +85,11 @@ def set_num_threads(n): C_grt_read_mod1d_from_file.restype = POINTER(c_GRT_MODEL1D) C_grt_read_mod1d_from_file.argtypes = [c_char_p, c_double, c_double, c_bool] +C_grt_set_mod1d_boundary = libgrt.grt_set_mod1d_boundary +"""设置模型边界条件并检查底界面""" +C_grt_set_mod1d_boundary.restype = None +C_grt_set_mod1d_boundary.argtypes = [POINTER(c_GRT_MODEL1D), c_int, c_int] + C_grt_free_mod1d = libgrt.grt_free_mod1d """释放C程序中申请的 GRT_MODEL1D 结构体内存""" C_grt_free_mod1d.restype = None diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index ce17a51..f34aa83 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -123,7 +123,10 @@ class c_GRT_MODEL1D(Structure): ('M_FA', c_RT_MATRIX), ('M_FB', c_RT_MATRIX), + ('topbound', c_int), ('M_top', c_RT_MATRIX), + ('botbound', c_int), + ('M_bot', c_RT_MATRIX), ('R_EV', CPLX * 4), ('R_EVL', CPLX), diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 253008d..23228d2 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -15,7 +15,7 @@ from obspy import read, Stream, Trace, UTCDateTime from scipy.fft import irfft, ifft from obspy.core import AttribDict -from typing import List, Dict, Union +from typing import List, Dict, Union, Literal import tempfile from time import time @@ -33,7 +33,9 @@ class PyModel1D: - def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False): + def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False, + topbound:Literal['free', 'rigid', 'halfspace']='free', + botbound:Literal['free', 'rigid', 'halfspace']='halfspace'): ''' Create 1D model instance, and insert the imaginary layer of source and receiver. @@ -41,12 +43,27 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b :param depsrc: source depth (km) :param deprcv: receiver depth (km) :param allowLiquid: whether liquid layers are allowed + :param topbound: boundary condition of the top layer + :param botbound: boundary condition of the bottom layer ''' self.depsrc:float = depsrc self.deprcv:float = deprcv self.c_mod1d:c_GRT_MODEL1D self.hasLiquid:bool = allowLiquid # 传入的模型是否有液体层 + self.topbound:str = topbound + self.botbound:str = botbound + + boundDct = { + 'free': 0, + 'rigid': 1, + 'halfspace': 2 + } + + if topbound not in boundDct: + raise ValueError(f"Unsupported topbound={topbound}.") + if botbound not in boundDct: + raise ValueError(f"Unsupported botbound={botbound}.") # 将modarr写入临时数组 with tempfile.NamedTemporaryFile(mode='w', delete=False) as tmpfile: @@ -57,10 +74,15 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b try: c_mod1d_ptr = C_grt_read_mod1d_from_file(tmp_path.encode("utf-8"), depsrc, deprcv, allowLiquid) self.c_mod1d = c_mod1d_ptr.contents # 这部分内存在C中申请，需由C函数释放。占用不多，这里跳过 + C_grt_set_mod1d_boundary(self.c_mod1d, boundDct[topbound], boundDct[botbound]) finally: if os.path.exists(tmp_path): os.unlink(tmp_path) + # 设置边界条件 + self.c_mod1d.topbound = boundDct[topbound] + self.c_mod1d.botbound = boundDct[botbound] + self.isrc = self.c_mod1d.isrc self.ircv = self.c_mod1d.ircv From 0a79bf50acb5f68f45760a2bb6afbcc7459705d4 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 9 Jan 2026 14:45:32 +0800 Subject: [PATCH 4/8] update --- docs/source/Formula/{bound.rst => boundary.rst} | 0 docs/source/Formula/index.rst | 4 ++-- 2 files changed, 2 insertions(+), 2 deletions(-) rename docs/source/Formula/{bound.rst => boundary.rst} (100%) diff --git a/docs/source/Formula/bound.rst b/docs/source/Formula/boundary.rst similarity index 100% rename from docs/source/Formula/bound.rst rename to docs/source/Formula/boundary.rst diff --git a/docs/source/Formula/index.rst b/docs/source/Formula/index.rst index ef9e673..d1667b1 100644 --- a/docs/source/Formula/index.rst +++ b/docs/source/Formula/index.rst @@ -22,7 +22,7 @@ - :doc:`RT_liquid` - :doc:`uiz` - :doc:`DS_zero` - - :doc:`bound` + - :doc:`boundary` **静力学震源参数** @@ -43,7 +43,7 @@ RT_liquid uiz DS_zero - bound + boundary static_uniform static_source \ No newline at end of file From bb0d8e94c099c51822ea8f10b07f1f911f223831 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 9 Jan 2026 16:59:38 +0800 Subject: [PATCH 5/8] update --- pygrt/C_extension/src/common/model.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index c380502..64378b3 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -435,7 +435,7 @@ void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOU mod1d->botbound = botbound; size_t n = mod1d->n; - bool atBottom = mod1d->isrc+1 == n || mod1d->ircv+1 == n; + bool atBottom = mod1d->isrc+1 >= n || mod1d->ircv+1 >= n; if (atBottom && botbound != GRT_BOUND_HALFSPACE){ GRTRaiseError( "When source/recevier is at the bottom, the bottom layer must be a halfspace. " From 29fce46316fdb0dcb524200cb2011a921288bead Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 9 Jan 2026 20:43:05 +0800 Subject: [PATCH 6/8] update --- .../C_extension/src/integral/kernel_template.c_ | 17 ++++++++++------- pygrt/pymod.py | 3 +-- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_ index 4438412..1448a52 100644 --- a/pygrt/C_extension/src/integral/kernel_template.c_ +++ b/pygrt/C_extension/src/integral/kernel_template.c_ @@ -84,8 +84,9 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) grt_reset_RT_matrix(M_bot); // 从顶到底进行矩阵递推, 公式(5.5.3) - for(size_t iy = 1; iy < mod1d->n - 1; ++iy){ // 因为n>=3, 故一定会进入该循环 - + size_t maxidx = (imax+1 == mod1d->n)? mod1d->n : mod1d->n - 1; + for(size_t iy = 1; iy < maxidx; ++iy) // 因为n>=3, 故一定会进入该循环 + { // 定义物理层内的反射透射系数矩阵 RT_MATRIX *M = &(RT_MATRIX){}; grt_reset_RT_matrix(M); @@ -157,11 +158,13 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; - // 递推RD_BL, 在模型定义阶段已检查震源/台站位于底层时的底层边界条件 - __grt_botbound_RD(mod1d); - if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; - grt_recursion_RD(M_BL, M_bot, M_BL); - if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + // 递推RD_BL, 在模型定义阶段已检查震源/台站的虚拟层成为底层时，只能是半无限空间 + if(imax+1 < mod1d->n){ // 此 if 仅减少不必要的计算 + __grt_botbound_RD(mod1d); + if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + grt_recursion_RD(M_BL, M_bot, M_BL); + if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + } BEFORE_RETURN: diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 23228d2..37398dc 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -182,7 +182,6 @@ def _get_grn_spectra( delayT0:float=0.0, delayV0:float=0.0, calc_upar:bool=False, - gf_source=['EX', 'VF', 'HF', 'DC'], statsfile:Union[str,None]=None, statsidxs:Union[np.ndarray,List[int],None]=None, print_runtime:bool=True @@ -461,7 +460,7 @@ def compute_grn( pygrnLst, pygrnLst_uiz, pygrnLst_uir = self._get_grn_spectra( distarr, nt, dt, upsampling_n, freqband, zeta, keepAllFreq, vmin_ref, keps, ampk, k0, Length, filonLength, safilonTol, filonCut, converg_method, - delayT0, delayV0, calc_upar, gf_source, + delayT0, delayV0, calc_upar, statsfile, statsidxs, print_runtime ) From c5a01f86dfc0477765c438477a385ef68fb88f47 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Sun, 11 Jan 2026 01:41:31 +0800 Subject: [PATCH 7/8] TEST: add ex16, boundary condition test --- docs/source/Gallery/ex16/ex16.rst | 50 ++ docs/source/Gallery/ex16/flip1.SVG | 523 +++++++++++++++++++++ docs/source/Gallery/ex16/infinite_space.py | 120 +++++ docs/source/Gallery/ex16/plot.py | 117 +++++ docs/source/Gallery/ex16/run.sh | 17 + docs/source/Gallery/gallery.rst | 2 +- test/greenfn/test_greenfn.py | 9 + test/greenfn/test_greenfn.sh | 6 + test/static_greenfn/test_static_greenfn.py | 8 + test/static_greenfn/test_static_greenfn.sh | 5 + 10 files changed, 856 insertions(+), 1 deletion(-) create mode 100644 docs/source/Gallery/ex16/ex16.rst create mode 100644 docs/source/Gallery/ex16/flip1.SVG create mode 100644 docs/source/Gallery/ex16/infinite_space.py create mode 100644 docs/source/Gallery/ex16/plot.py create mode 100644 docs/source/Gallery/ex16/run.sh 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/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 < Date: Sun, 11 Jan 2026 02:08:07 +0800 Subject: [PATCH 8/8] DOC: boundary condition --- docs/source/Advanced/boundary/boundary.rst | 29 + .../Advanced/boundary/plot_boundary.svg | 767 ++++++++++++++++++ docs/source/Module/explain_-Bbound.rst_ | 7 + docs/source/Module/greenfn.rst | 3 + docs/source/Module/static_greenfn.rst | 3 + docs/source/index.rst | 2 + 6 files changed, 811 insertions(+) create mode 100644 docs/source/Advanced/boundary/boundary.rst create mode 100644 docs/source/Advanced/boundary/plot_boundary.svg create mode 100644 docs/source/Module/explain_-Bbound.rst_ 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/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