From 6deaf5f69f2d6283b6f4d2321d60567bb747ff9f Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 7 Apr 2026 18:05:55 +0800 Subject: [PATCH] REFAC: add "grnspec.h/c" to handle the Green Functions spectra --- pygrt/C_extension/include/grt.h | 1 + pygrt/C_extension/include/grt/dynamic/grn.h | 36 +--- .../C_extension/include/grt/dynamic/grnspec.h | 62 ++++++ pygrt/C_extension/src/dynamic/grn.c | 69 +++---- pygrt/C_extension/src/dynamic/grnspec.c | 148 ++++++++++++++ pygrt/C_extension/src/dynamic/grt_greenfn.c | 185 +++++------------- pygrt/c_interfaces.py | 13 +- pygrt/c_structures.py | 28 ++- pygrt/pymod.py | 31 ++- 9 files changed, 341 insertions(+), 232 deletions(-) create mode 100644 pygrt/C_extension/include/grt/dynamic/grnspec.h create mode 100644 pygrt/C_extension/src/dynamic/grnspec.c diff --git a/pygrt/C_extension/include/grt.h b/pygrt/C_extension/include/grt.h index b7b5920..ee1acc2 100644 --- a/pygrt/C_extension/include/grt.h +++ b/pygrt/C_extension/include/grt.h @@ -41,6 +41,7 @@ #include "grt/dynamic/grn.h" +#include "grt/dynamic/grnspec.h" #include "grt/dynamic/layer.h" #include "grt/dynamic/signals.h" #include "grt/dynamic/source.h" diff --git a/pygrt/C_extension/include/grt/dynamic/grn.h b/pygrt/C_extension/include/grt/dynamic/grn.h index f6ded88..d7664e4 100755 --- a/pygrt/C_extension/include/grt/dynamic/grn.h +++ b/pygrt/C_extension/include/grt/dynamic/grn.h @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * 以下代码实现的是 广义反射透射系数矩阵+离散波数法 计算理论地震图，参考： + * 以下代码实现的是利用多种波数积分类方法计算理论地震图，参考： * * 1. 姚振兴, 谢小碧. 2022/03. 理论地震图及其应用（初稿）. * 2. Yao Z. X. and D. G. Harkrider. 1983. A generalized refelection-transmission coefficient @@ -15,43 +15,19 @@ #pragma once #include "grt/common/model.h" +#include "grt/dynamic/grnspec.h" #include "grt/integral/integ_method.h" /** * 积分计算Z, R, T三个分量格林函数的频谱的核心函数（被Python调用） * - * @param[in,out] mod1d `MODEL1D` 结构体指针 - * @param[in] nf1 开始计算频谱的频率索引值, 总范围在[nf1, nf2] - * @param[in] nf2 结束计算频谱的频率索引值, 总范围在[nf1, nf2] - * @param[in] freqs 所有频点的频率值（包括未计算的） - * @param[in] nr 震中距数量 - * @param[in] rs 震中距数组 - * @param[in] wI 虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$ - * @param[in] keepAllFreq 计算所有频点，不论频率多低 - * @param[in,out] Kmet 波数积分相关参数的结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in,out] Kmet 波数积分相关参数的结构体指针 + * @param[in,out] grn 计算的格林函数频谱 * @param[in] print_progressbar 是否打印进度条 - * @param[in] calc_upar 是否计算位移u的空间导数 - * @param[out] grn 不同震源不同阶数的格林函数的Z、R、T分量频谱结果 - * @param[out] grn_uiz 不同震源不同阶数的ui_z的Z、R、T分量频谱结果 - * @param[out] grn_uir 不同震源不同阶数的ui_r的Z、R、T分量频谱结果 - * - * @param[in] statsstr 积分过程输出目录 - * @param[in] nstatsidxs 输出积分过程的特定频点的个数 - * @param[in] statsidxs 特定频点的索引值 * */ -void grt_integ_grn_spec( - MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs, - size_t nr, real_t *rs, real_t wI, bool keepAllFreq, K_INTEG_METHOD *Kmet, - bool print_progressbar, bool calc_upar, - pcplxChnlGrid grn[nr], - pcplxChnlGrid grn_uiz[nr], - pcplxChnlGrid grn_uir[nr], - - const char *statsstr, // 积分过程输出 - size_t nstatsidxs, // 仅输出特定频点 - size_t *statsidxs -); +void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_METHOD *Kmet, GRNSPEC *grn, const bool print_progressbar); diff --git a/pygrt/C_extension/include/grt/dynamic/grnspec.h b/pygrt/C_extension/include/grt/dynamic/grnspec.h new file mode 100644 index 0000000..a7002ab --- /dev/null +++ b/pygrt/C_extension/include/grt/dynamic/grnspec.h @@ -0,0 +1,62 @@ +/** + * @file grnspec.h + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-04 + * + * 将格林函数频谱用一个结构体包裹 + * + * + */ + +#pragma once + +#include "grt/common/const.h" +#include "grt/common/myfftw.h" +#include "grt/common/sacio.h" + +/** 不同震中距、不同震源、不同分量的格林函数频谱 */ +typedef struct { + size_t nf; ///< 总频点数 + real_t *freqs; ///< 频率数组 freqs[nf] + size_t nf1; + size_t nf2; ///< 待计算的频率索引 + size_t nr; ///< 震中距数量 + real_t *rs; ///< 震中距数组 + real_t wI; ///< 虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$ + bool keepAllFreq; ///< 是否计算所有频点，不论频率多低 + bool calc_upar; ///< 是否计算位移u的空间导数 + + pcplxChnlGrid *u; ///< 不同震源不同阶数的格林函数的Z、R、T分量频谱结果 + pcplxChnlGrid *uiz; ///< 不同震源不同阶数的格林函数的Z、R、T分量对z偏导的频谱结果 + pcplxChnlGrid *uir; ///< 不同震源不同阶数的格林函数的Z、R、T分量对r偏导的频谱结果 + + char *statsstr; ///< 积分结果输出路径 + size_t nstatsidxs; ///< 仅输出特定频点的对应数量 + size_t *statsidxs; ///< 对应频点的频率索引 statsidxs[nstatsidxs] +} GRNSPEC; + +/** 申请 u, uiz, uir 的内存 */ +void grt_grnspec_allocate_u(GRNSPEC *grn); + +/** 释放 u, uiz, uir 的内存 */ +void grt_grnspec_free_u(GRNSPEC *grn); + +/** + * 将频谱 u, uiz, uir 变换到时域后以 SAC 格式保存到本地 + * + * @param[in] grn 格林函数频谱结构体 + * @param[in] travtPS 不同震源距的初至P、S到时 + * @param[in] begintimes 不同震中距的波形时移 + * @param[in] outputdirs 不同震中距的保存目录 + * @param[in] fh 控制反傅里叶变换的结构体 + * @param[in,out] sac SACTRACE 原型，在头段变量中记录了基本信息 + * @param[in] skipImagComps 跳过虚频率的补偿 + * @param[in] saveEX 保存爆炸源结果 + * @param[in] saveVF 保存垂直力源结果 + * @param[in] saveHF 保存水平力源结果 + * @param[in] saveDC 保存水平力源结果 + * + */ +void grt_grnspec_write_sac( + const GRNSPEC *grn, const real_t (*travtPS)[2], const real_t *begintimes, char **outputdirs, GRT_FFTW_HOLDER *fh, SACTRACE *sac, + const bool skipImagComps, const bool saveEX, const bool saveVF, const bool saveHF, const bool saveDC); diff --git a/pygrt/C_extension/src/dynamic/grn.c b/pygrt/C_extension/src/dynamic/grn.c index f98eaea..f3ba031 100755 --- a/pygrt/C_extension/src/dynamic/grn.c +++ b/pygrt/C_extension/src/dynamic/grn.c @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-07-24 * - * 以下代码实现的是 广义反射透射系数矩阵+离散波数法 计算理论地震图，参考： + * 以下代码实现的是利用多种波数积分类方法计算理论地震图，参考： * * 1. 姚振兴, 谢小碧. 2022/03. 理论地震图及其应用（初稿）. * 2. Yao Z. X. and D. G. Harkrider. 1983. A generalized refelection-transmission coefficient @@ -20,6 +20,7 @@ #include #include "grt/dynamic/grn.h" +#include "grt/dynamic/grnspec.h" #include "grt/integral/integ_method.h" #include "grt/common/const.h" @@ -29,48 +30,36 @@ /** - * 将计算好的复数形式的积分结果记录到GRN结构体中 + * 将计算好的复数形式的积分结果转化为频谱 * * @param[in] iw 当前频率索引值 * @param[in] nr 震中距个数 * @param[in] coef 统一系数 * @param[in] sumJ 积分结果 - * @param[out] grn 三分量频谱 + * @param[out] u 三分量频谱 */ -static void recordin_GRN( - size_t iw, size_t nr, cplx_t coef, cplxIntegGrid sumJ[nr], - pcplxChnlGrid grn[nr]) +static void recordin_GRN(size_t iw, size_t nr, cplx_t coef, cplxIntegGrid sumJ[nr], pcplxChnlGrid u[nr]) { // 局部变量，将某个频点的格林函数谱临时存放 - cplxChnlGrid *tmp_grn = (cplxChnlGrid *)calloc(nr, sizeof(*tmp_grn)); + cplxChnlGrid *tmp_u = (cplxChnlGrid *)calloc(nr, sizeof(*tmp_u)); for(size_t ir=0; irnf2 + 1, sizeof(int)); // 实际计算的频点数 - size_t nf_valid = nf2 - nf1 + 1; + size_t nf_valid = grn->nf2 - grn->nf1 + 1; - mod1d->omgref = PI2*freqs[nf2]; + mod1d->omgref = PI2*grn->freqs[grn->nf2]; // 频率omega循环 // schedule语句可以动态调度任务，最大程度地使用计算资源 #pragma omp parallel for schedule(guided) default(shared) - for(size_t iw=nf1; iw<=nf2; ++iw){ - real_t w = freqs[iw]*PI2; // 实频率 - cplx_t omega = w - wI*I; // 复数频率 omega = w - i*wI + for(size_t iw = grn->nf1; iw <= grn->nf2; ++iw){ + real_t w = grn->freqs[iw]*PI2; // 实频率 + cplx_t omega = w - grn->wI*I; // 复数频率 omega = w - i*wI // 如果在虚频率的帮助下，频率仍然距离原点太近， // 计算会有严重的数值问题，因此直接根据频率距离原点的距离， // 跳过该频率，没有必要再计算 - if( ! keepAllFreq ) + if( ! grn->keepAllFreq ) { real_t R = 0.1; // 完全经验性地设定，暂不必要暴露在用户可控层面 if(fabs(omega) < R){ #pragma omp critical { - GRTRaiseWarning("Skip low frequency (iw=%zu, freq=%.5e).", iw, freqs[iw]); + GRTRaiseWarning("Skip low frequency (iw=%zu, freq=%.5e).", iw, grn->freqs[iw]); nf_valid--; } if(nf_valid == 0) GRTRaiseError("NO VALID FREQUENCIES."); @@ -132,13 +121,13 @@ void grt_integ_grn_spec( grt_update_mod1d_state_omega(local_mstat, omega, false); // 是否要输出积分过程文件 - bool needfstats = (statsstr!=NULL && (grt_findElement_size_t(statsidxs, nstatsidxs, iw) >= 0)); + bool needfstats = (grn->statsstr!=NULL && (grt_findElement_size_t(grn->statsidxs, grn->nstatsidxs, iw) >= 0)); // 为当前频率创建波数积分记录文件 if(needfstats){ char *suffix = NULL; - GRT_SAFE_ASPRINTF(&suffix, "_%04zu_%.5e", iw, freqs[iw]); - grt_KMET_init_fstats(nr, rs, statsstr, suffix, local_Kmet); + GRT_SAFE_ASPRINTF(&suffix, "_%04zu_%.5e", iw, grn->freqs[iw]); + grt_KMET_init_fstats(grn->nr, grn->rs, grn->statsstr, suffix, local_Kmet); GRT_SAFE_FREE_PTR(suffix); } @@ -150,20 +139,20 @@ void grt_integ_grn_spec( // Wavenumber Integration // 波数积分上限 local_Kmet->kmax = hypot(local_Kmet->k0, local_Kmet->ampk * w / local_Kmet->vmin); - K_INTEG *Kint = grt_wavenumber_integral(local_mstat, nr, rs, local_Kmet, calc_upar, grt_kernel); + K_INTEG *Kint = grt_wavenumber_integral(local_mstat, grn->nr, grn->rs, local_Kmet, grn->calc_upar, grt_kernel); // 记录到格林函数结构体内 // 如果计算核函数过程中存在除零错误，则放弃该频率【通常在大震中距的低频段】 - recordin_GRN(iw, nr, coef, Kint->sumJ, grn); - if(calc_upar){ - recordin_GRN(iw, nr, coef, Kint->sumJz, grn_uiz); - recordin_GRN(iw, nr, coef, Kint->sumJr, grn_uir); + recordin_GRN(iw, grn->nr, coef, Kint->sumJ, grn->u); + if(grn->calc_upar){ + recordin_GRN(iw, grn->nr, coef, Kint->sumJz, grn->uiz); + recordin_GRN(iw, grn->nr, coef, Kint->sumJr, grn->uir); } // =================================================================================== freq_invstats[iw] = local_mstat->stats; - if(needfstats) grt_KMET_destroy_fstats(nr, local_Kmet); + if(needfstats) grt_KMET_destroy_fstats(grn->nr, local_Kmet); grt_free_mod1d_state(local_mstat); @@ -180,9 +169,9 @@ void grt_integ_grn_spec( } // END omega loop // 打印 freq_invstats - for(size_t iw=nf1; iw<=nf2; ++iw){ + for(size_t iw = grn->nf1; iw <= grn->nf2; ++iw){ if(freq_invstats[iw]==GRT_INVERSE_FAILURE){ - GRTRaiseWarning("iw=%zu, freq=%e(Hz), meet Zero Divison Error, results are filled with 0.\n", iw, freqs[iw]); + GRTRaiseWarning("iw=%zu, freq=%e(Hz), meet Zero Divison Error, results are filled with 0.\n", iw, grn->freqs[iw]); } } GRT_SAFE_FREE_PTR(freq_invstats); diff --git a/pygrt/C_extension/src/dynamic/grnspec.c b/pygrt/C_extension/src/dynamic/grnspec.c new file mode 100644 index 0000000..ce18278 --- /dev/null +++ b/pygrt/C_extension/src/dynamic/grnspec.c @@ -0,0 +1,148 @@ +/** + * @file grnspec.c + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-04 + * + * 将格林函数频谱用一个结构体包裹 + * + * + */ + +#include + +#include "grt/dynamic/grnspec.h" +#include "grt/common/myfftw.h" +#include "grt/common/sacio.h" + +void grt_grnspec_allocate_u(GRNSPEC *grn) +{ + grn->u = (pcplxChnlGrid *) calloc(grn->nr, sizeof(*grn->u)); + grn->uiz = (grn->calc_upar)? (pcplxChnlGrid *) calloc(grn->nr, sizeof(*grn->uiz)) : NULL; + grn->uir = (grn->calc_upar)? (pcplxChnlGrid *) calloc(grn->nr, sizeof(*grn->uir)) : NULL; + + for(size_t ir = 0; ir < grn->nr; ++ir){ + GRT_LOOP_ChnlGrid(im, c){ + grn->u[ir][im][c] = (cplx_t*) calloc(grn->nf, sizeof(cplx_t)); + if(grn->uiz) grn->uiz[ir][im][c] = (cplx_t*)calloc(grn->nf, sizeof(cplx_t)); + if(grn->uir) grn->uir[ir][im][c] = (cplx_t*)calloc(grn->nf, sizeof(cplx_t)); + } + } +} + +void grt_grnspec_free_u(GRNSPEC *grn) +{ + for(size_t ir = 0; ir < grn->nr; ++ir){ + GRT_LOOP_ChnlGrid(im, c){ + GRT_SAFE_FREE_PTR(grn->u[ir][im][c]); + if(grn->uiz) GRT_SAFE_FREE_PTR(grn->uiz[ir][im][c]); + if(grn->uir) GRT_SAFE_FREE_PTR(grn->uir[ir][im][c]); + } + } + GRT_SAFE_FREE_PTR(grn->u); + GRT_SAFE_FREE_PTR(grn->uiz); + GRT_SAFE_FREE_PTR(grn->uir); +} + + + +/** 将一条数据反变换回时间域再进行处理，保存到SAC文件 */ +static void write_one_to_sac( + const char *srcname, const char ch, GRT_FFTW_HOLDER *fh, const real_t wI, + SACTRACE *sac, const char *s_output_subdir, const char *s_prefix, + const int sgn, bool skipImagComps, const cplx_t *grncplx) +{ + snprintf(sac->hd.kcmpnm, sizeof(sac->hd.kcmpnm), "%s%s%c", s_prefix, srcname, ch); + + char *s_outpath = NULL; + GRT_SAFE_ASPRINTF(&s_outpath, "%s/%s.sac", s_output_subdir, sac->hd.kcmpnm); + + // 执行fft任务会修改数组，需重新置零 + grt_reset_fftw_holder_zero(fh); + + // 赋值复数，包括时移 + cplx_t cfac, ccoef; + cfac = exp(I*PI2*fh->df*sac->hd.b); + ccoef = sgn; + // 只赋值有效长度，其余均为0 + for(size_t i = 0; i < fh->nf_valid; ++i){ + fh->W_f[i] = grncplx[i] * ccoef; + ccoef *= cfac; + } + + + if(! fh->naive_inv){ + // 发起fft任务 + fftw_execute(fh->plan); + } else { + grt_naive_inverse_transform_double(fh); + } + + // 归一化，并处理虚频 + // 并转为 SAC 需要的单精度类型 + real_t fac, coef; + coef = fh->df; + fac = 1.0; + if (! skipImagComps) { + coef *= exp(sac->hd.b * wI); + fac = exp(wI*fh->dt); + } + for(size_t i = 0; i < fh->nt; ++i){ + sac->data[i] = fh->w_t[i] * coef; + coef *= fac; + } + + + // 以sac文件保存到本地 + grt_write_SACTRACE(s_outpath, sac); + + GRT_SAFE_FREE_PTR(s_outpath); +} + + +void grt_grnspec_write_sac( + const GRNSPEC *grn, const real_t (*travtPS)[2], const real_t *begintimes, char **outputdirs, GRT_FFTW_HOLDER *fh, SACTRACE *sac, + const bool skipImagComps, const bool saveEX, const bool saveVF, const bool saveHF, const bool saveDC) +{ + // 做反傅里叶变换，保存SAC文件 + for(size_t ir = 0; ir < grn->nr; ++ir){ + real_t dist = grn->rs[ir]; + sac->hd.dist = dist; + + // 文件保存子目录 + const char *dirpath = outputdirs[ir]; + GRTCheckMakeDir(dirpath); + + // 计算理论走时 + sac->hd.t0 = travtPS[ir][0]; + strcpy(sac->hd.kt0, "P"); + sac->hd.t1 = travtPS[ir][1]; + strcpy(sac->hd.kt1, "S"); + + // 时间延迟 + sac->hd.b = begintimes[ir]; + + GRT_LOOP_ChnlGrid(im, c){ + if(! saveEX && im==GRT_SRC_M_EX_INDEX) continue; + if(! saveVF && im==GRT_SRC_M_VF_INDEX) continue; + if(! saveHF && im==GRT_SRC_M_HF_INDEX) continue; + if(! saveDC && im>=GRT_SRC_M_DD_INDEX) continue; + + int modr = GRT_SRC_M_ORDERS[im]; + int sgn=1; // 用于反转Z分量 + + if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; // 跳过输出0阶的T分量 + + // Z分量反转 + sgn = (GRT_ZRT_CODES[c]=='Z') ? -1 : 1; + + char ch = GRT_ZRT_CODES[c]; + + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, grn->wI, sac, dirpath, "", sgn, skipImagComps, grn->u[ir][im][c]); + + if(grn->calc_upar){ + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, grn->wI, sac, dirpath, "z", sgn*(-1), skipImagComps, grn->uiz[ir][im][c]); + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, grn->wI, sac, dirpath, "r", sgn , skipImagComps, grn->uir[ir][im][c]); + } + } + } +} \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 4cb793c..449c16c 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -3,8 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2024-11-28 * - * 定义main函数，形成命令行式的用法（不使用python的entry_points，会牺牲性能） - * 计算不同震源的格林函数 + * 计算不同震源的层状介质全波解 * */ @@ -902,61 +901,6 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ -/** 将一条数据反变换回时间域再进行处理，保存到SAC文件 */ -static void write_one_to_sac( - const char *srcname, const char ch, GRT_FFTW_HOLDER *fh, const real_t wI, - SACTRACE *sac, const char *s_output_subdir, const char *s_prefix, - const int sgn, bool skipImagComps, const cplx_t *grncplx) -{ - snprintf(sac->hd.kcmpnm, sizeof(sac->hd.kcmpnm), "%s%s%c", s_prefix, srcname, ch); - - char *s_outpath = NULL; - GRT_SAFE_ASPRINTF(&s_outpath, "%s/%s.sac", s_output_subdir, sac->hd.kcmpnm); - - // 执行fft任务会修改数组，需重新置零 - grt_reset_fftw_holder_zero(fh); - - // 赋值复数，包括时移 - cplx_t cfac, ccoef; - cfac = exp(I*PI2*fh->df*sac->hd.b); - ccoef = sgn; - // 只赋值有效长度，其余均为0 - for(size_t i = 0; i < fh->nf_valid; ++i){ - fh->W_f[i] = grncplx[i] * ccoef; - ccoef *= cfac; - } - - - if(! fh->naive_inv){ - // 发起fft任务 - fftw_execute(fh->plan); - } else { - grt_naive_inverse_transform_double(fh); - } - - // 归一化，并处理虚频 - // 并转为 SAC 需要的单精度类型 - real_t fac, coef; - coef = fh->df; - fac = 1.0; - if (! skipImagComps) { - coef *= exp(sac->hd.b * wI); - fac = exp(wI*fh->dt); - } - for(size_t i = 0; i < fh->nt; ++i){ - sac->data[i] = fh->w_t[i] * coef; - coef *= fac; - } - - - // 以sac文件保存到本地 - grt_write_SACTRACE(s_outpath, sac); - - GRT_SAFE_FREE_PTR(s_outpath); -} - - - /** 子模块主函数 */ int greenfn_main(int argc, char **argv) { GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); @@ -1065,19 +1009,6 @@ int greenfn_main(int argc, char **argv) { GRTCheckMakeDir(Ctrl->S.s_statsdir); } - // 建立格林函数的complex数组 - pcplxChnlGrid *grn = (pcplxChnlGrid *) calloc(Ctrl->R.nr, sizeof(*grn)); - pcplxChnlGrid *grn_uiz = (Ctrl->e.active)? (pcplxChnlGrid *) calloc(Ctrl->R.nr, sizeof(*grn_uiz)) : NULL; - pcplxChnlGrid *grn_uir = (Ctrl->e.active)? (pcplxChnlGrid *) calloc(Ctrl->R.nr, sizeof(*grn_uir)) : NULL; - - for(size_t ir=0; irR.nr; ++ir){ - GRT_LOOP_ChnlGrid(im, c){ - grn[ir][im][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); - if(grn_uiz) grn_uiz[ir][im][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); - if(grn_uir) grn_uir[ir][im][c] = (cplx_t*)calloc(Ctrl->N.nf, sizeof(cplx_t)); - } - } - // 波数积分方法 K_INTEG_METHOD KMET = {0}; @@ -1106,29 +1037,37 @@ int greenfn_main(int argc, char **argv) { if(! Ctrl->s.active){ print_Ctrl(Ctrl); } + + // 格林函数频谱 + 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_integ_grn_spec( - mod1d, Ctrl->H.nf1, Ctrl->H.nf2, Ctrl->N.freqs, Ctrl->R.nr, Ctrl->R.rs, Ctrl->N.wI, Ctrl->N.keepAllFreq, - &KMET, !Ctrl->s.active, Ctrl->e.active, - grn, grn_uiz, grn_uir, - Ctrl->S.s_statsdir, Ctrl->S.nstatsidxs, Ctrl->S.statsidxs - ); + grt_integ_grn_spec(mod1d, &KMET, grn, !Ctrl->s.active); //============================================================================== // 使用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); - - real_t (* travtPS)[2] = (real_t (*)[2])calloc(Ctrl->R.nr, sizeof(real_t)*2); + 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); - char *s_output_dirprefx = NULL; - GRT_SAFE_ASPRINTF(&s_output_dirprefx, "%s/%s_%s_%s", Ctrl->O.s_output_dir, Ctrl->M.s_modelname, Ctrl->D.s_depsrc, Ctrl->D.s_deprcv); - - // 建立SAC头文件，包含必要的头变量 + // 建立 SAC 文件原型，包含必要的头变量 SACTRACE *sac = grt_new_SACTRACE(fh->dt, fh->nt, 0.0); // 发震时刻作为参考时刻 sac->hd.o = 0.0; @@ -1149,22 +1088,20 @@ int greenfn_main(int argc, char **argv) { sac->hd.user7 = mod1d->Vb[mod1d->isrc]; sac->hd.user8 = mod1d->Rho[mod1d->isrc]; - // 做反傅里叶变换，保存SAC文件 + // 为每个震中距设置对应的变量 + real_t (*travtPS)[2] = (real_t (*)[2])calloc(grn->nr, sizeof(real_t)*2); + real_t *begintimes = (real_t *)calloc(grn->nr, sizeof(real_t)); + char **outputdirs = (char **)calloc(grn->nr, sizeof(char *)); for(size_t ir = 0; ir < Ctrl->R.nr; ++ir){ real_t dist = Ctrl->R.rs[ir]; - sac->hd.dist = dist; - // 文件保存子目录 - char *s_output_subdir = NULL; - - GRT_SAFE_ASPRINTF(&s_output_subdir, "%s_%s", s_output_dirprefx, Ctrl->R.s_rs[ir]); - GRTCheckMakeDir(s_output_subdir); + outputdirs[ir] = NULL; + GRT_SAFE_ASPRINTF(&outputdirs[ir], "%s/%s_%s_%s_%s", + Ctrl->O.s_output_dir, Ctrl->M.s_modelname, Ctrl->D.s_depsrc, Ctrl->D.s_deprcv, Ctrl->R.s_rs[ir]); // 计算理论走时 - sac->hd.t0 = grt_compute_travt1d(mod1d->Thk, mod1d->Va, mod1d->n, mod1d->isrc, mod1d->ircv, dist); - strcpy(sac->hd.kt0, "P"); - sac->hd.t1 = grt_compute_travt1d(mod1d->Thk, mod1d->Vb, mod1d->n, mod1d->isrc, mod1d->ircv, dist); - strcpy(sac->hd.kt1, "S"); + travtPS[ir][0] = grt_compute_travt1d(mod1d->Thk, mod1d->Va, mod1d->n, mod1d->isrc, mod1d->ircv, dist); + travtPS[ir][1] = grt_compute_travt1d(mod1d->Thk, mod1d->Vb, mod1d->n, mod1d->isrc, mod1d->ircv, dist); // 时间延迟 { @@ -1173,45 +1110,18 @@ int greenfn_main(int argc, char **argv) { delayT = Ctrl->E.delayT0; if(Ctrl->E.delayV0 > 0.0) delayT += sqrt( GRT_SQUARE(dist) + GRT_SQUARE(Ctrl->D.deprcv - Ctrl->D.depsrc) ) / Ctrl->E.delayV0; } else { - delayT = Ctrl->E.delayT0 + sac->hd.t0; + delayT = Ctrl->E.delayT0 + travtPS[ir][0]; } // 修改SAC头段时间变量 - sac->hd.b = delayT; - } - - GRT_LOOP_ChnlGrid(im, c){ - if(! Ctrl->G.doEX && im==GRT_SRC_M_EX_INDEX) continue; - if(! Ctrl->G.doVF && im==GRT_SRC_M_VF_INDEX) continue; - if(! Ctrl->G.doHF && im==GRT_SRC_M_HF_INDEX) continue; - if(! Ctrl->G.doDC && im>=GRT_SRC_M_DD_INDEX) continue; - - int modr = GRT_SRC_M_ORDERS[im]; - int sgn=1; // 用于反转Z分量 - - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; // 跳过输出0阶的T分量 - - // Z分量反转 - sgn = (GRT_ZRT_CODES[c]=='Z') ? -1 : 1; - - char ch = GRT_ZRT_CODES[c]; - - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "", sgn, Ctrl->N.skipImagComps, grn[ir][im][c]); - - if(Ctrl->e.active){ - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "z", sgn*(-1), Ctrl->N.skipImagComps, grn_uiz[ir][im][c]); - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "r", sgn , Ctrl->N.skipImagComps, grn_uir[ir][im][c]); - } - } - - GRT_SAFE_FREE_PTR(s_output_subdir); - - // 记录走时 - if(travtPS != NULL){ - travtPS[ir][0] = sac->hd.t0; - travtPS[ir][1] = sac->hd.t1; + begintimes[ir] = delayT; } } + // 反傅里叶变换后保存到SAC文件 + grt_grnspec_write_sac( + grn, travtPS, begintimes, outputdirs, fh, sac, + Ctrl->N.skipImagComps, Ctrl->G.doEX, Ctrl->G.doVF, Ctrl->G.doHF, Ctrl->G.doDC); + // 输出警告：当震源位于液体层中时，仅允许计算爆炸源对应的格林函数 if(mod1d->isLiquid[mod1d->isrc]){ GRTRaiseWarning( @@ -1219,7 +1129,6 @@ int greenfn_main(int argc, char **argv) { "therefore only the Green's Funtions for the Explosion source will be computed.\n"); } - grt_free_SACTRACE(sac); // 打印走时 @@ -1235,19 +1144,13 @@ int greenfn_main(int argc, char **argv) { } // 释放内存 - for(size_t ir=0; irR.nr; ++ir){ - GRT_LOOP_ChnlGrid(im, c){ - GRT_SAFE_FREE_PTR(grn[ir][im][c]); - if(grn_uiz) GRT_SAFE_FREE_PTR(grn_uiz[ir][im][c]); - if(grn_uir) GRT_SAFE_FREE_PTR(grn_uir[ir][im][c]); - } - } - GRT_SAFE_FREE_PTR(grn); - GRT_SAFE_FREE_PTR(grn_uiz); - GRT_SAFE_FREE_PTR(grn_uir); + grt_grnspec_free_u(grn); GRT_SAFE_FREE_PTR(travtPS); - + GRT_SAFE_FREE_PTR(begintimes); + grt_free_SACTRACE(sac); grt_destroy_fftw_holder(fh); + for(size_t ir = 0; ir < grn->nr; ++ir) GRT_SAFE_FREE_PTR(outputdirs[ir]); + GRT_SAFE_FREE_PTR(outputdirs); free_Ctrl(Ctrl); return EXIT_SUCCESS; diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 9a6c8bc..dc09fef 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -29,18 +29,7 @@ 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), c_size_t, c_size_t, PREAL, - c_size_t, PREAL, REAL, c_bool, POINTER(c_K_INTEG_METHOD), - c_bool, c_bool, - POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM), - POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM), - POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM), - - c_char_p, - c_size_t, - POINTER(c_size_t) -] +C_grt_integ_grn_spec.argtypes = [POINTER(c_MODEL1D), POINTER(c_K_INTEG_METHOD), POINTER(c_GRNSPEC), c_bool] C_grt_integ_static_grn = libgrt.grt_integ_static_grn diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index 4d528a8..921fe13 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -33,7 +33,8 @@ "PCPLX", "c_MODEL1D", - "c_K_INTEG_METHOD" + "c_K_INTEG_METHOD", + "c_GRNSPEC" ] @@ -123,4 +124,29 @@ class c_K_INTEG_METHOD(Structure): ('fstats', c_void_p), ('ptam_fstatsnr', c_void_p), + ] + +class c_GRNSPEC(Structure): + """ + 和 C 结构体 GRNSPEC 作匹配 + """ + + _fields_ = [ + ('nf', c_size_t), + ('freqs', PREAL), + ('nf1', c_size_t), + ('nf2', c_size_t), + ('nr', c_size_t), + ('rs', PREAL), + ('wI', REAL), + ('keepAllFreq', c_bool), + ('calc_upar', c_bool), + + ('u', POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM)), + ('uiz', POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM)), + ('uir', POINTER((PCPLX*CHANNEL_NUM)*SRC_M_NUM)), + + ('statsstr', c_char_p), + ('nstatsidxs', c_size_t), + ('statsidxs', POINTER(c_size_t)), ] \ No newline at end of file diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 81809d7..ad34e2c 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -352,9 +352,8 @@ def _get_grn_spectra( print(f"statsfile_index=", statsidxs) - + # ==================================================================== KMET = c_K_INTEG_METHOD() - hs = max(abs(depsrc - deprcv), 1.0) KMET.k0 = k0 * np.pi / hs KMET.ampk = ampk @@ -373,6 +372,27 @@ def _get_grn_spectra( KMET.applyDCM = converg_method == 'DCM' KMET.applyPTAM = converg_method == 'PTAM' + # ==================================================================== + + + # ==================================================================== + 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库函数 #///////////////////////////////////////////////////////////////////////////////// @@ -381,12 +401,7 @@ def _get_grn_spectra( # 爆炸源 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, nf1, nf2, c_freqs, nrs, c_rs, wI, keepAllFreq, pointer(KMET), - print_runtime, calc_upar, - c_grnArr, c_grnArr_uiz, c_grnArr_uir, - c_statsfile, nstatsidxs, c_statsidxs - ) + C_grt_integ_grn_spec(self.c_mod1d, pointer(KMET), pointer(grn), print_runtime) #================================================================================= #/////////////////////////////////////////////////////////////////////////////////