Skip to content

Commit 6ab4fea

Browse files
Dengda98Dengda98
authored
DOC: add lamb1 (the first-kind Lamb problem) (#119)
1 parent cab9d59 commit 6ab4fea

File tree

9 files changed

+260
-0
lines changed

9 files changed

+260
-0
lines changed

‎docs/source/Lamb_problem/index.rst‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
Lamb 问题
2+
================
3+
4+
:Author: Zhu Dengda
5+
:Email: [email protected]
6+
7+
-----------------------------------------------------------
8+
9+
北京大学张海明教授已在其编著的 **《地震学中的 Lamb 问题》上下两册** 中对三类 Lamb 问题进行了详细论述，
10+
并分别在上下两册中给出了频域解和时域解，过程之详细令人叹为观止。这里不再对 Lamb 问题及其解法做过多介绍，
11+
详见张海明老师的书（强烈推荐！！）
12+
13+
**PyGRT** 程序包本身是在频域求解，这里我对照下册书也编程实现了 Lamb 问题的时域解（目前包括第一类 Lamb 问题，后续会再扩展），
14+
方便进行对比验证。以下简单示范用法以及绘图。
15+
16+
.. note::
17+
18+
程序输入的时间为无量纲时间，输出的位移为和阶跃函数卷积后的无量纲位移。
19+
20+
+ 无量纲时间 :math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}`
21+
+ 无量纲位移 :math:`\bar{\mathbf{G}}^H = \pi^2 \mu r \mathbf{G}^H` ，上标 :math:`H` 表示已和阶跃函数卷积。
22+
23+
24+
.. toctree::
25+
:maxdepth: 1
26+
27+
lamb1
28+

‎docs/source/Lamb_problem/lamb1.rst‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
第一类 Lamb 问题
2+
===================
3+
4+
:Author: Zhu Dengda
5+
:Email: [email protected]
6+
7+
-----------------------------------------------------------
8+
9+
第一类 Lamb 问题是指，在半空间模型中，源点和场点均位于地表，求解场点记录到的位移。
10+
11+
.. tabs::
12+
13+
.. group-tab:: C
14+
15+
C 程序 :command:`grt` 提供了模块 :doc:`/Module/lamb1` 求解第一类 Lamb 问题。
16+
17+
.. literalinclude:: run/run.sh
18+
:language: bash
19+
:start-after: BEGIN LAMB1
20+
:end-before: END LAMB1
21+
22+
使用重定向将结果保存到文件 *lamb1.txt* 中，其内容格式类似于
23+
24+
.. literalinclude:: run/head_lamb1
25+
:language: text
26+
27+
.. group-tab:: Python
28+
29+
Python 提供了函数 :func:`solve_lamb1() <pygrt.utils.solve_lamb1>` 求解第一类 Lamb 问题。
30+
31+
.. literalinclude:: run/lamb1_plot_time.py
32+
:language: python
33+
:start-after: BEGIN LAMB1
34+
:end-before: END LAMB1
35+
36+
最后绘制计算得到的格林函数。
37+
38+
:download:`lamb1_plot_time.py <run/lamb1_plot_time.py>`
39+
40+
.. image:: run/lamb1_time.png
41+
:align: center
42+
43+
44+
频域解和时域解的对比
45+
-------------------------------
46+
47+
对比观察可发现频域解在波形突变出有明显的 Gibbs 效应。
48+
49+
:download:`lamb1_plot_freq_time.py <run/lamb1_plot_freq_time.py>`
50+
51+
.. image:: run/lamb1_compare_freq_time.png
52+
:align: center
53+
54+

‎docs/source/Lamb_problem/run/lamb1_plot_freq_time.py‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
import pygrt
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
5+
plt.rcParams.update({
6+
"font.sans-serif": "Times New Roman",
7+
"mathtext.fontset": "cm"
8+
})
9+
10+
# 定义半无限空间模型
11+
Vp = 8 # km/s
12+
Vs = 4.62 # km/s
13+
Rho = 3.3 # g/cm^3
14+
15+
# 泊松比
16+
nu = 0.5 * (1 - 2*(Vs/Vp)**2) / (1 - (Vs/Vp)**2)
17+
18+
# 模型数组，半无限空间
19+
modarr = np.array([
20+
[0., Vp, Vs, Rho, 9e10, 9e10],
21+
])
22+
23+
# 剪切模量
24+
mu = Vs**2 * Rho
25+
26+
depsrc = 0.0 # 震源深度 km
27+
deprcv = 0.0 # 台站深度 km
28+
29+
rs = np.array([10]) # 震中距数组，km
30+
31+
nt = 1000 # 总点数，不要求2的幂次
32+
dt = 0.005 # 采样时间间隔(s)
33+
34+
idx = 0 # 震中距索引
35+
r = rs[idx]
36+
37+
t = np.arange(0, nt)*dt * Vs/r
38+
39+
40+
pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) # 整理好的模型对象
41+
# 计算格林函数频谱
42+
st = pymod.compute_grn(
43+
distarr=rs,
44+
nt=nt,
45+
dt=dt,
46+
)[0]
47+
48+
# 卷积阶跃函数
49+
pygrt.utils.stream_integral(st)
50+
51+
52+
# 时域解
53+
u = pygrt.utils.solve_lamb1(nu, t, 0).reshape(-1, 9)
54+
u = u[:, [0,2,4,6,8]]
55+
u.shape
56+
57+
# 频域解
58+
v = np.zeros_like(u)
59+
coef = np.pi * np.pi * mu * r
60+
v[:,0] = st.select(channel='HFR')[0].data * coef
61+
v[:,1] = st.select(channel='VFR')[0].data * coef
62+
v[:,2] = st.select(channel='HFT')[0].data * coef
63+
v[:,3] = st.select(channel='HFZ')[0].data * coef * (-1)
64+
v[:,4] = st.select(channel='VFZ')[0].data * coef * (-1)
65+
66+
fig, axs = plt.subplots(5, 2, figsize=(10, 10), sharex=True)
67+
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}$"]
68+
for i in range(5):
69+
ax = axs[i, 0]
70+
ax.plot(t, u[:,i])
71+
ax.set_ylim(-2, 2)
72+
ax.set_xlim(0, 2)
73+
ax.grid()
74+
ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12)
75+
76+
ax = axs[i, 1]
77+
ax.plot(t, v[:,i])
78+
ax.set_ylim(-2, 2)
79+
ax.set_xlim(0, 2)
80+
ax.grid()
81+
ax.text(0.05, 0.92, labels[i], transform=ax.transAxes, ha='left', va='top', fontsize=12)
82+
83+
axs[0,0].set_title("From Time-Domain")
84+
axs[0,1].set_title("From Frequency-Domain")
85+
86+
fig.savefig("lamb1_compare_freq_time.png", dpi=100, bbox_inches='tight')

‎docs/source/Lamb_problem/run/lamb1_plot_time.py‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
# BEGIN LAMB1
2+
import pygrt
3+
import numpy as np
4+
5+
ts = np.arange(0, 2, 1e-4)
6+
u = pygrt.utils.solve_lamb1(0.25, ts, 30)
7+
# END LAMB1
8+
9+
import matplotlib.pyplot as plt
10+
11+
fig, axs = plt.subplots(3, 3, figsize=(10, 5), sharex=True)
12+
for i in range(3):
13+
for j in range(3):
14+
ax = axs[i, j]
15+
ax.plot(ts, u[:, i, j])
16+
ax.set_xlim(0, 2)
17+
ax.set_ylim(-2, 2)
18+
19+
ax.text(0.1, 0.9, rf"$\bar{{G}}^H_{{{i+1}{j+1}}}$", transform=ax.transAxes, ha='left', va='top', fontsize=12)
20+
21+
fig.savefig("lamb1_time.png", dpi=100, bbox_inches='tight')

‎docs/source/Lamb_problem/run/run.sh‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
#!/bin/bash
2+
3+
# -----------------------------------------------------------------------------------
4+
# BEGIN LAMB1
5+
grt lamb1 -P0.25 -T0/2/1e-4 -A30 > lamb1.txt
6+
# END LAMB1
7+
# -----------------------------------------------------------------------------------
8+
9+
head -n 10 lamb1.txt > head_lamb1

‎docs/source/Module/index.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,7 @@
9595
- :doc:`sac2asc`
9696
- :doc:`ker2asc`
9797
- :doc:`travt`
98+
- :doc:`lamb1`
9899

99100

100101
.. toctree::
@@ -114,4 +115,5 @@
114115
sac2asc
115116
ker2asc
116117
travt
118+
lamb1
117119

‎docs/source/Module/lamb1.rst‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
.. include:: common_OPTs.rst_
2+
3+
4+
lamb1
5+
==============
6+
7+
:简介: 使用广义闭合解求解第一类 Lamb 问题
8+
9+
语法
10+
-----------
11+
12+
**grt lamb1**
13+
|-P|\ *nu*
14+
|-T|\ *t1/t2/dt*
15+
|-A|\ *azimuth*
16+
[ **-h** ]
17+
18+
19+
描述
20+
--------
21+
22+
第一类 Lamb1 问题指在半空间中，当源点和场点均位于地表时，求场点记录到的位移。
23+
**lamb1** 模块实现的理论基础来源于《地震学中的 Lamb 问题（下）》。
24+
结果为无量纲位移（距离物理解还需除 :math:`\pi^2 \mu r` ）和阶跃函数卷积后的结果，
25+
输出到标准输出。
26+
27+
必选选项
28+
----------
29+
30+
.. _-P:
31+
32+
**-P**\ *nu*
33+
半空间的泊松比 *nu*, 要求范围在 (0, 0.5) 。
34+
35+
.. _-T:
36+
37+
**-T**\ *t1/t2/dt*
38+
无量纲时间序列 :math:`\bar{t}` ，根据起点时刻 *t1*, 终点时刻 *t2* 和间隔 *dt* 定义，
39+
:math:`\bar{t} = \dfrac{t}{T_S}` ，其中 :math:`T_S = \dfrac{r}{\beta}` ，:math:`r` 为震中距，:math:`\beta` 为 S 波速度。
40+
41+
.. _-A:
42+
43+
**-A**\ *azimuth*
44+
方位角，单位为度。
45+
46+
47+
示例
48+
-------
49+
50+
详见 :doc:`/Lamb_problem/lamb1` 。

‎docs/source/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -120,6 +120,7 @@
120120
:hidden:
121121
:maxdepth: 1
122122

123+
Lamb_problem/index
123124
API/api
124125
changelog
125126
copyright

‎docs/source/run_all.sh‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,4 +60,13 @@ if [[ $1 == '7' || $1 == '' ]]; then
6060
python run.py
6161
python plot.py
6262
cd -
63+
fi
64+
65+
if [[ $1 == '8' || $1 == '' ]]; then
66+
cd Lamb_problem/run
67+
chmod +x *.sh
68+
./run.sh
69+
python lamb1_plot_time.py
70+
python lamb1_plot_freq_time.py
71+
cd -
6372
fi

0 commit comments

Comments
 (0)