From 136a81b748a6c34d9e96358f83193e3332f4cedf Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 6 Aug 2026 13:54:32 +0800 Subject: [PATCH] FIX: guard dk at zero rmax and enforce GRT_MIN_NK in dynamic GF --- pygrt/C_extension/include/grt/common/const.h | 2 +- pygrt/C_extension/src/dynamic/grn.c | 20 ++++++++++++++++--- pygrt/C_extension/src/dynamic/grt_greenfn.c | 7 ++++--- .../src/static/grt_static_greenfn.c | 7 ++++--- pygrt/C_extension/src/static/static_grn.c | 6 +++--- pygrt/pymod.py | 12 ++++++----- 6 files changed, 36 insertions(+), 18 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index 905c6b7f..f578392c 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -132,7 +132,7 @@ typedef double complex cplx_t; #define GRT_GTYPES_MAX 2 ///< 2, 所有震源根据是否使用格林函数导数分为两类 -#define GRT_MIN_STATIC_NK 500 ///< 波数积分中 k0 部分最少的点数，若不够则对应调整 dk +#define GRT_MIN_NK 500 ///< 波数积分最少采样点数；nk 不足时收紧 dk typedef cplx_t cplxChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef cplx_t* pcplxChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index efea10da..33afeec6 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -68,7 +68,6 @@ void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_PROCESS *Kproc, GRNSPEC *grn, co const real_t Rho = mod1d->Rho[mod1d->isrc]; // 震源区密度 const real_t fac = 1.0/(4.0*PI*Rho); - const real_t dk = Kproc->dk; // 进度条变量 int progress=0; @@ -112,8 +111,6 @@ void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_PROCESS *Kproc, GRNSPEC *grn, co real_t w = grn->freqs[iw]*PI2; // 实频率 cplx_t omega = w - grn->wI*I; // 复数频率 omega = w - i*wI - cplx_t coef = - dk*fac / GRT_SQUARE(omega); // 最终要乘上的系数 - K_INTEG_PROCESS *local_Kproc = NULL; #ifdef _OPENMP // 定义局部对象 @@ -192,7 +189,24 @@ void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_PROCESS *Kproc, GRNSPEC *grn, co grt_printprogressBar("Computing Green Functions: ", progress*100/nf_valid); } } + + // 若 nk 不够，适当调整 dk（与静态解一致） + if(nk < GRT_MIN_NK){ + real_t new_dk = local_Kproc->kmax / GRT_MIN_NK; + #pragma omp critical(grn_console) + { + if(print_log){ + printf("\r\033[K"); + GRTRaiseInfo("iw=%zu, freq=%.3e: To increase nk(%zu) to %d, adjust dk from %.3e to %.3e", + iw, w/PI2, nk, GRT_MIN_NK, local_Kproc->dk, new_dk); + } + } + local_Kproc->dk = new_dk; + } } + // coef 必须用最终 dk（可能因 nk 下限被收紧） + cplx_t coef = - local_Kproc->dk * fac / GRT_SQUARE(omega); + K_INTEG *Kint = grt_wavenumber_integral(local_mstat, grn->nr, grn->rs, local_Kproc, grn->calc_upar, grt_kernel); // 记录到格林函数结构体内 diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index ec46d947..5561df5e 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -909,11 +909,11 @@ int greenfn_main(int argc, char **argv) { } // 自动选择积分间隔，默认使用传统离散波数积分 - // 自动选择会给出很保守的值（较大的Length） + // 自动选择会给出很保守的值（较大的Length）；rmax=0 时保留默认值 if(Ctrl->L.Length == 0.0){ Ctrl->L.Length = 15.0; real_t jus = GRT_SQUARE(vmax*tmax) - GRT_SQUARE(Ctrl->D.deprcv - Ctrl->D.depsrc); - if(jus >= 0.0){ + if(jus >= 0.0 && !GRT_IS_ZERO(rmax)){ Ctrl->L.Length = GRT_MAX(1.0 + sqrt(jus)/rmax + 0.5, Ctrl->L.Length); // +0.5为保守值 } } @@ -975,7 +975,8 @@ int greenfn_main(int argc, char **argv) { KPROC.kcut = Ctrl->L.kcut / rmax; - KPROC.dk = PI2 / (Ctrl->L.Length * rmax); + // rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + KPROC.dk = PI2 / (Ctrl->L.Length * GRT_MAX(rmax, GRT_ZERO_DISTANCE)); KPROC.applyFIM = Ctrl->L.FIM.active; KPROC.filondk = (Ctrl->L.FIM.active) ? PI2 / (Ctrl->L.FIM.Length * rmax) : 0.0; diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index f48cedf1..e4f10790 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -616,11 +616,12 @@ int static_greenfn_main(int argc, char **argv){ KPROC.keps = (Ctrl->C.convmet != K_INTEG_CONVERG_AUTO)? 0.0 : Ctrl->K.keps; // 如果使用了显式收敛方法，则不使用keps进行收敛判断 // 最大震中距 - real_t rmax = Ctrl->rs[grt_findMax_real_t(Ctrl->rs, Ctrl->nr)]; - + real_t rmax = Ctrl->rs[grt_findMax_real_t(Ctrl->rs, Ctrl->nr)]; + KPROC.kcut = Ctrl->L.kcut / rmax; - KPROC.dk = PI2 / (Ctrl->L.Length * rmax); + // rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + KPROC.dk = PI2 / (Ctrl->L.Length * GRT_MAX(rmax, GRT_ZERO_DISTANCE)); KPROC.applyFIM = Ctrl->L.FIM.active; KPROC.filondk = (Ctrl->L.FIM.active) ? PI2 / (Ctrl->L.FIM.Length * rmax) : 0.0; diff --git a/pygrt/C_extension/src/static/static_grn.c b/pygrt/C_extension/src/static/static_grn.c index 2c8a97f6..eaf2c178 100644 --- a/pygrt/C_extension/src/static/static_grn.c +++ b/pygrt/C_extension/src/static/static_grn.c @@ -122,9 +122,9 @@ void grt_integ_static_grn( GRTRaiseInfo("kmax = %.3e, (DWM)nk = %zu, kref = %.3e, ncount = %zu", static_kmax, nk, Kproc->k0, ncount); // 若 nk 不够，适当调整 dk - if(nk < GRT_MIN_STATIC_NK){ - real_t new_dk = static_kmax / GRT_MIN_STATIC_NK; - GRTRaiseInfo("To increase nk(%zu) to %d, adjust dk from %.3e to %.3e", nk, GRT_MIN_STATIC_NK, Kproc->dk, new_dk); + if(nk < GRT_MIN_NK){ + real_t new_dk = static_kmax / GRT_MIN_NK; + GRTRaiseInfo("To increase nk(%zu) to %d, adjust dk from %.3e to %.3e", nk, GRT_MIN_NK, Kproc->dk, new_dk); Kproc->dk = new_dk; } diff --git a/pygrt/pymod.py b/pygrt/pymod.py index cc3394fa..ba35b5ee 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -275,11 +275,11 @@ def _get_grn_spectra( tmax += rmax/delayV0 # 设置波数积分间隔 - # 自动情况下给出保守值 + # 自动情况下给出保守值；rmax≈0 时保留默认 Length if Length == 0.0: Length = 15.0 jus = (self.vmax*tmax)**2 - (depsrc - deprcv)**2 - if jus >= 0.0: + if jus >= 0.0 and rmax > ZERO_DISTANCE: Length = 1.0 + np.sqrt(jus)/rmax + 0.5 # 0.5作保守值 if Length < 15.0: Length = 15.0 @@ -327,8 +327,9 @@ def _get_grn_spectra( KPROC.vmin = vmin_ref KPROC.kcut = filonCut / rmax - - KPROC.dk = 2.0*np.pi / (Length * rmax) + + # rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + KPROC.dk = 2.0*np.pi / (Length * max(rmax, ZERO_DISTANCE)) KPROC.applyFIM = filonLength > 0.0 KPROC.filondk = 2.0*np.pi / (filonLength * rmax) if filonLength > 0.0 else 0.0 @@ -675,7 +676,8 @@ def compute_static_grn( # 最大震中距 rmax = np.max(rs) KPROC.kcut = filonCut / rmax - KPROC.dk = 2.0*np.pi / (Length * rmax) + # rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + KPROC.dk = 2.0*np.pi / (Length * max(rmax, ZERO_DISTANCE)) KPROC.applyFIM = filonLength > 0.0 KPROC.filondk = 2.0*np.pi / (filonLength * rmax) if filonLength > 0.0 else 0.0