From 7dbf374dbaaa95b8c8c6cfc8cfcc83cada29634c Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 14 Nov 2025 08:51:43 +0800 Subject: [PATCH] DOC: add lamb1 (the first-kind Lamb problem) --- docs/source/Lamb_problem/index.rst | 28 ++++++ docs/source/Lamb_problem/lamb1.rst | 54 ++++++++++++ .../Lamb_problem/run/lamb1_plot_freq_time.py | 86 +++++++++++++++++++ .../Lamb_problem/run/lamb1_plot_time.py | 21 +++++ docs/source/Lamb_problem/run/run.sh | 9 ++ docs/source/Module/index.rst | 2 + docs/source/Module/lamb1.rst | 50 +++++++++++ docs/source/index.rst | 1 + docs/source/run_all.sh | 9 ++ 9 files changed, 260 insertions(+) create mode 100644 docs/source/Lamb_problem/index.rst create mode 100644 docs/source/Lamb_problem/lamb1.rst create mode 100644 docs/source/Lamb_problem/run/lamb1_plot_freq_time.py create mode 100644 docs/source/Lamb_problem/run/lamb1_plot_time.py create mode 100755 docs/source/Lamb_problem/run/run.sh create mode 100644 docs/source/Module/lamb1.rst diff --git a/docs/source/Lamb_problem/index.rst b/docs/source/Lamb_problem/index.rst new file mode 100644 index 0000000..d1f2d96 --- /dev/null +++ b/docs/source/Lamb_problem/index.rst @@ -0,0 +1,28 @@ +Lamb 问题 +================ + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +北京大学张海明教授已在其编著的 **《地震学中的 Lamb 问题》上下两册** 中对三类 Lamb 问题进行了详细论述， +并分别在上下两册中给出了频域解和时域解，过程之详细令人叹为观止。这里不再对 Lamb 问题及其解法做过多介绍， +详见张海明老师的书（强烈推荐！！） + +**PyGRT** 程序包本身是在频域求解，这里我对照下册书也编程实现了 Lamb 问题的时域解（目前包括第一类 Lamb 问题，后续会再扩展）， +方便进行对比验证。以下简单示范用法以及绘图。 + +.. note:: + + 程序输入的时间为无量纲时间，输出的位移为和阶跃函数卷积后的无量纲位移。 + + + 无量纲时间 :math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}` + + 无量纲位移 :math:`\bar{\mathbf{G}}^H = \pi^2 \mu r \mathbf{G}^H` ，上标 :math:`H` 表示已和阶跃函数卷积。 + + +.. toctree:: + :maxdepth: 1 + + lamb1 + diff --git a/docs/source/Lamb_problem/lamb1.rst b/docs/source/Lamb_problem/lamb1.rst new file mode 100644 index 0000000..0e1cda2 --- /dev/null +++ b/docs/source/Lamb_problem/lamb1.rst @@ -0,0 +1,54 @@ +第一类 Lamb 问题 +=================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +第一类 Lamb 问题是指，在半空间模型中，源点和场点均位于地表，求解场点记录到的位移。 + +.. tabs:: + + .. group-tab:: C + + C 程序 :command:`grt` 提供了模块 :doc:`/Module/lamb1` 求解第一类 Lamb 问题。 + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN LAMB1 + :end-before: END LAMB1 + + 使用重定向将结果保存到文件 *lamb1.txt* 中，其内容格式类似于 + + .. literalinclude:: run/head_lamb1 + :language: text + + .. group-tab:: Python + + Python 提供了函数 :func:`solve_lamb1() ` 求解第一类 Lamb 问题。 + + .. literalinclude:: run/lamb1_plot_time.py + :language: python + :start-after: BEGIN LAMB1 + :end-before: END LAMB1 + +最后绘制计算得到的格林函数。 + +:download:`lamb1_plot_time.py ` + +.. image:: run/lamb1_time.png + :align: center + + +频域解和时域解的对比 +------------------------------- + +对比观察可发现频域解在波形突变出有明显的 Gibbs 效应。 + +:download:`lamb1_plot_freq_time.py ` + +.. image:: run/lamb1_compare_freq_time.png + :align: center + + diff --git a/docs/source/Lamb_problem/run/lamb1_plot_freq_time.py b/docs/source/Lamb_problem/run/lamb1_plot_freq_time.py new file mode 100644 index 0000000..8811943 --- /dev/null +++ b/docs/source/Lamb_problem/run/lamb1_plot_freq_time.py @@ -0,0 +1,86 @@ +import pygrt +import numpy as np +import matplotlib.pyplot as plt + +plt.rcParams.update({ + "font.sans-serif": "Times New Roman", + "mathtext.fontset": "cm" +}) + +# 定义半无限空间模型 +Vp = 8 # km/s +Vs = 4.62 # km/s +Rho = 3.3 # g/cm^3 + +# 泊松比 +nu = 0.5 * (1 - 2*(Vs/Vp)**2) / (1 - (Vs/Vp)**2) + +# 模型数组，半无限空间 +modarr = np.array([ + [0., Vp, Vs, Rho, 9e10, 9e10], +]) + +# 剪切模量 +mu = Vs**2 * Rho + +depsrc = 0.0 # 震源深度 km +deprcv = 0.0 # 台站深度 km + +rs = np.array([10]) # 震中距数组，km + +nt = 1000 # 总点数，不要求2的幂次 +dt = 0.005 # 采样时间间隔(s) + +idx = 0 # 震中距索引 +r = rs[idx] + +t = np.arange(0, nt)*dt * Vs/r + + +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) # 整理好的模型对象 +# 计算格林函数频谱 +st = pymod.compute_grn( + distarr=rs, + nt=nt, + dt=dt, +)[0] + +# 卷积阶跃函数 +pygrt.utils.stream_integral(st) + + +# 时域解 +u = pygrt.utils.solve_lamb1(nu, t, 0).reshape(-1, 9) +u = u[:, [0,2,4,6,8]] +u.shape + +# 频域解 +v = np.zeros_like(u) +coef = np.pi * np.pi * mu * r +v[:,0] = st.select(channel='HFR')[0].data * coef +v[:,1] = st.select(channel='VFR')[0].data * coef +v[:,2] = st.select(channel='HFT')[0].data * coef +v[:,3] = st.select(channel='HFZ')[0].data * coef * (-1) +v[:,4] = st.select(channel='VFZ')[0].data * coef * (-1) + +fig, axs = plt.subplots(5, 2, figsize=(10, 10), sharex=True) +labels = [r"$\bar{{G}}^H_{11}$", r"$\bar{{G}}^H_{13}$", r"$\bar{{G}}^H_{22}$", r"$\bar{{G}}^H_{31}$", r"$\bar{{G}}^H_{33}$"] +for i in range(5): + ax = axs[i, 0] + ax.plot(t, u[:,i]) + ax.set_ylim(-2, 2) + ax.set_xlim(0, 2) + ax.grid() + ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12) + + ax = axs[i, 1] + ax.plot(t, v[:,i]) + ax.set_ylim(-2, 2) + ax.set_xlim(0, 2) + ax.grid() + ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12) + +axs[0,0].set_title("From Time-Domain") +axs[0,1].set_title("From Frequency-Domain") + +fig.savefig("lamb1_compare_freq_time.png", dpi=100, bbox_inches='tight') \ No newline at end of file diff --git a/docs/source/Lamb_problem/run/lamb1_plot_time.py b/docs/source/Lamb_problem/run/lamb1_plot_time.py new file mode 100644 index 0000000..71be7b2 --- /dev/null +++ b/docs/source/Lamb_problem/run/lamb1_plot_time.py @@ -0,0 +1,21 @@ +# BEGIN LAMB1 +import pygrt +import numpy as np + +ts = np.arange(0, 2, 1e-4) +u = pygrt.utils.solve_lamb1(0.25, ts, 30) +# END LAMB1 + +import matplotlib.pyplot as plt + +fig, axs = plt.subplots(3, 3, figsize=(10, 5), sharex=True) +for i in range(3): + for j in range(3): + ax = axs[i, j] + ax.plot(ts, u[:, i, j]) + ax.set_xlim(0, 2) + ax.set_ylim(-2, 2) + + ax.text(0.1, 0.9, rf"$\bar{{G}}^H_{{{i+1}{j+1}}}$", transform=ax.transAxes, ha='left', va='top', fontsize=12) + +fig.savefig("lamb1_time.png", dpi=100, bbox_inches='tight') \ No newline at end of file diff --git a/docs/source/Lamb_problem/run/run.sh b/docs/source/Lamb_problem/run/run.sh new file mode 100755 index 0000000..bac86e5 --- /dev/null +++ b/docs/source/Lamb_problem/run/run.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +# ----------------------------------------------------------------------------------- +# BEGIN LAMB1 +grt lamb1 -P0.25 -T0/2/1e-4 -A30 > lamb1.txt +# END LAMB1 +# ----------------------------------------------------------------------------------- + +head -n 10 lamb1.txt > head_lamb1 diff --git a/docs/source/Module/index.rst b/docs/source/Module/index.rst index 4e3c271..ccccd3d 100644 --- a/docs/source/Module/index.rst +++ b/docs/source/Module/index.rst @@ -95,6 +95,7 @@ - :doc:`sac2asc` - :doc:`ker2asc` - :doc:`travt` + - :doc:`lamb1` .. toctree:: @@ -114,4 +115,5 @@ sac2asc ker2asc travt + lamb1 diff --git a/docs/source/Module/lamb1.rst b/docs/source/Module/lamb1.rst new file mode 100644 index 0000000..85406d4 --- /dev/null +++ b/docs/source/Module/lamb1.rst @@ -0,0 +1,50 @@ +.. include:: common_OPTs.rst_ + + +lamb1 +============== + +:简介: 使用广义闭合解求解第一类 Lamb 问题 + +语法 +----------- + +**grt lamb1** +|-P|\ *nu* +|-T|\ *t1/t2/dt* +|-A|\ *azimuth* +[ **-h** ] + + +描述 +-------- + +第一类 Lamb1 问题指在半空间中，当源点和场点均位于地表时，求场点记录到的位移。 +**lamb1** 模块实现的理论基础来源于《地震学中的 Lamb 问题（下）》。 +结果为无量纲位移（距离物理解还需除 :math:`\pi^2 \mu r` ）和阶跃函数卷积后的结果， +输出到标准输出。 + +必选选项 +---------- + +.. _-P: + +**-P**\ *nu* + 半空间的泊松比 *nu*, 要求范围在 (0, 0.5) 。 + +.. _-T: + +**-T**\ *t1/t2/dt* + 无量纲时间序列 :math:`\bar{t}` ，根据起点时刻 *t1*, 终点时刻 *t2* 和间隔 *dt* 定义， + :math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}` ，:math:`r` 为震中距，:math:`\beta` 为 S 波速度。 + +.. _-A: + +**-A**\ *azimuth* + 方位角，单位为度。 + + +示例 +------- + +详见 :doc:`/Lamb_problem/lamb1` 。 \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst index 7e977d4..d650906 100755 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -120,6 +120,7 @@ :hidden: :maxdepth: 1 + Lamb_problem/index API/api changelog copyright diff --git a/docs/source/run_all.sh b/docs/source/run_all.sh index e978766..09d8e5f 100755 --- a/docs/source/run_all.sh +++ b/docs/source/run_all.sh @@ -60,4 +60,13 @@ if [[ $1 == '7' || $1 == '' ]]; then python run.py python plot.py cd - +fi + +if [[ $1 == '8' || $1 == '' ]]; then + cd Lamb_problem/run + chmod +x *.sh + ./run.sh + python lamb1_plot_time.py + python lamb1_plot_freq_time.py + cd - fi \ No newline at end of file