Skip to content

REFAC: add "grnspec.h/c" to handle the Green Functions spectra #242

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
Apr 7, 2026
Merged

REFAC: add "grnspec.h/c" to handle the Green Functions spectra #242

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
1 change: 1 addition & 0 deletions pygrt/C_extension/include/grt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@


#include "grt/dynamic/grn.h"
#include "grt/dynamic/grnspec.h"
#include "grt/dynamic/layer.h"
#include "grt/dynamic/signals.h"
#include "grt/dynamic/source.h"
Expand Down
36 changes: 6 additions & 30 deletions pygrt/C_extension/include/grt/dynamic/grn.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
* @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
* @date 2024-07-24
*
* 以下代码实现的是 广义反射透射系数矩阵+离散波数法 计算理论地震图，参考：
* 以下代码实现的是利用多种波数积分类方法计算理论地震图，参考：
*
* 1. 姚振兴, 谢小碧. 2022/03. 理论地震图及其应用（初稿）.
* 2. Yao Z. X. and D. G. Harkrider. 1983. A generalized refelection-transmission coefficient
Expand All @@ -15,43 +15,19 @@
#pragma once

#include "grt/common/model.h"
#include "grt/dynamic/grnspec.h"
#include "grt/integral/integ_method.h"

/**
* 积分计算Z, R, T三个分量格林函数的频谱的核心函数（被Python调用）
*
* @param[in,out] mod1d `MODEL1D` 结构体指针
* @param[in] nf1 开始计算频谱的频率索引值, 总范围在[nf1, nf2]
* @param[in] nf2 结束计算频谱的频率索引值, 总范围在[nf1, nf2]
* @param[in] freqs 所有频点的频率值（包括未计算的）
* @param[in] nr 震中距数量
* @param[in] rs 震中距数组
* @param[in] wI 虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$
* @param[in] keepAllFreq 计算所有频点，不论频率多低
* @param[in,out] Kmet 波数积分相关参数的结构体指针
* @param[in,out] mod1d `MODEL1D` 结构体指针
* @param[in,out] Kmet 波数积分相关参数的结构体指针
* @param[in,out] grn 计算的格林函数频谱
* @param[in] print_progressbar 是否打印进度条
* @param[in] calc_upar 是否计算位移u的空间导数
* @param[out] grn 不同震源不同阶数的格林函数的Z、R、T分量频谱结果
* @param[out] grn_uiz 不同震源不同阶数的ui_z的Z、R、T分量频谱结果
* @param[out] grn_uir 不同震源不同阶数的ui_r的Z、R、T分量频谱结果
*
* @param[in] statsstr 积分过程输出目录
* @param[in] nstatsidxs 输出积分过程的特定频点的个数
* @param[in] statsidxs 特定频点的索引值
*
*/
void grt_integ_grn_spec(
MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs,
size_t nr, real_t *rs, real_t wI, bool keepAllFreq, K_INTEG_METHOD *Kmet,
bool print_progressbar, bool calc_upar,
pcplxChnlGrid grn[nr],
pcplxChnlGrid grn_uiz[nr],
pcplxChnlGrid grn_uir[nr],

const char *statsstr, // 积分过程输出
size_t nstatsidxs, // 仅输出特定频点
size_t *statsidxs
);
void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_METHOD *Kmet, GRNSPEC *grn, const bool print_progressbar);



Expand Down
62 changes: 62 additions & 0 deletions pygrt/C_extension/include/grt/dynamic/grnspec.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
/**
* @file grnspec.h
* @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
* @date 2026-04
*
* 将格林函数频谱用一个结构体包裹
*
*
*/

#pragma once

#include "grt/common/const.h"
#include "grt/common/myfftw.h"
#include "grt/common/sacio.h"

/** 不同震中距、不同震源、不同分量的格林函数频谱 */
typedef struct {
size_t nf; ///< 总频点数
real_t *freqs; ///< 频率数组 freqs[nf]
size_t nf1;
size_t nf2; ///< 待计算的频率索引
size_t nr; ///< 震中距数量
real_t *rs; ///< 震中距数组
real_t wI; ///< 虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$
bool keepAllFreq; ///< 是否计算所有频点，不论频率多低
bool calc_upar; ///< 是否计算位移u的空间导数

pcplxChnlGrid *u; ///< 不同震源不同阶数的格林函数的Z、R、T分量频谱结果
pcplxChnlGrid *uiz; ///< 不同震源不同阶数的格林函数的Z、R、T分量对z偏导的频谱结果
pcplxChnlGrid *uir; ///< 不同震源不同阶数的格林函数的Z、R、T分量对r偏导的频谱结果

char *statsstr; ///< 积分结果输出路径
size_t nstatsidxs; ///< 仅输出特定频点的对应数量
size_t *statsidxs; ///< 对应频点的频率索引 statsidxs[nstatsidxs]
} GRNSPEC;

/** 申请 u, uiz, uir 的内存 */
void grt_grnspec_allocate_u(GRNSPEC *grn);

/** 释放 u, uiz, uir 的内存 */
void grt_grnspec_free_u(GRNSPEC *grn);

/**
* 将频谱 u, uiz, uir 变换到时域后以 SAC 格式保存到本地
*
* @param[in] grn 格林函数频谱结构体
* @param[in] travtPS 不同震源距的初至P、S到时
* @param[in] begintimes 不同震中距的波形时移
* @param[in] outputdirs 不同震中距的保存目录
* @param[in] fh 控制反傅里叶变换的结构体
* @param[in,out] sac SACTRACE 原型，在头段变量中记录了基本信息
* @param[in] skipImagComps 跳过虚频率的补偿
* @param[in] saveEX 保存爆炸源结果
* @param[in] saveVF 保存垂直力源结果
* @param[in] saveHF 保存水平力源结果
* @param[in] saveDC 保存水平力源结果
*
*/
void grt_grnspec_write_sac(
const GRNSPEC *grn, const real_t (*travtPS)[2], const real_t *begintimes, char **outputdirs, GRT_FFTW_HOLDER *fh, SACTRACE *sac,
const bool skipImagComps, const bool saveEX, const bool saveVF, const bool saveHF, const bool saveDC);
69 changes: 29 additions & 40 deletions pygrt/C_extension/src/dynamic/grn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
* @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
* @date 2024-07-24
*
* 以下代码实现的是 广义反射透射系数矩阵+离散波数法 计算理论地震图，参考：
* 以下代码实现的是利用多种波数积分类方法计算理论地震图，参考：
*
* 1. 姚振兴, 谢小碧. 2022/03. 理论地震图及其应用（初稿）.
* 2. Yao Z. X. and D. G. Harkrider. 1983. A generalized refelection-transmission coefficient
Expand All @@ -20,6 +20,7 @@
#include <sys/time.h>

#include "grt/dynamic/grn.h"
#include "grt/dynamic/grnspec.h"
#include "grt/integral/integ_method.h"

#include "grt/common/const.h"
Expand All @@ -29,48 +30,36 @@


/**
* 将计算好的复数形式的积分结果记录到GRN结构体中
* 将计算好的复数形式的积分结果转化为频谱
*
* @param[in] iw 当前频率索引值
* @param[in] nr 震中距个数
* @param[in] coef 统一系数
* @param[in] sumJ 积分结果
* @param[out] grn 三分量频谱
* @param[out] u 三分量频谱
*/
static void recordin_GRN(
size_t iw, size_t nr, cplx_t coef, cplxIntegGrid sumJ[nr],
pcplxChnlGrid grn[nr])
static void recordin_GRN(size_t iw, size_t nr, cplx_t coef, cplxIntegGrid sumJ[nr], pcplxChnlGrid u[nr])
{
// 局部变量，将某个频点的格林函数谱临时存放
cplxChnlGrid *tmp_grn = (cplxChnlGrid *)calloc(nr, sizeof(*tmp_grn));
cplxChnlGrid *tmp_u = (cplxChnlGrid *)calloc(nr, sizeof(*tmp_u));

for(size_t ir=0; ir<nr; ++ir){
grt_merge_Pk(sumJ[ir], tmp_grn[ir]);
grt_merge_Pk(sumJ[ir], tmp_u[ir]);

GRT_LOOP_ChnlGrid(im, c){
int modr = GRT_SRC_M_ORDERS[im];
if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue;

grn[ir][im][c][iw] = coef * tmp_grn[ir][im][c];
u[ir][im][c][iw] = coef * tmp_u[ir][im][c];
}
}

GRT_SAFE_FREE_PTR(tmp_grn);
GRT_SAFE_FREE_PTR(tmp_u);
}



void grt_integ_grn_spec(
MODEL1D *mod1d, size_t nf1, size_t nf2, real_t *freqs,
size_t nr, real_t *rs, real_t wI, bool keepAllFreq, K_INTEG_METHOD *Kmet,
bool print_progressbar, bool calc_upar,
pcplxChnlGrid grn[nr],
pcplxChnlGrid grn_uiz[nr],
pcplxChnlGrid grn_uir[nr],

const char *statsstr, // 积分结果输出
size_t nstatsidxs, // 仅输出特定频点
size_t *statsidxs)
void grt_integ_grn_spec(MODEL1D *mod1d, K_INTEG_METHOD *Kmet, GRNSPEC *grn, const bool print_progressbar)
{
// 程序运行开始时间
struct timeval begin_t;
Expand All @@ -84,30 +73,30 @@ void grt_integ_grn_spec(
int progress=0;

// 记录每个频率的计算中是否有除0错误
int *freq_invstats = (int *)calloc(nf2+1, sizeof(int));
int *freq_invstats = (int *)calloc(grn->nf2 + 1, sizeof(int));

// 实际计算的频点数
size_t nf_valid = nf2 - nf1 + 1;
size_t nf_valid = grn->nf2 - grn->nf1 + 1;

mod1d->omgref = PI2*freqs[nf2];
mod1d->omgref = PI2*grn->freqs[grn->nf2];

// 频率omega循环
// schedule语句可以动态调度任务，最大程度地使用计算资源
#pragma omp parallel for schedule(guided) default(shared)
for(size_t iw=nf1; iw<=nf2; ++iw){
real_t w = freqs[iw]*PI2; // 实频率
cplx_t omega = w - wI*I; // 复数频率 omega = w - i*wI
for(size_t iw = grn->nf1; iw <= grn->nf2; ++iw){
real_t w = grn->freqs[iw]*PI2; // 实频率
cplx_t omega = w - grn->wI*I; // 复数频率 omega = w - i*wI

// 如果在虚频率的帮助下，频率仍然距离原点太近，
// 计算会有严重的数值问题，因此直接根据频率距离原点的距离，
// 跳过该频率，没有必要再计算
if( ! keepAllFreq )
if( ! grn->keepAllFreq )
{
real_t R = 0.1; // 完全经验性地设定，暂不必要暴露在用户可控层面
if(fabs(omega) < R){
#pragma omp critical
{
GRTRaiseWarning("Skip low frequency (iw=%zu, freq=%.5e).", iw, freqs[iw]);
GRTRaiseWarning("Skip low frequency (iw=%zu, freq=%.5e).", iw, grn->freqs[iw]);
nf_valid--;
}
if(nf_valid == 0) GRTRaiseError("NO VALID FREQUENCIES.");
Expand All @@ -132,13 +121,13 @@ void grt_integ_grn_spec(
grt_update_mod1d_state_omega(local_mstat, omega, false);

// 是否要输出积分过程文件
bool needfstats = (statsstr!=NULL && (grt_findElement_size_t(statsidxs, nstatsidxs, iw) >= 0));
bool needfstats = (grn->statsstr!=NULL && (grt_findElement_size_t(grn->statsidxs, grn->nstatsidxs, iw) >= 0));

// 为当前频率创建波数积分记录文件
if(needfstats){
char *suffix = NULL;
GRT_SAFE_ASPRINTF(&suffix, "_%04zu_%.5e", iw, freqs[iw]);
grt_KMET_init_fstats(nr, rs, statsstr, suffix, local_Kmet);
GRT_SAFE_ASPRINTF(&suffix, "_%04zu_%.5e", iw, grn->freqs[iw]);
grt_KMET_init_fstats(grn->nr, grn->rs, grn->statsstr, suffix, local_Kmet);
GRT_SAFE_FREE_PTR(suffix);
}

Expand All @@ -150,20 +139,20 @@ void grt_integ_grn_spec(
// Wavenumber Integration
// 波数积分上限
local_Kmet->kmax = hypot(local_Kmet->k0, local_Kmet->ampk * w / local_Kmet->vmin);
K_INTEG *Kint = grt_wavenumber_integral(local_mstat, nr, rs, local_Kmet, calc_upar, grt_kernel);
K_INTEG *Kint = grt_wavenumber_integral(local_mstat, grn->nr, grn->rs, local_Kmet, grn->calc_upar, grt_kernel);

// 记录到格林函数结构体内
// 如果计算核函数过程中存在除零错误，则放弃该频率【通常在大震中距的低频段】
recordin_GRN(iw, nr, coef, Kint->sumJ, grn);
if(calc_upar){
recordin_GRN(iw, nr, coef, Kint->sumJz, grn_uiz);
recordin_GRN(iw, nr, coef, Kint->sumJr, grn_uir);
recordin_GRN(iw, grn->nr, coef, Kint->sumJ, grn->u);
if(grn->calc_upar){
recordin_GRN(iw, grn->nr, coef, Kint->sumJz, grn->uiz);
recordin_GRN(iw, grn->nr, coef, Kint->sumJr, grn->uir);
}
// ===================================================================================

freq_invstats[iw] = local_mstat->stats;

if(needfstats) grt_KMET_destroy_fstats(nr, local_Kmet);
if(needfstats) grt_KMET_destroy_fstats(grn->nr, local_Kmet);

grt_free_mod1d_state(local_mstat);

Expand All @@ -180,9 +169,9 @@ void grt_integ_grn_spec(
} // END omega loop

// 打印 freq_invstats
for(size_t iw=nf1; iw<=nf2; ++iw){
for(size_t iw = grn->nf1; iw <= grn->nf2; ++iw){
if(freq_invstats[iw]==GRT_INVERSE_FAILURE){
GRTRaiseWarning("iw=%zu, freq=%e(Hz), meet Zero Divison Error, results are filled with 0.\n", iw, freqs[iw]);
GRTRaiseWarning("iw=%zu, freq=%e(Hz), meet Zero Divison Error, results are filled with 0.\n", iw, grn->freqs[iw]);
}
}
GRT_SAFE_FREE_PTR(freq_invstats);
Expand Down
Loading
Loading