diff --git a/pygrt/C_extension/include/grt.h b/pygrt/C_extension/include/grt.h index c21ee300..47467f94 100644 --- a/pygrt/C_extension/include/grt.h +++ b/pygrt/C_extension/include/grt.h @@ -43,7 +43,6 @@ #include "grt/dynamic/grn.h" #include "grt/dynamic/grnspec.h" #include "grt/dynamic/layer.h" -#include "grt/dynamic/dyn_postprocess.h" #include "grt/dynamic/signals.h" #include "grt/dynamic/source.h" @@ -59,7 +58,6 @@ #include "grt/static/static_grn.h" #include "grt/static/static_layer.h" -#include "grt/static/static_postprocess.h" #include "grt/static/static_source.h" diff --git a/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h b/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h deleted file mode 100644 index 459d66e5..00000000 --- a/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h +++ /dev/null @@ -1,38 +0,0 @@ -/** - * @file dyn_postprocess.h - * @brief 动态位移偏导后处理(应力/应变/旋转) - * @date 2026-07 - */ - -#pragma once - -#include "grt/common/const.h" - -/** - * 由动态位移偏导合成应变张量。 - * 数组布局:u[分量][采样点]、upar[偏导方向][分量][采样点]、 - * res[第二分量][第一分量][采样点]。 - * - * ZRT 联络项:r≠0 用 u/r;r=0 改用 ∂_r u(upar[1][*]), - * 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。 - */ -void grt_compute_strain( - size_t npts, float dist, float *const u[GRT_CHANNEL_NUM], - float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE); - -/** 由动态位移偏导合成旋转张量。数组布局同 grt_compute_strain()。 */ -void grt_compute_rotation( - size_t npts, float dist, float *const u[GRT_CHANNEL_NUM], - float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE); - -/** - * 在频域由动态位移偏导合成应力张量。 - * 介质参数约定与命令行 stress 子模块的 SAC 头段一致。 - */ -void grt_compute_stress( - size_t npts, float dt, float dist, float va, float vb, float rho, - float Qainv, float Qbinv, float *const u[GRT_CHANNEL_NUM], - float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE); diff --git a/pygrt/C_extension/include/grt/dynamic/grn.h b/pygrt/C_extension/include/grt/dynamic/grn.h index 629c3150..dd4cbed9 100755 --- a/pygrt/C_extension/include/grt/dynamic/grn.h +++ b/pygrt/C_extension/include/grt/dynamic/grn.h @@ -19,7 +19,7 @@ #include "grt/integral/integ_process.h" /** - * 积分计算Z, R, T三个分量格林函数的频谱的核心函数(被Python调用) + * 积分计算Z, R, T三个分量格林函数的频谱的核心函数 * * @param[in,out] mod1d `MODEL1D` 结构体指针 * @param[in,out] Kproc 波数积分相关参数的结构体指针 diff --git a/pygrt/C_extension/include/grt/dynamic/signals.h b/pygrt/C_extension/include/grt/dynamic/signals.h index edf4bb71..3d7ecd4c 100644 --- a/pygrt/C_extension/include/grt/dynamic/signals.h +++ b/pygrt/C_extension/include/grt/dynamic/signals.h @@ -173,7 +173,7 @@ float * grt_get_ricker_wave(float dt, float f0, int *Nt); float * grt_get_custom_wave(int *Nt, const char *tfparams); /** - * 专用于在Python端释放C中申请的内存 + * 释放 C 侧 malloc 的一维数组 * * @param[out] pt 指针 */ diff --git a/pygrt/C_extension/include/grt/static/static_postprocess.h b/pygrt/C_extension/include/grt/static/static_postprocess.h deleted file mode 100644 index 68a2ed46..00000000 --- a/pygrt/C_extension/include/grt/static/static_postprocess.h +++ /dev/null @@ -1,45 +0,0 @@ -/** - * @file static_postprocess.h - * @brief 静态位移偏导后处理(应力/应变/旋转) - * @date 2026-07 - */ - -#pragma once - -#include "grt/common/const.h" - -/** - * 由静态位移偏导合成对称应力张量。 - * - * 数组布局:u[分量][点]、upar[偏导方向][分量][点]、 - * res[第二分量][第一分量][点]。仅写入 res 的上三角分量。 - * - * ZRT 联络项:r≠0 用 u/r;r=0 改用 ∂_r u(upar[1][*]), - * 与 syn 中轴点处 (1/r)∂_θ 有限部分配套。 - */ -void grt_static_compute_stress( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, - real_t *const u[GRT_CHANNEL_NUM], - real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - bool rot2ZNE, real_t mu, real_t lam); - -/** - * 由静态位移偏导合成对称应变张量。 - * 数组布局同 grt_static_compute_stress()。 - */ -void grt_static_compute_strain( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, - real_t *const u[GRT_CHANNEL_NUM], - real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE); - -/** - * 由静态位移偏导合成反对称旋转张量。 - * 数组布局同 grt_static_compute_stress()。仅写入 res 的非对角分量。 - */ -void grt_static_compute_rotation( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, - real_t *const u[GRT_CHANNEL_NUM], - real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE); diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 11166664..adf04d58 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -53,13 +53,8 @@ typedef struct { struct { bool active; size_t nt; - size_t nf; real_t dt; - real_t df; - real_t winT; ///< 时窗长度 - real_t zeta; ///< 虚频率系数, w <- w - zeta*PI/r* 1j - real_t wI; ///< 虚频率 zeta*PI/r - real_t *freqs; + real_t zeta; ///< 虚频率系数, w <- w - zeta*PI/T * 1j size_t upsample_n; ///< 升采样倍数 bool keepAllFreq; ///< 计算所有频率,不论频率多低 bool skipImagComps; ///< 跳过虚频率的补偿 @@ -74,8 +69,6 @@ typedef struct { bool active; real_t freq1; real_t freq2; - size_t nf1; - size_t nf2; } H; /** 波数积分间隔以及方法 */ struct { @@ -163,9 +156,6 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ GRT_SAFE_FREE_PTR(Ctrl->D.s_depsrc); GRT_SAFE_FREE_PTR(Ctrl->D.s_deprcv); - // N - GRT_SAFE_FREE_PTR(Ctrl->N.freqs); - // O GRT_SAFE_FREE_PTR(Ctrl->O.s_output_dir); @@ -851,11 +841,11 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** - * 频谱积分前的共享准备:自动 Length、wI、频点网格、nf1/nf2、填充 KPROC 与 GRNSPEC 元数据。 - * grn->freqs 由本函数 malloc,调用方负责释放;grn->rs 借用调用方指针。 - * 不分配 u/uiz/uir,不处理 stats 路径。 + * 频谱积分前准备:自动 Length、wI、频点网格、nf1/nf2、填充 KPROC 与 GRNSPEC 元数据 + * grn->freqs 由本函数 malloc,调用方负责释放;grn->rs 借用调用方指针 + * 不分配 u/uiz/uir,不处理 stats 路径 */ -void grt_prepare_grn_spec( +static void prepare_grn_spec( MODEL1D *mod1d, size_t nr, real_t *rs, size_t nt, real_t dt, real_t zeta, bool keepAllFreq, @@ -989,7 +979,7 @@ int greenfn_main(int argc, char **argv) { K_INTEG_PROCESS KPROC = {0}; GRNSPEC grn_storage = {0}; GRNSPEC *grn = &grn_storage; - grt_prepare_grn_spec( + prepare_grn_spec( mod1d, Ctrl->R.nr, Ctrl->R.rs, Ctrl->N.nt, Ctrl->N.dt, Ctrl->N.zeta, Ctrl->N.keepAllFreq, @@ -1004,19 +994,12 @@ int greenfn_main(int argc, char **argv) { 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; + real_t df = 1.0 / (Ctrl->N.nt * Ctrl->N.dt); // 如果只传入了 -S, 未指定索引,则默认所有频率索引 if(Ctrl->S.active && Ctrl->S.statsidxs == NULL){ // 另外两个字符相关的指针仍指向 NULL - Ctrl->S.nstatsidxs = Ctrl->N.nf; + Ctrl->S.nstatsidxs = grn->nf; Ctrl->S.statsidxs = (size_t*)realloc(Ctrl->S.statsidxs, sizeof(size_t)*(Ctrl->S.nstatsidxs)); for(size_t i=0; i < Ctrl->S.nstatsidxs; ++i){ Ctrl->S.statsidxs[i] = i; @@ -1047,7 +1030,7 @@ int greenfn_main(int argc, char **argv) { // 使用fftw3做反傅里叶变换,并保存到 SAC // 其中考虑了升采样倍数 - GRT_FFTW_HOLDER *fh = grt_create_fftw_holder_C2R_1D(Ctrl->N.nt*Ctrl->N.upsample_n, Ctrl->N.dt/Ctrl->N.upsample_n, Ctrl->N.nf, Ctrl->N.df); + GRT_FFTW_HOLDER *fh = grt_create_fftw_holder_C2R_1D(Ctrl->N.nt*Ctrl->N.upsample_n, Ctrl->N.dt/Ctrl->N.upsample_n, grn->nf, df); // 建立 SAC 文件原型,包含必要的头变量 SACTRACE *sac = grt_new_SACTRACE(fh->dt, fh->nt, 0.0); @@ -1058,7 +1041,7 @@ int greenfn_main(int argc, char **argv) { sac->hd.evdp = Ctrl->D.depsrc; // km sac->hd.stel = (-1.0)*Ctrl->D.deprcv*1e3; // m // 写入虚频率 - sac->hd.user0 = Ctrl->N.wI; + sac->hd.user0 = grn->wI; // 写入接受点的Vp,Vs,rho sac->hd.user1 = mod1d->Va[mod1d->ircv]; sac->hd.user2 = mod1d->Vb[mod1d->ircv]; @@ -1112,6 +1095,7 @@ int greenfn_main(int argc, char **argv) { // 释放内存 grt_grnspec_free_u(grn); + GRT_SAFE_FREE_PTR(grn->freqs); GRT_SAFE_FREE_PTR(travtPS); GRT_SAFE_FREE_PTR(begintimes); grt_free_SACTRACE(sac); diff --git a/pygrt/C_extension/src/dynamic/grt_rotation.c b/pygrt/C_extension/src/dynamic/grt_rotation.c index d983e579..c4186408 100644 --- a/pygrt/C_extension/src/dynamic/grt_rotation.c +++ b/pygrt/C_extension/src/dynamic/grt_rotation.c @@ -51,7 +51,8 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ GRTCheckOptionSet(argc > 1); } -void grt_compute_rotation( +/** 由位移偏导合成旋转张量 */ +static void compute_rotation( size_t npts, float dist, float *const u[GRT_CHANNEL_NUM], float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE) @@ -131,7 +132,7 @@ int rotation_main(int argc, char **argv){ res[c2][c] = calloc(npts, sizeof(*res[c2][c])); } } - grt_compute_rotation(npts, dist, u, upar, res, rot2ZNE); + compute_rotation(npts, dist, u, upar, res, rot2ZNE); // 写出3个分量 for(int i1=0; i1<2; ++i1){ diff --git a/pygrt/C_extension/src/dynamic/grt_strain.c b/pygrt/C_extension/src/dynamic/grt_strain.c index 2d098c85..11074c4a 100644 --- a/pygrt/C_extension/src/dynamic/grt_strain.c +++ b/pygrt/C_extension/src/dynamic/grt_strain.c @@ -49,7 +49,8 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ GRTCheckOptionSet(argc > 1); } -void grt_compute_strain( +/** 由位移偏导合成应变张量 */ +static void compute_strain( size_t npts, float dist, float *const u[GRT_CHANNEL_NUM], float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE) @@ -133,7 +134,7 @@ int strain_main(int argc, char **argv){ res[c2][c] = calloc(npts, sizeof(*res[c2][c])); } } - grt_compute_strain(npts, dist, u, upar, res, rot2ZNE); + compute_strain(npts, dist, u, upar, res, rot2ZNE); // 写出6个分量 for(int i1=0; i1<3; ++i1){ diff --git a/pygrt/C_extension/src/dynamic/grt_stress.c b/pygrt/C_extension/src/dynamic/grt_stress.c index e1a40388..63c4ea05 100644 --- a/pygrt/C_extension/src/dynamic/grt_stress.c +++ b/pygrt/C_extension/src/dynamic/grt_stress.c @@ -51,7 +51,8 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ GRTCheckOptionSet(argc > 1); } -void grt_compute_stress( +/** 在频域由位移偏导合成应力张量 */ +static void compute_stress( size_t npts, float dt, float dist, float va, float vb, float rho, float Qainv, float Qbinv, float *const u[GRT_CHANNEL_NUM], float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], @@ -201,7 +202,7 @@ int stress_main(int argc, char **argv){ res[c2][c] = calloc(npts, sizeof(*res[c2][c])); } } - grt_compute_stress(npts, dt, dist, va, vb, rho, Qainv, Qbinv, u, upar, res, rot2ZNE); + compute_stress(npts, dt, dist, va, vb, rho, Qainv, Qbinv, u, upar, res, rot2ZNE); // 写出6个分量 for(int i1=0; i1<3; ++i1){ diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index a705e8b6..bcf146c0 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -531,18 +531,18 @@ static void syn_accum_from_gf( /** - * 由动态格林函数合成三分量地震图(及可选空间偏导)。 + * 由动态格林函数合成三分量地震图(及可选空间偏导) * * 数组布局:gf[震源][分量][采样点]、syn[分量][采样点]、 * syn_upar[偏导方向][分量][采样点]。gf_uiz/gf_uir 在 calc_upar=false - * 时可传 NULL;单个分量指针为 NULL 时跳过该道。 + * 时可传 NULL;单个分量指针为 NULL 时跳过该道 * * r=0 时强制 *azrad=0(e_r→N、e_θ→E)并告警; - * 并用 ∂_r 格林函数合成 (1/r)∂_θ 的有限部分(见函数内注释)。 + * 并用 ∂_r 格林函数合成 (1/r)∂_θ 的有限部分(见函数内注释) * * @param[in,out] azrad 方位角(弧度);r=0 时写回 0 */ -void grt_syn_from_gf( +static void syn_from_gf( size_t npts, float dist, const pfloatChnlGrid gf, const pfloatChnlGrid gf_uiz, const pfloatChnlGrid gf_uir, GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, real_t *azrad, @@ -768,7 +768,7 @@ int syn_main(int argc, char **argv) } } - grt_syn_from_gf( + syn_from_gf( (size_t)npts, Ctrl->dist, gf, calc_upar ? gf_uiz : NULL, calc_upar ? gf_uir : NULL, Ctrl->computeType, Ctrl->S.M0, Ctrl->VpVs_ratio, &Ctrl->A.azrad, Ctrl->mchn, diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index a8d1659f..a2d5aa26 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -567,10 +567,10 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** - * 静态积分前准备:默认 Length,填充 K_INTEG_PROCESS。 - * 不分配输出缓冲,不处理 stats 路径。 + * 静态积分前准备:默认 Length,填充 K_INTEG_PROCESS + * 不分配输出缓冲,不处理 stats 路径 */ -void grt_prepare_static_grn( +static void prepare_static_grn( MODEL1D *mod1d, size_t nr, real_t *rs, real_t Length, @@ -649,7 +649,7 @@ int static_greenfn_main(int argc, char **argv){ realChnlGrid *grn_uir = (Ctrl->e.active)? (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; K_INTEG_PROCESS KPROC = {0}; - grt_prepare_static_grn( + prepare_static_grn( mod1d, Ctrl->nr, Ctrl->rs, Ctrl->L.Length, diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 0485753a..50461383 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -49,7 +49,8 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 暂不支持设置其它参数 } -void grt_static_compute_rotation( +/** 由静态位移偏导合成旋转张量 */ +static void compute_rotation( size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, real_t *const u[GRT_CHANNEL_NUM], real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], @@ -186,7 +187,7 @@ int static_rotation_main(int argc, char **argv){ } } - grt_static_compute_rotation(nnorth, neast, norths, easts, u, upar, res, rot2ZNE); + compute_rotation(nnorth, neast, norths, easts, u, upar, res, rot2ZNE); // 写入 nc 文件 for(int c=0; c