Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 0 additions & 13 deletions docs/source/Tutorial/dynamic/integ_converg.rst

This file was deleted.

2 changes: 1 addition & 1 deletion docs/source/Tutorial/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@
static/static_syn
dynamic/strain_stress
static/static_strain_stress
dynamic/integ_converg
integ_converg/integ_converg



Expand Down
205 changes: 205 additions & 0 deletions docs/source/Tutorial/integ_converg/integ_converg.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
积分收敛性
===================

:Author: Zhu Dengda
:Email: zhudengda@mail.iggcas.ac.cn

-----------------------------------------------------------

通过输出核函数文件,观察源点和场点深度接近时,积分收敛性的变化,以及峰谷平均法的作用。 **具体积分表达式以及分类详见** :ref:`gfunc_rst`。

核函数文件
---------------
该文件记录了波数积分过程中不同波数对应的核函数值。以下用法展示了如何输出该文件。

.. tabs::

.. tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN DGRN
:end-before: END DGRN

输出的核函数文件路径会在 :rst:dir:`GRN_grtstats/milrow_{depsrc}_{deprcv}/` 路径下。

.. tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN DGRN
:end-before: END DGRN

输出的核函数文件路径会在自定义路径下。


C和Python导出的核函数文件是一致的,底层调用的是相同的函数。文件名称格式为 ``K_{iw}_{freq}``,其中 ``{iw}`` 表示频率索引值, ``{freq}`` 表示对应频率(Hz)。文件为自定义的二进制文件, **强烈建议使用Python进行读取及后续处理**。这里还是给出两种读取方法。

.. tabs::

.. tab:: C

:command:`grt.k2a` 程序可将单个核函数文件转为文本格式。

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN grt.k2a
:end-before: END grt.k2a

输出的文件如下,

.. literalinclude:: run/stats_head
:language: text

后续你可以选择习惯的方式读取和处理。

.. tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN read statsfile
:end-before: END read statsfile

其中除了波数 ``k`` 外,每条结果的命名格式均为 ``{srcType}_{q/w/v}{m}``,与 :ref:`gfunc_rst` 部分介绍的积分公式中的核函数 :math:`q_m, w_m, v_m` 保持一致。


可视化
-------------
以下将使用Python进行图件绘制。 **在Python函数中指定震源类型、阶数、积分类型,可自动绘制核函数、被积函数和积分值随波数的变化**,其中积分类型对应 :ref:`gfunc_rst` 部分介绍的4种类型。

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN plot stats
:end-before: END plot stats

.. image:: run/DC_20.png
:align: center


可以通过指定 ``RorI=False`` 参数来指定绘制虚部,传入 ``RorI=2`` 表示实虚部都绘制。

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN plot stats RI
:end-before: END plot stats RI

.. image:: run/DC_20_RI.png
:align: center


从图中可以看到,当波数增加到一定范围以上,积分收敛良好,无振荡现象。



当震源和场点深度接近/相等时,积分收敛速度变慢
----------------------------------------------

假设其它参数不变,我们调整 **震源深度为0.1km**,再计算格林函数,读取核函数文件,绘制图像。

.. tabs::

.. tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN DEPSRC 0.0 DGRN
:end-before: END DEPSRC 0.0 DGRN

.. tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN DEPSRC 0.0 DGRN
:end-before: END DEPSRC 0.0 DGRN


.. image:: run/DC_20_0.0_RI.png
:align: center

从图中可以清晰地看到,相比震源深度2km时,积分收敛速度明显变慢,积分值振荡,这要求增加波数积分上限,但这必然降低计算效率。

你可以尝试 **当震源和场点位于同一深度**,收敛振荡现象也很明显。要注意的是这和是否在地表无关,即使震源和场点都在地下,深度接近时积分收敛也会变慢。


峰谷平均法
-------------------------
具体原理很简洁易懂,从以下图像中你也能了解大概。详见 :ref:`(Zhang et al., 2003) <zhang_2003>` :ref:`(张海明, 2021) <zhang_book_2021>`。

在以上示例中,当震源和场点深度接近时,你会注意到输出的核函数文件增加了,这是因为程序内部设定 **当震源和场点深度差小于1km时,自动使用峰谷平均法(PTAM)**。具体流程为,程序中在进行完离散波数积分后,继续增加k(不同震中距使用不同dk),使用PTAM寻找足够数量的波峰和波谷(内部设定为36个),再对这些波峰波谷取缩减序列 :math:`M_i \leftarrow 0.5\times(M_i + M_{i+1})` ,得到估计的积分收敛值。取缩减序列的C代码如下,

.. code-block:: C

for(int n=n1; n>1; --n){
for(int i=0; i<n-1; ++i){
arr[i] = 0.5*(arr[i] + arr[i+1]);
}
}


在 ``K_{iw}_{freq}`` 文件同级目录下,程序把 **PTAM过程中的核函数以及积分峰谷位置分为两个文件** 保存在 ``PTAM_{ir}_{dist}/`` 目录下( ``{ir}`` 为震中距索引, ``{dist}`` 为震中距),其中 ``PTAM_{ir}_{dist}/K_{iw}_{freq}`` 为核函数文件(格式不变), ``PTAM_{ir}_{dist}/PTAM_{iw}_{freq}`` 为峰谷文件,其中记录积分值的峰谷。

.. note::

:command:`grt.k2a` 程序也支持将 ``PTAM_{ir}_{dist}/PTAM_{iw}_{freq}`` 文件转为文本格式,

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN grt.k2a ptam
:end-before: END grt.k2a ptam

输出的文件如下,

.. literalinclude:: run/ptam_stats_head
:language: text

记录了不同震源、不同积分类型的峰谷位置。


故为了绘制完整的波数积分+PTAM过程,涉及三个文件。不过Python函数已经做了优化,由于程序输出的目录结构固定,文件命名格式固定,只需传入 ``PTAM_{ir}_{dist}/PTAM_{iw}_{freq}`` 文件,Python函数会自动读取对应的三个文件。

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN plot ptam
:end-before: END plot ptam

.. image:: run/DC_20_0.0_ptam_RI.png
:align: center


红十字为选取的波峰波谷,再取缩减序列即可得到估计的积分收敛值。



静态解的积分收敛性
-------------------------
以上部分是以动态解为例,静态解从积分类型、收敛特征、文件格式、绘图完全类似,只是不再有频率索引值。

假设震源深度0.1km,场点位于地表,场点仅定义一个点(2,2)做示例,这里直接给出脚本。

.. tabs::

.. tab:: C

.. literalinclude:: run/run.sh
:language: bash
:start-after: BEGIN SGRN
:end-before: END SGRN

.. tab:: Python

.. literalinclude:: run/run.py
:language: python
:start-after: BEGIN SGRN
:end-before: END SGRN


+ 只使用离散波数积分

.. image:: run/DC_20_0.1_static.png
:align: center

+ 使用峰谷平均法

.. image:: run/DC_20_0.1_ptam_static.png
:align: center

Binary file added docs/source/Tutorial/integ_converg/run/DC_20.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
9 changes: 9 additions & 0 deletions docs/source/Tutorial/integ_converg/run/milrow
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
0.2 3.4 1.7 2.3 9e10 9e10
0.6 3.7 1.9 2.4 9e10 9e10
0.5 4.2 2.1 2.4 9e10 9e10
0.5 4.6 2.3 2.5 9e10 9e10
0.7 4.9 2.8 2.6 9e10 9e10
0.5 5.1 2.9 2.7 9e10 9e10
6.0 5.9 3.3 2.7 9e10 9e10
1. 6.9 4.0 2.8 9e10 9e10
0.0 8.2 4.7 3.2 9e10 9e10
11 changes: 11 additions & 0 deletions docs/source/Tutorial/integ_converg/run/ptam_stats_head
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
# sum_EXP_00_k sum_EXP_00 sum_VF_00_k sum_VF_00 sum_DC_00_k sum_DC_00 sum_EXP_02_k sum_EXP_02 sum_VF_02_k sum_VF_02 sum_DC_02_k sum_DC_02 sum_HF_10_k sum_HF_10 sum_DC_10_k sum_DC_10 sum_HF_11_k sum_HF_11 sum_DC_11_k sum_DC_11 sum_HF_12_k sum_HF_12 sum_DC_12_k sum_DC_12 sum_HF_13_k sum_HF_13 sum_DC_13_k sum_DC_13 sum_DC_20_k sum_DC_20 sum_DC_21_k sum_DC_21 sum_DC_22_k sum_DC_22 sum_DC_23_k sum_DC_23
2.67622465e+01 6.66555330e+02-6.52355489e+03J 2.67622568e+01 -6.71825014e+01+1.37589835e+03J 2.67622465e+01 -2.66622132e+03+2.60942196e+04J 2.66053630e+01 1.35857785e+04+7.55803841e+02J 2.66053633e+01 -1.97736289e+03-2.98716114e+01J 2.66053630e+01 -5.43431139e+04-3.02321536e+03J 2.66053604e+01 -6.25363762e+02-3.80092044e+02J 2.65219514e+01 2.42729886e-11-5.31849265e-11J 2.67622542e+01 -1.03466177e+01+9.34738740e+00J 2.65154089e+01 6.92762231e-13-1.64735342e-13J 2.67622568e+01 6.71825014e+01-1.37589835e+03J 2.66205772e+01 1.83164486e-10+3.18094236e-12J 2.66053605e+01 1.38503321e+03+8.64167945e+02J 2.65098060e+01 0.00000000e+00+0.00000000e+00J 2.67622465e+01 1.33311066e+03-1.30471098e+04J 2.66045993e+01 -4.06156039e+02-2.45114486e+02J 2.66046025e+01 2.71917459e+04+1.23617601e+03J 2.67622465e+01 -1.80515555e+04+2.77535386e+04J
2.70764072e+01 2.52519538e+03-6.56643813e+03J 2.70764171e+01 -1.32379593e+02+1.37827970e+03J 2.70764072e+01 -1.01007815e+04+2.62657525e+04J 2.69195212e+01 1.26971000e+04+7.88692689e+02J 2.69195216e+01 -1.79715336e+03-3.51777813e+01J 2.69195212e+01 -5.07883999e+04-3.15477075e+03J 2.69195189e+01 -4.86255714e+02-3.83331419e+02J 2.66699769e+01 2.62038907e-11-5.32231922e-11J 2.70764150e+01 -1.00454716e+01+9.34041115e+00J 2.66063826e+01 6.77250543e-13-1.64858007e-13J 2.70764171e+01 1.32379593e+02-1.37827970e+03J 2.66898801e+01 1.86562483e-10+3.02839093e-12J 2.69195190e+01 1.60593887e+03+8.59028414e+02J 2.65883458e+01 0.00000000e+00+0.00000000e+00J 2.70764072e+01 5.05039077e+03-1.31328763e+04J 2.69187667e+01 -4.22517373e+02-2.44734408e+02J 2.69187695e+01 2.54141606e+04+1.30196624e+03J 2.70764073e+01 -1.21489937e+04+2.76171490e+04J
2.73905679e+01 6.73937049e+02-6.52446155e+03J 2.73905774e+01 -6.93843871e+01+1.37603618e+03J 2.73905679e+01 -2.69574820e+03+2.60978462e+04J 2.72336795e+01 1.35651038e+04+7.57392855e+02J 2.72336799e+01 -1.97297785e+03-3.01233966e+01J 2.72336795e+01 -5.42604152e+04-3.02957142e+03J 2.72336774e+01 -6.23162563e+02-3.80200555e+02J 2.67559430e+01 3.41902954e-11-5.32000344e-11J 2.73905757e+01 -1.03383271e+01+9.34715594e+00J 2.66773218e+01 6.90513962e-13-1.65285048e-13J 2.73905774e+01 6.93843871e+01-1.37603618e+03J 2.69233162e+01 2.01682508e-10+3.07569034e-12J 2.72336775e+01 1.38857526e+03+8.64017751e+02J 2.66668856e+01 0.00000000e+00+0.00000000e+00J 2.73905679e+01 1.34787410e+03-1.30489231e+04J 2.72329339e+01 -4.06424282e+02-2.45106147e+02J 2.72329362e+01 2.71503815e+04+1.23935525e+03J 2.73905680e+01 -1.80264891e+04+2.77511940e+04J
2.77047285e+01 2.51888998e+03-6.56562927e+03J 2.77047378e+01 -1.30349223e+02+1.37815606e+03J 2.77047285e+01 -1.00755599e+04+2.62625171e+04J 2.75478378e+01 1.27159979e+04+7.87267798e+02J 2.75478383e+01 -1.80122353e+03-3.49512546e+01J 2.75478378e+01 -5.08639916e+04-3.14907119e+03J 2.75478359e+01 -4.88323031e+02-3.83233078e+02J 2.69858115e+01 1.64622008e-11-5.31797570e-11J 2.77047363e+01 -1.00533787e+01+9.34063305e+00J 2.69057304e+01 7.38538662e-13-1.65293921e-13J 2.77047378e+01 1.30349223e+02-1.37815606e+03J 2.70777608e+01 1.97480539e-10+3.19880445e-12J 2.75478361e+01 1.60258550e+03+8.59167341e+02J 2.67454255e+01 0.00000000e+00+0.00000000e+00J 2.77047285e+01 5.03777996e+03-1.31312585e+04J 2.75471010e+01 -4.22260093e+02-2.44742526e+02J 2.75471029e+01 2.54519700e+04+1.29911541e+03J 2.77047287e+01 -1.21712548e+04+2.76192913e+04J
2.80188892e+01 6.79285408e+02-6.52518567e+03J 2.80188982e+01 -7.12624018e+01+1.37614758e+03J 2.80188892e+01 -2.71714163e+03+2.61007427e+04J 2.78619961e+01 1.35477811e+04+7.58676203e+02J 2.78619967e+01 -1.96918808e+03-3.03281499e+01J 2.78619961e+01 -5.41911242e+04-3.03470481e+03J 2.78619945e+01 -6.21215273e+02-3.80290105e+02J 2.71012124e+01 3.30180096e-11-5.30607034e-11J 2.80188970e+01 -1.03307817e+01+9.34694370e+00J 2.69926585e+01 6.94850926e-13-1.64914019e-13J 2.80188982e+01 7.12624018e+01-1.37614758e+03J 2.71405501e+01 1.98200699e-10+3.14969160e-12J 2.78619946e+01 1.39175621e+03+8.63888863e+02J 2.68239653e+01 0.00000000e+00+0.00000000e+00J 2.80188892e+01 1.35857082e+03-1.30503713e+04J 2.78612679e+01 -4.06671081e+02-2.45098278e+02J 2.78612694e+01 2.71157239e+04+1.24192286e+03J 2.80188894e+01 -1.80067770e+04+2.77492322e+04J
2.83330498e+01 2.51439617e+03-6.56497902e+03J 2.83330585e+01 -1.28607138e+02+1.37805527e+03J 2.83330498e+01 -1.00575847e+04+2.62599161e+04J 2.81761545e+01 1.27319178e+04+7.86107252e+02J 2.81761552e+01 -1.80476247e+03-3.47653851e+01J 2.81761545e+01 -5.09276711e+04-3.14442901e+03J 2.81761531e+01 -4.90162208e+02-3.83151174e+02J 2.71755016e+01 2.84552368e-11-5.29861787e-11J 2.83330576e+01 -1.00605833e+01+9.34083574e+00J 2.70775518e+01 7.18655762e-13-1.64877890e-13J 2.83330585e+01 1.28607138e+02-1.37805527e+03J 2.72057560e+01 1.95856080e-10+3.27609856e-12J 2.81761532e+01 1.59956258e+03+8.59287249e+02J 2.69025051e+01 0.00000000e+00+0.00000000e+00J 2.83330498e+01 5.02879234e+03-1.31299580e+04J 2.81754346e+01 -4.22023282e+02-2.44750129e+02J 2.81754358e+01 2.54838206e+04+1.29679353e+03J 2.83330500e+01 -1.21886443e+04+2.76210913e+04J
2.86472104e+01 6.83013242e+02-6.52577117e+03J 2.86472189e+01 -7.28827184e+01+1.37623912e+03J 2.86472104e+01 -2.73205297e+03+2.61030847e+04J 2.84903130e+01 1.35331158e+04+7.59729619e+02J 2.84903137e+01 -1.96587449e+03-3.04975494e+01J 2.84903130e+01 -5.41324632e+04-3.03891848e+03J 2.84903117e+01 -6.19473898e+02-3.80365318e+02J 2.73193464e+01 3.17829035e-11-5.29670298e-11J 2.86472181e+01 -1.03238981e+01+9.34675032e+00J 2.71675681e+01 7.25092543e-13-1.64995602e-13J 2.86472189e+01 7.28827184e+01-1.37623912e+03J 2.73247322e+01 2.02381773e-10+3.32872257e-12J 2.84903119e+01 1.39463388e+03+8.63777020e+02J 2.69810449e+01 0.00000000e+00+0.00000000e+00J 2.86472104e+01 1.36602648e+03-1.30515423e+04J 2.84896011e+01 -4.06898398e+02-2.45090950e+02J 2.84896021e+01 2.70863837e+04+1.24403039e+03J 2.86472106e+01 -1.79915090e+04+2.77475776e+04J
2.89613710e+01 2.51135736e+03-6.56445052e+03J 2.89613792e+01 -1.27096316e+02+1.37797185e+03J 2.89613710e+01 -1.00454294e+04+2.62578021e+04J 2.88044714e+01 1.27454566e+04+7.85147773e+02J 2.88044722e+01 -1.80787280e+03-3.46104256e+01J 2.88044714e+01 -5.09818265e+04-3.14059109e+03J 2.88044703e+01 -4.91814753e+02-3.83081845e+02J 2.75044851e+01 4.05272699e-11-5.30701150e-11J 2.89613787e+01 -1.00671646e+01+9.34102012e+00J 2.73243991e+01 7.48357530e-13-1.64621873e-13J 2.89613792e+01 1.27096316e+02-1.37797185e+03J 2.75113883e+01 1.92670301e-10+3.35712528e-12J 2.88044705e+01 1.59681879e+03+8.59391825e+02J 2.70595847e+01 0.00000000e+00+0.00000000e+00J 2.89613710e+01 5.02271471e+03-1.31289010e+04J 2.88037675e+01 -4.21804977e+02-2.44757180e+02J 2.88037682e+01 2.55109071e+04+1.29487395e+03J 2.89613713e+01 -1.22019695e+04+2.76226148e+04J
2.92755316e+01 6.85430224e+02-6.52624925e+03J 2.92755395e+01 -7.42947329e+01+1.37631540e+03J 2.92755316e+01 -2.74172089e+03+2.61049970e+04J 2.91186299e+01 1.35205919e+04+7.60606336e+02J 2.91186308e+01 -1.96294815e+03-3.06397871e+01J 2.91186299e+01 -5.40823677e+04-3.04242534e+03J 2.91186290e+01 -6.17902335e+02-3.80429446e+02J 2.75540770e+01 4.00135248e-11-5.30596745e-11J 2.92755392e+01 -1.03176014e+01+9.34657458e+00J 2.74979689e+01 7.12410051e-13-1.64503655e-13J 2.92755395e+01 7.42947329e+01-1.37631540e+03J 2.75966853e+01 1.96164020e-10+3.36759361e-12J 2.91186292e+01 1.39725398e+03+8.63679015e+02J 2.71381245e+01 0.00000000e+00+0.00000000e+00J 2.92755316e+01 1.37086045e+03-1.30524985e+04J 2.91179337e+01 -4.07108158e+02-2.45084173e+02J 2.91179342e+01 2.70613279e+04+1.24578436e+03J 2.92755318e+01 -1.79799669e+04+2.77461729e+04J
...
Loading