diff --git a/docs/source/Advanced/integ_converg/run_dcm/plot_depth_kernel.py b/docs/source/Advanced/integ_converg/run_dcm/plot_depth_kernel.py index bedf5b2b..a0ca10c1 100644 --- a/docs/source/Advanced/integ_converg/run_dcm/plot_depth_kernel.py +++ b/docs/source/Advanced/integ_converg/run_dcm/plot_depth_kernel.py @@ -13,10 +13,9 @@ def gray_mapping(x): return (0.7 - x / 0.7) ** 0.8 def _plot_one(ax:Axes, pattern:str, ktype:str): - norm = 0.0 - yLst = [] - evdpLst = [] + plotLst = [] maxnorm = 0 + kmax = 0.0 for path in glob.glob(pattern): data = pygrt.utils.read_statsfile(path) @@ -26,16 +25,16 @@ def _plot_one(ax:Axes, pattern:str, ktype:str): Farr = data[ktype] norm = np.linalg.norm(np.real(Farr), np.inf) - yLst.append(np.real(Farr)) - evdpLst.append(evdp) + plotLst.append((karr, np.real(Farr), evdp)) maxnorm = max(maxnorm, norm) + kmax = max(kmax, karr[-1]) ax.axhline(y=0, color='k', ls='--', linewidth=1) - for y, evdp in zip(yLst, evdpLst): + for karr, y, evdp in plotLst: ax.plot(karr, y/maxnorm * 1e8, c=str(gray_mapping(evdp)), zorder=int(evdp*1e6)) - ax.set_xlim(0, karr[-1]) + ax.set_xlim(0, kmax * 0.4) ax.grid(linewidth=0.4) diff --git a/docs/source/Advanced/integ_converg/run_ptam/run.py b/docs/source/Advanced/integ_converg/run_ptam/run.py index 881915d9..520711c5 100755 --- a/docs/source/Advanced/integ_converg/run_ptam/run.py +++ b/docs/source/Advanced/integ_converg/run_ptam/run.py @@ -12,7 +12,7 @@ distarr = [5,8,10] # 设置 converg_method='PTAM' 进行收敛 stgrnLst = pymod.compute_grn( - distarr=distarr, nt=500, dt=0.02, converg_method='PTAM', + distarr=distarr, nt=500, dt=0.02, converg_method='PTAM', k0=2, ampk=1.2, use_kmax_ref=True, statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100] ) # END DEPSRC 0.0 DGRN @@ -46,7 +46,7 @@ xarr = np.array([2.0]) yarr = np.array([2.0]) -static_grn = pymod.compute_static_grn(xarr, yarr, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}") +static_grn = pymod.compute_static_grn(xarr, yarr, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}", k0=3, use_kmax_ref=True) ir = 0 statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"static_pygrtstats_{depsrc}_{deprcv}/PTAM_{ir:04d}_*/PTAM") diff --git a/docs/source/Advanced/k_integ/kmax.rst b/docs/source/Advanced/k_integ/kmax.rst index 8fe56183..24451147 100644 --- a/docs/source/Advanced/k_integ/kmax.rst +++ b/docs/source/Advanced/k_integ/kmax.rst @@ -10,22 +10,38 @@ P_m(\omega) = \Delta k \sum_{j=0}^{\infty} F_m(k_j,\omega)J_m(k_j r)k_j -其中 :math:`\Delta k = 2\pi/L, k_j=j\Delta k`，:math:`L` 为特征长度，即 :math:`L` 控制波数积分的积分间隔，默认根据 :doc:`/Tutorial/dynamic/gfunc` 部分介绍的约束条件自动确定。 +其中 :math:`\Delta k = 2\pi/L, k_j=j\Delta k`，:math:`L` 为特征长度，即 :math:`L` 控制波数积分的积分间隔， +默认根据 :doc:`/Tutorial/dynamic/gfunc` 部分介绍的约束条件自动确定。 -程序会取一个较大值 :math:`k_{\text{max}}` 作为波数积分的上限，经验公式为 +程序首先根据经验公式确定波数积分搜索区间的上界 :math:`k_{\text{max,ref}}` 。对动态全波解， .. math:: - k_{\text{max}} = \sqrt{k_0^2 + \left(s \cdot \dfrac{\omega}{v_{\text{min}}}\right)^2} + k_{\text{max,ref}} = \sqrt{\left(k_0 \cdot \dfrac{\pi}{h_s}\right)^2 + \left(s \cdot \dfrac{\omega}{v_{\text{min}}}\right)^2} 其中 -+ :math:`k_0` 为零频的波数积分上限，默认为 :math:`\dfrac{5\pi}{h_s}`， :math:`h_s` 为震源和场点的深度差（绝对值），限制最小为1km，若实际深度差小于该值，则自动使用直接收敛法（DCM）； ++ :math:`k_0` 为零频项的系数，默认为 50，程序内部使用 :math:`k_0 \cdot \pi / h_s` ； + :math:`h_s=\max(|z_s-z_r|, 0.1)` km 为震源和场点的深度差； + :math:`\omega` 为角频率； -+ :math:`v_{\text{min}}` 为参考最小速度，默认取自模型中的最小速度，且限制在0.1km/s以上。 -+ :math:`s` 为放大系数，默认为1.15； ++ :math:`v_{\text{min}}` 为参考最小速度，默认取自模型中的最小速度， + 且限制在 0.1 km/s 以上； ++ :math:`s` 为放大系数（``ampk``），默认为 2.0。 -程序支持自动判断积分收敛 (|yao1983|) 。当所有积分满足如下表达式时，自动退出波数循环（若达到 :math:`k_{\text{max}}` 则强制退出 ）， +对静态解，:math:`k_{\text{max,ref}} = k_0 \cdot \pi / h_s` 。 + +默认情况下，程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索 +实际积分上限 :math:`k_{\text{max}}` 。 +同深度时判断核函数是否逼近常数，异深度时判断振幅是否衰减至 0 。 +若指定 ``use_kmax_ref=True`` （C 模块 **+f**），则直接使用 :math:`k_{\text{max,ref}}` 作为 +:math:`k_{\text{max}}` 。 + +若振幅搜索达到 :math:`k_{\text{max,ref}}` 仍未满足收敛判据， +程序在默认（Auto）模式下将自动启用直接收敛法（DCM）处理积分收敛。 +当震源与场点完全同深度时，也会自动使用 DCM 。 + +程序还支持提前判断积分收敛 (|yao1983|) 。 +当所有积分满足如下表达式时，自动退出波数循环（若达到 :math:`k_{\text{max}}` 则强制退出 ）。 .. math:: @@ -46,13 +62,18 @@ .. group-tab:: Python - :func:`compute_grn() ` 函数支持以下可选参数来控制波数积分，具体说明详见API。 + :func:`compute_grn() ` 函数支持以下可选参数来控制波数积分， + 具体说明详见API。 - + ``k0:float``, 对应公式中的 :math:`k_0` 的系数，默认为5 - + ``ampk:float``, 对应公式中的 :math:`s` ，默认为1.15 - + ``keps:float`` 对应公式中的 :math:`\epsilon`，默认为-1（不使用） + + ``k0:float``, 对应公式中零频项的系数 :math:`k_0` ，默认为 50 + + ``ampk:float``, 对应公式中的 :math:`s` ，默认为 2.0 + + ``keps:float`` 对应公式中的 :math:`\epsilon`，默认为 -1（不使用） + + ``use_kmax_ref:bool`` 为 True 时直接使用 :math:`k_{\text{max,ref}}` 作为 + 积分上限 - :func:`compute_static_grn() ` 函数支持以下可选参数来控制波数积分，参数与上面对应，具体说明详见API。 + :func:`compute_static_grn() ` 函数支持以下可选参数来控制波数积分， + 参数与上面对应，具体说明详见API。 + ``k0:float`` + ``keps:float`` + + ``use_kmax_ref:bool`` diff --git a/docs/source/Advanced/kernel_old/run/run.py b/docs/source/Advanced/kernel_old/run/run.py index 2029c1d8..91f5a8cb 100755 --- a/docs/source/Advanced/kernel_old/run/run.py +++ b/docs/source/Advanced/kernel_old/run/run.py @@ -12,7 +12,7 @@ # 不指定statsidx，默认输出全部频率点的积分过程文件 # vmin_ref 显式给定参考速度（用于定义波数积分上限），避免使用PTAM # Length 给定波数积分间隔dk -_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, k0_is_fixed=True, converg_method='none', statsfile="pygrtstats") +_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, use_kmax_ref=True, converg_method='none', statsfile="pygrtstats") # END GRN # ----------------------------------------------------------------- diff --git a/docs/source/Gallery/ex17/run.sh b/docs/source/Gallery/ex17/run.sh index 2716b952..4bc61776 100644 --- a/docs/source/Gallery/ex17/run.sh +++ b/docs/source/Gallery/ex17/run.sh @@ -13,8 +13,8 @@ cat > halfspace_Q <} \cdot \left(\dfrac{\omega}{v_{\text{min}}}\right)^2} - - + **+k**\ *k0* - 控制零频的积分上限 [5.0]，其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 1.0)` 。 - + **+s**\ *ampk* - 放大倍数 [1.15] 。 - + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，详见 - Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。 - + **+v**\ *vmin* - 参考最小速度，默认 :math:`\max{(\min\limits_{i} (\alpha_i \cup \beta_i), 0.1)}` 。 + k_{\text{max,ref}} = \sqrt{ \left(k_0 \cdot \dfrac{\pi}{\Delta h}\right)^2 + \textit{}^2 \cdot \left(\dfrac{\omega}{v_{\text{min}}}\right)^2 } + + 程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索实际积分上限 :math:`k_{\text{max}}` 。 + 若搜索达到 :math:`k_{\text{max,ref}}` 仍未收敛，或震源与场点完全同深度时， + 默认模式下将自动启用 DCM 。 + + + **+k**\ *k0* - 零频项系数 [50.0]， + 其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 0.1)` 。 + + **+f** - 直接使用 :math:`k_{\text{max,ref}}` 作为积分上限， + 不进行振幅搜索。 + + **+s**\ *ampk* - 放大倍数 [2.0] 。 + + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]， + 详见 Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。 + + **+v**\ *vmin* - 参考最小速度， + 默认 :math:`\max{(\min\limits_{i} (\alpha_i \cup \beta_i), 0.1)}` 。 .. include:: explain_-Cconverg.rst_ diff --git a/docs/source/Module/static_greenfn.rst b/docs/source/Module/static_greenfn.rst index 4bffe7b1..cbd7a1fc 100644 --- a/docs/source/Module/static_greenfn.rst +++ b/docs/source/Module/static_greenfn.rst @@ -22,7 +22,7 @@ static_greenfn [ |-B|\ **f|F|r|R|h|H** ] [ |-L|\ *length*\ [**+l**\ *Flength*][**+a**\ *Ftol*][**+o**\ *offset*] ] [ |-C|\ **d|p|n** ] -[ |-K|\ [**+k**\ *k0*][**+e**\ *keps*] ] +[ |-K|\ [**+k**\ *k0*][**+f**][**+e**\ *keps*] ] [ |-S| ] [ **-e** ] [ **-h** ] @@ -66,12 +66,19 @@ static_greenfn .. _-K: -**-K**\ [**+k**\ *k0*][**+e**\ *keps*] - 控制波数积分上限 :math:`k_0 \cdot \dfrac{\pi}{\Delta h}` +**-K**\ [**+k**\ *k0*][**+f**][**+e**\ *keps*] + 控制波数积分搜索区间的上界 :math:`k_{\text{max,ref}} = k_0 \cdot \dfrac{\pi}{\Delta h}` - + **+k**\ *k0* - 控制零频的积分上限 [5.0]，其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 1.0)` 。 - + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，详见 - Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。 + 程序在 :math:`[\Delta k, k_{\text{max,ref}}]` 内基于核函数振幅搜索实际积分上限 :math:`k_{\text{max}}` 。 + 若搜索达到 :math:`k_{\text{max,ref}}` 仍未收敛，或震源与场点完全同深度时， + 默认模式下将自动启用 DCM 。 + + + **+k**\ *k0* - 零频项系数 [50.0]， + 其中深度差 :math:`\Delta h = \max(|z_s - z_r|, 0.1)` 。 + + **+f** - 直接使用 :math:`k_{\text{max,ref}}` 作为积分上限， + 不进行振幅搜索。 + + **+e**\ *keps* - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]， + 详见 Yao and Harkrider (1983) 和 :doc:`/Advanced/k_integ/kmax` 。 .. include:: explain_-Cconverg.rst_ diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index ba2b690e..6da83b5d 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -55,7 +55,7 @@ typedef double complex cplx_t; #define DEG1 0.017453292519943295 ///< \f$ \frac{\pi}{180} \f$ #define GOLDEN_RATIO 0.6180339887498949 ///< \f$ \frac{\sqrt{5}-1}{2} \f$ -#define GRT_MIN_DEPTH_GAP_SRC_RCV 1.0 ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限，以及判断是否需要其它收敛方法） +#define GRT_MIN_DEPTH_GAP_SRC_RCV 0.1 ///< 震源和台站的最小深度差（不做绝对限制，仅用于 kmax_ref 中 hs 的下限） #define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数 #define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值 diff --git a/pygrt/C_extension/include/grt/integral/integ_process.h b/pygrt/C_extension/include/grt/integral/integ_process.h index 47cec29c..7e0733ac 100644 --- a/pygrt/C_extension/include/grt/integral/integ_process.h +++ b/pygrt/C_extension/include/grt/integral/integ_process.h @@ -34,15 +34,19 @@ typedef enum { // 描述不同波数积分方法的结构体 typedef struct { - real_t k0; ///< 波数积分的上限 \f$ \tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2} \f$ ，k循环必须退出, hs=max(震源和台站深度差,1.0) - bool k0_is_fixed; ///< 固定 k0，默认在程序中自动调整 k0 - real_t ampk; ///< 影响波数k积分上限的系数 - real_t keps; ///< 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 - real_t vmin; ///< 参考最小速度，用于定义波数积分的上限 + real_t k0; ///< 用户参数 k0 经 \f$ \pi/hs \f$ 缩放后的零频项， + ///< 与 \f$ ampk*\omega/vmin \f$ 共同确定搜索上界 kmax_ref； + ///< hs=max(震源和台站深度差, 0.1) + bool use_kmax_ref; ///< 为 true 时直接将 kmax_ref 作为积分上限，不进行基于振幅的搜索 + real_t ampk; ///< 影响 kmax_ref 中频率相关项的系数，默认 2.0 + real_t keps; ///< 波数积分的收敛条件，要求在某震中距下所有格林函数都收敛； + ///< 为负数代表不提前判断收敛 + real_t vmin; ///< 参考最小速度，用于定义 kmax_ref real_t kcut; ///< 波数积分和Filon积分的分割点 - real_t kmax; ///< 全局波数最大值，程序运行中会随频率变动 + real_t kmax; ///< 实际波数积分上限，默认在 [dk, kmax_ref] 内由振幅搜索确定， + ///< 程序运行中会随频率变动 real_t dk; ///< DWM 的波数积分间隔 diff --git a/pygrt/C_extension/include/grt/integral/kmax.h b/pygrt/C_extension/include/grt/integral/kmax.h index 2ee2fe01..ef16e570 100644 --- a/pygrt/C_extension/include/grt/integral/kmax.h +++ b/pygrt/C_extension/include/grt/integral/kmax.h @@ -18,14 +18,14 @@ * * 搜索在 log(k) 上等间距推进（自适应几何因子），目标约 40 步覆盖 [kmax_init, kmax_ref] * 在 kmax_low 之前不判断收敛 - * 同深度使用当前核函数与前一采样点的差值 dF 相对于搜索过程中的最大振幅 Fmax 判断逼近常数； - * 异深度用振幅相对峰值衰减判断逼近 0 + * 同深度使用当前核函数与前一采样点的差值 dF 相对于搜索过程中的最大振幅 Fmax + * 判断逼近常数；异深度用振幅相对峰值衰减判断逼近 0 * * @param[in,out] mstat 已设置频率的模型状态 * @param[in] kerfunc 待检查的核函数 * @param[in] kmax_init 扫描的初始波数 * @param[in] kmax_low 允许判断收敛的最低波数 - * @param[in] kmax_ref 最大上限 + * @param[in] kmax_ref 搜索区间的上界（经验公式） * @param[out] Ncount 估计过程中计算核函数的次数，可为 NULL * @return 估计的 kmax */ diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index 163622a1..efea10da 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -150,9 +150,9 @@ void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_PROCESS *Kproc, GRNSPEC *grn, co // =================================================================================== // Wavenumber Integration - // 每个频率均直接根据动态核函数估计积分上限 + // 每个频率根据 kmax_ref 与核函数振幅搜索实际积分上限 real_t kmax_ref = hypot(local_Kproc->k0, local_Kproc->ampk * w / local_Kproc->vmin); - if(local_Kproc->k0_is_fixed){ + if(local_Kproc->use_kmax_ref){ local_Kproc->kmax = kmax_ref; size_t nk = floor(local_Kproc->kmax / local_Kproc->dk) + 1; #pragma omp critical(grn_console) diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 6907db91..7cd95166 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -18,7 +18,7 @@ #define GRT_GREENFN_H_FREQ2 -1.0 #define GRT_GREENFN_K_VMIN 0.1 #define GRT_GREENFN_K_K0 50.0 -#define GRT_GREENFN_K_AMPK 1.15 +#define GRT_GREENFN_K_AMPK 2.0 #define GRT_GREENFN_G_EX true #define GRT_GREENFN_G_VF true #define GRT_GREENFN_G_HF true @@ -103,7 +103,7 @@ typedef struct { real_t ampk; real_t k0; real_t vmin; - bool k0_is_fixed; + bool use_kmax_ref; } K; /** 时间延迟 */ struct { @@ -303,11 +303,11 @@ printf("\n" " [0, k*] and [k*, kmax], in which k*=/rmax,\n" " the former uses DWM and the latter uses FIM/SAFIM.\n" "\n" -" -Cd|p|n Set convergence method,\n" +" -Cd|p|n Set global convergence method.\n" " + d: Direct Convergence Method (DCM).\n" " + p: Peak-Trough Averaging Method (PTAM).\n" " + n: None.\n" -" Default use -Cd when fabs(depsrc-deprcv) <= %.1f.\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( +" DCM may still be applied internally when needed.\n" "\n" " -E[p][/]\n" " Introduce the time shift in results. The times series \n" @@ -323,22 +323,23 @@ printf("\n" " If -E is not used, then the first time sample will be the origin time.\n" "\n" " -K[+k][+f][+s][+e][+v]\n" -" Define the wavenumber integration upper bound\n" -" sqrt( ^2 + (*w/)^2 ),\n" -" : maximum upper bound residual k for 0 frequency, \n" +" Define the reference upper bound of wavenumber integration\n" +" kmax_ref = sqrt( (*PI/hs)^2 + (*w/)^2 ),\n" +" : coefficient for the zero-frequency term in kmax_ref,\n" " default is %.1f, and multiply PI/hs in program, \n", GRT_GREENFN_K_K0); printf( " where hs = max(fabs(depsrc-deprcv), %.1f).\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( -" The program will choose the proper residual in [0, k0].\n" -" If k0 is not enough, convergence method will be applied.\n" -" If use +f, directly set k0 as the residual.\n" +" The program searches kmax in [dk, kmax_ref] based on\n" +" kernel amplitude. If the search reaches kmax_ref without\n" +" convergence, or source and receiver are at the same depth,\n" +" DCM will be applied (in Auto mode).\n" +" If use +f, directly set kmax to kmax_ref.\n" " : amplification factor, default is %.2f.\n", GRT_GREENFN_K_AMPK); printf( " : a threshold for break wavenumber \n" " integration in advance. See \n" " (Yao and Harkrider, 1983) for details.\n" " Default 0.0 not use.\n" " : Minimum velocity (km/s) for reference. This\n" -" is designed to define the upper bound \n" -" of wavenumber integration.\n" +" is designed to define kmax_ref.\n" " There are 2 cases:\n" " + (default) not set or set 0.0.\n"); printf( " will be the minimum velocity\n" @@ -667,7 +668,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ break; case 'f': - Ctrl->K.k0_is_fixed = true; + Ctrl->K.use_kmax_ref = true; break; default: @@ -957,7 +958,7 @@ int greenfn_main(int argc, char **argv) { { real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); KPROC.k0 = Ctrl->K.k0 * PI / hs; - KPROC.k0_is_fixed = Ctrl->K.k0_is_fixed; + KPROC.use_kmax_ref = Ctrl->K.use_kmax_ref; KPROC.ampk = Ctrl->K.ampk; KPROC.keps = (Ctrl->C.convmet != K_INTEG_CONVERG_AUTO)? 0.0 : Ctrl->K.keps; // 如果使用了显式收敛方法，则不使用keps进行收敛判断 KPROC.vmin = Ctrl->K.vmin; diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index e2785d1e..22c7a08c 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -61,7 +61,7 @@ typedef struct { bool active; real_t keps; real_t k0; - bool k0_is_fixed; + bool use_kmax_ref; } K; /** 波数积分过程的核函数文件 */ struct { @@ -203,11 +203,11 @@ printf("\n" " [0, k*] and [k*, kmax], in which k*=/rmax,\n" " the former uses DWM and the latter uses FIM/SAFIM.\n" "\n" -" -Cd|p|n Set convergence method,\n" +" -Cd|p|n Set global convergence method.\n" " + d: Direct Convergence Method (DCM).\n" " + p: Peak-Trough Averaging Method (PTAM).\n" " + n: None.\n" -" Default use -Cd when fabs(depsrc-deprcv) <= %.1f.\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( +" DCM may still be applied internally when needed.\n" "\n" " -Bf|F|r|R|h|H\n" " Boundary condition of top layer (lowercase) and\n" @@ -217,12 +217,15 @@ printf("\n" " h|H: Halfspace.\n" "\n" " -K[+k][+f][+e]\n" -" Define the wavenumber integration upperbound\n" -" : maximum upper bound, default is %.1f, and multiply PI/hs \n", GRT_GREENFN_K_K0); printf( +" Define the reference upper bound of wavenumber integration\n" +" kmax_ref = * PI / hs,\n" +" : coefficient, default is %.1f, and multiply PI/hs \n", GRT_GREENFN_K_K0); printf( " in program, where hs = max(fabs(depsrc-deprcv), %.1f).\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( -" The program will choose the proper upper bound in [0, k0].\n" -" If k0 is not enough, convergence method will be applied.\n" -" If use +f, directly set k0 as the upper bound.\n" +" The program searches kmax in [dk, kmax_ref] based on\n" +" kernel amplitude. If the search reaches kmax_ref without\n" +" convergence, or source and receiver are at the same depth,\n" +" DCM will be applied (in Auto mode).\n" +" If use +f, directly set kmax to kmax_ref.\n" " : a threshold for break wavenumber \n" " integration in advance. See \n" " (Yao and Harkrider, 1983) for details.\n" @@ -408,7 +411,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ break; case 'f': - Ctrl->K.k0_is_fixed = true; + Ctrl->K.use_kmax_ref = true; break; default: @@ -599,7 +602,7 @@ int static_greenfn_main(int argc, char **argv){ { real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); KPROC.k0 = Ctrl->K.k0 * PI / hs; - KPROC.k0_is_fixed = Ctrl->K.k0_is_fixed; + KPROC.use_kmax_ref = Ctrl->K.use_kmax_ref; KPROC.keps = (Ctrl->C.convmet != K_INTEG_CONVERG_AUTO)? 0.0 : Ctrl->K.keps; // 如果使用了显式收敛方法，则不使用keps进行收敛判断 // 最大震中距 diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index 7db410a9..2c8a97f6 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -105,8 +105,8 @@ void grt_integ_static_grn( // =================================================================================== // Wavenumber Integration - // 波数积分上限 - if(Kproc->k0_is_fixed){ + // 波数积分上限：在 [dk, kmax_ref] 内基于核函数振幅搜索 + if(Kproc->use_kmax_ref){ Kproc->kmax = Kproc->k0; size_t nk = floor(Kproc->kmax / Kproc->dk) + 1; GRTRaiseInfo("kmax = %.3e, nk = %zu", Kproc->kmax, nk); diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index 83246b39..158eda07 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -23,6 +23,7 @@ "ZNEchs", "qwvchs", "MECHANISM_NUM", + "MIN_DEPTH_GAP_SRC_RCV", "NPCT_REAL_TYPE", "NPCT_CMPLX_TYPE", @@ -50,6 +51,7 @@ ZNEchs = ['Z', 'N', 'E'] qwvchs = ['q', 'w', 'v'] MECHANISM_NUM = 6 +MIN_DEPTH_GAP_SRC_RCV = 0.1 NPCT_REAL_TYPE = 'f8' NPCT_CMPLX_TYPE = 'c16' @@ -109,11 +111,14 @@ class c_MODEL1D(Structure): class c_K_INTEG_PROCESS(Structure): """ 和C结构体 K_INTEG_PROCESS 作匹配 + + 波数积分参数。k0、ampk、vmin 共同确定搜索上界 kmax_ref； + 默认在 [dk, kmax_ref] 内基于核函数振幅搜索实际积分上限 kmax。 """ _fields_ = [ ('k0', REAL), - ('k0_is_fixed', c_bool), + ('use_kmax_ref', c_bool), ('ampk', REAL), ('keps', REAL), ('vmin', REAL), diff --git a/pygrt/pymod.py b/pygrt/pymod.py index ae40f8ff..b1046beb 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -175,9 +175,9 @@ def _get_grn_spectra( keepAllFreq:bool=False, vmin_ref:float=0.0, keps:float=-1.0, - ampk:float=1.15, + ampk:float=2.0, k0:float=50.0, - k0_is_fixed:bool=False, + use_kmax_ref:bool=False, Length:float=0.0, filonLength:float=0.0, safilonTol:float=0.0, @@ -318,9 +318,9 @@ def _get_grn_spectra( # ==================================================================== KPROC = c_K_INTEG_PROCESS() - hs = max(abs(depsrc - deprcv), 1.0) + hs = max(abs(depsrc - deprcv), MIN_DEPTH_GAP_SRC_RCV) KPROC.k0 = k0 * np.pi / hs - KPROC.k0_is_fixed = k0_is_fixed + KPROC.use_kmax_ref = use_kmax_ref KPROC.ampk = ampk KPROC.keps = keps if converg_method.upper() != 'AUTO' else 0.0 KPROC.vmin = vmin_ref @@ -462,9 +462,9 @@ def compute_grn( keepAllFreq:bool=False, vmin_ref:float=0.0, keps:float=-1.0, - ampk:float=1.15, + ampk:float=2.0, k0:float=50.0, - k0_is_fixed:bool=False, + use_kmax_ref:bool=False, Length:float=0.0, filonLength:float=0.0, safilonTol:float=0.0, @@ -492,12 +492,18 @@ def compute_grn( :math:`\tilde{\omega} = \omega - j*w_I, w_I = \zeta*\pi/T, T=nt*dt` . see Bouchon (1981) and 张海明 (2021) for more details and tests. :param keepAllFreq: calculate all frequency points, no matter how low the frequency is - :param vmin_ref: minimum reference velocity (km/s). the default vmin=max(minimum velocity, 0.1), used to define the upper bound of k integral + :param vmin_ref: minimum reference velocity (km/s). + the default vmin=max(minimum velocity, 0.1), used to define kmax_ref :param keps: automatic convergence condition, see Yao and Harkrider (1983) for more details. negative value denotes not use. - :param ampk: The factor that affect the upper bound of the k integral, see below. - :param k0: k0 used to define the maximum offset of upper bound :math:`\tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2}` , hs=max(abs(depsrc-deprcv),1.0) - :param k0_is_fixed: directly use k0, rather than choosing a proper offset in [0, k0] + :param ampk: amplification factor in kmax_ref, see below. + :param k0: coefficient in kmax_ref + :math:`k_{\text{max,ref}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*\omega/vmin_{ref})^2}` , + hs=max(abs(depsrc-deprcv),0.1). + The actual kmax is searched in [dk, kmax_ref] based on kernel amplitude; + if the search reaches kmax_ref without convergence, + or source and receiver are at the same depth, DCM is applied in Auto mode. + :param use_kmax_ref: directly use kmax_ref as kmax, without amplitude search :param Length: integration step `dk=2\pi / (L*rmax)`, see Bouchon (1981) and 张海明 (2021) for the criterion, default set automatically. :param filonLength: integration step of Fixed-Interval Filon's Integration Method :param safilonTol: precision of Self-Adaptive Filon's Integration Method @@ -524,7 +530,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, k0_is_fixed, Length, filonLength, safilonTol, filonCut, converg_method, + vmin_ref, keps, ampk, k0, use_kmax_ref, Length, filonLength, safilonTol, filonCut, converg_method, delayT0, delayV0, calc_upar, statsfile, statsidxs, print_log ) @@ -545,7 +551,7 @@ def compute_static_grn( distarr:Union[np.ndarray,List[float],float,None]=None, keps:float=-1.0, k0:float=50.0, - k0_is_fixed:bool=False, + use_kmax_ref:bool=False, Length:float=15.0, filonLength:float=0.0, safilonTol:float=0.0, @@ -565,8 +571,12 @@ def compute_static_grn( :param distarr: equal to "xarr=[0.0], yarr=distarr" :param keps: automatic convergence condition, see (Yao and Harkrider (1983) for more details. negative value denotes not use. - :param k0: k0 used to define the maximum offset of upper bound :math:`\tilde{k_{max}}=(k_{0}*\pi/hs)^2`, hs=max(abs(depsrc-deprcv),1.0) - :param k0_is_fixed: directly use k0, rather than choosing a proper offset in [0, k0] + :param k0: coefficient in kmax_ref :math:`k_{\text{max,ref}}=k_{0}*\pi/hs`, + hs=max(abs(depsrc-deprcv),0.1). + The actual kmax is searched in [dk, kmax_ref] based on kernel amplitude; + if the search reaches kmax_ref without convergence, + or source and receiver are at the same depth, DCM is applied in Auto mode. + :param use_kmax_ref: directly use kmax_ref as kmax, without amplitude search :param Length: integration step `dk=2\pi / (L*rmax)`, default L=15 :param filonLength: integration step of Fixed-Interval Filon's Integration Method :param safilonTol: precision of Self-Adaptive Filon's Integration Method @@ -652,9 +662,9 @@ def compute_static_grn( # ==================================================================== KPROC = c_K_INTEG_PROCESS() - hs = max(abs(depsrc - deprcv), 1.0) + hs = max(abs(depsrc - deprcv), MIN_DEPTH_GAP_SRC_RCV) KPROC.k0 = k0 * np.pi / hs - KPROC.k0_is_fixed = k0_is_fixed + KPROC.use_kmax_ref = use_kmax_ref KPROC.keps = keps if converg_method.upper() != 'AUTO' else 0.0 # 最大震中距