Skip to content

REFAC: move Generalized R/T matrix to struct GRT_MODEL1D #146

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Nov 27, 2025
Merged

REFAC: move Generalized R/T matrix to struct GRT_MODEL1D #146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions pygrt/C_extension/include/grt/common/RT_matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 上传反射系数矩阵
Expand Down
2 changes: 1 addition & 1 deletion pygrt/C_extension/include/grt/common/kernel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]);


Expand Down
90 changes: 45 additions & 45 deletions pygrt/C_extension/include/grt/common/kernel_template.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <stdio.h>
#include <complex.h>
Expand Down Expand Up @@ -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
Expand All @@ -167,33 +174,27 @@ 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
size_t imin, imax; // 相对浅层深层层位
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; iy<mod1d->n; ++iy){ // 因为n>=3, 故一定会进入该循环
Expand Down Expand Up @@ -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;
Expand All @@ -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，加入了自由表面的效应
Expand All @@ -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; i<GRT_SRC_M_NUM; ++i){
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, 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; i<GRT_SRC_M_NUM; ++i){
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, 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);
Expand All @@ -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; i<GRT_SRC_M_NUM; ++i){
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, 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; i<GRT_SRC_M_NUM; ++i){
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, 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]);
}
}

Expand Down
20 changes: 20 additions & 0 deletions pygrt/C_extension/include/grt/common/model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include <complex.h>
#include <stdbool.h>
#include "grt/common/const.h"
#include "grt/common/RT_matrix.h"


/** 1D 模型结构体，包括多个水平层，以及复数形式的弹性参数 */
Expand Down Expand Up @@ -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;


Expand Down
30 changes: 11 additions & 19 deletions pygrt/C_extension/include/grt/dynamic/layer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);


/**
Expand Down
9 changes: 4 additions & 5 deletions pygrt/C_extension/include/grt/dynamic/source.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
29 changes: 10 additions & 19 deletions pygrt/C_extension/include/grt/static/static_layer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);


/**
Expand Down
Loading