From a379ec7d0a64b644685f18e625a4f133d256c29f Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 4 Aug 2026 20:24:07 +0800 Subject: [PATCH] REFAC: expose dynamic GF synthesis through C API --- pygrt/C_extension/include/grt.h | 1 - pygrt/C_extension/include/grt/common/const.h | 1 + pygrt/C_extension/include/grt/dynamic/syn.h | 25 -- pygrt/C_extension/src/dynamic/grt_syn.c | 426 +++++++++++++------ pygrt/C_extension/src/dynamic/syn.c | 58 --- pygrt/c_interfaces.py | 14 + pygrt/utils.py | 124 ++++-- 7 files changed, 388 insertions(+), 261 deletions(-) delete mode 100644 pygrt/C_extension/include/grt/dynamic/syn.h delete mode 100644 pygrt/C_extension/src/dynamic/syn.c diff --git a/pygrt/C_extension/include/grt.h b/pygrt/C_extension/include/grt.h index 7c4a2895..c21ee300 100644 --- a/pygrt/C_extension/include/grt.h +++ b/pygrt/C_extension/include/grt.h @@ -46,7 +46,6 @@ #include "grt/dynamic/dyn_postprocess.h" #include "grt/dynamic/signals.h" #include "grt/dynamic/source.h" -#include "grt/dynamic/syn.h" #include "grt/modal/eigenfn.h" diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index 6da83b5d..cb6742e0 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -138,6 +138,7 @@ typedef cplx_t cplxChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef cplx_t* pcplxChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef real_t realChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef real_t* prealChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; +typedef float* pfloatChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef int intChnlGrid[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]; typedef cplx_t cplxIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM]; diff --git a/pygrt/C_extension/include/grt/dynamic/syn.h b/pygrt/C_extension/include/grt/dynamic/syn.h deleted file mode 100644 index 05c4805d..00000000 --- a/pygrt/C_extension/include/grt/dynamic/syn.h +++ /dev/null @@ -1,25 +0,0 @@ -/** - * @file syn.h - * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) - * @date 2025-12 - * - * 将合成地震图部分单独用一个源文件和若干函数来管理 - * - */ - -#pragma once - -#include "grt/common/const.h" -#include "grt/common/sacio.h" -#include "grt/common/radiation.h" - -/** - * 根据已有系数（方向因子）合成理论地震图 - * - * @param[in] srcRadi 不同震源不同分量的系数 - * @param[in] computeType 要计算的震源类型 - * @param[in] dirpath 格林函数所在目录路径 - * @param[in] prefix 格林函数文件名前缀 - * @param[out] synsac 三分量 SACTRACE - */ -void grt_syn(const realChnlGrid srcRadi, const GRT_SYN_TYPE computeType, const char *dirpath, const char *prefix, SACTRACE *synsac[GRT_CHANNEL_NUM]); \ No newline at end of file diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index 7a666cc7..933c34be 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -487,198 +487,362 @@ static void save_to_sac(GRT_MODULE_CTRL *Ctrl, const char *pfx, const char ch, S } -static void data_zrt2zne(SACTRACE *synsac[3], SACTRACE *synparsac[3][3], real_t azrad) +/** 判断该震源类型是否参与合成 */ +static bool syn_need_src(GRT_SYN_TYPE computeType, int im) { - real_t dblsyn[3] = {0}; - real_t dblupar[3][3] = {0}; + if (computeType == GRT_SYN_EX) { + return im == GRT_SRC_M_EX_INDEX; + } else if (computeType == GRT_SYN_SF) { + return im == GRT_SRC_M_VF_INDEX || im == GRT_SRC_M_HF_INDEX; + } else if (computeType == GRT_SYN_DC) { + return im >= GRT_SRC_M_DD_INDEX; + } else if (computeType == GRT_SYN_TS || computeType == GRT_SYN_MT) { + return im >= GRT_SRC_M_DD_INDEX || im == GRT_SRC_M_EX_INDEX; + } + return false; +} - bool doupar = (synparsac[0][0]!=NULL); - float dist = synsac[0]->hd.dist; - int nt = synsac[0]->hd.npts; +/** 线性叠加：out += coef * gf，跳过零系数；非零系数时 gf 不可为 NULL */ +static void syn_accum_from_gf( + size_t npts, const pfloatChnlGrid gf, + const realChnlGrid srcRadi, float *const out[GRT_CHANNEL_NUM]) +{ + GRT_LOOP_ChnlGrid(im, c) { + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + + const real_t coef = srcRadi[im][c]; + if(coef == 0.0) continue; + if(gf[im][c] == NULL){ + GRTRaiseError("Missing Green function for %s%c.\n", + GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + } - // 对每一个时间点 - for(int n = 0; n < nt; ++n){ - // 复制数据，以调用函数 - for(int i1=0; i1<3; ++i1){ - dblsyn[i1] = synsac[i1]->data[n]; - for(int i2=0; i2<3; ++i2){ - if(doupar) dblupar[i1][i2] = synparsac[i1][i2]->data[n]; + float *dst = out[c]; + const float *src = gf[im][c]; + for(size_t n = 0; n < npts; ++n){ + dst[n] += (float)(src[n] * coef); + } + } +} + + +/** + * 由动态格林函数合成三分量地震图（及可选空间偏导）。 + * + * 数组布局：gf[震源][分量][采样点]、syn[分量][采样点]、 + * syn_upar[偏导方向][分量][采样点]。gf_uiz/gf_uir 在 calc_upar=false + * 时可传 NULL；单个分量指针为 NULL 时跳过该道。 + * + * @param[in] azrad 方位角（弧度） + */ +void grt_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, + const real_t mchn[GRT_MECHANISM_NUM], + bool rot2ZNE, bool calc_upar, + float *const syn[GRT_CHANNEL_NUM], float *const syn_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + const real_t az = azrad; + + int calcUTypes = calc_upar ? 4 : 1; + realChnlGrid srcRadi = {0}; + + // 清零输出 + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + memset(syn[c], 0, npts * sizeof(float)); + if(calc_upar){ + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + memset(syn_upar[c][c2], 0, npts * sizeof(float)); + } + } + } + + for(int ityp = 0; ityp < calcUTypes; ++ityp){ + real_t upar_scale = 1.0; + // 求位移空间导数时，需调整比例系数（1e-5: km→cm） + // ZRT 协变导数拆两步：此处算 (1/r)∂_θ u，后处理再补 ±u/r。 + if(ityp > 0){ + switch (GRT_ZRT_CODES[ityp-1]){ + case 'Z': case 'R': + upar_scale = 1e-5; + break; + case 'T': + upar_scale = 1e-5 / dist; + break; + default: + break; } } - if(doupar) { - grt_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5); // 1e5 km 转为 cm + float *const (*up)[GRT_CHANNEL_NUM] = gf; + if(ityp == 1){ + up = gf_uiz; + } else if(ityp == 2){ + up = gf_uir; + } + + memset(srcRadi, 0, sizeof(srcRadi)); + grt_set_source_radiation(srcRadi, computeType, (ityp == 3), M0, upar_scale, VpVs_ratio, az, mchn); + + float *out_ptrs[GRT_CHANNEL_NUM]; + if(ityp == 0){ + for(int c = 0; c < GRT_CHANNEL_NUM; ++c) out_ptrs[c] = syn[c]; } else { - grt_rot_zxy2zrt_vec(-azrad, dblsyn); + for(int c = 0; c < GRT_CHANNEL_NUM; ++c) out_ptrs[c] = syn_upar[ityp-1][c]; } + syn_accum_from_gf(npts, up, srcRadi, out_ptrs); + } - // 将结果写入原数组 - for(int i1=0; i1<3; ++i1){ - synsac[i1]->data[n] = dblsyn[i1]; - for(int i2=0; i2<3; ++i2){ - if(doupar) synparsac[i1][i2]->data[n] = dblupar[i1][i2]; + // 是否转到 ZNE + if(rot2ZNE){ + real_t dblsyn[3] = {0}; + real_t dblupar[3][3] = {0}; + for(size_t n = 0; n < npts; ++n){ + for(int i1 = 0; i1 < GRT_CHANNEL_NUM; ++i1){ + dblsyn[i1] = syn[i1][n]; + if(calc_upar){ + for(int i2 = 0; i2 < GRT_CHANNEL_NUM; ++i2){ + dblupar[i1][i2] = syn_upar[i1][i2][n]; + } + } + } + if(calc_upar){ + grt_rot_zrt2zxy_upar(az, dblsyn, dblupar, dist * 1e5); // 1e5 km→cm + } else { + grt_rot_zxy2zrt_vec(-az, dblsyn); + } + for(int i1 = 0; i1 < GRT_CHANNEL_NUM; ++i1){ + syn[i1][n] = (float)dblsyn[i1]; + if(calc_upar){ + for(int i2 = 0; i2 < GRT_CHANNEL_NUM; ++i2){ + syn_upar[i1][i2][n] = (float)dblupar[i1][i2]; + } + } } } } } +/** 读取一道格林函数 SAC；不存在则返回 NULL（允许 m=0 的 T） */ +static SACTRACE *syn_load_one_gf(const char *dirpath, const char *prefix, int im, int c) +{ + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') return NULL; + + char *grnpath = NULL; + GRT_SAFE_ASPRINTF(&grnpath, "%s/%s%s%c.sac", dirpath, prefix, GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + if(access(grnpath, F_OK) != 0){ + GRT_SAFE_FREE_PTR(grnpath); + return NULL; + } + SACTRACE *sac = grt_read_SACTRACE(grnpath, false); + GRT_SAFE_FREE_PTR(grnpath); + return sac; +} -/** 子模块主函数 */ -int syn_main(int argc, char **argv){ - GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); - - getopt_from_command(Ctrl, argc, argv); - // 输出分量格式，即是否需要旋转到ZNE - bool rot2ZNE = Ctrl->N.active; +/** 对一道合成结果做时间函数卷积 / 积分 / 微分 */ +static void syn_postprocess_trace(SACTRACE *sac, SACTRACE *tfsac, int int_times, int dif_times) +{ + float dt = sac->hd.delta; + int nt = sac->hd.npts; - // 根据参数设置，选择分量名 - const char *chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; + if(tfsac != NULL){ + float wI = sac->hd.user0; + float fac = 1.0f; + float dfac = expf(-wI * dt); + for(int n = 0; n < nt; ++n){ + sac->data[n] *= fac; + if(n < tfsac->hd.npts) tfsac->data[n] *= fac; + fac *= dfac; + } - SACTRACE *synsac[GRT_CHANNEL_NUM] = {0}; - SACTRACE *synparsac[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = {0}; - SACTRACE **sacs = NULL; - SACTRACE *tfsac = NULL; + float *convarr = (float *)calloc(nt, sizeof(float)); + grt_oaconvolve(sac->data, nt, tfsac->data, tfsac->hd.npts, convarr, nt, true); + fac = 1.0f; + dfac = expf(wI * dt); + for(int n = 0; n < nt; ++n){ + sac->data[n] = convarr[n] * fac * dt; + if(n < tfsac->hd.npts) tfsac->data[n] *= fac; + fac *= dfac; + } + GRT_SAFE_FREE_PTR(convarr); + } - real_t upar_scale=1.0; + for(int i = 0; i < int_times; ++i){ + grt_trap_integral(sac->data, nt, dt); + } + for(int i = 0; i < dif_times; ++i){ + grt_differential(sac->data, nt, dt); + } +} - // 计算和位移相关量的种类（1-位移，2-ui_z，3-ui_r，4-ui_t） - int calcUTypes = (Ctrl->e.active)? 4 : 1; - for(int ityp=0; itypdist; - break; - - default: - break; + bool rot2ZNE = Ctrl->N.active; + bool calc_upar = Ctrl->e.active; + const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; + + // 读入格林函数到内存 + SACTRACE *gf_sac[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = {{0}}; + SACTRACE *gf_uiz_sac[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = {{0}}; + SACTRACE *gf_uir_sac[GRT_SRC_M_NUM][GRT_CHANNEL_NUM] = {{0}}; + pfloatChnlGrid gf = {{0}}; + pfloatChnlGrid gf_uiz = {{0}}; + pfloatChnlGrid gf_uir = {{0}}; + + SACTRACE *tmpl = NULL; + GRT_LOOP_ChnlGrid(im, c) { + if(!syn_need_src(Ctrl->computeType, im)) continue; + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + + gf_sac[im][c] = syn_load_one_gf(Ctrl->G.s_grnpath, "", im, c); + if(gf_sac[im][c] == NULL){ + GRTRaiseError("Failed to read Green function %s%s%c.sac\n", + "", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); } - - if (ityp == 0) { - sacs = &synsac[0]; - } else { - sacs = &synparsac[ityp-1][0]; + gf[im][c] = gf_sac[im][c]->data; + if(tmpl == NULL) tmpl = gf_sac[im][c]; + + if(calc_upar){ + gf_uiz_sac[im][c] = syn_load_one_gf(Ctrl->G.s_grnpath, "z", im, c); + gf_uir_sac[im][c] = syn_load_one_gf(Ctrl->G.s_grnpath, "r", im, c); + if(gf_uiz_sac[im][c] == NULL || gf_uir_sac[im][c] == NULL){ + GRTRaiseError("Failed to read Green function derivatives for %s%c\n", + GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + } + gf_uiz[im][c] = gf_uiz_sac[im][c]->data; + gf_uir[im][c] = gf_uir_sac[im][c]->data; } + } + if(tmpl == NULL){ + GRTRaiseError("No Green functions loaded.\n"); + } - // 重新计算方向因子 - grt_set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, (ityp==3), Ctrl->S.M0, upar_scale, Ctrl->VpVs_ratio, Ctrl->A.azrad, Ctrl->mchn); + int npts = tmpl->hd.npts; + float dt = tmpl->hd.delta; - // 合成地震图 - if (ityp==0 || ityp==3) { - grt_syn(Ctrl->srcRadi, Ctrl->computeType, Ctrl->G.s_grnpath, "", sacs); - } else { - char prefix[] = {tolower(GRT_ZRT_CODES[ityp-1]), '\0'}; - grt_syn(Ctrl->srcRadi, Ctrl->computeType, Ctrl->G.s_grnpath, prefix, sacs); - } - - // 首次读取获得时间函数 - if(Ctrl->D.active && tfsac == NULL){ - int tfnt; - float *tfarr = grt_get_time_function(&tfnt, sacs[0]->hd.delta, Ctrl->D.tftype, Ctrl->D.tfparams); - if(tfarr == NULL){ - GRTRaiseError("get time function error.\n"); + // 分配合成结果（与格林函数同头段，数据清零） + SACTRACE *synsac[GRT_CHANNEL_NUM] = {0}; + SACTRACE *synparsac[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = {{0}}; + float *syn[GRT_CHANNEL_NUM]; + float *syn_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = {{0}}; + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + synsac[c] = grt_copy_SACTRACE(tmpl, true); + syn[c] = synsac[c]->data; + if(calc_upar){ + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + synparsac[c][c2] = grt_copy_SACTRACE(tmpl, true); + syn_upar[c][c2] = synparsac[c][c2]->data; } - tfsac = grt_new_SACTRACE(sacs[0]->hd.delta, tfnt, 0.0); - memcpy(tfsac->data, tfarr, sizeof(float)*tfnt); - GRT_SAFE_FREE_PTR(tfarr); - } - - for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ - float dt = sacs[0]->hd.delta; - int nt = sacs[0]->hd.npts; - - // 时域循环卷积 - if(tfsac != NULL){ - float wI; - float fac, dfac; - // 虚频率幅值压制 - wI = sacs[0]->hd.user0; - fac = 1.0; - dfac = expf(- wI*dt); - for(int n = 0; n < nt; ++n){ - sacs[c]->data[n] *= fac; - if (n < tfsac->hd.npts) tfsac->data[n] *= fac; - fac *= dfac; - } + } + } - float *convarr = (float *)calloc(nt, sizeof(float)); - grt_oaconvolve(sacs[c]->data, nt, tfsac->data, tfsac->hd.npts, convarr, nt, true); - // 虚频率振幅恢复 - fac = 1.0; - dfac = expf(wI*dt); - for(int n = 0; n < nt; ++n){ - sacs[c]->data[n] = convarr[n] * fac * dt; // dt是连续卷积的系数 - if (n < tfsac->hd.npts) tfsac->data[n] *= fac; - fac *= dfac; - } - GRT_SAFE_FREE_PTR(convarr); - } + grt_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, + false, calc_upar, syn, syn_upar); - // 时域积分或求导 - for(int i=0; iI.int_times; ++i){ - grt_trap_integral(sacs[c]->data, nt, dt); - } - for(int i=0; iJ.dif_times; ++i){ - grt_differential(sacs[c]->data, nt, dt); + // 时间函数 / 积分 / 微分（在旋转前，与历史行为一致） + SACTRACE *tfsac = NULL; + if(Ctrl->D.active){ + int tfnt; + float *tfarr = grt_get_time_function(&tfnt, dt, Ctrl->D.tftype, Ctrl->D.tfparams); + if(tfarr == NULL){ + GRTRaiseError("get time function error.\n"); + } + tfsac = grt_new_SACTRACE(dt, tfnt, 0.0); + memcpy(tfsac->data, tfarr, sizeof(float) * tfnt); + GRT_SAFE_FREE_PTR(tfarr); + } + + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + syn_postprocess_trace(synsac[c], tfsac, Ctrl->I.int_times, Ctrl->J.dif_times); + if(calc_upar){ + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + syn_postprocess_trace(synparsac[c][c2], tfsac, Ctrl->I.int_times, Ctrl->J.dif_times); } } } - // 是否需要旋转 + // 旋转到 ZNE（在时域处理后） if(rot2ZNE){ - data_zrt2zne(synsac, synparsac, Ctrl->A.azrad); + real_t dblsyn[3] = {0}; + real_t dblupar[3][3] = {0}; + for(int n = 0; n < npts; ++n){ + for(int i1 = 0; i1 < GRT_CHANNEL_NUM; ++i1){ + dblsyn[i1] = synsac[i1]->data[n]; + if(calc_upar){ + for(int i2 = 0; i2 < GRT_CHANNEL_NUM; ++i2){ + dblupar[i1][i2] = synparsac[i1][i2]->data[n]; + } + } + } + if(calc_upar){ + grt_rot_zrt2zxy_upar(Ctrl->A.azrad, dblsyn, dblupar, Ctrl->dist * 1e5); + } else { + grt_rot_zxy2zrt_vec(-Ctrl->A.azrad, dblsyn); + } + for(int i1 = 0; i1 < GRT_CHANNEL_NUM; ++i1){ + synsac[i1]->data[n] = (float)dblsyn[i1]; + if(calc_upar){ + for(int i2 = 0; i2 < GRT_CHANNEL_NUM; ++i2){ + synparsac[i1][i2]->data[n] = (float)dblupar[i1][i2]; + } + } + } + } } - // 保存到SAC文件 - for(int i1=0; i1e.active){ - for(int i2=0; i2O.s_output_dir); grt_write_SACTRACE(buffer, tfsac); GRT_SAFE_FREE_PTR(buffer); } - - if(! Ctrl->s.active) { + + if(!Ctrl->s.active){ GRTRaiseInfo("Under \"%s\"", Ctrl->O.s_output_dir); GRTRaiseInfo("Synthetic Seismograms of %-13s source done.", srcTypeFullName[Ctrl->computeType]); if(tfsac != NULL) GRTRaiseInfo("Time Function saved."); } if(tfsac != NULL) grt_free_SACTRACE(tfsac); - for(int i=0; i<3; ++i){ + for(int i = 0; i < GRT_CHANNEL_NUM; ++i){ grt_free_SACTRACE(synsac[i]); - for(int j=0; j<3; ++j){ - grt_free_SACTRACE(synparsac[i][j]); + for(int j = 0; j < GRT_CHANNEL_NUM; ++j){ + if(synparsac[i][j] != NULL) grt_free_SACTRACE(synparsac[i][j]); } } + GRT_LOOP_ChnlGrid(im, c) { + if(gf_sac[im][c] != NULL) grt_free_SACTRACE(gf_sac[im][c]); + if(gf_uiz_sac[im][c] != NULL) grt_free_SACTRACE(gf_uiz_sac[im][c]); + if(gf_uir_sac[im][c] != NULL) grt_free_SACTRACE(gf_uir_sac[im][c]); + } free_Ctrl(Ctrl); return EXIT_SUCCESS; diff --git a/pygrt/C_extension/src/dynamic/syn.c b/pygrt/C_extension/src/dynamic/syn.c deleted file mode 100644 index af0aabde..00000000 --- a/pygrt/C_extension/src/dynamic/syn.c +++ /dev/null @@ -1,58 +0,0 @@ -/** - * @file syn.c - * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) - * @date 2025-12 - * - * 将合成地震图部分单独用一个源文件和若干函数来管理 - * - */ - -#include -#include -#include "grt/dynamic/syn.h" - -void grt_syn(const realChnlGrid srcRadi, const GRT_SYN_TYPE computeType, const char *dirpath, const char *prefix, SACTRACE *synsac[GRT_CHANNEL_NUM]) -{ - GRT_LOOP_ChnlGrid(im, c) { - int modr = GRT_SRC_M_ORDERS[im]; - if(modr == 0 && GRT_QWV_CODES[c] == 'v') continue; - - if (computeType == GRT_SYN_EX) { - if (im != GRT_SRC_M_EX_INDEX) continue; - } else if (computeType == GRT_SYN_SF) { - if (im != GRT_SRC_M_VF_INDEX && im != GRT_SRC_M_HF_INDEX) continue; - } else if (computeType == GRT_SYN_DC) { - if (im < GRT_SRC_M_DD_INDEX) continue; - } else if (computeType == GRT_SYN_TS || computeType == GRT_SYN_MT) { - if (im < GRT_SRC_M_DD_INDEX && im != GRT_SRC_M_EX_INDEX) continue; - } else { - GRTRaiseError("Not Supported."); - } - - const real_t coef = srcRadi[im][c]; - const char ch = GRT_ZRT_CODES[c]; - - // 读取格林函数 - SACTRACE *grnsac = NULL; - { - char *grnpath = NULL; - GRT_SAFE_ASPRINTF(&grnpath, "%s/%s%s%c.sac", dirpath, prefix, GRT_SRC_M_NAME_ABBR[im], ch); - grnsac = grt_read_SACTRACE(grnpath, false); - GRT_SAFE_FREE_PTR(grnpath); - } - - // 如果是第一次读取，先初始化结果指针 - if (synsac[0] == NULL) { - for(int ii = 0; ii < GRT_CHANNEL_NUM; ++ii){ - synsac[ii] = grt_copy_SACTRACE(grnsac, true); - } - } - - // 线性叠加 - for(int n = 0; n < grnsac->hd.npts; ++n){ - synsac[c]->data[n] += grnsac->data[n] * coef; - } - - grt_free_SACTRACE(grnsac); - } -} \ No newline at end of file diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 1fbfa045..1cff66d9 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -66,6 +66,20 @@ REAL_CHNL_PTRS = PREAL * CHANNEL_NUM REAL_CHNL_PTR_MAT = REAL_CHNL_PTRS * CHANNEL_NUM +C_grt_syn_from_gf = libgrt.grt_syn_from_gf +"""由动态格林函数合成三分量地震图（及可选空间偏导）""" +C_grt_syn_from_gf.restype = None +C_grt_syn_from_gf.argtypes = [ + c_size_t, c_float, + POINTER(FPOINTER * CHANNEL_NUM), + POINTER(FPOINTER * CHANNEL_NUM), + POINTER(FPOINTER * CHANNEL_NUM), + c_int, REAL, REAL, REAL, REAL * MECHANISM_NUM, + c_bool, c_bool, + FPOINTER * CHANNEL_NUM, POINTER(FPOINTER * CHANNEL_NUM), +] + + C_grt_compute_stress = libgrt.grt_compute_stress """由动态位移偏导合成应力张量""" C_grt_compute_stress.restype = None diff --git a/pygrt/utils.py b/pygrt/utils.py index a047ad9b..bb5daacb 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -90,65 +90,97 @@ def _gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:fl :param ZNE: 是否以ZNE分量输出? """ - chs = ZRTchs - sacin_prefixes = ["", "z", "r", ""] # 输入通道名 - sacout_prefixes = ["", "z", "r", "t"] # 输出通道名 - srcName = ["EX", "VF", "HF", "DD", "DS", "SS"] allchs = [tr.stats.channel for tr in st] # 为张裂计算 Vp/Vs src_va = st[0].stats.sac['user6'] src_vb = st[0].stats.sac['user7'] - VpVs_ratio = src_va / src_vb + VpVs_ratio = float(src_va / src_vb) + + azrad = float(np.deg2rad(az)) + dist = float(st[0].stats.sac['dist']) + npts = int(st[0].stats.npts) baz = 180 + az if baz > 360: baz -= 360 - azrad = np.deg2rad(az) - - calcUTypes = 4 if calc_upar else 1 + FPtrs = FPOINTER * CHANNEL_NUM + FGrid = FPtrs * SRC_M_NUM + FMat = FPtrs * CHANNEL_NUM + + def _trace_data(channel:str): + if channel not in allchs: + raise ValueError(f"Failed, channel=\"{channel}\" not exists.") + return np.ascontiguousarray(st.select(channel=channel)[0].data, dtype=np.float32) + + # 格林函数数组布局：gf[震源][分量][采样点] + gf_hold = [[None] * CHANNEL_NUM for _ in range(SRC_M_NUM)] + gf_uiz_hold = [[None] * CHANNEL_NUM for _ in range(SRC_M_NUM)] + gf_uir_hold = [[None] * CHANNEL_NUM for _ in range(SRC_M_NUM)] + for im, src_name in enumerate(SRC_M_NAME_ABBR): + for ic, ch in enumerate(ZRTchs): + # m=0 无 T 分量 + if SRC_M_ORDERS[im] == 0 and ch == 'T': + continue + channel = f'{src_name}{ch}' + # 仅加载存在的道（未用到的震源可缺省；C 侧对非零系数会校验） + if channel not in allchs: + continue + gf_hold[im][ic] = _trace_data(channel) + if calc_upar: + gf_uiz_hold[im][ic] = _trace_data(f'z{src_name}{ch}') + gf_uir_hold[im][ic] = _trace_data(f'r{src_name}{ch}') + + def _src_chnl_ptrs(hold): + return FGrid(*( + FPtrs(*( + (arr.ctypes.data_as(FPOINTER) if arr is not None else FPOINTER()) + for arr in row + )) + for row in hold + )) + + gf_ptrs = _src_chnl_ptrs(gf_hold) + gf_uiz_ptrs = _src_chnl_ptrs(gf_uiz_hold) if calc_upar else None + gf_uir_ptrs = _src_chnl_ptrs(gf_uir_hold) if calc_upar else None + + synarr = np.zeros((CHANNEL_NUM, npts), dtype=np.float32) + syn_upar_arr = np.zeros((CHANNEL_NUM, CHANNEL_NUM, npts), dtype=np.float32) + syn_ptrs = FPtrs(*(synarr[c].ctypes.data_as(FPOINTER) for c in range(CHANNEL_NUM))) + syn_upar_ptrs = FMat(*( + FPtrs(*(syn_upar_arr[d, c].ctypes.data_as(FPOINTER) for c in range(CHANNEL_NUM))) + for d in range(CHANNEL_NUM) + )) + + # ======================================================================== + # 调用 C 函数 + mchn = _set_source_mechanism(compute_type, **kwargs) + C_grt_syn_from_gf( + npts, dist, + gf_ptrs, gf_uiz_ptrs, gf_uir_ptrs, + compute_type.value, M0, VpVs_ratio, azrad, npct.as_ctypes(mchn), + ZNE, calc_upar, + syn_ptrs, syn_upar_ptrs, + ) + # ======================================================================== + out_chs = ZNEchs if ZNE else ZRTchs stall = Stream() + for c, ch in enumerate(out_chs): + tr:Trace = st[0].copy() + tr.data = synarr[c].copy() + tr.stats.channel = kcmpnm = f'{ch}' + __check_trace_attr_sac(tr, az=az, baz=baz, kcmpnm=kcmpnm) + stall.append(tr) + if calc_upar: + for d, dch in enumerate(out_chs): + tr = st[0].copy() + tr.data = syn_upar_arr[d, c].copy() + tr.stats.channel = kcmpnm = f'{dch.lower()}{ch}' + __check_trace_attr_sac(tr, az=az, baz=baz, kcmpnm=kcmpnm) + stall.append(tr) - dist = st[0].stats.sac['dist'] - upar_scale:float = 1.0 - for ityp in range(calcUTypes): - if ityp > 0: - upar_scale = 1e-5 - if ityp == 3: - upar_scale /= dist - - srcRadi = _set_source_radi(ityp==3, upar_scale, compute_type, M0, azrad, VpVs_ratio=VpVs_ratio, **kwargs) - - inpref = sacin_prefixes[ityp] - outpref = sacout_prefixes[ityp] - - for c in range(CHANNEL_NUM): - ch = chs[c] - tr:Trace = st[0].copy() - tr.data[:] = 0.0 - tr.stats.channel = kcmpnm = f'{outpref}{ch}' - __check_trace_attr_sac(tr, az=az, baz=baz, kcmpnm=kcmpnm) - for k in range(SRC_M_NUM): - coef = srcRadi[k, c] - if coef==0.0: - continue - - # 读入数据 - channel = f'{inpref}{srcName[k]}{ch}' - if channel not in allchs: - raise ValueError(f"Failed, channel=\"{channel}\" not exists.") - - tr0 = st.select(channel=channel)[0].copy() - - tr.data += coef*tr0.data - - stall.append(tr) - - if ZNE: - stall = _data_zrt2zne(stall) - return stall