Skip to content

REFAC: define the data type of the QWV-related arrays #164

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

REFAC: define the data type of the QWV-related arrays #164

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: 3 additions & 0 deletions pygrt/C_extension/include/grt/common/const.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,9 @@ typedef double complex cplx_t;

#define GRT_GTYPES_MAX 2 ///< 2, 所有震源根据是否使用格林函数导数分为两类

typedef cplx_t QWVgrid[GRT_SRC_M_NUM][GRT_QWV_NUM];
typedef cplx_t INTEGgrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];

/** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */
enum {
GRT_SRC_M_EX_INDEX = 0,
Expand Down
6 changes: 3 additions & 3 deletions pygrt/C_extension/include/grt/common/integral.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
*/
void grt_int_Pk(
real_t k, real_t r,
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV,
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);

Expand Down Expand Up @@ -60,7 +60,7 @@ void grt_merge_Pk(
*/
void grt_int_Pk_filon(
real_t k, real_t r, bool iscos,
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV,
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);

Expand All @@ -79,6 +79,6 @@ void grt_int_Pk_filon(
*/
void grt_int_Pk_sa_filon(
const real_t k3[3], real_t r,
const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV3[3],
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
3 changes: 1 addition & 2 deletions pygrt/C_extension/include/grt/common/iostats.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,7 @@
* @note 文件记录的值均为波数积分的中间结果，与最终的结果还差一系列的系数，
* 记录其值主要用于参考其变化趋势。
*/
void grt_write_stats(
FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]);
void grt_write_stats(FILE *f0, real_t k, const QWVgrid QWV);


/**
Expand Down
20 changes: 10 additions & 10 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 @@ -16,8 +16,8 @@
* 计算核函数的函数指针，动态与静态的接口一致
*/
typedef void (*GRT_KernelFunc) (
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_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz);


/**
Expand Down Expand Up @@ -89,26 +89,26 @@ typedef void (*GRT_KernelFunc) (
*
*/
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]);
GRT_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz);

/** 构建广义反射透射系数矩阵。作为 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]);
GRT_MODEL1D *mod1d, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz);

/** 静态解的核函数 */
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]);
GRT_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz);

/** 静态广义反射透射系数矩阵 */
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]);
GRT_MODEL1D *mod1d, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz);
5 changes: 3 additions & 2 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 @@ -67,8 +67,9 @@ typedef struct {
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$
/* 震源处的震源系数 \f$ P_m, SV_m, SH_m */
QWVgrid src_coefD;
QWVgrid src_coefU;

} GRT_MODEL1D;

Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/src/common/dwm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,8 @@ real_t grt_discrete_integ(
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM] = {0};

// 不同震源不同阶数的核函数 F(k, w)
cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM] = {0};
cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM] = {0};
QWVgrid QWV = {0};
QWVgrid QWV_uiz = {0};

real_t k = 0.0;
size_t ik = 0;
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/src/common/fim.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ real_t grt_linear_filon_integ(
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM];

// 不同震源不同阶数的核函数 F(k, w)
cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM];
cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM];
QWVgrid QWV = {0};
QWVgrid QWV_uiz = {0};

real_t k=k0;
size_t ik=0;
Expand Down
12 changes: 6 additions & 6 deletions pygrt/C_extension/src/common/integral.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
void grt_int_Pk(
real_t k, real_t r,
// F(ki,w)， 第一个维度表示不同震源，不同阶数，第二个维度3代表三类系数qm,wm,vm
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV,
// F(ki,w)Jm(ki*r)ki，第一个维度表示不同震源，不同阶数，第二个维度代表4种类型的F(k,w)Jm(kr)k的类型
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
Expand Down Expand Up @@ -72,7 +72,7 @@ void grt_int_Pk(

void grt_int_Pk_filon(
real_t k, real_t r, bool iscos,
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV,
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
{
Expand Down Expand Up @@ -161,7 +161,7 @@ static cplx_t interg_quad_cos(

void grt_int_Pk_sa_filon(
const real_t k3[3], real_t r,
const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
const QWVgrid QWV3[3],
bool calc_uir,
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
{
Expand All @@ -173,9 +173,9 @@ void grt_int_Pk_sa_filon(

// 对sqrt(k)*F(k,w)进行二次曲线拟合，再计算 (a*k^2 + b*k + c) * cos(kr - (2m+1)/4) 的积分
// 拟合二次函数的参数
cplx_t quad_a[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
cplx_t quad_b[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
cplx_t quad_c[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
QWVgrid quad_a={0};
QWVgrid quad_b={0};
QWVgrid quad_c={0};
for(int im=0; im<GRT_SRC_M_NUM; ++im){
int modr = GRT_SRC_M_ORDERS[im];
for(int c=0; c<GRT_QWV_NUM; ++c){
Expand Down
3 changes: 1 addition & 2 deletions pygrt/C_extension/src/common/iostats.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,7 @@



void grt_write_stats(
FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM])
void grt_write_stats(FILE *f0, real_t k, const QWVgrid QWV)
{
fwrite(&k, sizeof(real_t), 1, f0);

Expand Down
39 changes: 20 additions & 19 deletions pygrt/C_extension/src/common/kernel.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,35 +40,36 @@
* @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 震源系数，\f$ P_m, SV_m，SH_m\f$ ，维度2表示下行波(p=0)和上行波(p=1)
* @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$
* @param[in] coefD 下行震源系数，\f$ P_m, SV_m，SH_m\f$
* @param[in] coefU 上行震源系数，\f$ P_m, SV_m，SH_m\f$
* @param[out] QWV 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$
*/
inline GCC_ALWAYS_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[GRT_QWV_NUM][2], cplx_t qwv[GRT_QWV_NUM])
const QWVgrid coefD, const QWVgrid coefU, QWVgrid QWV)
{
cplx_t qw0[2], qw1[2], v0;
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[2][0] + coef[2][1]);
} else {
grt_cmat2x1_mul(R2, coefU, qw0);
qw0[0] += coefD[0]; qw0[1] += coefD[1];
v0 = RL1 * (coef[2][0] + RL2*coef[2][1]);
for(int i = 0; i < GRT_SRC_M_NUM; ++i)
{
if(ircvup){
grt_cmat2x1_mul(R2, coefD[i], QWV[i]);
QWV[i][0] += coefU[i][0];
QWV[i][1] += coefU[i][1];
QWV[i][2] = RL1 * (RL2*coefD[i][2] + coefU[i][2]);
} else {
grt_cmat2x1_mul(R2, coefU[i], QWV[i]);
QWV[i][0] += coefD[i][0];
QWV[i][1] += coefD[i][1];
QWV[i][2] = RL1 * (coefD[i][2] + RL2*coefU[i][2]);
}
grt_cmat2x1_mul(R1, QWV[i], QWV[i]);
}
grt_cmat2x1_mul(R1, qw0, qw1);

qwv[0] = qw1[0];
qwv[1] = qw1[1];
qwv[2] = v0;
}




// ================ 动态解 =====================
#define __DYNAMIC_KERNEL__
#include "kernel_template.c_"
Expand Down
24 changes: 8 additions & 16 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 @@ -154,8 +154,8 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k)


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])
GRT_MODEL1D *mod1d, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz)
{
// 初始化qwv为0
memset(QWV, 0, sizeof(cplx_t)*GRT_SRC_M_NUM*GRT_QWV_NUM);
Expand Down Expand Up @@ -198,9 +198,7 @@ void __BUILD_QWV__(
tmpRL = M_FB->invTL / (1.0 - M_BL->RDL * M_FB->RUL);
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, mod1d->src_coef[i], QWV[i]);
}
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_BL->RD, M_BL->RDL, mod1d->src_coefD, mod1d->src_coefU, QWV);


if(calc_uiz){
Expand All @@ -209,9 +207,7 @@ void __BUILD_QWV__(
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, mod1d->src_coef[i], QWV_uiz[i]);
}
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_BL->RD, M_BL->RDL, mod1d->src_coefD, mod1d->src_coefU, QWV_uiz);
}
}
else { // A震源 B接收
Expand All @@ -236,9 +232,7 @@ void __BUILD_QWV__(
tmpRL = M_AL->invTL / (1.0 - M_FA->RUL * M_AL->RDL);
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, mod1d->src_coef[i], QWV[i]);
}
grt_construct_qwv(ircvup, tmp2x2, tmpRL2, M_FA->RU, M_FA->RUL, mod1d->src_coefD, mod1d->src_coefU, QWV);


if(calc_uiz){
Expand All @@ -247,9 +241,7 @@ void __BUILD_QWV__(
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, mod1d->src_coef[i], QWV_uiz[i]);
}
grt_construct_qwv(ircvup, tmp2x2_uiz, tmpRL2, M_FA->RU, M_FA->RUL, mod1d->src_coefD, mod1d->src_coefU, QWV_uiz);
}

} // END if
Expand All @@ -271,8 +263,8 @@ void __BUILD_QWV__(


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_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
bool calc_uiz, QWVgrid QWV_uiz)
{
__GRT_MATRIX_FUNC__(mod1d, k);
__BUILD_QWV__(mod1d, QWV, calc_uiz, QWV_uiz);
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/src/common/ptam.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,8 +132,8 @@ void grt_PTA_method(
real_t k=0.0;

// 不同震源不同阶数的核函数 F(k, w)
cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM];
cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM];
QWVgrid QWV = {0};
QWVgrid QWV_uiz = {0};

// 使用宏函数，方便定义
#define __CALLOC_ARRAY(VAR, TYP, __ARR) \
Expand Down
12 changes: 6 additions & 6 deletions pygrt/C_extension/src/common/safim.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,8 @@ typedef enum {
// 区间结构体
typedef struct {
real_t k3[3];
cplx_t F3[3][GRT_SRC_M_NUM][GRT_QWV_NUM];
cplx_t F3_uiz[3][GRT_SRC_M_NUM][GRT_QWV_NUM];
QWVgrid F3[3];
QWVgrid F3_uiz[3];
} KInterval;

// 区间栈结构体
Expand Down Expand Up @@ -104,7 +104,7 @@ static KInterval stack_pop(KIntervalStack *stack) {
static cplx_t simpson(const KInterval *item_pt, int im, int iqwv, bool isuiz, SIMPSON_INTV stats) {
cplx_t Fint = 0.0;
real_t klen = item_pt->k3[2] - item_pt->k3[0];
const cplx_t (*F3)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? item_pt->F3_uiz : item_pt->F3;
const QWVgrid *F3 = (isuiz)? item_pt->F3_uiz : item_pt->F3;

// 使用F(k)*sqrt(k)来衡量积分值，这可以平衡后续计算F(k)*Jm(kr)k积分时的系数
real_t sk[3];
Expand Down Expand Up @@ -136,7 +136,7 @@ static cplx_t simpson(const KInterval *item_pt, int im, int iqwv, bool isuiz, SI
}

/** 比较QWV的最大绝对值 */
static void get_maxabsQWV(const cplx_t F[GRT_SRC_M_NUM][GRT_QWV_NUM], real_t maxabsF[GRT_GTYPES_MAX]){
static void get_maxabsQWV(const QWVgrid F, real_t maxabsF[GRT_GTYPES_MAX]){
real_t tmp;
for(int i=0; i<GRT_SRC_M_NUM; ++i){
for(int c=0; c<GRT_QWV_NUM; ++c){
Expand All @@ -158,8 +158,8 @@ static bool check_fit(
cplx_t S11, S12, S21, S22;

// 核函数
const cplx_t (*F3L)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvL->F3_uiz : ptKitvL->F3;
const cplx_t (*F3R)[GRT_SRC_M_NUM][GRT_QWV_NUM] = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3;
const QWVgrid *F3L = (isuiz)? ptKitvL->F3_uiz : ptKitvL->F3;
const QWVgrid *F3R = (isuiz)? ptKitvR->F3_uiz : ptKitvR->F3;

// 取近似积分 \int_k1^k2 k^0.5 dk
real_t kcoef13 = RTWOTHIRD*( ptKitv->k3[2]*sqrt(ptKitv->k3[2]) - ptKitv->k3[0]*sqrt(ptKitv->k3[0]) );
Expand Down
Loading