FEAT: set the max frequency as the attenuation reference frequency

docs/source/Gallery/ex17/ex17.rst

@@ -0,0 +1,28 @@
+:author: 朱邓达
+:date: 2025-01-23
+
+(17) 衰减对波形的影响
+-------------------------------------
+
+下载示例： :download:`ex17.tar.gz`
+
+以半空间模型为例，比较衰减参数对波形的影响。
+半空间 :math:`V_P=8.0` km/s，:math:`V_S=4.62` km/s，:math:`\rho=3.3` g/cm^3，
+:math:`Q_P=80` ， :math:`Q_S=20` 。
+
+.. note::
+
+    + 在强衰减模型中，核函数在极点附近曲线很平滑而非尖锐变化，
+      此时对于源点和场点深度接近的极端情况，积分上限截止之后的部分对于积分的贡献逐渐明显，
+      因此收敛方法（如DCM，PTAM）需要更高的积分上限才能给出稳定结果。
+
+    + 自 v0.15.0 版本起，对于衰减参数的处理，程序不再使用常见默认的 1 Hz 为参考频率，
+      而使用当前计算的最高频作为参考频率，这样计算出来的波形才符合因果关系
+      （之前的结果会发现部分波形出现在到时之前，也就是这个原因）。
+      对应公式详见 :c:func:`grt_attenuation_law` 函数。
+
+
+.. figure:: compare.svg
+    :align: center
+
+    蓝虚线为 P 波到时，红虚线为 S 波到时

docs/source/Gallery/ex17/plot.py

@@ -0,0 +1,46 @@
+import numpy as np
+from obspy import *
+import matplotlib.pyplot as plt
+
+import pygrt
+
+st = read("GRN/*/HF?.sac")
+st_Q = read("GRN_Q/*/HF?.sac")
+
+nt = st[0].stats.npts
+dt = st[0].stats.delta
+
+for _ in range(1):
+    st.integrate()
+    st_Q.integrate()
+
+
+
+ts0 = 5 / 4.62
+t = np.arange(nt)*dt / ts0
+
+fig, axs = plt.subplots(2, 3, figsize=(12, 6), sharex=True, gridspec_kw=dict(hspace=0))
+
+def _plot(ax, st, channel):
+    tr = st.select(channel=channel)[0]
+    tp = tr.stats.sac['t0'] / ts0
+    ts = 1.0 
+    ax.plot(t, tr.data, c='k', lw=1.0)
+    ylim = ax.get_ylim()
+    ax.vlines(tp, *ylim, colors='b', ls='--', lw=0.5)
+    ax.vlines(ts, *ylim, colors='r', ls='--', lw=0.5)
+    ax.text(0.05, 0.9, channel, transform=ax.transAxes, bbox=dict(fc='w'), ha='left')
+    ax.set_xmargin(0)
+    ax.ticklabel_format(axis='y', style='sci', scilimits=[0,0])
+
+_plot(axs[0,0], st,   'HFZ')
+_plot(axs[1,0], st_Q, 'HFZ')
+_plot(axs[0,1], st,   'HFR')
+_plot(axs[1,1], st_Q, 'HFR')
+_plot(axs[0,2], st,   'HFT')
+_plot(axs[1,2], st_Q, 'HFT')
+
+axs[0,0].set_ylabel("No attenuation")
+axs[1,0].set_ylabel(r"$Q_P=80, Q_S=20$")
+
+fig.savefig("compare.svg", bbox_inches='tight')

docs/source/Gallery/ex17/run.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+set -euo pipefail
+
+rm -rf GRN* *.svg *.tar.gz
+
+cat > halfspace <<EOF
+0.0   8.0   4.62    3.3
+EOF
+
+cat > halfspace_Q <<EOF
+0.0   8.0   4.62    3.3    80    20
+EOF
+
+# 强衰减介质需要更高的积分上限
+K="-K+k60+s2"
+grt greenfn -Mhalfspace   -N2000/0.001+a -D0/0  -R5  -OGRN     $K -Cp
+grt greenfn -Mhalfspace_Q -N2000/0.001+a -D0/0  -R5  -OGRN_Q   $K -Cp
+python plot.py
+
+cp compare.svg cover.svg
+
+ex=$(basename $(pwd))
+cd .. && tar -czvf ${ex}.tar.gz ${ex} && mv ${ex}.tar.gz ${ex} && cd -
+
+rm -rf halfspace* GRN* 
+
+
+

docs/source/Gallery/gallery.rst

@@ -14,7 +14,7 @@
 
     .. jinja::
 
-        {% for i in range(1, 17) %}
+        {% for i in range(1, 18) %}
         {% set i = '%02d' % i %}
         .. grid-item-card:: :doc:`ex{{i}}/ex{{i}}`
             :padding: 1

pygrt/C_extension/include/grt/common/attenuation.h

@@ -19,13 +19,9 @@
  * 其中虚数部分的正负号和书中不同，是因为书中使用的傅里叶变换的e指数符号和我们通常使用的相反。
  * 
  * @param[in]    Qinv     1/Q
- * @param[in]    omega    复数频率\f$ \tilde{\omega} =\omega - i\zeta \f$ 
+ * @param[in]    omgref   参考圆频率；原公式默认使用 1Hz 为参考频率
+ * @param[in]    omega    复数圆频率\f$ \tilde{\omega} =\omega - i\zeta \f$ 
  * 
  * @return atncoef 系数因子，作用在 \f$ k=\omega / c(\omega)\f$的计算
  */
-cplx_t grt_attenuation_law(real_t Qinv, cplx_t omega);
-
-/**
- * attenuation_law函数在python中被调用的版本，长度2的数组分别表示复数的实部和虚部
- */
-void grt_py_attenuation_law(real_t Qinv, real_t omg[2], real_t atte[2]);
+cplx_t grt_attenuation_law(real_t Qinv, cplx_t omgref, cplx_t omega);

pygrt/C_extension/include/grt/common/model.h

@@ -24,6 +24,7 @@ typedef struct {
     bool ircvup; ///< 接收点位于浅层, ircv < isrc
     bool io_depth; ///< 读取的模型首列为每层顶界面深度
 
+    cplx_t omgref;   ///< 参考圆频率, 用于计算衰减
     cplx_t omega;   ///< 圆频率
     real_t k;   ///< 波数
     cplx_t c_phase;   ///< 当前相速度

pygrt/C_extension/src/common/attenuation.c

@@ -12,15 +12,6 @@
 
 
 
-cplx_t grt_attenuation_law(real_t Qinv, cplx_t omega){
-    return 1.0 + Qinv/PI * log(omega/PI2) + 0.5*Qinv*I;
-    // return 1.0;
-}
-
-void grt_py_attenuation_law(real_t Qinv, real_t omg[2], real_t atte[2]){
-    // 用于在python中调用attenuation_law
-    cplx_t omega = omg[0] + I*omg[1];
-    cplx_t atte0 = grt_attenuation_law(Qinv, omega);
-    atte[0] = creal(atte0);
-    atte[1] = cimag(atte0);
+cplx_t grt_attenuation_law(real_t Qinv, cplx_t omgref, cplx_t omega){
+    return 1.0 + Qinv/PI * log(omega/omgref) + 0.5*Qinv*I;
 }

pygrt/C_extension/src/common/model.c

@@ -121,6 +121,7 @@ void grt_free_mod1d(GRT_MODEL1D *mod1d){
 GRT_MODEL1D * grt_init_mod1d(size_t n){
     GRT_MODEL1D *mod1d = (GRT_MODEL1D *)calloc(1, sizeof(GRT_MODEL1D));
     mod1d->n = n;
+    mod1d->omgref = PI2*1.0;
 
     #define X(P, T)  mod1d->P = (T*)calloc(n, sizeof(T));
         GRT_FOR_EACH_MODEL_QUANTITY_ARRAY
@@ -158,8 +159,8 @@ void grt_attenuate_mod1d(GRT_MODEL1D *mod1d, cplx_t omega){
 
         // 圆频率实部为负数表明不考虑模型的 Q 值属性
         // 在读入模型后需要需要运行一次本函数以填充弹性模量，见 grt_read_mod1d_from_file 函数
-        atna = (creal(omega) >= 0.0 && mod1d->Qainv[i] > 0.0)? grt_attenuation_law(mod1d->Qainv[i], omega) : 1.0;
-        atnb = (creal(omega) >= 0.0 && mod1d->Qbinv[i] > 0.0)? grt_attenuation_law(mod1d->Qbinv[i], omega) : 1.0;
+        atna = (creal(omega) >= 0.0 && mod1d->Qainv[i] > 0.0)? grt_attenuation_law(mod1d->Qainv[i], mod1d->omgref, omega) : 1.0;
+        atnb = (creal(omega) >= 0.0 && mod1d->Qbinv[i] > 0.0)? grt_attenuation_law(mod1d->Qbinv[i], mod1d->omgref, omega) : 1.0;
 
         mod1d->atna[i] = atna;
         mod1d->atnb[i] = atnb;

pygrt/C_extension/src/dynamic/grn.c

@@ -90,6 +90,8 @@ void grt_integ_grn_spec(
     // 实际计算的频点数
     size_t nf_valid = nf2 - nf1 + 1;
 
+    mod1d->omgref = PI2*freqs[nf2];
+
     // 频率omega循环
     // schedule语句可以动态调度任务，最大程度地使用计算资源
     #pragma omp parallel for schedule(guided) default(shared) 

pygrt/C_extension/src/dynamic/grt_kernel.c

@@ -379,6 +379,8 @@ int kernel_main(int argc, char **argv)
     const real_t Rho = mod1d->Rho[mod1d->isrc]; // 震源区密度
     const real_t fac = 1.0/(4.0*PI*Rho);
 
+    mod1d->omgref = PI2*Ctrl->F.freqs[Ctrl->F.nf-1];
+
     // 频率循环
     #pragma omp parallel for schedule(guided) default(shared) 
     for(size_t iw = 0; iw < Ctrl->F.nf; ++iw)

pygrt/C_extension/src/dynamic/grt_stress.c

@@ -125,8 +125,8 @@ int stress_main(int argc, char **argv){
         freq = (i==0) ? 0.01f : df*i; // 计算衰减因子不能为0频
         w = PI2 * freq;
         fftwf_complex atta, attb;
-        atta = grt_attenuation_law(Qainv, w);
-        attb = grt_attenuation_law(Qbinv, w);
+        atta = grt_attenuation_law(Qainv, PI2*(nf-1)*df, w);
+        attb = grt_attenuation_law(Qbinv, PI2*(nf-1)*df, w);
         // 乘上1e10，转为dyne/(cm^2)
         mus[i] = vb*vb*attb*attb*rho*1e10;
         lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i];

pygrt/c_interfaces.py

@@ -141,10 +141,10 @@ def set_num_threads(n):
 # -------------------------------------------------------------------
 #                      C函数定义的衰减函数
 # -------------------------------------------------------------------
-C_grt_py_attenuation_law = libgrt.grt_py_attenuation_law
+C_grt_attenuation_law = libgrt.grt_attenuation_law
 """品质因子Q 对 波速的影响"""
-C_grt_py_attenuation_law.restype = None
-C_grt_py_attenuation_law.argtypes = [REAL, DPOINTER, DPOINTER]  # double, double[2], double[2]
+C_grt_attenuation_law.restype = CPLX
+C_grt_attenuation_law.argtypes = [REAL, CPLX, CPLX]  # double, cplx, cplx
 
 
 # -------------------------------------------------------------------

pygrt/c_structures.py

@@ -29,6 +29,7 @@
 
     "REAL",
     "PREAL",
+    "CPLX",
     "PCPLX",
 
     "c_GRT_MODEL1D",
@@ -51,9 +52,14 @@
 NPCT_CMPLX_TYPE = 'c16'
 
 
+class CPLX(Structure):
+    _fields_ = [
+        ('real', c_double),
+        ('imag', c_double),
+    ]
+
 
 REAL = c_double
-CPLX = REAL*2
 PREAL = POINTER(REAL)
 PCPLX = POINTER(CPLX)
 
@@ -91,6 +97,7 @@ class c_GRT_MODEL1D(Structure):
         ('ircvup', c_bool),
         ('io_depth', c_bool),
 
+        ('omgref', CPLX),
         ('omega', CPLX),
         ('k', REAL),
         ('c_phase', CPLX),

pygrt/utils.py

@@ -737,17 +737,18 @@ def _compute_stress(st_syn:Stream):
     # 计算不同频率下的拉梅系数
     mus = np.zeros((nf,), dtype='c16')
     lams = np.zeros((nf,), dtype='c16')
-    omega = (REAL*2)(0.0, 0.0)
-    atte = (REAL*2)(0.0, 0.0)
+    omega = CPLX(0.0, 0.0)
+    atte = CPLX(0.0, 0.0)
+    omgref = CPLX(2.0*np.pi*(nf-1)*df, 0.0)
     for i in range(nf):
         freq = 0.01 if i==0 else df*i 
         w = 2.0*np.pi*freq 
-        omega[0] = w
-        C_grt_py_attenuation_law(Qbinv, omega, atte)
-        attb = atte[0] + atte[1]*1j
+        omega.real = w
+        atte = C_grt_attenuation_law(Qbinv, omgref, omega)
+        attb = atte.real + atte.imag*1j
         mus[i] = vb*vb*attb*attb*rho*1e10
-        C_grt_py_attenuation_law(Qainv, omega, atte)
-        atta = atte[0] + atte[1]*1j
+        atte = C_grt_attenuation_law(Qainv, omgref, omega)
+        atta = atte.real + atte.imag*1j
         lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i]
 
     del omega, atte