From eec774fe205a133cc1cfcefdcb44d4cb0edb340f Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 7 Aug 2026 10:11:47 +0800 Subject: [PATCH] REFAC: share GF preparation in C for dynamic/static Python paths --- pygrt/C_extension/src/dynamic/grt_greenfn.c | 247 +++++++++------- .../src/static/grt_static_greenfn.c | 87 ++++-- pygrt/c_interfaces.py | 31 ++ pygrt/c_structures.py | 4 - pygrt/pymod.py | 266 ++++++------------ 5 files changed, 312 insertions(+), 323 deletions(-) diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 5561df5e..11166664 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -846,17 +846,122 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ fclose(fp); GRT_SAFE_FREE_PTR(dummy); +} + + + +/** + * 频谱积分前的共享准备：自动 Length、wI、频点网格、nf1/nf2、填充 KPROC 与 GRNSPEC 元数据。 + * grn->freqs 由本函数 malloc，调用方负责释放；grn->rs 借用调用方指针。 + * 不分配 u/uiz/uir，不处理 stats 路径。 + */ +void grt_prepare_grn_spec( + MODEL1D *mod1d, + size_t nr, real_t *rs, + size_t nt, real_t dt, real_t zeta, bool keepAllFreq, + real_t freq1, real_t freq2, + real_t Length, + real_t filonLength, real_t safilonTol, real_t filonCut, + real_t k0, real_t ampk, real_t keps, real_t vmin_ref, bool use_kmax_ref, + int convmet, + real_t delayT0, real_t delayV0, bool refFirstP, + bool calc_upar, + K_INTEG_PROCESS *Kproc, + GRNSPEC *grn) +{ // FIM/SAFIM 面向远震中距，公式含 1/r、1/√r，不能用于 r=0（含 kcut 分段） - if(Ctrl->L.FIM.active || Ctrl->L.SAFIM.active){ - for(size_t ir=0; irR.nr; ++ir){ - if(GRT_IS_ZERO(Ctrl->R.rs[ir])){ - GRTBadOptionError(L, "FIM/SAFIM cannot be used with zero epicentral distance."); + if(filonLength > 0.0 || safilonTol > 0.0){ + for(size_t ir = 0; ir < nr; ++ir){ + if(GRT_IS_ZERO(rs[ir])){ + GRTRaiseError("FIM/SAFIM cannot be used with zero epicentral distance."); } } } -} + real_t vmin, vmax; + grt_get_mod1d_vmin_vmax(mod1d, &vmin, &vmax); + // 参考最小速度 + if(vmin_ref == 0.0){ + vmin_ref = GRT_MAX(vmin, GRT_GREENFN_K_VMIN); + } + + real_t winT = nt * dt; + real_t rmax = rs[grt_findMax_real_t(rs, nr)]; + + // 时窗最大截止时刻（用于自动 Length） + real_t tmax = 0.0; + if(!refFirstP){ + tmax = delayT0 + winT; + if(delayV0 > 0.0) tmax += rmax / delayV0; + } else { + real_t maxP = grt_compute_travt1d(mod1d->Thk, mod1d->Va, mod1d->n, mod1d->isrc, mod1d->ircv, rmax); + tmax = delayT0 + maxP + winT; + } + + // 自动选择积分间隔；rmax=0 时保留默认值 + if(Length == 0.0){ + Length = 15.0; + real_t jus = GRT_SQUARE(vmax * tmax) - GRT_SQUARE(mod1d->deprcv - mod1d->depsrc); + if(jus >= 0.0 && !GRT_IS_ZERO(rmax)){ + Length = GRT_MAX(1.0 + sqrt(jus) / rmax + 0.5, Length); // +0.5为保守值 + } + } + + real_t wI = zeta * PI / winT; + + size_t nf = nt / 2 + 1; + real_t df = 1.0 / winT; + real_t *freqs = (real_t *)malloc(nf * sizeof(real_t)); + for(size_t i = 0; i < nf; ++i){ + freqs[i] = i * df; + } + + size_t nf1 = 0; + size_t nf2 = nf - 1; + if(freq1 > 0.0){ + nf1 = GRT_MIN(ceil(freq1 / df), nf - 1); + } + if(freq2 > 0.0){ + nf2 = GRT_MIN(floor(freq2 / df), nf - 1); + } + nf2 = GRT_MAX(nf1, nf2); + + memset(Kproc, 0, sizeof(*Kproc)); + { + real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); + Kproc->k0 = k0 * PI / hs; + Kproc->use_kmax_ref = use_kmax_ref; + Kproc->ampk = ampk; + // 显式收敛方法时不使用 keps + Kproc->keps = (convmet != K_INTEG_CONVERG_AUTO) ? 0.0 : keps; + Kproc->vmin = vmin_ref; + + Kproc->kcut = filonCut / rmax; + + // rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + Kproc->dk = PI2 / (Length * GRT_MAX(rmax, GRT_ZERO_DISTANCE)); + + Kproc->applyFIM = filonLength > 0.0; + Kproc->filondk = (filonLength > 0.0) ? PI2 / (filonLength * rmax) : 0.0; + + Kproc->applySAFIM = safilonTol > 0.0; + Kproc->sa_tol = safilonTol; + + Kproc->cvgmet = convmet; + } + + memset(grn, 0, sizeof(*grn)); + grn->nf = nf; + grn->freqs = freqs; + grn->nf1 = nf1; + grn->nf2 = nf2; + grn->nr = nr; + grn->rs = rs; + grn->wI = wI; + grn->keepAllFreq = keepAllFreq; + grn->calc_upar = calc_upar; +} /** 子模块主函数 */ @@ -881,53 +986,32 @@ int greenfn_main(int argc, char **argv) { Ctrl->G.doHF = Ctrl->G.doVF = Ctrl->G.doDC = false; } - // 最大最小速度 - real_t vmin, vmax; - grt_get_mod1d_vmin_vmax(mod1d, &vmin, &vmax); - - // 参考最小速度 - if(Ctrl->K.vmin == 0.0){ - Ctrl->K.vmin = GRT_MAX(vmin, GRT_GREENFN_K_VMIN); - } - - // 时窗长度 - Ctrl->N.winT = Ctrl->N.nt*Ctrl->N.dt; - - // 最大震中距 - real_t rmax = Ctrl->R.rs[grt_findMax_real_t(Ctrl->R.rs, Ctrl->R.nr)]; - - // 时窗最大截止时刻 - real_t tmax = 0.0; - { - if (! Ctrl->E.refFirstP){ - tmax = Ctrl->E.delayT0 + Ctrl->N.winT; - if(Ctrl->E.delayV0 > 0.0) tmax += rmax/Ctrl->E.delayV0; - } else { - real_t maxP = grt_compute_travt1d(mod1d->Thk, mod1d->Va, mod1d->n, mod1d->isrc, mod1d->ircv, rmax); - tmax = Ctrl->E.delayT0 + maxP + Ctrl->N.winT; - } - } - - // 自动选择积分间隔，默认使用传统离散波数积分 - // 自动选择会给出很保守的值（较大的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 && !GRT_IS_ZERO(rmax)){ - Ctrl->L.Length = GRT_MAX(1.0 + sqrt(jus)/rmax + 0.5, Ctrl->L.Length); // +0.5为保守值 - } - } - - // 虚频率 - Ctrl->N.wI = Ctrl->N.zeta*PI/Ctrl->N.winT; - - // 定义要计算的频率、时窗等 - Ctrl->N.nf = Ctrl->N.nt/2 + 1; - Ctrl->N.df = 1.0/Ctrl->N.winT; - Ctrl->N.freqs = (real_t*)malloc(Ctrl->N.nf*sizeof(real_t)); - for(size_t i=0; iN.nf; ++i){ - Ctrl->N.freqs[i] = i*Ctrl->N.df; - } + K_INTEG_PROCESS KPROC = {0}; + GRNSPEC grn_storage = {0}; + GRNSPEC *grn = &grn_storage; + grt_prepare_grn_spec( + mod1d, + Ctrl->R.nr, Ctrl->R.rs, + Ctrl->N.nt, Ctrl->N.dt, Ctrl->N.zeta, Ctrl->N.keepAllFreq, + Ctrl->H.freq1, Ctrl->H.freq2, + Ctrl->L.Length, + Ctrl->L.FIM.active ? Ctrl->L.FIM.Length : 0.0, + Ctrl->L.SAFIM.active ? Ctrl->L.SAFIM.tol : 0.0, + Ctrl->L.kcut, + Ctrl->K.k0, Ctrl->K.ampk, Ctrl->K.keps, Ctrl->K.vmin, Ctrl->K.use_kmax_ref, + Ctrl->C.convmet, + Ctrl->E.delayT0, Ctrl->E.delayV0, Ctrl->E.refFirstP, + Ctrl->e.active, + &KPROC, grn); + + // 写回 Ctrl，供后续 IFFT/SAC/free_Ctrl 使用 + Ctrl->N.winT = Ctrl->N.nt * Ctrl->N.dt; + Ctrl->N.nf = grn->nf; + Ctrl->N.df = 1.0 / Ctrl->N.winT; + Ctrl->N.freqs = grn->freqs; + Ctrl->N.wI = grn->wI; + Ctrl->H.nf1 = grn->nf1; + Ctrl->H.nf2 = grn->nf2; // 如果只传入了 -S, 未指定索引，则默认所有频率索引 if(Ctrl->S.active && Ctrl->S.statsidxs == NULL){ @@ -940,17 +1024,6 @@ int greenfn_main(int argc, char **argv) { Ctrl->S.s_raw = strdup("(all)"); } - // 自定义频段 - Ctrl->H.nf1 = 0; - Ctrl->H.nf2 = Ctrl->N.nf-1; - if(Ctrl->H.freq1 > 0.0){ - Ctrl->H.nf1 = GRT_MIN(ceil(Ctrl->H.freq1/Ctrl->N.df), Ctrl->N.nf-1); - } - if(Ctrl->H.freq2 > 0.0){ - Ctrl->H.nf2 = GRT_MIN(floor(Ctrl->H.freq2/Ctrl->N.df), Ctrl->N.nf-1); - } - Ctrl->H.nf2 = GRT_MAX(Ctrl->H.nf1, Ctrl->H.nf2); - // 波数积分中间文件输出目录 if(Ctrl->S.active){ Ctrl->S.s_statsdir = NULL; @@ -962,50 +1035,10 @@ int greenfn_main(int argc, char **argv) { GRTCheckMakeDir(Ctrl->S.s_statsdir); } - - // 波数积分方法 - K_INTEG_PROCESS KPROC = {0}; - { - real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); - KPROC.k0 = Ctrl->K.k0 * PI / hs; - 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; - - KPROC.kcut = Ctrl->L.kcut / 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; - - KPROC.applySAFIM = Ctrl->L.SAFIM.active; - KPROC.sa_tol = Ctrl->L.SAFIM.tol; - - KPROC.cvgmet = Ctrl->C.convmet; - } - - // 格林函数频谱 - GRNSPEC *grn = &(GRNSPEC){0}; - { - grn->nf = Ctrl->N.nf; - grn->freqs = Ctrl->N.freqs; - grn->nf1 = Ctrl->H.nf1; - grn->nf2 = Ctrl->H.nf2; - grn->nr = Ctrl->R.nr; - grn->rs = Ctrl->R.rs; - grn->wI = Ctrl->N.wI; - grn->keepAllFreq = Ctrl->N.keepAllFreq; - grn->calc_upar = Ctrl->e.active; - - grt_grnspec_allocate_u(grn); - grn->statsstr = Ctrl->S.s_statsdir; - grn->nstatsidxs = Ctrl->S.nstatsidxs; - grn->statsidxs = Ctrl->S.statsidxs; - } - + grt_grnspec_allocate_u(grn); + grn->statsstr = Ctrl->S.s_statsdir; + grn->nstatsidxs = Ctrl->S.nstatsidxs; + grn->statsidxs = Ctrl->S.statsidxs; //============================================================================== // 计算格林函数 diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index e4f10790..0b967213 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -560,17 +560,58 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } } +} + + +/** + * 静态积分前准备：默认 Length，填充 K_INTEG_PROCESS。 + * 不分配输出缓冲，不处理 stats 路径。 + */ +void grt_prepare_static_grn( + MODEL1D *mod1d, + size_t nr, real_t *rs, + real_t Length, + real_t filonLength, real_t safilonTol, real_t filonCut, + real_t k0, real_t keps, bool use_kmax_ref, + int convmet, + K_INTEG_PROCESS *Kproc) +{ // FIM/SAFIM 面向远震中距，公式含 1/r、1/√r，不能用于 r=0（含 kcut 分段） - if(Ctrl->L.FIM.active || Ctrl->L.SAFIM.active){ - for(size_t ir=0; irnr; ++ir){ - if(GRT_IS_ZERO(Ctrl->rs[ir])){ - GRTBadOptionError(L, "FIM/SAFIM cannot be used with zero epicentral distance."); + if(filonLength > 0.0 || safilonTol > 0.0){ + for(size_t ir = 0; ir < nr; ++ir){ + if(GRT_IS_ZERO(rs[ir])){ + GRTRaiseError("FIM/SAFIM cannot be used with zero epicentral distance."); } } } -} + if(Length == 0.0){ + Length = GRT_GREENFN_L_LENGTH; + } + + memset(Kproc, 0, sizeof(*Kproc)); + { + real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); + Kproc->k0 = k0 * PI / hs; + Kproc->use_kmax_ref = use_kmax_ref; + // 显式收敛方法时不使用 keps + Kproc->keps = (convmet != K_INTEG_CONVERG_AUTO) ? 0.0 : keps; + + real_t rmax = rs[grt_findMax_real_t(rs, nr)]; + Kproc->kcut = filonCut / rmax; + + // rmax=0 时用阈值防止除零，此时 dk 偏大，后续由 GRT_MIN_NK 收紧 + Kproc->dk = PI2 / (Length * GRT_MAX(rmax, GRT_ZERO_DISTANCE)); + Kproc->applyFIM = filonLength > 0.0; + Kproc->filondk = (filonLength > 0.0) ? PI2 / (filonLength * rmax) : 0.0; + + Kproc->applySAFIM = safilonTol > 0.0; + Kproc->sa_tol = safilonTol; + + Kproc->cvgmet = convmet; + } +} /** 子模块主函数 */ @@ -588,9 +629,6 @@ int static_greenfn_main(int argc, char **argv){ // 边界条件 grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); - - // 设置积分间隔默认值 - if(Ctrl->L.Length == 0.0) Ctrl->L.Length = GRT_GREENFN_L_LENGTH; // 波数积分输出目录 if(Ctrl->S.active){ @@ -607,30 +645,17 @@ int static_greenfn_main(int argc, char **argv){ realChnlGrid *grn_uiz = (Ctrl->e.active)? (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; realChnlGrid *grn_uir = (Ctrl->e.active)? (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; - // 波数积分方法 K_INTEG_PROCESS KPROC = {0}; - { - real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); - KPROC.k0 = Ctrl->K.k0 * PI / hs; - KPROC.use_kmax_ref = Ctrl->K.use_kmax_ref; - 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)]; - - KPROC.kcut = Ctrl->L.kcut / 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; - - KPROC.applySAFIM = Ctrl->L.SAFIM.active; - KPROC.sa_tol = Ctrl->L.SAFIM.tol; - - KPROC.cvgmet = Ctrl->C.convmet; - } + grt_prepare_static_grn( + mod1d, + Ctrl->nr, Ctrl->rs, + Ctrl->L.Length, + Ctrl->L.FIM.active ? Ctrl->L.FIM.Length : 0.0, + Ctrl->L.SAFIM.active ? Ctrl->L.SAFIM.tol : 0.0, + Ctrl->L.kcut, + Ctrl->K.k0, Ctrl->K.keps, Ctrl->K.use_kmax_ref, + Ctrl->C.convmet, + &KPROC); //============================================================================== // 计算静态格林函数 diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 87e4d06c..e85aed33 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -26,11 +26,42 @@ """libgrt库""" +C_grt_prepare_grn_spec = libgrt.grt_prepare_grn_spec +"""动态格林函数：频谱积分前准备（Length/wI/freqs/KPROC/GRNSPEC 元数据）""" +C_grt_prepare_grn_spec.restype = None +C_grt_prepare_grn_spec.argtypes = [ + POINTER(c_MODEL1D), + c_size_t, PREAL, + c_size_t, REAL, REAL, c_bool, + REAL, REAL, + REAL, + REAL, REAL, REAL, + REAL, REAL, REAL, REAL, c_bool, + c_int, + REAL, REAL, c_bool, + c_bool, + POINTER(c_K_INTEG_PROCESS), + POINTER(c_GRNSPEC), +] + C_grt_integ_grn_spec = libgrt.grt_integ_grn_spec """C库中计算格林函数的主函数 integ_grn_spec, 详见C API同名函数""" C_grt_integ_grn_spec.argtypes = [POINTER(c_MODEL1D), POINTER(c_K_INTEG_PROCESS), POINTER(c_GRNSPEC), c_bool] +C_grt_prepare_static_grn = libgrt.grt_prepare_static_grn +"""静态格林函数：积分前准备（Length 默认 + KPROC）""" +C_grt_prepare_static_grn.restype = None +C_grt_prepare_static_grn.argtypes = [ + POINTER(c_MODEL1D), + c_size_t, PREAL, + REAL, + REAL, REAL, REAL, + REAL, REAL, c_bool, + c_int, + POINTER(c_K_INTEG_PROCESS), +] + C_grt_integ_static_grn = libgrt.grt_integ_static_grn """计算静态格林函数""" C_grt_integ_static_grn.restype = None diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index e1129ce6..23151c26 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -23,8 +23,6 @@ "ZNEchs", "qwvchs", "MECHANISM_NUM", - "MIN_DEPTH_GAP_SRC_RCV", - "ZERO_DISTANCE", "NPCT_REAL_TYPE", "NPCT_CMPLX_TYPE", @@ -52,8 +50,6 @@ ZNEchs = ['Z', 'N', 'E'] qwvchs = ['q', 'w', 'v'] MECHANISM_NUM = 6 -MIN_DEPTH_GAP_SRC_RCV = 0.1 -ZERO_DISTANCE = 1e-8 # 与 C 侧 GRT_ZERO_DISTANCE 一致 NPCT_REAL_TYPE = 'f8' NPCT_CMPLX_TYPE = 'c16' diff --git a/pygrt/pymod.py b/pygrt/pymod.py index ba35b5ee..f491cff6 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -17,11 +17,12 @@ from obspy.core import AttribDict from typing import List, Dict, Union, Literal import tempfile +import os from time import time from copy import deepcopy -from ctypes import Array, pointer +from ctypes import Array, pointer, c_char_p from ctypes import _Pointer from .c_interfaces import * from .c_structures import * @@ -190,10 +191,7 @@ def _get_grn_spectra( statsidxs:Union[np.ndarray,List[int],None]=None, print_log:bool=True ): - - depsrc = self.depsrc - deprcv = self.deprcv - + # 仅做最基本的正负号等检查；物理预处理交给 C grt_prepare_grn_spec if np.any(distarr < 0): raise ValueError(f"distarr < 0") if nt < 0: @@ -214,162 +212,92 @@ def _get_grn_spectra( if filonCut < 0.0: raise ValueError(f"filonCut ({filonCut}) < 0") if safilonTol < 0.0: - raise ValueError(f"filonCut ({safilonTol}) < 0") + raise ValueError(f"safilonTol ({safilonTol}) < 0") # 只能设置一种filon积分方法 if safilonTol > 0.0 and filonLength > 0.0: raise ValueError(f"You should only set one of filonLength and safilonTol.") - # FIM/SAFIM 面向远震中距，不能用于 r=0（含 filonCut 分段） - if (filonLength > 0.0 or safilonTol > 0.0) and np.any(np.asarray(distarr) <= ZERO_DISTANCE): - raise ValueError("FIM/SAFIM cannot be used with zero epicentral distance.") - - nf = nt//2+1 - df = 1/(nt*dt) - fnyq = 1/(2*dt) - # 确定频带范围 f1, f2 = freqband if f1 >= f2 and f1 >= 0 and f2 >= 0: raise ValueError(f"freqband f1({f1}) >= f2({f2})") - - if f1 < 0: - f1 = 0 - if f2 < 0: - f2 = fnyq+df - - f1 = max(0, f1) - f2 = min(f2, fnyq + df) - nf1 = min(int(np.ceil(f1/df)), nf-1) - nf2 = min(int(np.floor(f2/df)), nf-1) - if nf2 < nf1: - nf2 = nf1 - - # 所有频点 - freqs = (np.arange(0, nf)*df).astype(NPCT_REAL_TYPE) - - # 虚频率 - wI = zeta * np.pi/(nt*dt) + distarr = np.asarray(distarr, dtype=NPCT_REAL_TYPE) nrs = len(distarr) - for ir in range(nrs): - if(distarr[ir] < 0.0): - raise ValueError(f"r({distarr[ir]}) < 0") - - # 最大震中距 - rmax = np.max(distarr) - - # 转为C类型 - c_freqs = npct.as_ctypes(freqs) - c_rs = npct.as_ctypes(np.array(distarr).astype(NPCT_REAL_TYPE) ) + c_rs = npct.as_ctypes(distarr) - # 参考最小速度 - if vmin_ref == 0.0: - vmin_ref = max(self.vmin, 0.1) - - # 时窗长度 - winT = nt*dt - - # 时窗最大截止时刻 - tmax = delayT0 + winT - if delayV0 > 0.0: - 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 and rmax > ZERO_DISTANCE: - Length = 1.0 + np.sqrt(jus)/rmax + 0.5 # 0.5作保守值 - if Length < 15.0: - Length = 15.0 - - # 初始化格林函数 - pygrnLst, c_grnArr = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, '') - - pygrnLst_uiz = [] - c_grnArr_uiz = None - pygrnLst_uir = [] - c_grnArr_uir = None - if calc_upar: - pygrnLst_uiz, c_grnArr_uiz = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, 'z') - pygrnLst_uir, c_grnArr_uir = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, 'r') - - - c_statsfile = None - if statsfile is not None: - os.makedirs(statsfile, exist_ok=True) - c_statsfile = c_char_p(statsfile.encode('utf-8')) - - nstatsidxs = 0 - if statsidxs is None: - statsidxs = np.arange(nf) - - statsidxs = np.array(statsidxs) - # 不能有负数 - if np.any(statsidxs < 0): - raise ValueError("negative value in statsidxs is not supported.") - - c_statsidxs = npct.as_ctypes(np.array(statsidxs).astype(np.uint64)) # size_t - nstatsidxs = len(statsidxs) - else: - c_statsfile = c_statsidxs = None - nstatsidxs = 0 - - - # ==================================================================== KPROC = c_K_INTEG_PROCESS() - hs = max(abs(depsrc - deprcv), MIN_DEPTH_GAP_SRC_RCV) - KPROC.k0 = k0 * np.pi / hs - 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 - - KPROC.kcut = filonCut / 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 - - KPROC.applySAFIM = safilonTol > 0.0 - KPROC.sa_tol = safilonTol - - KPROC.cvgmet = K_INTEG_CVGMET_DICT[converg_method.upper()] - # ==================================================================== - - - # ==================================================================== grn = c_GRNSPEC() - grn.nf = nf - grn.freqs = c_freqs - grn.nf1 = nf1 - grn.nf2 = nf2 - grn.nr = nrs - grn.rs = c_rs - grn.wI = wI - grn.keepAllFreq = keepAllFreq - grn.calc_upar = calc_upar - grn.u = c_grnArr - grn.uiz = c_grnArr_uiz - grn.uir = c_grnArr_uir - grn.statsstr = c_statsfile - grn.nstatsidxs = nstatsidxs - grn.statsidxs = c_statsidxs - # ==================================================================== + C_grt_prepare_grn_spec( + self.c_mod1d, + nrs, c_rs, + nt, dt, zeta, keepAllFreq, + float(f1), float(f2), + Length, + filonLength, safilonTol, filonCut, + k0, ampk, keps, vmin_ref, use_kmax_ref, + K_INTEG_CVGMET_DICT[converg_method.upper()], + delayT0, delayV0, False, + calc_upar, + pointer(KPROC), pointer(grn), + ) - # 运行C库函数 - #///////////////////////////////////////////////////////////////////////////////// - # 计算得到的格林函数的单位： - # 单力源 HF[ZRT],VF[ZR] 1e-15 cm/dyne - # 爆炸源 EX[ZR] 1e-20 cm/(dyne*cm) - # 剪切源 DD[ZR],DS[ZRT],SS[ZRT] 1e-20 cm/(dyne*cm) - #================================================================================= - C_grt_integ_grn_spec(self.c_mod1d, pointer(KPROC), pointer(grn), print_log) - #================================================================================= - #///////////////////////////////////////////////////////////////////////////////// + try: + nf = grn.nf + freqs = npct.as_array(grn.freqs, shape=(nf,)).copy() + freqs.flags.writeable = False + wI = float(grn.wI) + + # 初始化格林函数（缓冲仍由 Python 持有） + pygrnLst, c_grnArr = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, '') + + pygrnLst_uiz = [] + c_grnArr_uiz = None + pygrnLst_uir = [] + c_grnArr_uir = None + if calc_upar: + pygrnLst_uiz, c_grnArr_uiz = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, 'z') + pygrnLst_uir, c_grnArr_uir = self._init_grn(distarr, nt, dt, upsampling_n, freqs, wI, 'r') + + c_statsfile = None + if statsfile is not None: + os.makedirs(statsfile, exist_ok=True) + c_statsfile = c_char_p(statsfile.encode('utf-8')) + + if statsidxs is None: + statsidxs = np.arange(nf) + + statsidxs = np.array(statsidxs) + if np.any(statsidxs < 0): + raise ValueError("negative value in statsidxs is not supported.") + + c_statsidxs = npct.as_ctypes(np.array(statsidxs).astype(np.uint64)) # size_t + nstatsidxs = len(statsidxs) + else: + c_statsidxs = None + nstatsidxs = 0 + + grn.u = c_grnArr + grn.uiz = c_grnArr_uiz + grn.uir = c_grnArr_uir + grn.statsstr = c_statsfile + grn.nstatsidxs = nstatsidxs + grn.statsidxs = c_statsidxs + + # 运行C库函数 + #///////////////////////////////////////////////////////////////////////////////// + # 计算得到的格林函数的单位： + # 单力源 HF[ZRT],VF[ZR] 1e-15 cm/dyne + # 爆炸源 EX[ZR] 1e-20 cm/(dyne*cm) + # 剪切源 DD[ZR],DS[ZRT],SS[ZRT] 1e-20 cm/(dyne*cm) + #================================================================================= + C_grt_integ_grn_spec(self.c_mod1d, pointer(KPROC), pointer(grn), print_log) + #================================================================================= + #///////////////////////////////////////////////////////////////////////////////// + finally: + if grn.freqs: + C_grt_free(grn.freqs) + grn.freqs = None return pygrnLst, pygrnLst_uiz, pygrnLst_uir @@ -610,9 +538,6 @@ def compute_static_grn( if safilonTol > 0.0 and filonLength > 0.0: raise ValueError(f"You should only set one of filonLength and safilonTol.") - depsrc = self.depsrc - deprcv = self.deprcv - if distarr is not None: if isinstance(distarr, float) or isinstance(distarr, int): distarr = np.array([distarr*1.0]) @@ -643,14 +568,6 @@ def compute_static_grn( rs[ix + iy*nx] = np.hypot(xarr[ix], yarr[iy]) c_rs = npct.as_ctypes(rs) - # FIM/SAFIM 面向远震中距，不能用于 r=0（含 filonCut 分段） - if (filonLength > 0.0 or safilonTol > 0.0) and np.any(rs <= ZERO_DISTANCE): - raise ValueError("FIM/SAFIM cannot be used with zero epicentral distance.") - - # 设置波数积分间隔 - if Length == 0.0: - Length = 15.0 - # 积分状态文件 c_statsfile = None if statsfile is not None: @@ -664,31 +581,18 @@ def compute_static_grn( if not calc_upar: c_pygrn_uiz = c_pygrn_uir = None - - # ==================================================================== + # 仅做最基本的正负号等检查；物理预处理交给 C grt_prepare_static_grn KPROC = c_K_INTEG_PROCESS() - hs = max(abs(depsrc - deprcv), MIN_DEPTH_GAP_SRC_RCV) - KPROC.k0 = k0 * np.pi / hs - KPROC.use_kmax_ref = use_kmax_ref - KPROC.keps = keps if converg_method.upper() != 'AUTO' else 0.0 - - # 最大震中距 - rmax = np.max(rs) - KPROC.kcut = filonCut / 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 - - KPROC.applySAFIM = safilonTol > 0.0 - KPROC.sa_tol = safilonTol - - KPROC.cvgmet = K_INTEG_CVGMET_DICT[converg_method.upper()] - # ==================================================================== - - + C_grt_prepare_static_grn( + self.c_mod1d, + nr, c_rs, + Length, + filonLength, safilonTol, filonCut, + k0, keps, use_kmax_ref, + K_INTEG_CVGMET_DICT[converg_method.upper()], + pointer(KPROC), + ) # 运行C库函数 #/////////////////////////////////////////////////////////////////////////////////