Skip to content

Commit 65a5454

Browse files
Dengda98Dengda98
authored
REFAC: define the data type of the QWV-related arrays (#164)
1 parent 41c6686 commit 65a5454

File tree

15 files changed

+107
-110
lines changed

15 files changed

+107
-110
lines changed

‎pygrt/C_extension/include/grt/common/const.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -123,6 +123,9 @@ typedef double complex cplx_t;
123123

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

126+
typedef cplx_t QWVgrid[GRT_SRC_M_NUM][GRT_QWV_NUM];
127+
typedef cplx_t INTEGgrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];
128+
126129
/** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */
127130
enum {
128131
GRT_SRC_M_EX_INDEX = 0,

‎pygrt/C_extension/include/grt/common/integral.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@
2626
*/
2727
void grt_int_Pk(
2828
real_t k, real_t r,
29-
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
29+
const QWVgrid QWV,
3030
bool calc_uir,
3131
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
3232

@@ -60,7 +60,7 @@ void grt_merge_Pk(
6060
*/
6161
void grt_int_Pk_filon(
6262
real_t k, real_t r, bool iscos,
63-
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
63+
const QWVgrid QWV,
6464
bool calc_uir,
6565
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
6666

@@ -79,6 +79,6 @@ void grt_int_Pk_filon(
7979
*/
8080
void grt_int_Pk_sa_filon(
8181
const real_t k3[3], real_t r,
82-
const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
82+
const QWVgrid QWV3[3],
8383
bool calc_uir,
8484
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);

‎pygrt/C_extension/include/grt/common/iostats.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,7 @@
2626
* @note 文件记录的值均为波数积分的中间结果，与最终的结果还差一系列的系数，
2727
* 记录其值主要用于参考其变化趋势。
2828
*/
29-
void grt_write_stats(
30-
FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]);
29+
void grt_write_stats(FILE *f0, real_t k, const QWVgrid QWV);
3130

3231

3332
/**

‎pygrt/C_extension/include/grt/common/kernel.h‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,8 @@
1616
* 计算核函数的函数指针，动态与静态的接口一致
1717
*/
1818
typedef void (*GRT_KernelFunc) (
19-
GRT_MODEL1D *mod1d, const real_t k, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
20-
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
19+
GRT_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
20+
bool calc_uiz, QWVgrid QWV_uiz);
2121

2222

2323
/**
@@ -89,26 +89,26 @@ typedef void (*GRT_KernelFunc) (
8989
*
9090
*/
9191
void grt_kernel(
92-
GRT_MODEL1D *mod1d, const real_t k, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
93-
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
92+
GRT_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
93+
bool calc_uiz, QWVgrid QWV_uiz);
9494

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

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

103103
/** 静态解的核函数 */
104104
void grt_static_kernel(
105-
GRT_MODEL1D *mod1d, const real_t k, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
106-
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
105+
GRT_MODEL1D *mod1d, const real_t k, QWVgrid QWV,
106+
bool calc_uiz, QWVgrid QWV_uiz);
107107

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

111111
/** 静态 QWV */
112112
void grt_static_GRT_build_QWV(
113-
GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
114-
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
113+
GRT_MODEL1D *mod1d, QWVgrid QWV,
114+
bool calc_uiz, QWVgrid QWV_uiz);

‎pygrt/C_extension/include/grt/common/model.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -67,8 +67,9 @@ typedef struct {
6767
cplx_t uiz_R_EV[2][2];
6868
cplx_t uiz_R_EVL;
6969

70-
/* 震源处的震源系数 */
71-
cplx_t src_coef[GRT_SRC_M_NUM][GRT_QWV_NUM][2]; ///< 震源系数 \f$ P_m, SV_m, SH_m \f$
70+
/* 震源处的震源系数 \f$ P_m, SV_m, SH_m */
71+
QWVgrid src_coefD;
72+
QWVgrid src_coefU;
7273

7374
} GRT_MODEL1D;
7475

‎pygrt/C_extension/src/common/dwm.c‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,8 @@ real_t grt_discrete_integ(
3636
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM] = {0};
3737

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

4242
real_t k = 0.0;
4343
size_t ik = 0;

‎pygrt/C_extension/src/common/fim.c‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,8 +39,8 @@ real_t grt_linear_filon_integ(
3939
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM];
4040

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

4545
real_t k=k0;
4646
size_t ik=0;

‎pygrt/C_extension/src/common/integral.c‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222
void grt_int_Pk(
2323
real_t k, real_t r,
2424
// F(ki,w)， 第一个维度表示不同震源，不同阶数，第二个维度3代表三类系数qm,wm,vm
25-
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
25+
const QWVgrid QWV,
2626
// F(ki,w)Jm(ki*r)ki，第一个维度表示不同震源，不同阶数，第二个维度代表4种类型的F(k,w)Jm(kr)k的类型
2727
bool calc_uir,
2828
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
@@ -72,7 +72,7 @@ void grt_int_Pk(
7272

7373
void grt_int_Pk_filon(
7474
real_t k, real_t r, bool iscos,
75-
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
75+
const QWVgrid QWV,
7676
bool calc_uir,
7777
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
7878
{
@@ -161,7 +161,7 @@ static cplx_t interg_quad_cos(
161161

162162
void grt_int_Pk_sa_filon(
163163
const real_t k3[3], real_t r,
164-
const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
164+
const QWVgrid QWV3[3],
165165
bool calc_uir,
166166
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM])
167167
{
@@ -173,9 +173,9 @@ void grt_int_Pk_sa_filon(
173173

174174
// 对sqrt(k)*F(k,w)进行二次曲线拟合，再计算 (a*k^2 + b*k + c) * cos(kr - (2m+1)/4) 的积分
175175
// 拟合二次函数的参数
176-
cplx_t quad_a[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
177-
cplx_t quad_b[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
178-
cplx_t quad_c[GRT_SRC_M_NUM][GRT_QWV_NUM]={0};
176+
QWVgrid quad_a={0};
177+
QWVgrid quad_b={0};
178+
QWVgrid quad_c={0};
179179
for(int im=0; im<GRT_SRC_M_NUM; ++im){
180180
int modr = GRT_SRC_M_ORDERS[im];
181181
for(int c=0; c<GRT_QWV_NUM; ++c){

‎pygrt/C_extension/src/common/iostats.c‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,7 @@
1717

1818

1919

20-
void grt_write_stats(
21-
FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM])
20+
void grt_write_stats(FILE *f0, real_t k, const QWVgrid QWV)
2221
{
2322
fwrite(&k, sizeof(real_t), 1, f0);
2423

‎pygrt/C_extension/src/common/kernel.c‎

Lines changed: 20 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -40,35 +40,36 @@
4040
* @param[in] RL1 SH波， \f$ R_1\f$
4141
* @param[in] R2 P-SV波，\f$\mathbf{R_2}\f$矩阵
4242
* @param[in] RL2 SH波， \f$ R_2\f$
43-
* @param[in] coef 震源系数，\f$ P_m, SV_m，SH_m\f$ ，维度2表示下行波(p=0)和上行波(p=1)
44-
* @param[out] qwv 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$
43+
* @param[in] coefD 下行震源系数，\f$ P_m, SV_m，SH_m\f$
44+
* @param[in] coefU 上行震源系数，\f$ P_m, SV_m，SH_m\f$
45+
* @param[out] QWV 最终通过矩阵传播计算出的在台站位置的\f$ q_m,w_m,v_m\f$
4546
*/
4647
inline GCC_ALWAYS_INLINE void grt_construct_qwv(
4748
bool ircvup,
4849
const cplx_t R1[2][2], cplx_t RL1,
4950
const cplx_t R2[2][2], cplx_t RL2,
50-
const cplx_t coef[GRT_QWV_NUM][2], cplx_t qwv[GRT_QWV_NUM])
51+
const QWVgrid coefD, const QWVgrid coefU, QWVgrid QWV)
5152
{
52-
cplx_t qw0[2], qw1[2], v0;
53-
cplx_t coefD[2] = {coef[0][0], coef[1][0]};
54-
cplx_t coefU[2] = {coef[0][1], coef[1][1]};
55-
if(ircvup){
56-
grt_cmat2x1_mul(R2, coefD, qw0);
57-
qw0[0] += coefU[0]; qw0[1] += coefU[1];
58-
v0 = RL1 * (RL2*coef[2][0] + coef[2][1]);
59-
} else {
60-
grt_cmat2x1_mul(R2, coefU, qw0);
61-
qw0[0] += coefD[0]; qw0[1] += coefD[1];
62-
v0 = RL1 * (coef[2][0] + RL2*coef[2][1]);
53+
for(int i = 0; i < GRT_SRC_M_NUM; ++i)
54+
{
55+
if(ircvup){
56+
grt_cmat2x1_mul(R2, coefD[i], QWV[i]);
57+
QWV[i][0] += coefU[i][0];
58+
QWV[i][1] += coefU[i][1];
59+
QWV[i][2] = RL1 * (RL2*coefD[i][2] + coefU[i][2]);
60+
} else {
61+
grt_cmat2x1_mul(R2, coefU[i], QWV[i]);
62+
QWV[i][0] += coefD[i][0];
63+
QWV[i][1] += coefD[i][1];
64+
QWV[i][2] = RL1 * (coefD[i][2] + RL2*coefU[i][2]);
65+
}
66+
grt_cmat2x1_mul(R1, QWV[i], QWV[i]);
6367
}
64-
grt_cmat2x1_mul(R1, qw0, qw1);
65-
66-
qwv[0] = qw1[0];
67-
qwv[1] = qw1[1];
68-
qwv[2] = v0;
6968
}
7069

7170

71+
72+
7273
// ================ 动态解 =====================
7374
#define __DYNAMIC_KERNEL__
7475
#include "kernel_template.c_"

0 commit comments

Comments
 (0)