Skip to content

REFAC: split some kernel-related functions into P-SV and SH parts #228

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
Jan 30, 2026
Merged

REFAC: split some kernel-related functions into P-SV and SH parts #228

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
2 changes: 2 additions & 0 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 @@ -39,6 +39,8 @@ typedef struct {
* @param[out] M RT_MATRIX 指针
*/
void grt_reset_RT_matrix(RT_MATRIX *M);
void grt_reset_RT_matrix_PSV(RT_MATRIX *M);
void grt_reset_RT_matrix_SH(RT_MATRIX *M);

/**
* 合并 recursion_RD_PSV(SH) ，仅计算RD/RDL
Expand Down
8 changes: 6 additions & 2 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 @@ -24,7 +24,8 @@
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_top
*
*/
void grt_topbound_RU(GRT_MODEL1D *mod1d);
void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d);
void grt_topbound_RU_SH(GRT_MODEL1D *mod1d);

/**
* 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14)
Expand All @@ -33,7 +34,8 @@ void grt_topbound_RU(GRT_MODEL1D *mod1d);
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot
*
*/
void grt_botbound_RD(GRT_MODEL1D *mod1d);
void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d);
void grt_botbound_RD_SH(GRT_MODEL1D *mod1d);


/**
Expand Down Expand Up @@ -125,6 +127,8 @@ void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M
*
*/
void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);
void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);
void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M);

/**
* 为虚拟层的广义 R/T 矩阵添加时间延迟因子
Expand Down
6 changes: 4 additions & 2 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 @@ -25,15 +25,17 @@
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_top
*
*/
void grt_static_topbound_RU(GRT_MODEL1D *mod1d);
void grt_static_topbound_RU_PSV(GRT_MODEL1D *mod1d);
void grt_static_topbound_RU_SH(GRT_MODEL1D *mod1d);

/**
* 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12)
*
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot
*
*/
void grt_static_botbound_RD(GRT_MODEL1D *mod1d);
void grt_static_botbound_RD_PSV(GRT_MODEL1D *mod1d);
void grt_static_botbound_RD_SH(GRT_MODEL1D *mod1d);

/**
* 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35)
Expand Down
34 changes: 20 additions & 14 deletions pygrt/C_extension/src/common/RT_matrix.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,20 +22,26 @@

void grt_reset_RT_matrix(RT_MATRIX *M)
{
RT_MATRIX *M0 = &(RT_MATRIX){
.RD = GRT_INIT_ZERO_2x2_MATRIX,
.RU = GRT_INIT_ZERO_2x2_MATRIX,
.TD = GRT_INIT_IDENTITY_2x2_MATRIX,
.TU = GRT_INIT_IDENTITY_2x2_MATRIX,
.RDL = 0.0,
.RUL = 0.0,
.TDL = 1.0,
.TUL = 1.0,
.invT = GRT_INIT_ZERO_2x2_MATRIX,
.invTL = 0.0,
.stats = GRT_INVERSE_SUCCESS
};
memcpy(M, M0, sizeof(*M));
memset(M, 0, sizeof(*M));
M->TD[0][0] = M->TD[1][1] =
M->TU[0][0] = M->TU[1][1] =
M->TDL = M->TUL = 1.0;
M->stats = GRT_INVERSE_SUCCESS;
}

void grt_reset_RT_matrix_PSV(RT_MATRIX *M)
{
memset(M, 0, sizeof(*M));
M->TD[0][0] = M->TD[1][1] =
M->TU[0][0] = M->TU[1][1] = 1.0;
M->stats = GRT_INVERSE_SUCCESS;
}

void grt_reset_RT_matrix_SH(RT_MATRIX *M)
{
M->RDL = M->RUL = M->invTL = 0.0;
M->TDL = M->TUL = 1.0;
M->stats = GRT_INVERSE_SUCCESS;
}

void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M)
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/src/common/model.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -384,8 +384,8 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea
nlay++;
}

mod1d->isrc = isrc;
mod1d->ircv = ircv;
mod1d->isrc = (depsrc >= 0.0)? isrc : nlay; // 使用明显越界的索引，表明没有插入震源层或台站层
mod1d->ircv = (deprcv >= 0.0)? ircv : nlay;
mod1d->ircvup = ircvup;
mod1d->n = nlay;
mod1d->depsrc = depsrc;
Expand Down
129 changes: 117 additions & 12 deletions pygrt/C_extension/src/dynamic/layer.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
T N##2 = mod1d->N[iy];\


static int __freeBound_R(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML)
static int __freeBound_R_PSV(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2])
{
if( ! isLiquid){
// 固体表面
Expand All @@ -47,18 +47,29 @@ static int __freeBound_R(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, re
M[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta * adsgn;
M[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta * adsgn;
M[1][1] = M[0][0];
*ML = 1.0;
}
else {
// 液体表面
M[0][0] = -1.0;
M[1][1] = M[0][1] = M[1][0] = 0.0;
}
return GRT_INVERSE_SUCCESS;
}

static int __freeBound_R_SH(bool isLiquid, cplx_t *ML)
{
if( ! isLiquid){
// 固体表面
*ML = 1.0;
}
else {
// 液体表面
*ML = 0.0;
}
return GRT_INVERSE_SUCCESS;
}

static int __rigidBound_R(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML)
static int __rigidBound_R_PSV(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2])
{
if( ! isLiquid){
// 固体表面
Expand All @@ -73,31 +84,60 @@ static int __rigidBound_R(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, c
M[0][1] = 2.0 * xb0 * Delta * adsgn;
M[1][0] = 2.0 * xa0 * Delta * adsgn;
M[1][1] = M[0][0];
*ML = - 1.0;
}
else {
// 液体表面
M[0][0] = 1.0;
M[1][1] = M[0][1] = M[1][0] = 0.0;
}
return GRT_INVERSE_SUCCESS;
}

static int __rigidBound_R_SH(bool isLiquid, cplx_t *ML)
{
if( ! isLiquid){
// 固体表面
*ML = - 1.0;
}
else {
// 液体表面
*ML = 0.0;
}
return GRT_INVERSE_SUCCESS;
}

void grt_topbound_RU(GRT_MODEL1D *mod1d)
void grt_topbound_RU_PSV(GRT_MODEL1D *mod1d)
{
cplx_t cbcb = mod1d->cbcb[0];
cplx_t xa = mod1d->xa[0];
cplx_t xb = mod1d->xb[0];
bool isLiquid = mod1d->isLiquid[0];
if(mod1d->topbound == GRT_BOUND_FREE){
mod1d->M_top.stats = __freeBound_R(isLiquid, cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL);
mod1d->M_top.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, 1.0, mod1d->M_top.RU);
}
else if(mod1d->topbound == GRT_BOUND_RIGID){
mod1d->M_top.stats = __rigidBound_R(isLiquid, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL);
mod1d->M_top.stats = __rigidBound_R_PSV(isLiquid, xa, xb, 1.0, mod1d->M_top.RU);
}
else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4);
mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
}
else{
GRTRaiseError("Wrong execution.");
}
// RU 不需要时延
}

void grt_topbound_RU_SH(GRT_MODEL1D *mod1d)
{
bool isLiquid = mod1d->isLiquid[0];
if(mod1d->topbound == GRT_BOUND_FREE){
mod1d->M_top.stats = __freeBound_R_SH(isLiquid, &mod1d->M_top.RUL);
}
else if(mod1d->topbound == GRT_BOUND_RIGID){
mod1d->M_top.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_top.RUL);
}
else if(mod1d->topbound == GRT_BOUND_HALFSPACE){
mod1d->M_top.RUL = 0.0;
mod1d->M_top.stats = GRT_INVERSE_SUCCESS;
}
Expand All @@ -107,22 +147,21 @@ void grt_topbound_RU(GRT_MODEL1D *mod1d)
// RU 不需要时延
}

void grt_botbound_RD(GRT_MODEL1D *mod1d)
void grt_botbound_RD_PSV(GRT_MODEL1D *mod1d)
{
size_t nlay = mod1d->n;
cplx_t cbcb = mod1d->cbcb[nlay-2];
cplx_t xa = mod1d->xa[nlay-2];
cplx_t xb = mod1d->xb[nlay-2];
cplx_t isLiquid = mod1d->isLiquid[nlay-2];
bool isLiquid = mod1d->isLiquid[nlay-2];
if(mod1d->botbound == GRT_BOUND_FREE){
mod1d->M_bot.stats = __freeBound_R(isLiquid, cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
mod1d->M_bot.stats = __freeBound_R_PSV(isLiquid, cbcb, xa, xb, -1.0, mod1d->M_bot.RD);
}
else if(mod1d->botbound == GRT_BOUND_RIGID){
mod1d->M_bot.stats = __rigidBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL);
mod1d->M_bot.stats = __rigidBound_R_PSV(isLiquid, xa, xb, -1.0, mod1d->M_bot.RD);
}
else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot);
grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
}
else{
GRTRaiseError("Wrong execution.");
Expand All @@ -141,6 +180,33 @@ void grt_botbound_RD(GRT_MODEL1D *mod1d)

mod1d->M_bot.RD[0][0] *= ex2a; mod1d->M_bot.RD[0][1] *= exab;
mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b;
}

void grt_botbound_RD_SH(GRT_MODEL1D *mod1d)
{
size_t nlay = mod1d->n;
cplx_t xb = mod1d->xb[nlay-2];
bool isLiquid = mod1d->isLiquid[nlay-2];
if(mod1d->botbound == GRT_BOUND_FREE){
mod1d->M_bot.stats = __freeBound_R_SH(isLiquid, &mod1d->M_bot.RDL);
}
else if(mod1d->botbound == GRT_BOUND_RIGID){
mod1d->M_bot.stats = __rigidBound_R_SH(isLiquid, &mod1d->M_bot.RDL);
}
else if(mod1d->botbound == GRT_BOUND_HALFSPACE){
grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot);
}
else{
GRTRaiseError("Wrong execution.");
}

// 时延 RD
real_t k = mod1d->k;
real_t thk = mod1d->Thk[nlay-2];

cplx_t exb, ex2b;
exb = exp(- k*thk*xb);
ex2b = exb * exb;

mod1d->M_bot.RDL *= ex2b;
}
Expand Down Expand Up @@ -558,6 +624,45 @@ void grt_delay_RT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M
M->TUL *= exb;
}

void grt_delay_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
{
real_t thk = mod1d->Thk[iy-1];
cplx_t xa1 = mod1d->xa[iy-1];
cplx_t xb1 = mod1d->xb[iy-1];
real_t k = mod1d->k;

cplx_t exa, exb, ex2a, ex2b, exab;
exa = exp(- k*thk*xa1);
exb = exp(- k*thk*xb1);
ex2a = exa * exa;
ex2b = exb * exb;
exab = exa * exb;

M->RD[0][0] *= ex2a; M->RD[0][1] *= exab;
M->RD[1][0] *= exab; M->RD[1][1] *= ex2b;

M->TD[0][0] *= exa; M->TD[0][1] *= exb;
M->TD[1][0] *= exa; M->TD[1][1] *= exb;

M->TU[0][0] *= exa; M->TU[0][1] *= exa;
M->TU[1][0] *= exb; M->TU[1][1] *= exb;
}

void grt_delay_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
{
real_t thk = mod1d->Thk[iy-1];
cplx_t xb1 = mod1d->xb[iy-1];
real_t k = mod1d->k;

cplx_t exb, ex2b;
exb = exp(- k*thk*xb1);
ex2b = exb * exb;

M->RDL *= ex2b;
M->TDL *= exb;
M->TUL *= exb;
}



void grt_delay_GRT_matrix(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M)
Expand Down
12 changes: 8 additions & 4 deletions pygrt/C_extension/src/integral/kernel_template.c_
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,10 @@
#define __grt_RT_matrix_SH CONCAT_PREFIX(RT_matrix_SH)
#define __grt_delay_RT_matrix CONCAT_PREFIX(delay_RT_matrix)
#define __grt_source_coef CONCAT_PREFIX(source_coef)
#define __grt_topbound_RU CONCAT_PREFIX(topbound_RU)
#define __grt_botbound_RD CONCAT_PREFIX(botbound_RD)
#define __grt_topbound_RU_PSV CONCAT_PREFIX(topbound_RU_PSV)
#define __grt_topbound_RU_SH CONCAT_PREFIX(topbound_RU_SH)
#define __grt_botbound_RD_PSV CONCAT_PREFIX(botbound_RD_PSV)
#define __grt_botbound_RD_SH CONCAT_PREFIX(botbound_RD_SH)
#define __grt_wave2qwv_REV_PSV CONCAT_PREFIX(wave2qwv_REV_PSV)
#define __grt_wave2qwv_REV_SH CONCAT_PREFIX(wave2qwv_REV_SH)
#define __grt_wave2qwv_z_REV_PSV CONCAT_PREFIX(wave2qwv_z_REV_PSV)
Expand Down Expand Up @@ -153,14 +155,16 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k)
//===================================================================================

// 递推RU_FA
__grt_topbound_RU(mod1d);
__grt_topbound_RU_PSV(mod1d);
__grt_topbound_RU_SH(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;

// 递推RD_BL, 在模型定义阶段已检查震源/台站的虚拟层成为底层时，只能是半无限空间
if(imax+1 < mod1d->n){ // 此 if 仅减少不必要的计算
__grt_botbound_RD(mod1d);
__grt_botbound_RD_PSV(mod1d);
__grt_botbound_RD_SH(mod1d);
if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN;
grt_recursion_RD(M_BL, M_bot, M_BL);
if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN;
Expand Down
Loading