Skip to content

FIX: guard dk at zero rmax and enforce GRT_MIN_NK in dynamic GF #282

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Aug 6, 2026
Merged

FIX: guard dk at zero rmax and enforce GRT_MIN_NK in dynamic GF #282

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion pygrt/C_extension/include/grt/common/const.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
20 changes: 17 additions & 3 deletions pygrt/C_extension/src/dynamic/grn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
// 定义局部对象
Expand Down Expand Up @@ -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);

// 记录到格林函数结构体内
Expand Down
7 changes: 4 additions & 3 deletions pygrt/C_extension/src/dynamic/grt_greenfn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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为保守值
}
}
Expand Down Expand Up @@ -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;
Expand Down
7 changes: 4 additions & 3 deletions pygrt/C_extension/src/static/grt_static_greenfn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
6 changes: 3 additions & 3 deletions pygrt/C_extension/src/static/static_grn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down
12 changes: 7 additions & 5 deletions pygrt/pymod.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading