Skip to content

Support liquid layer #34

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
merged 7 commits into from
Jun 24, 2025
Merged

Support liquid layer #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
6aa1ae2
old/wrong liquid version for backup
Dengda98 Jun 17, 2025
553692f
fix formula
Dengda98 Jun 17, 2025
7ef7099
fix formula
Dengda98 Jun 23, 2025
316333d
fix formula
Dengda98 Jun 24, 2025
e0c0017
update `get_pymod_vmin_vmax` doc
Dengda98 Jun 24, 2025
45002f8
add `allowLiquid` argument in function `read_pymod_from_file`
Dengda98 Jun 24, 2025
7a7bc09
update
Dengda98 Jun 24, 2025
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
2 changes: 2 additions & 0 deletions pygrt/C_extension/include/common/const.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,8 @@ typedef int MYINT; ///< 整数
#define MIN_DEPTH_GAP_SRC_RCV 1.0 ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限）
#define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数

#define GRT_SWAP(type, a, b) { type temp = a; a = b; b = temp; } ///< 交换两个变量的值

#define GRT_STRING_FMT "%18s" ///< 字符串输出格式
#define GRT_REAL_FMT "%18.8e" ///< 浮点数输出格式
#define GRT_CMPLX_FMT "%18.8e%-+14.8ej" ///< 复数输出格式
Expand Down
14 changes: 13 additions & 1 deletion pygrt/C_extension/include/common/model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,8 +143,20 @@ PYMODEL1D * init_pymod(MYINT n);
* @param[in] modelpath 模型文件路径
* @param[in] depsrc 震源深度
* @param[in] deprcv 接收深度
* @param[in] allowLiquid 是否允许液体层
*
* @return `PYMODEL1D` 结构体指针
*
*/
PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv);
PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid);


/**
* 计算PYMODEL1D结构体中的最大最小速度（非零值）
*
* @param pymod (in)`PYMODEL1D` 结构体指针
* @param vmin (out)最小速度
* @param vmax (out)最大速度
*
*/
void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax);
25 changes: 24 additions & 1 deletion pygrt/C_extension/include/dynamic/layer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ void calc_uiz_R_EV(
* @param[in] kbkb2 下层的S波水平波数的平方 \f$ k_b^2=(\frac{\omega}{V_b})^2 \f$
* @param[in] mu2 下层的剪切模量
* @param[in] thk 上层层厚
* @param[in] omega 角频率
* @param[in] k 波数
* @param[out] RD P-SV 下传反射系数矩阵
* @param[out] RDL SH 下传反射系数
Expand All @@ -104,11 +105,33 @@ void calc_uiz_R_EV(
void calc_RT_2x2(
MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1,
MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2,
MYREAL thk, MYREAL k,
MYREAL thk, MYCOMPLEX omega, MYREAL k,
MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL,
MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats);

/** 固-固 界面，函数参数与 calc_RT_2x2 函数一致 */
void calc_RT_ss_2x2(
MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1,
MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2,
MYREAL thk, MYCOMPLEX omega, MYREAL k,
MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL,
MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats);

/** 液-液 界面，函数参数与 calc_RT_2x2 函数一致 */
void calc_RT_ll_2x2(
MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1,
MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2,
MYREAL thk, MYCOMPLEX omega, MYREAL k,
MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL,
MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats);

/** 液-固 界面，函数参数与 calc_RT_2x2 函数一致 */
void calc_RT_ls_2x2(
MYREAL Rho1, MYCOMPLEX xa1, MYCOMPLEX xb1, MYCOMPLEX kbkb1, MYCOMPLEX mu1,
MYREAL Rho2, MYCOMPLEX xa2, MYCOMPLEX xb2, MYCOMPLEX kbkb2, MYCOMPLEX mu2,
MYREAL thk, MYCOMPLEX omega, MYREAL k,
MYCOMPLEX RD[2][2], MYCOMPLEX *RDL, MYCOMPLEX RU[2][2], MYCOMPLEX *RUL,
MYCOMPLEX TD[2][2], MYCOMPLEX *TDL, MYCOMPLEX TU[2][2], MYCOMPLEX *TUL, MYINT *stats);

/**
* 【未使用，仅用于代码测试】
Expand Down
27 changes: 23 additions & 4 deletions pygrt/C_extension/src/common/model.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,7 @@ void update_mod1d_omega(MODEL1D *mod1d, MYCOMPLEX omega){
atnb = attenuation_law(lay->Qbinv, omega);

ka0 = omega/(Va0*atna);
kb0 = omega/(Vb0*atnb);
kb0 = (Vb0>RZERO)? omega/(Vb0*atnb) : CZERO;
lay->kaka = ka0*ka0;
lay->kbkb = kb0*kb0;

Expand Down Expand Up @@ -255,7 +255,7 @@ void realloc_pymod(PYMODEL1D *pymod, MYINT n){
}


PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv){
PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, double depsrc, double deprcv, bool allowLiquid){
FILE *fp;
if((fp = fopen(modelpath, "r")) == NULL){
fprintf(stderr, "[%s] " BOLD_RED "Model file open error.\n" DEFAULT_RESTORE, command);
Expand Down Expand Up @@ -303,8 +303,13 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou
return NULL;
};

if(va <= 0.0 || vb <= 0.0 || rho <= 0.0 || qa <= 0.0 || qb <= 0.0){
fprintf(stderr, "[%s] " BOLD_RED "In line %d, nonpositive value is not supported.\n" DEFAULT_RESTORE, command, iline);
if(va <= 0.0 || rho <= 0.0 || qa <= 0.0 || qb <= 0.0){
fprintf(stderr, "[%s] " BOLD_RED "In model file, line %d, nonpositive value is not supported.\n" DEFAULT_RESTORE, command, iline);
return NULL;
}

if(!allowLiquid && vb <=0){
fprintf(stderr, "[%s] " BOLD_RED "In model file, line %d, nonpositive Vs is not supported.\n" DEFAULT_RESTORE, command, iline);
return NULL;
}

Expand Down Expand Up @@ -391,4 +396,18 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou
free(modarr);

return pymod;
}


void get_pymod_vmin_vmax(const PYMODEL1D *pymod, double *vmin, double *vmax){
*vmin = __DBL_MAX__;
*vmax = RZERO;
const MYREAL *Va = pymod->Va;
const MYREAL *Vb = pymod->Vb;
for(MYINT i=0; i<pymod->n; ++i){
if(Va[i] < *vmin) *vmin = Va[i];
if(Va[i] > *vmax) *vmax = Va[i];
if(Vb[i] < *vmin && Vb[i] > RZERO) *vmin = Vb[i];
if(Vb[i] > *vmax && Vb[i] > RZERO) *vmax = Vb[i];
}
}
9 changes: 2 additions & 7 deletions pygrt/C_extension/src/dynamic/grt_main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -991,17 +991,12 @@ int main(int argc, char **argv) {
getopt_from_command(argc, argv);

// 读入模型文件
if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv)) ==NULL){
if((pymod = read_pymod_from_file(command, s_modelpath, depsrc, deprcv, true)) ==NULL){
exit(EXIT_FAILURE);
}

// 最大最小速度
vmax = pymod->Va[findMinMax_MYREAL(pymod->Va, pymod->n, true)];
vmin = pymod->Vb[findMinMax_MYREAL(pymod->Vb, pymod->n, false)];
if(vmin > vmax) {
double tmp;
tmp = vmin; vmin = vmax; vmax = tmp;
}
get_pymod_vmin_vmax(pymod, &vmin, &vmax);

// 参考最小速度
if(vmin_ref == 0.0){
Expand Down
Loading