Skip to content

REFAC: assemble source coefficients functions of P-SV and SH #140

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: assemble source coefficients functions of P-SV and SH #140

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/kernel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,7 @@ typedef void (*GRT_KernelFunc) (
#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)
#define __grt_source_coef_PSV CONCAT_PREFIX(source_coef_PSV)
#define __grt_source_coef_SH CONCAT_PREFIX(source_coef_SH)
#define __grt_source_coef CONCAT_PREFIX(source_coef)
#define __grt_topfree_RU CONCAT_PREFIX(topfree_RU)
#define __grt_wave2qwv_REV_PSV CONCAT_PREFIX(wave2qwv_REV_PSV)
#define __grt_wave2qwv_REV_SH CONCAT_PREFIX(wave2qwv_REV_SH)
Expand Down
28 changes: 12 additions & 16 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 @@ -58,28 +58,26 @@
* @param[in] RL1 SH波， \f$ R_1\f$
* @param[in] R2 P-SV波，\f$\mathbf{R_2}\f$矩阵
* @param[in] RL2 SH波， \f$ R_2\f$
* @param[in] coef_PSV P-SV 波震源系数，\f$ P_m, SV_m\f$ ，维度2表示下行波(p=0)和上行波(p=1)
* @param[in] coef_SH SH 波震源系数，\f$ SH_m \f$ ，维度2表示下行波(p=0)和上行波(p=1)
* @param[in] coef 震源系数，\f$ P_m, SV_m，SH_m\f$ ，维度2表示下行波(p=0)和上行波(p=1)
* @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$
*/
inline void grt_construct_qwv(
bool ircvup,
const cplx_t R1[2][2], cplx_t RL1,
const cplx_t R2[2][2], cplx_t RL2,
const cplx_t coef_PSV[GRT_QWV_NUM-1][2], const cplx_t coef_SH[2],
cplx_t qwv[GRT_QWV_NUM])
const cplx_t coef[GRT_QWV_NUM][2], cplx_t qwv[GRT_QWV_NUM])
{
cplx_t qw0[2], qw1[2], v0;
cplx_t coefD[2] = {coef_PSV[0][0], coef_PSV[1][0]};
cplx_t coefU[2] = {coef_PSV[0][1], coef_PSV[1][1]};
cplx_t coefD[2] = {coef[0][0], coef[1][0]};
cplx_t coefU[2] = {coef[0][1], coef[1][1]};
if(ircvup){
grt_cmat2x1_mul(R2, coefD, qw0);
qw0[0] += coefU[0]; qw0[1] += coefU[1];
v0 = RL1 * (RL2*coef_SH[0] + coef_SH[1]);
v0 = RL1 * (RL2*coef[2][0] + coef[2][1]);
} else {
grt_cmat2x1_mul(R2, coefU, qw0);
qw0[0] += coefD[0]; qw0[1] += coefD[1];
v0 = RL1 * (coef_SH[0] + RL2*coef_SH[1]);
v0 = RL1 * (coef[2][0] + RL2*coef[2][1]);
}
grt_cmat2x1_mul(R1, qw0, qw1);

Expand Down Expand Up @@ -256,10 +254,8 @@ static void __KERNEL_FUNC__(


// 计算震源系数
cplx_t src_coef_PSV[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2] = {0};
cplx_t src_coef_SH[GRT_SRC_M_NUM][2] = {0};
__grt_source_coef_PSV(mod1d, src_coef_PSV);
__grt_source_coef_SH(mod1d, src_coef_SH);
cplx_t src_coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2] = {0};
__grt_source_coef(mod1d, src_coef);

// 临时中转矩阵 (temperary)
cplx_t tmpR2[2][2], tmp2x2[2][2], tmpRL, tmp2x2_uiz[2][2], tmpRL2;
Expand Down Expand Up @@ -295,7 +291,7 @@ static void __KERNEL_FUNC__(
tmpRL2 = 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_PSV[i], src_coef_SH[i], QWV[i]);
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, M_BL->RDL, src_coef[i], QWV[i]);
}


Expand All @@ -306,7 +302,7 @@ static void __KERNEL_FUNC__(
tmpRL2 = 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_PSV[i], src_coef_SH[i], QWV_uiz[i]);
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, M_BL->RDL, src_coef[i], QWV_uiz[i]);
}
}
}
Expand Down Expand Up @@ -334,7 +330,7 @@ static void __KERNEL_FUNC__(
tmpRL2 = 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_PSV[i], src_coef_SH[i], QWV[i]);
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, M_FA->RUL, src_coef[i], QWV[i]);
}


Expand All @@ -345,7 +341,7 @@ static void __KERNEL_FUNC__(
tmpRL2 = 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_PSV[i], src_coef_SH[i], QWV_uiz[i]);
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, M_FA->RUL, src_coef[i], QWV_uiz[i]);
}
}

Expand Down
7 changes: 5 additions & 2 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 @@ -24,7 +24,10 @@
* @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$
*
*/
void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]);
void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);

/* P-SV 波的震源系数 */
void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);

/* SH 波的震源系数 */
void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]);
void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);
8 changes: 6 additions & 2 deletions pygrt/C_extension/include/grt/static/static_source.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,11 @@
* @param[in] mod1d 模型结构体指针
* @param[out] coef 震源系数 \f$ P_m, SV_m, SH_m \f$
*/
void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2]);
void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);


/* P-SV 波的静态震源系数 */
void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);

/* SH 波的静态震源系数 */
void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2]);
void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]);
94 changes: 45 additions & 49 deletions pygrt/C_extension/src/dynamic/source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,93 +11,89 @@

#include <stdio.h>
#include <complex.h>
#include <string.h>

#include "grt/dynamic/source.h"
#include "grt/common/model.h"
#include "grt/common/matrix.h"
#include "grt/common/prtdbg.h"


void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2])
inline void _source_PSV(
const cplx_t xa, const cplx_t caca,
const cplx_t xb, const cplx_t cbcb, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
// 先全部赋0
for(int i=0; i<GRT_SRC_M_NUM; ++i){
for(int j=0; j<GRT_QWV_NUM-1; ++j){
for(int p=0; p<2; ++p){
coef[i][j][p] = 0.0;
}
}
}

size_t isrc = mod1d->isrc;
cplx_t xa = mod1d->xa[isrc];
cplx_t caca = mod1d->caca[isrc];
cplx_t xb = mod1d->xb[isrc];
cplx_t cbcb = mod1d->cbcb[isrc];
real_t k = mod1d->k;

cplx_t tmp;


// 爆炸源， 通过(4.9.8)的矩张量源公式，提取各向同性的量(M11+M22+M33)，-a+k^2/a -> ka^2/a
coef[0][0][0] = tmp = (caca / xa) * k; coef[0][0][1] = tmp;
coef[0][0][0] = tmp = (caca / xa) * k; coef[0][0][1] = tmp;

// 垂直力源 (4.6.15)
coef[1][0][0] = tmp = -1.0; coef[1][0][1] = - tmp;
coef[1][0][0] = tmp = -1.0; coef[1][0][1] = - tmp;
coef[1][1][0] = tmp = -1.0 / xb; coef[1][1][1] = tmp;

// 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向
// 推导可发现，r方向的力形成P,SV波, theta方向的力形成SH波
// 方向性因子包含水平力方向与震源台站连线方向的夹角
// 水平力源 (4.6.21,26)
coef[2][0][0] = tmp = -1.0 / xa; coef[2][0][1] = tmp;
coef[2][1][0] = tmp = -1.0; coef[2][1][1] = - tmp;
coef[2][1][0] = tmp = -1.0; coef[2][1][1] = - tmp;

// 剪切位错 (4.8.34)
// m=0
coef[3][0][0] = tmp = ((2.0*caca - 3.0) / xa) * k; coef[3][0][1] = tmp;
coef[3][1][0] = tmp = -3.0*k; coef[3][1][1] = - tmp;
coef[3][1][0] = tmp = -3.0*k; coef[3][1][1] = - tmp;
// m=1
coef[4][0][0] = tmp = 2.0*k; coef[4][0][1] = - tmp;
coef[4][1][0] = tmp = ((2.0 - cbcb) / xb) * k; coef[4][1][1] = tmp;

// m=2
coef[5][0][0] = tmp = - (1.0 / xa) * k; coef[5][0][1] = tmp;
coef[5][0][0] = tmp = - (1.0 / xa) * k; coef[5][0][1] = tmp;
coef[5][1][0] = tmp = - k; coef[5][1][1] = - tmp;

}

inline void _source_SH(const cplx_t xb, const cplx_t cbcb, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
cplx_t tmp;

// 水平力源 (4.6.21,26)
coef[2][2][0] = tmp = cbcb / xb; coef[2][2][1] = tmp;

// 剪切位错 (4.8.34)
// m=1
coef[4][2][0] = tmp = - cbcb * k; coef[4][2][1] = - tmp;

// m=2
coef[5][2][0] = tmp = (cbcb / xb) * k; coef[5][2][1] = tmp;
}


void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2])
void grt_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
// 先全部赋0
for(int i=0; i<GRT_SRC_M_NUM; ++i){
for(int p=0; p<2; ++p){
coef[i][p] = 0.0;
}
}
memset(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);
}


void grt_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
size_t isrc = mod1d->isrc;
cplx_t xa = mod1d->xa[isrc];
cplx_t caca = mod1d->caca[isrc];
cplx_t xb = mod1d->xb[isrc];
cplx_t cbcb = mod1d->cbcb[isrc];
real_t k = mod1d->k;

// cplx_t b = k*src_xb;
cplx_t tmp;

// 水平力源 (4.6.21,26), 这里可以把x1,x2方向的力转到r,theta方向
// 推导可发现，r方向的力形成P,SV波, theta方向的力形成SH波
// 方向性因子包含水平力方向与震源台站连线方向的夹角
coef[2][0] = tmp = cbcb / xb; coef[2][1] = tmp;
_source_PSV(xa, caca, xb, cbcb, k, coef);
}

// 剪切位错 (4.8.34)
// m=1
coef[4][0] = tmp = - cbcb * k; coef[4][1] = - tmp;

// m=2
coef[5][0] = tmp = (cbcb / xb) * k; coef[5][1] = tmp;
void grt_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
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);
}




1 change: 1 addition & 0 deletions pygrt/C_extension/src/static/static_layer.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ void grt_static_wave2qwv_REV_PSV(const GRT_MODEL1D *mod1d, const cplx_t R[2][2],

void grt_static_wave2qwv_REV_SH(const GRT_MODEL1D *mod1d, cplx_t RL, cplx_t *R_EVL)
{
(void)mod1d; // 暂停该参数 -Wunused-parameter 警告
*R_EVL = (1.0 + (RL));
}

Expand Down
69 changes: 38 additions & 31 deletions pygrt/C_extension/src/static/static_source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,25 +12,13 @@

#include <stdio.h>
#include <complex.h>
#include <string.h>

#include "grt/static/static_source.h"
#include "grt/common/const.h"

void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM-1][2])
{
// 先全部赋0
for(int i=0; i<GRT_SRC_M_NUM; ++i){
for(int j=0; j<GRT_QWV_NUM-1; ++j){
for(int p=0; p<2; ++p){
coef[i][j][p] = 0.0;
}
}
}

size_t isrc = mod1d->isrc;
cplx_t delta = mod1d->delta[isrc];
real_t k = mod1d->k;

inline void _source_PSV(const cplx_t delta, const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
cplx_t tmp;
cplx_t A = 1.0+delta;

Expand All @@ -39,46 +27,65 @@ void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_

// 垂直力源
coef[1][0][0] = tmp = -1.0/(2.0*A*k); coef[1][0][1] = - tmp;
coef[1][1][0] = tmp; coef[1][1][1] = - tmp;
coef[1][1][0] = tmp; coef[1][1][1] = - tmp;

// 水平力源
coef[2][0][0] = tmp = 1.0/(2.0*A*k); coef[2][0][1] = tmp;
coef[2][0][0] = tmp = 1.0/(2.0*A*k); coef[2][0][1] = tmp;
coef[2][1][0] = - tmp; coef[2][1][1] = - tmp;

// 剪切位错
// m=0
coef[3][0][0] = tmp = (-1.0+4.0*delta)/(2.0*A); coef[3][0][1] = tmp;
coef[3][1][0] = tmp = -3.0/(2.0*A); coef[3][1][1] = tmp;
// m=1
coef[4][0][0] = tmp = -delta/A; coef[4][0][1] = -tmp;
coef[4][0][0] = tmp = -delta/A; coef[4][0][1] = -tmp;
coef[4][1][0] = tmp = 1.0/A; coef[4][1][1] = -tmp;
// m=2
coef[5][0][0] = tmp = 1.0/(2.0*A); coef[5][0][1] = tmp;
coef[5][1][0] = tmp = -1.0/(2.0*A); coef[5][1][1] = tmp;
}


void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][2])
inline void _source_SH(const real_t k, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
real_t k = mod1d->k;

// 先全部赋0
for(int i=0; i<GRT_SRC_M_NUM; ++i){
for(int p=0; p<2; ++p){
coef[i][p] = 0.0;
}
}

cplx_t tmp;

// 水平力源
coef[2][0] = tmp = -1.0/k; coef[2][1] = tmp;
coef[2][2][0] = tmp = - 1.0/k; coef[2][2][1] = tmp;

// 剪切位错
// m=1
coef[4][0] = tmp = 1.0; coef[4][1] = -tmp;
coef[4][2][0] = tmp = 1.0; coef[4][2][1] = -tmp;
// m=2
coef[5][0] = tmp = -1.0; coef[5][1] = tmp;
coef[5][2][0] = tmp = - 1.0; coef[5][2][1] = tmp;
}


void grt_static_source_coef(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
// 先全部赋0
memset(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);
}


void grt_static_source_coef_PSV(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
size_t isrc = mod1d->isrc;
cplx_t delta = mod1d->delta[isrc];
real_t k = mod1d->k;

_source_PSV(delta, k, coef);
}


void grt_static_source_coef_SH(const GRT_MODEL1D *mod1d, cplx_t coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2])
{
real_t k = mod1d->k;

_source_SH(k, coef);
}