From 56f3266330c7d2f3fb9f6af8dfd6373708d4ffbf Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 27 Nov 2025 11:12:35 +0800 Subject: [PATCH] REFAC: move Generalized R/T matrix to struct `GRT_MODEL1D` --- .../include/grt/common/RT_matrix.h | 3 +- pygrt/C_extension/include/grt/common/kernel.h | 2 +- .../include/grt/common/kernel_template.h | 90 +++++++++---------- pygrt/C_extension/include/grt/common/model.h | 20 +++++ pygrt/C_extension/include/grt/dynamic/layer.h | 30 +++---- .../C_extension/include/grt/dynamic/source.h | 9 +- .../include/grt/static/static_layer.h | 29 +++--- .../include/grt/static/static_source.h | 9 +- pygrt/C_extension/src/common/dwm.c | 3 +- pygrt/C_extension/src/common/fim.c | 6 +- pygrt/C_extension/src/common/model.c | 22 ++--- pygrt/C_extension/src/common/ptam.c | 3 +- pygrt/C_extension/src/common/safim.c | 10 +-- pygrt/C_extension/src/dynamic/layer.c | 58 ++++++------ pygrt/C_extension/src/dynamic/source.c | 16 ++-- pygrt/C_extension/src/static/static_layer.c | 45 +++++----- pygrt/C_extension/src/static/static_source.c | 16 ++-- pygrt/c_structures.py | 71 +++++++-------- 18 files changed, 219 insertions(+), 223 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h index b4904b37..d92b01c8 100644 --- a/pygrt/C_extension/include/grt/common/RT_matrix.h +++ b/pygrt/C_extension/include/grt/common/RT_matrix.h @@ -13,8 +13,7 @@ #include "grt/common/const.h" - -/** 应该把 PSV 和 SH 的小矩阵统一放在一个结构体中管理 */ +/* 使用结构体管理上行和下行的 R/T 矩阵 */ typedef struct { cplx_t RD[2][2]; ///< P-SV 下传反射系数矩阵 cplx_t RU[2][2]; ///< P-SV 上传反射系数矩阵 diff --git a/pygrt/C_extension/include/grt/common/kernel.h b/pygrt/C_extension/include/grt/common/kernel.h index e62b2304..945b6638 100644 --- a/pygrt/C_extension/include/grt/common/kernel.h +++ b/pygrt/C_extension/include/grt/common/kernel.h @@ -17,7 +17,7 @@ * 计算核函数的函数指针，动态与静态的接口一致 */ typedef void (*GRT_KernelFunc) ( - const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + GRT_MODEL1D *mod1d, const real_t k, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]); diff --git a/pygrt/C_extension/include/grt/common/kernel_template.h b/pygrt/C_extension/include/grt/common/kernel_template.h index 6b9da811..37f4519d 100644 --- a/pygrt/C_extension/include/grt/common/kernel_template.h +++ b/pygrt/C_extension/include/grt/common/kernel_template.h @@ -14,9 +14,15 @@ // 注意这里不要使用 #pragma once + + // 只有定义了特定宏才可以被 include #if defined(__DYNAMIC_KERNEL__) || defined(__STATIC_KERNEL__) +// 不可以同时定义 +#if defined(__DYNAMIC_KERNEL__) && defined(__STATIC_KERNEL__) +#error "Cannot define both __DYNAMIC_KERNEL__ and __STATIC_KERNEL__ simultaneously" +#endif #include #include @@ -149,14 +155,15 @@ inline void grt_construct_qwv( * 即空间划分为FA,AB,BL, 计算这三个广义层的系数矩阵，再讨论震源层和接收层的深浅， * 计算相应的矩阵。 * - * @param[in] mod1d `MODEL1D` 结构体指针 + * @param[in,out] mod1d `MODEL1D` 结构体指针 + * @param[in] k 波数 * @param[out] QWV 不同震源，不同阶数的核函数 \f$ q_m, w_m, v_m \f$ * @param[in] calc_uiz 是否计算ui_z（位移u对坐标z的偏导） * @param[out] QWV_uiz 不同震源，不同阶数的核函数对z的偏导 \f$ \frac{\partial q_m}{\partial z}, \frac{\partial w_m}{\partial z}, \frac{\partial v_m}{\partial z} \f$ * */ static void __KERNEL_FUNC__( - const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], + GRT_MODEL1D *mod1d, const real_t k, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM], bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]) { // 初始化qwv为0 @@ -167,6 +174,12 @@ static void __KERNEL_FUNC__( } } + mod1d->k = k; + // 为动态解计算xa,xb,caca,cbcb + #if defined(__DYNAMIC_KERNEL__) + grt_mod1d_xa_xb(mod1d, k); + #endif + bool ircvup = mod1d->ircvup; size_t isrc = mod1d->isrc; // 震源所在虚拟层位, isrc>=1 size_t ircv = mod1d->ircv; // 接收点所在虚拟层位, ircv>=1, ircv != isrc @@ -174,26 +187,14 @@ static void __KERNEL_FUNC__( imin = GRT_MIN(isrc, ircv); imax = GRT_MAX(isrc, ircv); - // 初始化广义反射透射系数矩阵 - // BL - grt_init_RT_matrix(M_BL); - // AL - grt_init_RT_matrix(M_AL); - // RS - grt_init_RT_matrix(M_RS); - // FA (实际先计算ZA，再递推到FA) - grt_init_RT_matrix(M_FA); - // FB (实际先计算ZB，再递推到FB) - grt_init_RT_matrix(M_FB); - - // 自由表面的反射系数 - grt_init_RT_matrix(M_top); - - // 接收点处的接收矩阵(转为位移u的(B_m, C_m, P_m)系分量) - cplx_t R_EV[2][2], R_EVL; - - // 接收点处的接收矩阵(转为位移导数ui_z的(B_m, C_m, P_m)系分量) - cplx_t uiz_R_EV[2][2], uiz_R_EVL; + // 广义层的反射透射系数 + RT_MATRIX *M_BL = &mod1d->M_BL; + RT_MATRIX *M_AL = &mod1d->M_AL; + RT_MATRIX *M_RS = &mod1d->M_RS; + RT_MATRIX *M_FA = &mod1d->M_FA; + RT_MATRIX *M_FB = &mod1d->M_FB; + // 自由表面 + RT_MATRIX *M_top = &mod1d->M_top; // 从顶到底进行矩阵递推, 公式(5.5.3) for(size_t iy=1; iyn; ++iy){ // 因为n>=3, 故一定会进入该循环 @@ -254,14 +255,13 @@ static void __KERNEL_FUNC__( // 计算震源系数 - cplx_t src_coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2] = {0}; - __grt_source_coef(mod1d, src_coef); + __grt_source_coef(mod1d); // 临时中转矩阵 (temperary) cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2; // 递推RU_FA - __grt_topfree_RU(mod1d, M_top); + __grt_topfree_RU(mod1d); if(M_top->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; @@ -270,8 +270,8 @@ static void __KERNEL_FUNC__( if(ircvup){ // A接收 B震源 // 计算R_EV - __grt_wave2qwv_REV_PSV(mod1d, M_FA->RU, R_EV); - __grt_wave2qwv_REV_SH(mod1d, M_FA->RUL, &R_EVL); + __grt_wave2qwv_REV_PSV(mod1d); + __grt_wave2qwv_REV_SH(mod1d); // 递推RU_FS grt_recursion_RU(M_FA, M_RS, M_FB); // 已从ZR变为FR，加入了自由表面的效应 @@ -285,31 +285,31 @@ static void __KERNEL_FUNC__( if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - grt_cmat2x2_mul(R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(mod1d->R_EV, tmp2x2, tmp2x2); tmpRL = M_FB->invTL / (1.0 - M_BL->RDL * M_FB->RUL); - tmpRL2 = R_EVL * tmpRL; + tmpRL2 = mod1d->R_EVL * tmpRL; for(int i=0; iRD, M_BL->RDL, src_coef[i], QWV[i]); + grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, M_BL->RDL, mod1d->src_coef[i], QWV[i]); } if(calc_uiz){ - __grt_wave2qwv_z_REV_PSV(mod1d, M_FA->RU, uiz_R_EV); - __grt_wave2qwv_z_REV_SH(mod1d, M_FA->RUL, &uiz_R_EVL); - grt_cmat2x2_mul(uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); - tmpRL2 = uiz_R_EVL * tmpRL; + __grt_wave2qwv_z_REV_PSV(mod1d); + __grt_wave2qwv_z_REV_SH(mod1d); + grt_cmat2x2_mul(mod1d->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL2 = mod1d->uiz_R_EVL * tmpRL; for(int i=0; iRD, M_BL->RDL, src_coef[i], QWV_uiz[i]); + grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, M_BL->RDL, mod1d->src_coef[i], QWV_uiz[i]); } } } else { // A震源 B接收 // 计算R_EV - __grt_wave2qwv_REV_PSV(mod1d, M_BL->RD, R_EV); - __grt_wave2qwv_REV_SH(mod1d, M_BL->RDL, &R_EVL); + __grt_wave2qwv_REV_PSV(mod1d); + __grt_wave2qwv_REV_SH(mod1d); // 递推RD_SL grt_recursion_RD(M_RS, M_BL, M_AL); @@ -323,23 +323,23 @@ static void __KERNEL_FUNC__( if(calc_uiz) grt_cmat2x2_assign(tmp2x2, tmp2x2_uiz); // 为后续计算空间导数备份 - grt_cmat2x2_mul(R_EV, tmp2x2, tmp2x2); + grt_cmat2x2_mul(mod1d->R_EV, tmp2x2, tmp2x2); tmpRL = M_AL->invTL / (1.0 - M_FA->RUL * M_AL->RDL); - tmpRL2 = R_EVL * tmpRL; + tmpRL2 = mod1d->R_EVL * tmpRL; for(int i=0; iRU, M_FA->RUL, src_coef[i], QWV[i]); + grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, M_FA->RUL, mod1d->src_coef[i], QWV[i]); } if(calc_uiz){ - __grt_wave2qwv_z_REV_PSV(mod1d, M_BL->RD, uiz_R_EV); - __grt_wave2qwv_z_REV_SH(mod1d, M_BL->RDL, &uiz_R_EVL); - grt_cmat2x2_mul(uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); - tmpRL2 = uiz_R_EVL * tmpRL; + __grt_wave2qwv_z_REV_PSV(mod1d); + __grt_wave2qwv_z_REV_SH(mod1d); + grt_cmat2x2_mul(mod1d->uiz_R_EV, tmp2x2_uiz, tmp2x2_uiz); + tmpRL2 = mod1d->uiz_R_EVL * tmpRL; for(int i=0; iRU, M_FA->RUL, src_coef[i], QWV_uiz[i]); + grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, M_FA->RUL, mod1d->src_coef[i], QWV_uiz[i]); } } diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 2e539ed2..dc168532 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -11,6 +11,7 @@ #include #include #include "grt/common/const.h" +#include "grt/common/RT_matrix.h" /** 1D 模型结构体，包括多个水平层，以及复数形式的弹性参数 */ @@ -50,6 +51,25 @@ typedef struct { /* 状态变量，非 0 为异常值 */ int stats; + /* 根据震源和台站划分的广义 R/T 系数 */ + RT_MATRIX M_AL; + RT_MATRIX M_BL; + RT_MATRIX M_RS; + RT_MATRIX M_FA; + RT_MATRIX M_FB; + + /* 自由表面的反射系数矩阵，仅 RU, RUL 有用 */ + RT_MATRIX M_top; + + /* 接收点处的接收矩阵 (转为位移u和位移导数uiz的(B_m, C_m, P_m)系分量) */ + cplx_t R_EV[2][2]; + cplx_t R_EVL; + cplx_t uiz_R_EV[2][2]; + cplx_t uiz_R_EVL; + + /* 震源处的震源系数 */ + cplx_t src_coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]; ///< 震源系数 \f$ P_m, SV_m, SH_m \f$ + } GRT_MODEL1D; diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index b2f6669a..0ff8b7c5 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -20,56 +20,48 @@ /** * 计算自由表面的反射系数，公式(5.3.10-14) * - * @param[in] mod1d 模型结构体指针 - * @param[out] M R/T矩阵，仅填充RU + * + * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop * */ -void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M); +void grt_topfree_RU(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 P-SV 波接收矩阵，将波场转为位移，公式(5.2.19) + (5.7.7,25) * - * @param[in] mod1d 模型结构体指针 - * @param[in] R P-SV波场 - * @param[out] R_EV P-SV接收函数矩阵 + * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV * */ -void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, cplx_t R[2][2], cplx_t R_EV[2][2]); +void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 SH 波接收矩阵，将波场转为位移，公式(5.2.19) + (5.7.7,25) * - * @param[in] mod1d 模型结构体指针 - * @param[in] RL SH波场 - * @param[out] R_EVL SH接收函数值 + * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL * */ -void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); +void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 P-SV 波接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in] mod1d 模型结构体指针 - * @param[in] R P-SV波场 - * @param[out] R_EV P-SV接收函数矩阵 + * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); +void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 SH 波接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in] mod1d 模型结构体指针 - * @param[in] RL SH波场 - * @param[out] R_EVL SH接收函数值 + * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); +void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); /** diff --git a/pygrt/C_extension/include/grt/dynamic/source.h b/pygrt/C_extension/include/grt/dynamic/source.h index f0ff9520..f2d471d9 100755 --- a/pygrt/C_extension/include/grt/dynamic/source.h +++ b/pygrt/C_extension/include/grt/dynamic/source.h @@ -20,14 +20,13 @@ * 数组形状代表在[i][j][p]时表示i类震源的 * P(j=0),SV(j=1)的震源系数(分别对应q,w)，且分为下行波(p=0)和上行波(p=1). * - * @param[in] mod1d 模型结构体指针 - * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ + * @param[in,out] mod1d 模型结构体指针，结果保存在 src_coef * */ -void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); +void grt_source_coef(GRT_MODEL1D *mod1d); /* P-SV 波的震源系数 */ -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); +void grt_source_coef_PSV(GRT_MODEL1D *mod1d); /* SH 波的震源系数 */ -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); +void grt_source_coef_SH(GRT_MODEL1D *mod1d); diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index afd69132..9115697f 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -22,54 +22,45 @@ /** * 计算自由表面的静态反射系数，公式(6.3.12) * - * @param[in] mod1d 模型结构体指针 - * @param[out] R_tilt R/T矩阵，仅填充RU + * @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop * */ -void grt_static_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M); +void grt_static_topfree_RU(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35) * - * @param[in] mod1d 模型结构体指针 - * @param[in] R P-SV波场 - * @param[out] R_EV P-SV接收函数矩阵 + * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EV * */ -void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); +void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d); /** * 计算接收点位置的 SH 静态接收矩阵，将波场转为位移，公式(6.3.37) * - * @param[in] mod1d 模型结构体指针 - * @param[in] RL SH波场 - * @param[out] R_EVL SH接收函数值 + * @param[in,out] mod1d 模型结构体指针，结果保存在 R_EVL * */ -void grt_static_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); +void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 P-SV 静态接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in] mod1d 模型结构体指针 - * @param[in] R P-SV波场 - * @param[out] R_EV P-SV接收函数矩阵 + * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EV * */ -void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]); +void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d); /** * 计算接收点位置的ui_z的 SH 静态接收矩阵，即将波场转为ui_z。 * 公式本质是推导ui_z关于q_m, w_m, v_m的连接矩阵（就是应力推导过程的一部分） * - * @param[in] mod1d 模型结构体指针 - * @param[in] RL SH波场 - * @param[out] R_EVL SH接收函数值 + * @param[in,out] mod1d 模型结构体指针，结果保存在 uiz_R_EVL * */ -void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL); +void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d); /** diff --git a/pygrt/C_extension/include/grt/static/static_source.h b/pygrt/C_extension/include/grt/static/static_source.h index ef591cd6..13172ca5 100644 --- a/pygrt/C_extension/include/grt/static/static_source.h +++ b/pygrt/C_extension/include/grt/static/static_source.h @@ -19,14 +19,13 @@ * 数组形状代表在[i][j][p]时表示i类震源的 * P(j=0),SV(j=1)的震源系数(分别对应q,w)，且分为下行波(p=0)和上行波(p=1). * - * @param[in] mod1d 模型结构体指针 - * @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$ + * @param[in,out] mod1d 模型结构体指针，结果保存在 src_coef */ -void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); +void grt_static_source_coef(GRT_MODEL1D *mod1d); /* P-SV 波的静态震源系数 */ -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); +void grt_static_source_coef_PSV(GRT_MODEL1D *mod1d); /* SH 波的静态震源系数 */ -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]); \ No newline at end of file +void grt_static_source_coef_SH(GRT_MODEL1D *mod1d); \ No newline at end of file diff --git a/pygrt/C_extension/src/common/dwm.c b/pygrt/C_extension/src/common/dwm.c index 5b517b91..20a88d36 100644 --- a/pygrt/C_extension/src/common/dwm.c +++ b/pygrt/C_extension/src/common/dwm.c @@ -55,8 +55,7 @@ real_t grt_discrete_integ( k += dk; // 计算核函数 F(k, w) - grt_mod1d_xa_xb(mod1d, k); - kerfunc(mod1d, QWV, calc_upar, QWV_uiz); + kerfunc(mod1d, k, QWV, calc_upar, QWV_uiz); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 记录积分核函数 diff --git a/pygrt/C_extension/src/common/fim.c b/pygrt/C_extension/src/common/fim.c index 3d9f6de5..4f0834dc 100755 --- a/pygrt/C_extension/src/common/fim.c +++ b/pygrt/C_extension/src/common/fim.c @@ -60,8 +60,7 @@ real_t grt_linear_filon_integ( k += dk; // 计算核函数 F(k, w) - grt_mod1d_xa_xb(mod1d, k); - kerfunc(mod1d, QWV, calc_upar, QWV_uiz); + kerfunc(mod1d, k, QWV, calc_upar, QWV_uiz); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 记录积分结果 @@ -175,8 +174,7 @@ real_t grt_linear_filon_integ( } // 计算核函数 F(k, w) - grt_mod1d_xa_xb(mod1d, k0N); - kerfunc(mod1d, QWV, calc_upar, QWV_uiz); + kerfunc(mod1d, k0N, QWV, calc_upar, QWV_uiz); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; for(size_t ir=0; irn); - size_t n = mod1d1->n; - mod1d2->n = mod1d1->n; - mod1d2->depsrc = mod1d1->depsrc; - mod1d2->deprcv = mod1d1->deprcv; - mod1d2->isrc = mod1d1->isrc; - mod1d2->ircv = mod1d1->ircv; - mod1d2->ircvup = mod1d1->ircvup; - mod1d2->io_depth = mod1d1->io_depth; + GRT_MODEL1D *mod1d2 = (GRT_MODEL1D *)calloc(1, sizeof(GRT_MODEL1D)); - mod1d2->omega = mod1d1->omega; - mod1d2->k = mod1d1->k; - mod1d2->c_phase = mod1d1->c_phase; + // 先直接赋值，实现浅拷贝 + *mod1d2 = *mod1d1; - mod1d2->stats = mod1d1->stats; + // 对指针部分再重新申请内存并赋值，实现深拷贝 + size_t n = mod1d1->n; + #define X(P, T) \ + mod1d2->P = (T*)calloc(n, sizeof(T));\ + memcpy(mod1d2->P, mod1d1->P, sizeof(T)*n);\ - #define X(P, T) memcpy(mod1d2->P, mod1d1->P, sizeof(T)*n); GRT_FOR_EACH_MODEL_QUANTITY_ARRAY #undef X diff --git a/pygrt/C_extension/src/common/ptam.c b/pygrt/C_extension/src/common/ptam.c index e01100d8..cbabf528 100755 --- a/pygrt/C_extension/src/common/ptam.c +++ b/pygrt/C_extension/src/common/ptam.c @@ -217,8 +217,7 @@ void grt_PTA_method( k += dk; // 计算核函数 F(k, w) - grt_mod1d_xa_xb(mod1d, k); - kerfunc(mod1d, QWV, calc_upar, QWV_uiz); + kerfunc(mod1d, k, QWV, calc_upar, QWV_uiz); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; // 记录核函数 diff --git a/pygrt/C_extension/src/common/safim.c b/pygrt/C_extension/src/common/safim.c index 3bcfb159..89906245 100644 --- a/pygrt/C_extension/src/common/safim.c +++ b/pygrt/C_extension/src/common/safim.c @@ -323,8 +323,7 @@ real_t grt_sa_filon_integ( .F3_uiz = {{{0}}} }; for(int i=0; i<3; ++i) { - grt_mod1d_xa_xb(mod1d, Kitv.k3[i]); - kerfunc(mod1d, Kitv.F3[i], calc_upar, Kitv.F3_uiz[i]); + kerfunc(mod1d, Kitv.k3[i], Kitv.F3[i], calc_upar, Kitv.F3_uiz[i]); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; } stack_push(&stack, Kitv); @@ -371,12 +370,11 @@ real_t grt_sa_filon_integ( memcpy(Kitv_right.F3_uiz[0], Kitv.F3_uiz[1], sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM); memcpy(Kitv_right.F3_uiz[2], Kitv.F3_uiz[2], sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM); } - grt_mod1d_xa_xb(mod1d, Kitv_left.k3[1]); - kerfunc(mod1d, Kitv_left.F3[1], calc_upar, Kitv_left.F3_uiz[1]); + + kerfunc(mod1d, Kitv_left.k3[1], Kitv_left.F3[1], calc_upar, Kitv_left.F3_uiz[1]); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; - grt_mod1d_xa_xb(mod1d, Kitv_right.k3[1]); - kerfunc(mod1d, Kitv_right.F3[1], calc_upar, Kitv_right.F3_uiz[1]); + kerfunc(mod1d, Kitv_right.k3[1], Kitv_right.F3[1], calc_upar, Kitv_right.F3_uiz[1]); if(mod1d->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index 62f20264..af6d915f 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -29,7 +29,7 @@ T N##2 = mod1d->N[iy];\ -void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) +void grt_topfree_RU(GRT_MODEL1D *mod1d) { cplx_t cbcb0 = mod1d->cbcb[0]; if(cbcb0 != 0.0){ @@ -44,26 +44,26 @@ void grt_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 Delta = -1.0 + xa0*xb0 + cbcb0 - cbcb02; if(Delta == 0.0){ - M->stats = GRT_INVERSE_FAILURE; + mod1d->M_top.stats = GRT_INVERSE_FAILURE; return; } Delta = 1.0 / Delta; - M->RU[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; - M->RU[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta; - M->RU[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta; - M->RU[1][1] = M->RU[0][0]; - M->RUL = 1.0; + mod1d->M_top.RU[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; + mod1d->M_top.RU[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta; + mod1d->M_top.RU[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta; + mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][0]; + mod1d->M_top.RUL = 1.0; } else { // 液体表面 - M->RU[0][0] = -1.0; - M->RU[1][1] = M->RU[0][1] = M->RU[1][0] = 0.0; - M->RUL = 0.0; + mod1d->M_top.RU[0][0] = -1.0; + mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][1] = mod1d->M_top.RU[1][0] = 0.0; + mod1d->M_top.RUL = 0.0; } } -void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, cplx_t R[2][2], cplx_t R_EV[2][2]) +void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; @@ -88,16 +88,16 @@ void grt_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, cplx_t R[2][2], cplx_t R_EV[ // 公式(5.7.7,25) if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, R, R_EV); - grt_cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->R_EV); + grt_cmat2x2_add(D11, mod1d->R_EV, mod1d->R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, R, R_EV); - grt_cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->R_EV); + grt_cmat2x2_add(D12, mod1d->R_EV, mod1d->R_EV); } } -void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) +void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xb = mod1d->xb[ircv]; @@ -106,15 +106,19 @@ void grt_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) if(xb != 1.0){ // 位于固体层 // 公式(5.2.19) - *R_EVL = (1.0 + (RL))*k; + if(mod1d->ircvup){// 震源更深 + mod1d->R_EVL = (1.0 + mod1d->M_FA.RUL)*k; + } else { + mod1d->R_EVL = (1.0 + mod1d->M_BL.RDL)*k; + } } else { // 位于液体层 - *R_EVL = 0.0; + mod1d->R_EVL = 0.0; } } -void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]) +void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; @@ -138,16 +142,16 @@ void grt_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx // 公式(5.7.7,25) if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, R, R_EV); - grt_cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->uiz_R_EV); + grt_cmat2x2_add(D11, mod1d->uiz_R_EV, mod1d->uiz_R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, R, R_EV); - grt_cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->uiz_R_EV); + grt_cmat2x2_add(D12, mod1d->uiz_R_EV, mod1d->uiz_R_EV); } } -void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) +void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xb = mod1d->xb[ircv]; @@ -160,13 +164,13 @@ void grt_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) if(xb != 1.0){ // 位于固体层 if(mod1d->ircvup){// 震源更深 - *R_EVL = (1.0 - (RL))*bk; + mod1d->uiz_R_EVL = (1.0 - mod1d->M_FA.RUL)*bk; } else { // 接收点更深 - *R_EVL = (RL - 1.0)*bk; + mod1d->uiz_R_EVL = (mod1d->M_BL.RDL - 1.0)*bk; } } else { // 位于液体层 - *R_EVL = 0.0; + mod1d->uiz_R_EVL = 0.0; } } diff --git a/pygrt/C_extension/src/dynamic/source.c b/pygrt/C_extension/src/dynamic/source.c index 4d96f068..531db0bf 100755 --- a/pygrt/C_extension/src/dynamic/source.c +++ b/pygrt/C_extension/src/dynamic/source.c @@ -63,17 +63,17 @@ inline void _source_SH(const cplx_t xb, const cplx_t cbcb, const real_t k, cplx_ } -void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_source_coef(GRT_MODEL1D *mod1d) { // 先全部赋0 - memset(coef, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM*2); + memset(mod1d->src_coef, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM*2); - grt_source_coef_PSV(mod1d, coef); - grt_source_coef_SH(mod1d, coef); + grt_source_coef_PSV(mod1d); + grt_source_coef_SH(mod1d); } -void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_source_coef_PSV(GRT_MODEL1D *mod1d) { size_t isrc = mod1d->isrc; cplx_t xa = mod1d->xa[isrc]; @@ -82,18 +82,18 @@ void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GR cplx_t cbcb = mod1d->cbcb[isrc]; real_t k = mod1d->k; - _source_PSV(xa, caca, xb, cbcb, k, coef); + _source_PSV(xa, caca, xb, cbcb, k, mod1d->src_coef); } -void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_source_coef_SH(GRT_MODEL1D *mod1d) { size_t isrc = mod1d->isrc; cplx_t xb = mod1d->xb[isrc]; cplx_t cbcb = mod1d->cbcb[isrc]; real_t k = mod1d->k; - _source_SH(xb, cbcb, k, coef); + _source_SH(xb, cbcb, k, mod1d->src_coef); } diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index a8779c63..687f9272 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -25,38 +25,41 @@ T N##2 = mod1d->N[iy];\ -void grt_static_topfree_RU(const GRT_MODEL1D *mod1d, RT_MATRIX *M) +void grt_static_topfree_RU(GRT_MODEL1D *mod1d) { cplx_t delta1 = mod1d->delta[0]; // 公式(6.3.12) - M->RU[0][0] = M->RU[1][1] = 0.0; - M->RU[0][1] = -delta1; - M->RU[1][0] = -1.0/delta1; - M->RUL = 1.0; + mod1d->M_top.RU[0][0] = mod1d->M_top.RU[1][1] = 0.0; + mod1d->M_top.RU[0][1] = -delta1; + mod1d->M_top.RU[1][0] = -1.0/delta1; + mod1d->M_top.RUL = 1.0; } -void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]) +void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) { cplx_t D11[2][2] = {{1.0, -1.0}, {1.0, 1.0}}; cplx_t D12[2][2] = {{1.0, -1.0}, {-1.0, -1.0}}; // 公式(6.3.35,37) if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, R, R_EV); - grt_cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->R_EV); + grt_cmat2x2_add(D11, mod1d->R_EV, mod1d->R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, R, R_EV); - grt_cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->R_EV); + grt_cmat2x2_add(D12, mod1d->R_EV, mod1d->R_EV); } } -void grt_static_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) +void grt_static_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) { - (void)mod1d; // 暂停该参数 -Wunused-parameter 警告 - *R_EVL = (1.0 + (RL)); + if(mod1d->ircvup){// 震源更深 + mod1d->R_EVL = 1.0 + mod1d->M_FA.RUL; + } else { + mod1d->R_EVL = 1.0 + mod1d->M_BL.RDL; + } } -void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2], cplx_t R_EV[2][2]) +void grt_static_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) { real_t k = mod1d->k; size_t ircv = mod1d->ircv; @@ -67,22 +70,22 @@ void grt_static_wave2qwv_z_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2 cplx_t D11[2][2] = {{k, -k-kd2}, {k, k-kd2}}; cplx_t D12[2][2] = {{-k, k+kd2}, {k, k-kd2}}; if(mod1d->ircvup){// 震源更深 - grt_cmat2x2_mul(D12, R, R_EV); - grt_cmat2x2_add(D11, R_EV, R_EV); + grt_cmat2x2_mul(D12, mod1d->M_FA.RU, mod1d->uiz_R_EV); + grt_cmat2x2_add(D11, mod1d->uiz_R_EV, mod1d->uiz_R_EV); } else { // 接收点更深 - grt_cmat2x2_mul(D11, R, R_EV); - grt_cmat2x2_add(D12, R_EV, R_EV); + grt_cmat2x2_mul(D11, mod1d->M_BL.RD, mod1d->uiz_R_EV); + grt_cmat2x2_add(D12, mod1d->uiz_R_EV, mod1d->uiz_R_EV); } } -void grt_static_wave2qwv_z_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL) +void grt_static_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) { real_t k = mod1d->k; // 新推导公式 if(mod1d->ircvup){// 震源更深 - *R_EVL = (1.0 - (RL))*k; + mod1d->uiz_R_EVL = (1.0 - mod1d->M_FA.RUL)*k; } else { // 接收点更深 - *R_EVL = (RL - 1.0)*k; + mod1d->uiz_R_EVL = (mod1d->M_BL.RDL - 1.0)*k; } } diff --git a/pygrt/C_extension/src/static/static_source.c b/pygrt/C_extension/src/static/static_source.c index a1136b72..f6bb30e9 100644 --- a/pygrt/C_extension/src/static/static_source.c +++ b/pygrt/C_extension/src/static/static_source.c @@ -61,31 +61,31 @@ inline void _source_SH(const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2 } -void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_static_source_coef(GRT_MODEL1D *mod1d) { // 先全部赋0 - memset(coef, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM*2); + memset(mod1d->src_coef, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM*2); - grt_static_source_coef_PSV(mod1d, coef); - grt_static_source_coef_SH(mod1d, coef); + grt_static_source_coef_PSV(mod1d); + grt_static_source_coef_SH(mod1d); } -void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_static_source_coef_PSV(GRT_MODEL1D *mod1d) { size_t isrc = mod1d->isrc; cplx_t delta = mod1d->delta[isrc]; real_t k = mod1d->k; - _source_PSV(delta, k, coef); + _source_PSV(delta, k, mod1d->src_coef); } -void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]) +void grt_static_source_coef_SH(GRT_MODEL1D *mod1d) { real_t k = mod1d->k; - _source_SH(k, coef); + _source_SH(k, mod1d->src_coef); } diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index 191dbb5d..71f19dec 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -55,45 +55,31 @@ PREAL = POINTER(REAL) PCPLX = POINTER(CPLX) -class c_GRT_MODEL1D(Structure): +class c_RT_MATRIX(Structure): + """ + 和C结构体 RT_MATRIX 匹配 """ - 和C结构体 GRT_MODEL1D 作匹配 - :field n: 层数 - :filed depsrc: 震源深度 km - :filed deprcv: 接收点深度 km - :field isrc: 震源所在层位 - :field ircv: 台站所在层位 - :field ircvup: 台站层位是否高于震源 - :field io_depth: 模型读入的第一列是否为每层顶界面深度 - - :field omega: 圆频率 - :field omega: 波数 - :field c_phase: 相速度 - - :field thk: 数组, 每层层厚(km) - :field dep: 数组, 每层顶界面深度(km) - :field Va: 数组, 每层P波速度(km/s) - :field Vb: 数组, 每层S波速度(km/s) - :field Rho: 数组, 每层密度(g/cm^3) - :field Qa: 数组, 每层P波品质因子Q_P - :field Qb: 数组, 每层S波品质因子Q_S - :field Qainv: - :field Qbinv: - - :field mu: - :field lambda: - :field delta: - :field atna: - :field atnb: - :field xa: - :field xb: - :field caca: - :field cbcb: - - :field stats: + _fields_ = [ + ('RD', CPLX * 4), + ('RU', CPLX * 4), + ('TD', CPLX * 4), + ('TU', CPLX * 4), + ('RDL', CPLX), + ('RUL', CPLX), + ('TDL', CPLX), + ('TUL', CPLX), + ('invT', CPLX * 4), + ('invTL', CPLX), + ('stats', c_int) + ] + +class c_GRT_MODEL1D(Structure): """ + 和C结构体 GRT_MODEL1D 作匹配 + """ + _fields_ = [ ('n', c_size_t), ("depsrc", REAL), @@ -128,4 +114,19 @@ class c_GRT_MODEL1D(Structure): ('cbcb', PCPLX), ('stats', c_int), + + ('M_AL', c_RT_MATRIX), + ('M_BL', c_RT_MATRIX), + ('M_RS', c_RT_MATRIX), + ('M_FA', c_RT_MATRIX), + ('M_FB', c_RT_MATRIX), + + ('M_top', c_RT_MATRIX), + + ('R_EV', CPLX * 4), + ('R_EVL', CPLX), + ('uiz_R_EV', CPLX * 4), + ('uiz_R_EVL', CPLX), + + ('src_coef', CPLX * SRC_M_NUM * QWV_NUM * 2) ]