Skip to content

FEAT: split the original kernel function into functions GRT_matrix and GRT_build_QWV #156

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
Dec 2, 2025
Merged

FEAT: split the original kernel function into functions GRT_matrix and GRT_build_QWV #156

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
16 changes: 16 additions & 0 deletions 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 @@ -92,7 +92,23 @@ void grt_kernel(
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]);

/** 构建广义反射透射系数矩阵。作为 kernel 函数中的第一部分 */
void grt_GRT_matrix(GRT_MODEL1D *mod1d, const real_t k);

/** 从广义 R/T 矩阵出发，计算每个震源对应的核函数 QWV。 作为 kernel 函数中的第二部分 */
void grt_GRT_build_QWV(
GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);

/** 静态解的核函数 */
void grt_static_kernel(
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]);

/** 静态广义反射透射系数矩阵 */
void grt_static_GRT_matrix(GRT_MODEL1D *mod1d, const real_t k);

/** 静态 QWV */
void grt_static_GRT_build_QWV(
GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
52 changes: 38 additions & 14 deletions pygrt/C_extension/src/common/kernel_template.c_
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@
#endif

#define __KERNEL_FUNC__ CONCAT_PREFIX(kernel)
#define __GRT_MATRIX_FUNC__ CONCAT_PREFIX(GRT_matrix)
#define __BUILD_QWV__ CONCAT_PREFIX(GRT_build_QWV)
#define __grt_RT_matrix_PSV CONCAT_PREFIX(RT_matrix_PSV)
#define __grt_RT_matrix_SH CONCAT_PREFIX(RT_matrix_SH)
#define __grt_delay_RT_matrix CONCAT_PREFIX(delay_RT_matrix)
Expand All @@ -52,14 +54,8 @@
#define __grt_wave2qwv_z_REV_SH CONCAT_PREFIX(wave2qwv_z_REV_SH)


void __KERNEL_FUNC__(
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])
void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k)
{
// 初始化qwv为0
memset(QWV, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM);
if(calc_uiz) memset(QWV_uiz, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM);

mod1d->k = k;
// 为动态解计算xa,xb,caca,cbcb
#if defined(__DYNAMIC_KERNEL__)
Expand Down Expand Up @@ -149,20 +145,38 @@ void __KERNEL_FUNC__(
} // END for loop
//===================================================================================

// return;
// 递推RU_FA
__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;

BEFORE_RETURN:
return;
}


void __BUILD_QWV__(
GRT_MODEL1D *mod1d, 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
memset(QWV, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM);
if(calc_uiz) memset(QWV_uiz, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM);

bool ircvup = mod1d->ircvup;
// 广义层的反射透射系数
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;

// 计算震源系数
__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);
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;

// 根据震源和台站相对位置，计算最终的系数
if(ircvup){ // A接收 B震源
Expand Down Expand Up @@ -256,12 +270,22 @@ void __KERNEL_FUNC__(
if(calc_uiz) QWV_uiz[GRT_SRC_M_DS_INDEX][c] = 0.0;
}
}
}


void __KERNEL_FUNC__(
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])
{
__GRT_MATRIX_FUNC__(mod1d, k);
__BUILD_QWV__(mod1d, QWV, calc_uiz, QWV_uiz);
}


#undef CONCAT_PREFIX
#undef __KERNEL_FUNC__
#undef __GRT_MATRIX_FUNC__
#undef __BUILD_QWV__
#undef __grt_RT_matrix_PSV
#undef __grt_RT_matrix_SH
#undef __grt_delay_RT_matrix
Expand Down