Skip to content

REFAC: rearrange kernel template #155

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

REFAC: rearrange kernel template #155

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
108 changes: 75 additions & 33 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 @@ -11,7 +11,6 @@
#pragma once

#include "grt/common/model.h"
#include "grt/common/matrix.h"

/**
* 计算核函数的函数指针，动态与静态的接口一致
Expand All @@ -21,36 +20,79 @@ typedef void (*GRT_KernelFunc) (
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);


#if defined(__DYNAMIC_KERNEL__) || defined(__STATIC_KERNEL__)

// ================ 动态解 =====================
#ifdef __DYNAMIC_KERNEL__
#define CONCAT_PREFIX(base) grt_##base

#include "grt/dynamic/layer.h"
#include "grt/dynamic/source.h"
#endif

// ================ 静态解 =====================
#ifdef __STATIC_KERNEL__
#define CONCAT_PREFIX(base) grt_static_##base

#include "grt/static/static_layer.h"
#include "grt/static/static_source.h"
#endif


#define __KERNEL_FUNC__ CONCAT_PREFIX(kernel)
#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 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)
#define __grt_wave2qwv_z_REV_PSV CONCAT_PREFIX(wave2qwv_z_REV_PSV)
#define __grt_wave2qwv_z_REV_SH CONCAT_PREFIX(wave2qwv_z_REV_SH)

#include "grt/common/kernel_template.h"
/**
* kernel函数根据(5.5.3)式递推计算广义反射透射矩阵， 再根据公式得到
*
* 1.EX 爆炸源， (P0)
* 2.VF 垂直力源, (P0, SV0)
* 3.HF 水平力源, (P1, SV1, SH1)
* 4.DC 剪切源, (P0, SV0), (P1, SV1, SH1), (P2, SV2, SH2)
*
* 的 \f$ q_m, w_m, v_m \f$ 系数(\f$ m=0,1,2 \f$),
*
* eg. DC_qwv[i][j]表示 \f$ m=i \f$ 阶时的 \f$ q_m(j=0), w_m(j=1), v_m(j=2) \f$ 系数
*
* 在递推得到广义反射透射矩阵后，计算位移系数的公式本质是类似的，但根据震源和接受点的相对位置，
* 空间划分为多个层，公式也使用不同的矩阵，具体为
*
*
* \f[
* \begin{array}{c}
* \\\\ \hline
* \hspace{5cm}\text{Free Surface(自由表面)}\hspace{5cm} \\\\
* \text{...} \\\\ \hline
* \text{Source/Receiver interface(震源/接收虚界面) (A面)} \\\\
* \text{...} \\\\ \hline
* \text{Receiver/Source interface(接收/震源虚界面) (B面)} \\\\
* \text{...} \\\\ \hline
* \text{Lower interface(底界面)} \\\\
* \text{...} \\
* \text{(无穷深)} \\
* \text{...} \\
*
*
* \end{array}
* \f]
*
* 界面之间构成一个广义层，每个层都对应2个反射系数矩阵RD/RU和2个透射系数矩阵TD/TU,
* 根据公式的整理结果，但实际需要的矩阵为：
*
* | 广义层 | **台站在震源上方** | **台站在震源下方** |
* |----------|-------------------|-------------------|
* | FS (震源 <-> 表面) | RU | RD, RU, TD, TU |
* | FR (接收 <-> 表面) | RD, RU, TD, TU | / |
* | RS (震源 <-> 接收) | RD, RU, TD, TU | RD, RU, TD, TU |
* | SL (震源 <-> 底面) | RD | RD |
* | RL (接收 <-> 底面) | / | RD |
*
*
*
*
* @note 关于与自由表面相关的系数矩阵要注意，FS表示(z1, zR+)之间的效应，但通常我们
* 定义KP表示(zK+, zP+)之间的效应，所以这里F表示已经加入了自由表面的作用，
* 对应的我们使用ZR表示(z1+, zR+)的效应，FR和ZR也满足类似的递推关系。
* @note 从公式推导上，例如RD_RS，描述的是(zR+, zS-)的效应，但由于我们假定
* 震源位于介质层内，则z=zS并不是介质的物理分界面，此时 \f$ D_{j-1}^{-1} * D_j = I \f$，
* 故在程序可更方便的编写。（这个在静态情况下不成立，不能以此优化）
* @note 接收点位于自由表面的情况 不再单独考虑，合并在接受点浅于震源的情况
*
*
* 为了尽量减少冗余的计算，且保证程序的可读性，可将震源层和接收层抽象为A,B层，
* 即空间划分为FA,AB,BL, 计算这三个广义层的系数矩阵，再讨论震源层和接收层的深浅，
* 计算相应的矩阵。
*
* @param[in,out] mod1d `MODEL1D` 结构体指针
* @param[in] k 波数
* @param[out] QWV 不同震源，不同阶数的核函数 \f$ q_m, w_m, v_m \f$
* @param[in] calc_uiz 是否计算ui_z（位移u对坐标z的偏导）
* @param[out] QWV_uiz 不同震源，不同阶数的核函数对z的偏导 \f$ \frac{\partial q_m}{\partial z}, \frac{\partial w_m}{\partial z}, \frac{\partial v_m}{\partial z} \f$
*
*/
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]);

#endif
/** 静态解的核函数 */
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]);
1 change: 1 addition & 0 deletions pygrt/C_extension/include/grt/common/matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

#pragma once

#include <stdio.h>
#include "grt/common/const.h"

#define GRT_INIT_ZERO_2x2_MATRIX {{0, 0}, {0, 0}} ///< 初始化复数0矩阵
Expand Down
80 changes: 80 additions & 0 deletions pygrt/C_extension/src/common/kernel.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
/**
* @file kernel.c
* @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
* @date 2025-12
*
* 动态或静态下计算核函数
*
*/

#include <stdbool.h>
#include "grt/common/matrix.h"
#include "grt/common/kernel.h"

/**
* 最终公式(5.7.12,13,26,27)简化为 (P-SV波) :
* + 当台站在震源上方时：
*
* \f[
* \begin{pmatrix} q_m \\ w_m \end{pmatrix} = \mathbf{R_1}
* \left[
* \mathbf{R_2} \begin{pmatrix} P_m^+ \\ SV_m^+ \end{pmatrix}
* + \begin{pmatrix} P_m^- \\ SV_m^- \end{pmatrix}
* \right]
* \f]
*
* + 当台站在震源下方时：
*
* \f[
* \begin{pmatrix} q_m \\ w_m \end{pmatrix} = \mathbf{R_1}
* \left[
* \begin{pmatrix} P_m^+ \\ SV_m^+ \end{pmatrix}
* + \mathbf{R_2} \begin{pmatrix} P_m^- \\ SV_m^- \end{pmatrix}
* \right]
* \f]
*
* SH波类似，但是是标量形式。
*
* @param[in] ircvup 接收层是否浅于震源层
* @param[in] R1 P-SV波，\f$\mathbf{R_1}\f$矩阵
* @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$
*/
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])
{
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]);
}
grt_cmat2x1_mul(R1, qw0, qw1);

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


// ================ 动态解 =====================
#define __DYNAMIC_KERNEL__
#include "kernel_template.c_"
#undef __DYNAMIC_KERNEL__

// ================ 静态解 =====================
#define __STATIC_KERNEL__
#include "kernel_template.c_"
#undef __STATIC_KERNEL__
Loading