From 464944c9c1825e4617d809e8d19cd4b63ddb2fbc Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Tue, 4 Aug 2026 13:00:21 +0800 Subject: [PATCH] REFAC: extract stress/strain/rotation into C APIs and call them from Python --- pygrt/C_extension/include/grt.h | 2 + .../include/grt/dynamic/dyn_postprocess.h | 35 ++ .../include/grt/static/static_postprocess.h | 42 +++ pygrt/C_extension/src/dynamic/grt_rotation.c | 82 +++-- pygrt/C_extension/src/dynamic/grt_strain.c | 92 ++++-- pygrt/C_extension/src/dynamic/grt_stress.c | 227 +++++++------ .../src/static/grt_static_rotation.c | 56 ++-- .../src/static/grt_static_strain.c | 63 ++-- .../src/static/grt_static_stress.c | 77 ++--- pygrt/c_interfaces.py | 59 +++- pygrt/utils.py | 312 +++++------------- 11 files changed, 551 insertions(+), 496 deletions(-) create mode 100644 pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h create mode 100644 pygrt/C_extension/include/grt/static/static_postprocess.h diff --git a/pygrt/C_extension/include/grt.h b/pygrt/C_extension/include/grt.h index cca67a87..7c4a2895 100644 --- a/pygrt/C_extension/include/grt.h +++ b/pygrt/C_extension/include/grt.h @@ -43,6 +43,7 @@ #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" #include "grt/dynamic/syn.h" @@ -59,6 +60,7 @@ #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 new file mode 100644 index 00000000..db06c83b --- /dev/null +++ b/pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h @@ -0,0 +1,35 @@ +/** + * @file dyn_postprocess.h + * @brief 动态位移偏导后处理（应力/应变/旋转） + * @date 2026-07 + */ + +#pragma once + +#include "grt/common/const.h" + +/** + * 由动态位移偏导合成应变张量。 + * 数组布局：u[分量][采样点]、upar[偏导方向][分量][采样点]、 + * res[第二分量][第一分量][采样点]。 + */ +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/static/static_postprocess.h b/pygrt/C_extension/include/grt/static/static_postprocess.h new file mode 100644 index 00000000..36b83838 --- /dev/null +++ b/pygrt/C_extension/include/grt/static/static_postprocess.h @@ -0,0 +1,42 @@ +/** + * @file static_postprocess.h + * @brief 静态位移偏导后处理（应力/应变/旋转） + * @date 2026-07 + */ + +#pragma once + +#include "grt/common/const.h" + +/** + * 由静态位移偏导合成对称应力张量。 + * + * 数组布局：u[分量][点]、upar[偏导方向][分量][点]、 + * res[第二分量][第一分量][点]。仅写入 res 的上三角分量。 + */ +void grt_static_compute_stress( + size_t nx, size_t ny, const real_t *xs, const real_t *ys, + 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 nx, size_t ny, const real_t *xs, const real_t *ys, + 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 nx, size_t ny, const real_t *xs, const real_t *ys, + 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_rotation.c b/pygrt/C_extension/src/dynamic/grt_rotation.c index 472d1717..c218df9d 100644 --- a/pygrt/C_extension/src/dynamic/grt_rotation.c +++ b/pygrt/C_extension/src/dynamic/grt_rotation.c @@ -51,6 +51,29 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ GRTCheckOptionSet(argc > 1); } +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) +{ + const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; + + for(size_t i=0; is_synpath, chs[c]); + insac = grt_read_SACTRACE(s_filepath, false); + u[c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + for(int c2=0; c2s_synpath, tolower(chs[c2]), chs[c]); + insac = grt_read_SACTRACE(s_filepath, false); + upar[c2][c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + res[c2][c] = calloc(npts, sizeof(*res[c2][c])); + } + } + grt_compute_rotation(npts, dist, u, upar, res, rot2ZNE); + + // 写出3个分量 for(int i1=0; i1<2; ++i1){ c1 = chs[i1]; for(int i2=i1+1; i2<3; ++i2){ c2 = chs[i2]; - - // 读取数据 u_{i,j} - GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c2), c1); - insac = grt_read_SACTRACE(s_filepath, false); - - // 累加 - for(int i=0; idata[i] += insac->data[i]; - - // 读取数据 u_{j,i} - GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c1), c2); - insac = grt_read_SACTRACE(s_filepath, false); - - // 累加 - for(int i=0; idata[i] = (outsac->data[i] - insac->data[i]) * 0.5f; - - // 特殊情况需加上协变导数，1e-5是因为km->cm - if(c1=='R' && c2=='T'){ - // 读取数据 u_T - GRT_SAFE_ASPRINTF(&s_filepath, "%s/T.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - for(int i=0; idata[i] -= 0.5f * insac->data[i] / dist * 1e-5; - } - - // 保存到SAC + memcpy(outsac->data, res[i2][i1], sizeof(*outsac->data)*npts); sprintf(outsac->hd.kcmpnm, "%c%c", c1, c2); GRT_SAFE_ASPRINTF(&s_filepath, "%s/rotation_%c%c.sac", Ctrl->s_synpath, c1, c2); grt_write_SACTRACE(s_filepath, outsac); - - // 置零 - for(int i=0; idata[i] = 0.0f; } } - grt_free_SACTRACE(insac); + for(int c=0; c 1); } +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) +{ + const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; + + for(size_t i=0; is_synpath, chs[c]); + insac = grt_read_SACTRACE(s_filepath, false); + u[c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + for(int c2=0; c2s_synpath, tolower(chs[c2]), chs[c]); + insac = grt_read_SACTRACE(s_filepath, false); + upar[c2][c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + res[c2][c] = calloc(npts, sizeof(*res[c2][c])); + } + } + grt_compute_strain(npts, dist, u, upar, res, rot2ZNE); + + // 写出6个分量 for(int i1=0; i1<3; ++i1){ c1 = chs[i1]; for(int i2=i1; i2<3; ++i2){ c2 = chs[i2]; - - // 读取数据 u_{i,j} - GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c2), c1); - insac = grt_read_SACTRACE(s_filepath, false); - - // 累加 - for(int i=0; idata[i] += insac->data[i]; - - // 读取数据 u_{j,i} - GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c1), c2); - insac = grt_read_SACTRACE(s_filepath, false); - - // 累加 - for(int i=0; idata[i] = (outsac->data[i] + insac->data[i]) * 0.5f; - - // 特殊情况需加上协变导数，1e-5是因为km->cm - if(c1=='R' && c2=='T'){ - // 读取数据 u_T - GRT_SAFE_ASPRINTF(&s_filepath, "%s/T.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - for(int i=0; idata[i] -= 0.5f * insac->data[i] / dist * 1e-5; - } - else if(c1=='T' && c2=='T'){ - // 读取数据 u_R - GRT_SAFE_ASPRINTF(&s_filepath, "%s/R.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - for(int i=0; idata[i] += insac->data[i] / dist * 1e-5; - } - - // 保存到SAC + memcpy(outsac->data, res[i2][i1], sizeof(*outsac->data)*npts); sprintf(outsac->hd.kcmpnm, "%c%c", c1, c2); GRT_SAFE_ASPRINTF(&s_filepath, "%s/strain_%c%c.sac", Ctrl->s_synpath, c1, c2); grt_write_SACTRACE(s_filepath, outsac); - - // 置零 - for(int i=0; idata[i] = 0.0f; } } - grt_free_SACTRACE(insac); + for(int c=0; c 1); } +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) +{ + const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; + size_t nf = npts/2 + 1; + float df = 1.0f/(npts*dt); + fftwf_complex *lam_ukk = fftwf_malloc(sizeof(*lam_ukk)*nf); + fftwf_complex *lams = fftwf_malloc(sizeof(*lams)*nf); + fftwf_complex *mus = fftwf_malloc(sizeof(*mus)*nf); + GRT_FFTWF_HOLDER *fwd = grt_create_fftwf_holder_R2C_1D(npts, dt, nf, df); + GRT_FFTWF_HOLDER *inv = grt_create_fftwf_holder_C2R_1D(npts, dt, nf, df); + + memset(lam_ukk, 0, sizeof(*lam_ukk)*nf); + for(size_t i=0; iw_t, upar[c][c], sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i]; + } + + // ZRT 联络项 u/r（1e-5: km→cm）；时域先算好再 FFT，避免频域再分支缩放 + float *ur_over_r = (float *)malloc(sizeof(float)*npts); + float *ut_over_r = (float *)malloc(sizeof(float)*npts); + for(size_t i=0; iw_t, ur_over_r, sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i]; + } + for(size_t i=0; iw_t, upar[c2][c], sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i] += fwd->W_f[i]; + + memcpy(fwd->w_t, upar[c][c2], sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i] = (inv->W_f[i] + fwd->W_f[i]) * mus[i]; + if(c == c2){ + for(size_t i=0; iW_f[i] += lam_ukk[i]; + } + if(chs[c]=='R' && chs[c2]=='T'){ + memcpy(fwd->w_t, ut_over_r, sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i] -= mus[i]*fwd->W_f[i]; + } + else if(chs[c]=='T' && chs[c2]=='T'){ + memcpy(fwd->w_t, ur_over_r, sizeof(float)*npts); + fftwf_execute(fwd->plan); + for(size_t i=0; iW_f[i] += 2.0f*mus[i]*fwd->W_f[i]; + } + fftwf_execute(inv->plan); + for(size_t i=0; iw_t[i]/npts; + grt_reset_fftwf_holder_zero(inv); + } + } + + GRT_SAFE_FREE_PTR(ur_over_r); + GRT_SAFE_FREE_PTR(ut_over_r); + + grt_destroy_fftwf_holder(fwd); + grt_destroy_fftwf_holder(inv); + fftwf_free(lam_ukk); + fftwf_free(lams); + fftwf_free(mus); +} + int stress_main(int argc, char **argv){ GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); @@ -87,8 +171,6 @@ int stress_main(int argc, char **argv){ int npts = insac->hd.npts; float dt = insac->hd.delta; float dist = insac->hd.dist; - float df = 1.0/(npts*dt); - int nf = npts/2 + 1; float va = insac->hd.user1; float vb = insac->hd.user2; float rho = insac->hd.user3; @@ -100,135 +182,48 @@ int stress_main(int argc, char **argv){ SACTRACE *outsac = grt_copy_SACTRACE(insac, true); grt_free_SACTRACE(insac); - // 申请内存 - // lamda * 体积应变 - fftwf_complex *lam_ukk = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); - // 不同频率的lambda和mu - fftwf_complex *lams = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); - fftwf_complex *mus = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); - // 分配FFTW - GRT_FFTWF_HOLDER *fwd_fftw_holder = grt_create_fftwf_holder_R2C_1D(npts, dt, nf, df); - GRT_FFTWF_HOLDER *inv_fftw_holder = grt_create_fftwf_holder_C2R_1D(npts, dt, nf, df); - // 初始化 - memset(lam_ukk, 0, sizeof(fftwf_complex)*nf); - memset(lams, 0, sizeof(fftwf_complex)*nf); - memset(mus, 0, sizeof(fftwf_complex)*nf); - // 计算不同频率下的拉梅系数 - for(int i=0; is_synpath, tolower(c1), c1); + float *u[GRT_CHANNEL_NUM]; + float *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + float *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + for(int c=0; cs_synpath, chs[c]); insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - - // 累加 - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i]; - } - // 加上协变导数 - if(!rot2ZNE){ - GRT_SAFE_ASPRINTF(&s_filepath, "%s/R.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i]/dist*1e-5; + u[c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + for(int c2=0; c2s_synpath, tolower(chs[c2]), chs[c]); + insac = grt_read_SACTRACE(s_filepath, false); + upar[c2][c] = insac->data; + insac->data = NULL; + grt_free_SACTRACE(insac); + res[c2][c] = calloc(npts, sizeof(*res[c2][c])); + } } + grt_compute_stress(npts, dt, dist, va, vb, rho, Qainv, Qbinv, u, upar, res, rot2ZNE); - // 乘上lambda系数 - for(int i=0; is_synpath, tolower(c2), c1); - insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - - // 累加 - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i] += fwd_fftw_holder->W_f[i]; - - // 读取数据 u_{j,i} - GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c1), c2); - insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - - // 累加 - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i] = (inv_fftw_holder->W_f[i] + fwd_fftw_holder->W_f[i]) * mus[i]; - - // 对于对角线分量，需加上lambda * u_kk - if(c1 == c2){ - for(int i=0; iW_f[i] += lam_ukk[i]; - } - - // 特殊情况需加上协变导数，1e-5是因为km->cm - if(c1=='R' && c2=='T'){ - // 读取数据 u_T - GRT_SAFE_ASPRINTF(&s_filepath, "%s/T.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i] -= mus[i] * fwd_fftw_holder->W_f[i] / dist * 1e-5; - } - else if(c1=='T' && c2=='T'){ - // 读取数据 u_R - GRT_SAFE_ASPRINTF(&s_filepath, "%s/R.sac", Ctrl->s_synpath); - insac = grt_read_SACTRACE(s_filepath, false); - memcpy(fwd_fftw_holder->w_t, insac->data, sizeof(float)*npts); - fftwf_execute(fwd_fftw_holder->plan); - for(int i=0; iW_f[i] += 2.0f * mus[i] * fwd_fftw_holder->W_f[i] / dist * 1e-5; - } - - // 保存到SAC - fftwf_execute(inv_fftw_holder->plan); - for(int i=0; iw_t[i] /= npts; - memcpy(outsac->data, inv_fftw_holder->w_t, sizeof(float)*npts); + memcpy(outsac->data, res[i2][i1], sizeof(*outsac->data)*npts); sprintf(outsac->hd.kcmpnm, "%c%c", c1, c2); GRT_SAFE_ASPRINTF(&s_filepath, "%s/stress_%c%c.sac", Ctrl->s_synpath, c1, c2); grt_write_SACTRACE(s_filepath, outsac); - - // 置零 - grt_reset_fftwf_holder_zero(inv_fftw_holder); } } - - grt_destroy_fftwf_holder(fwd_fftw_holder); - grt_destroy_fftwf_holder(inv_fftw_holder); - - GRT_SAFE_FFTW_FREE_PTR(lam_ukk, f); - GRT_SAFE_FFTW_FREE_PTR(lams, f); - GRT_SAFE_FFTW_FREE_PTR(mus, f); - - grt_free_SACTRACE(insac); + for(int c=0; ccm - if(chs[c]=='R' && chs[c2]=='T'){ - val -= 0.5 * u[2][ir] / dist * 1e-5; - } - - res[c2][c][ir] = val; - } - } - } - } + grt_static_compute_rotation(nx, ny, xs, ys, u, upar, res, rot2ZNE); // 写入 nc 文件 for(int c=0; ccm - if(chs[c]=='R' && chs[c2]=='T'){ - val -= 0.5 * u[2][ir] / dist * 1e-5; - } - else if(chs[c]=='T' && chs[c2]=='T'){ - val += u[1][ir] / dist * 1e-5; - } - - res[c2][c][ir] = val; - } - } - } - } + grt_static_compute_strain(nx, ny, xs, ys, u, upar, res, rot2ZNE); // 写入 nc 文件 for(int c=0; ccm - if(chs[c]=='R' && chs[c2]=='T'){ - val -= rcv_mu * u[2][ir] / dist * 1e-5; - } - else if(chs[c]=='T' && chs[c2]=='T'){ - val += 2.0 * rcv_mu * u[1][ir] / dist * 1e-5; - } - - res[c2][c][ir] = val; - } - } - } - } + grt_static_compute_stress(nx, ny, xs, ys, u, upar, res, rot2ZNE, rcv_mu, rcv_lam); // 写入 nc 文件 for(int c=0; c 0: chs = ZNEchs + npts = st_syn[0].stats.npts dist = st_syn[0].stats.sac['dist'] + u, upar, u_ptrs, upar_ptrs = _prepare_dynamic_postprocess_arrays(st_syn, chs, npts) + resarr, res_ptrs = _prepare_dynamic_postprocess_result(npts) + if Type == 'strain': + C_grt_compute_strain(npts, dist, u_ptrs, upar_ptrs, res_ptrs, chs == ZNEchs) + else: + C_grt_compute_rotation(npts, dist, u_ptrs, upar_ptrs, res_ptrs, chs == ZNEchs) - # ---------------------------------------------------------------------------------- - # 循环6/3个分量 stres = Stream() for i1 in range(i1_end): c1 = chs[i1] for i2 in range(i1+i2_offset, 3): c2 = chs[i2] - - channel = f"{c2.lower()}{c1}" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - tr = st[0].copy() - - channel = f"{c1.lower()}{c2}" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - tr.data = (tr.data + sgn*st[0].data) * 0.5 - - # 特殊情况加上协变导数 - if c1=='R' and c2=='T': - channel = f"T" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - tr.data -= 0.5*st[0].data / dist * 1e-5 - - elif c1=='T' and c2=='T': - channel = f"R" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - tr.data += st[0].data / dist * 1e-5 - - # 修改通道名 + tr = st_syn.select(channel=f"{c2.lower()}{c1}")[0].copy() + tr.data = resarr[i2, i1] tr.stats.channel = tr.stats.sac['kcmpnm'] = f"{c1}{c2}" - stres.append(tr) return stres +def _prepare_dynamic_postprocess_arrays(st_syn:Stream, chs:List[str], npts:int): + """收集动态位移/偏导数组，并构造 ctypes 通道指针表。""" + def data(channel:str): + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + if st[0].stats.npts != npts: + raise ValueError("All dynamic traces must have the same number of samples.") + return np.ascontiguousarray(st[0].data, dtype=np.float32) + + u = [data(c) for c in chs] + upar = [[data(f"{d.lower()}{c}") for c in chs] for d in chs] + u_ptrs = FLOAT_CHNL_PTRS(*(arr.ctypes.data_as(FPOINTER) for arr in u)) + upar_ptrs = FLOAT_CHNL_PTR_MAT(*( + FLOAT_CHNL_PTRS(*(arr.ctypes.data_as(FPOINTER) for arr in row)) for row in upar)) + return u, upar, u_ptrs, upar_ptrs + + +def _prepare_dynamic_postprocess_result(npts:int): + """分配动态后处理结果数组及其 ctypes 通道指针表。""" + resarr = np.zeros((CHANNEL_NUM, CHANNEL_NUM, npts), dtype=np.float32) + res_ptrs = FLOAT_CHNL_PTR_MAT(*( + FLOAT_CHNL_PTRS(*(resarr[c2, c1].ctypes.data_as(FPOINTER) for c1 in range(CHANNEL_NUM))) + for c2 in range(CHANNEL_NUM))) + return resarr, res_ptrs + + +def _prepare_static_postprocess_arrays(syn:dict, chs:List[str]): + """整理静态后处理所需的连续内存数组与 ctypes 通道指针表。""" + xarr = np.ascontiguousarray(syn['_xarr'], dtype=np.float64) + yarr = np.ascontiguousarray(syn['_yarr'], dtype=np.float64) + if xarr.ndim != 1 or yarr.ndim != 1: + raise ValueError("'_xarr' and '_yarr' must be one-dimensional arrays.") + + u = [np.ascontiguousarray(syn[c], dtype=np.float64) for c in chs] + upar = [ + [np.ascontiguousarray(syn[f"{d.lower()}{c}"], dtype=np.float64) for c in chs] + for d in chs + ] + expected_shape = (len(xarr), len(yarr)) + if any(arr.shape != expected_shape for arr in u) or any( + arr.shape != expected_shape for row in upar for arr in row): + raise ValueError("Static displacement and derivative arrays must match '_xarr'/'_yarr'.") + + u_ptrs = REAL_CHNL_PTRS(*(arr.ctypes.data_as(PREAL) for arr in u)) + upar_ptrs = REAL_CHNL_PTR_MAT(*( + REAL_CHNL_PTRS(*(arr.ctypes.data_as(PREAL) for arr in row)) for row in upar)) + resarr = np.zeros((CHANNEL_NUM, CHANNEL_NUM, *expected_shape), dtype=np.float64) + res_ptrs = REAL_CHNL_PTR_MAT(*( + REAL_CHNL_PTRS(*(resarr[c2, c1].ctypes.data_as(PREAL) for c1 in range(CHANNEL_NUM))) + for c2 in range(CHANNEL_NUM))) + return xarr, yarr, u, upar, u_ptrs, upar_ptrs, resarr, res_ptrs + + def _compute_static_strain_rotation(syn:dict, Type:str): r""" Compute static strain/rotation tensor from synthetic spatial derivatives. @@ -604,11 +633,9 @@ def _compute_static_strain_rotation(syn:dict, Type:str): """ if Type == 'strain': - sgn = 1 i1_end = 3 i2_offset = 0 elif Type == 'rotation': - sgn = -1 i1_end = 2 i2_offset = 1 else: @@ -620,8 +647,8 @@ def _compute_static_strain_rotation(syn:dict, Type:str): if f"nN" in syn.keys(): chs = ZNEchs - xarr:np.ndarray = syn['_xarr'] - yarr:np.ndarray = syn['_yarr'] + xarr, yarr, u, upar, u_ptrs, upar_ptrs, resarr, res_ptrs = \ + _prepare_static_postprocess_arrays(syn, chs) # 结果字典 resDct = {} @@ -632,48 +659,18 @@ def _compute_static_strain_rotation(syn:dict, Type:str): continue resDct[k] = deepcopy(syn[k]) - # 6/3个分量建立数组 - for i1 in range(i1_end): - c1 = chs[i1] - for i2 in range(i1+i2_offset, 3): - c2 = chs[i2] - channel = f"{c1}{c2}" - resDct[channel] = np.zeros((len(xarr), len(yarr)), dtype='f8') - - - for iy in range(len(yarr)): - for ix in range(len(xarr)): - # 震中距 - dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) - - # ---------------------------------------------------------------------------------- - # 循环6/3个分量 - for i1 in range(i1_end): - c1 = chs[i1] - for i2 in range(i1+i2_offset, 3): - c2 = chs[i2] - - channel = f"{c2.lower()}{c1}" - v12 = syn[channel][ix, iy] - - channel = f"{c1.lower()}{c2}" - v21 = syn[channel][ix, iy] - - val = 0.5*(v12 + sgn*v21) - - # 特殊情况加上协变导数 - if c1=='R' and c2=='T': - channel = f"T" - v0 = syn[channel][ix, iy] - val -= 0.5*v0 / dist * 1e-5 - - elif c1=='T' and c2=='T': - channel = f"R" - v0 = syn[channel][ix, iy] - val += v0 / dist * 1e-5 + if Type == 'strain': + C_grt_static_compute_strain( + len(xarr), len(yarr), xarr.ctypes.data_as(PREAL), yarr.ctypes.data_as(PREAL), + u_ptrs, upar_ptrs, res_ptrs, chs == ZNEchs) + else: + C_grt_static_compute_rotation( + len(xarr), len(yarr), xarr.ctypes.data_as(PREAL), yarr.ctypes.data_as(PREAL), + u_ptrs, upar_ptrs, res_ptrs, chs == ZNEchs) - channel = f"{c1}{c2}" - resDct[channel][ix, iy] = val + for i1 in range(i1_end): + for i2 in range(i1+i2_offset, CHANNEL_NUM): + resDct[f"{chs[i1]}{chs[i2]}"] = resarr[i2, i1] return resDct @@ -736,99 +733,27 @@ def _compute_stress(st_syn:Stream): nt = st_syn[0].stats.npts dt = st_syn[0].stats.delta dist = st_syn[0].stats.sac['dist'] - df = 1.0/(nt*dt) - nf = nt//2 + 1 va = st_syn[0].stats.sac['user1'] vb = st_syn[0].stats.sac['user2'] rho = st_syn[0].stats.sac['user3'] Qainv = st_syn[0].stats.sac['user4'] Qbinv = st_syn[0].stats.sac['user5'] - # 计算不同频率下的拉梅系数 - mus = np.zeros((nf,), dtype='c16') - lams = np.zeros((nf,), dtype='c16') - omega = CPLX(0.0, 0.0) - atte = CPLX(0.0, 0.0) - omgref = CPLX(2.0*np.pi*(nf-1)*df, 0.0) - for i in range(nf): - freq = 0.01 if i==0 else df*i - w = 2.0*np.pi*freq - omega.real = w - atte = C_grt_attenuation_law(Qbinv, omgref, omega) - attb = atte.real + atte.imag*1j - mus[i] = vb*vb*attb*attb*rho*1e10 - atte = C_grt_attenuation_law(Qainv, omgref, omega) - atta = atte.real + atte.imag*1j - lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i] - - del omega, atte - - # ---------------------------------------------------------------------------------- - # 先计算体积应变u_kk = u_11 + u22 + u33 和 lamda的乘积 - lam_ukk = np.zeros((nf,), dtype='c16') - for i in range(3): - c = chs[i] - channel = f"{c.lower()}{c}" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - lam_ukk[:] += rfft(st[0].data, nt) - - # 加上协变导数 - if not rot2ZNE: - channel = f"R" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - lam_ukk[:] += rfft(st[0].data, nt) / dist * 1e-5 + u, upar, u_ptrs, upar_ptrs = _prepare_dynamic_postprocess_arrays(st_syn, chs, nt) + resarr, res_ptrs = _prepare_dynamic_postprocess_result(nt) + C_grt_compute_stress( + nt, dt, dist, va, vb, rho, Qainv, Qbinv, + u_ptrs, upar_ptrs, res_ptrs, rot2ZNE) - lam_ukk[:] *= lams - - # ---------------------------------------------------------------------------------- - # 循环6个分量 stres = Stream() for i1 in range(3): c1 = chs[i1] for i2 in range(i1, 3): c2 = chs[i2] - channel = f"{c2.lower()}{c1}" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - tr = st[0].copy() - fftarr = np.zeros((nf,), dtype='c16') - - channel = f"{c1.lower()}{c2}" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - fftarr[:] = rfft(tr.data + st[0].data, nt) * mus - - # 对于对角线分量，需加上lambda * u_kk - if c1==c2: - fftarr[:] += lam_ukk - - # 特殊情况加上协变导数 - if c1=='R' and c2=='T': - channel = f"T" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - fftarr[:] -= mus*rfft(st[0].data, nt) / dist * 1e-5 - - elif c1=='T' and c2=='T': - channel = f"R" - st = st_syn.select(channel=channel) - if len(st) == 0: - raise NameError(f"{channel} not exists.") - fftarr[:] += 2.0*mus*rfft(st[0].data, nt) / dist * 1e-5 - - # 修改通道名 + tr = st_syn.select(channel=f"{c2.lower()}{c1}")[0].copy() + tr.data = resarr[i2, i1] tr.stats.channel = tr.stats.sac['kcmpnm'] = f"{c1}{c2}" - - tr.data = irfft(fftarr, nt) - stres.append(tr) return stres @@ -852,8 +777,8 @@ def _compute_static_stress(syn:dict): chs = ZNEchs rot2ZNE = True - xarr:np.ndarray = syn['_xarr'] - yarr:np.ndarray = syn['_yarr'] + xarr, yarr, u, upar, u_ptrs, upar_ptrs, resarr, res_ptrs = \ + _prepare_static_postprocess_arrays(syn, chs) va = syn['_rcv_va'] vb = syn['_rcv_vb'] rho = syn['_rcv_rho'] @@ -869,67 +794,13 @@ def _compute_static_stress(syn:dict): continue resDct[k] = deepcopy(syn[k]) - # 6个分量建立数组 - for i1 in range(3): - c1 = chs[i1] - for i2 in range(i1, 3): - c2 = chs[i2] - channel = f"{c1}{c2}" - resDct[channel] = np.zeros((len(xarr), len(yarr)), dtype='f8') - - - for iy in range(len(yarr)): - for ix in range(len(xarr)): - # 震中距 - dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) - - # ---------------------------------------------------------------------------------- - # 先计算体积应变u_kk = u_11 + u22 + u33 和 lamda的乘积 - lam_ukk = 0.0 - for i in range(3): - c = chs[i] - channel = f"{c.lower()}{c}" - lam_ukk += syn[channel][ix, iy] - - # 加上协变导数 - if not rot2ZNE: - channel = f"R" - lam_ukk += syn[channel][ix, iy] / dist * 1e-5 - - lam_ukk *= lam - - # ---------------------------------------------------------------------------------- - # 循环6个分量 - for i1 in range(3): - c1 = chs[i1] - for i2 in range(i1, 3): - c2 = chs[i2] - - channel = f"{c2.lower()}{c1}" - v12 = syn[channel][ix, iy] - - channel = f"{c1.lower()}{c2}" - v21 = syn[channel][ix, iy] + C_grt_static_compute_stress( + len(xarr), len(yarr), xarr.ctypes.data_as(PREAL), yarr.ctypes.data_as(PREAL), + u_ptrs, upar_ptrs, res_ptrs, rot2ZNE, mu, lam) - val = mu*(v12 + v21) - - # 对于对角线分量，需加上lambda * u_kk - if c1==c2: - val += lam_ukk - - # 特殊情况加上协变导数 - if c1=='R' and c2=='T': - channel = f"T" - v0 = syn[channel][ix, iy] - val -= mu*v0 / dist * 1e-5 - - elif c1=='T' and c2=='T': - channel = f"R" - v0 = syn[channel][ix, iy] - val += 2.0*mu*v0 / dist * 1e-5 - - channel = f"{c1}{c2}" - resDct[channel][ix, iy] = val + for i1 in range(CHANNEL_NUM): + for i2 in range(i1, CHANNEL_NUM): + resDct[f"{chs[i1]}{chs[i2]}"] = resarr[i2, i1] return resDct @@ -1602,4 +1473,3 @@ def solve_lamb1(nu:float, ts:np.ndarray, azimuth:float): C_grt_solve_lamb1(nu, npct.as_ctypes(ts), nt, azimuth, npct.as_ctypes(u.ravel())) return u -