Skip to content

FEAT: support Free, Rigid and Halfspace boundary condition #212

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 3 commits into from
Jan 10, 2026
Merged

FEAT: support Free, Rigid and Halfspace boundary condition #212

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fc9744a
support boundary condition setting
Dengda98 Jan 8, 2026
bb0d8e9
update
Dengda98 Jan 9, 2026
29fce46
update
Dengda98 Jan 9, 2026
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
22 changes: 6 additions & 16 deletions pygrt/C_extension/include/grt/common/RT_matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,22 +33,12 @@ typedef struct {
} RT_MATRIX;


/** 初始化 R/T 矩阵 */
#define grt_init_RT_matrix(M) \
RT_MATRIX *M = &(RT_MATRIX){\
.RD = GRT_INIT_ZERO_2x2_MATRIX,\
.RU = GRT_INIT_ZERO_2x2_MATRIX,\
.TD = GRT_INIT_IDENTITY_2x2_MATRIX,\
.TU = GRT_INIT_IDENTITY_2x2_MATRIX,\
.RDL = 0.0,\
.RUL = 0.0,\
.TDL = 1.0,\
.TUL = 1.0,\
.invT = GRT_INIT_ZERO_2x2_MATRIX,\
.invTL = 0.0,\
.stats = GRT_INVERSE_SUCCESS \
};\

/**
* 重新设置 R/T 矩阵为默认值
*
* @param[out] M RT_MATRIX 指针
*/
void grt_reset_RT_matrix(RT_MATRIX *M);

/**
* 合并 recursion_RD_PSV(SH) ，仅计算RD/RDL
Expand Down
6 changes: 6 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 @@ -141,6 +141,12 @@ typedef size_t sizeIntegGrid[GRT_SRC_M_NUM][GRT_INTEG_NUM];
for(int im = 0; im < GRT_SRC_M_NUM; ++im) \
for(int v = 0; v < GRT_INTEG_NUM; ++v)

/** 边界条件 */
typedef enum {
GRT_BOUND_FREE = 0,
GRT_BOUND_RIGID,
GRT_BOUND_HALFSPACE,
} GRT_BOUND_TYPE;

/** 不同震源类型在大小为 GRT_SRC_M_NUM 的数组中的索引 */
enum {
Expand Down
15 changes: 14 additions & 1 deletion 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 @@ -58,8 +58,13 @@ typedef struct {
RT_MATRIX M_FA;
RT_MATRIX M_FB;

/* 自由表面的反射系数矩阵，仅 RU, RUL 有用 */

/* 顶界面的反射系数矩阵，仅 RU, RUL 有用 */
GRT_BOUND_TYPE topbound;
RT_MATRIX M_top;
/* 底界面的反射系数矩阵，仅 RD, RDL 有用 */
GRT_BOUND_TYPE botbound;
RT_MATRIX M_bot;

/* 接收点处的接收矩阵 (转为位移u和位移导数uiz的(B_m, C_m, P_m)系分量) */
cplx_t R_EV[2][2];
Expand Down Expand Up @@ -146,6 +151,14 @@ void grt_realloc_mod1d(GRT_MODEL1D *mod1d, size_t n);
*/
GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid);

/**
* 设置模型的边界条件，并对底界面做检查
*
* @param[in,out] mod1d `MODEL1D` 结构体指针
* @param[in] topbound 顶层边界条件
* @param[in] botbound 底层边界条件
*/
void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound);

/**
* 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数
Expand Down
15 changes: 12 additions & 3 deletions pygrt/C_extension/include/grt/dynamic/layer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,22 @@
#include "grt/common/RT_matrix.h"

/**
* 计算自由表面的反射系数，公式(5.3.10-14)
* 计算不同边界条件下顶界面的反射系数RU，其中自由表面的公式见(5.3.10-14)
*
*
* @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_top
*
*/
void grt_topfree_RU(GRT_MODEL1D *mod1d);
void grt_topbound_RU(GRT_MODEL1D *mod1d);

/**
* 计算不同边界条件下底界面的反射系数RD，其中自由表面的公式见(5.3.10-14)
*
*
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot
*
*/
void grt_botbound_RD(GRT_MODEL1D *mod1d);


/**
Expand Down
13 changes: 10 additions & 3 deletions pygrt/C_extension/include/grt/static/static_layer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,20 @@
#include "grt/common/model.h"

/**
* 计算自由表面的静态反射系数，公式(6.3.12)
* 计算不同边界条件下顶界面的静态反射系数RU，其中自由表面的公式见(6.3.12)
*
* @param[in,out] mod1d 模型结构体指针，结果保存在 Mtop
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_top
*
*/
void grt_static_topfree_RU(GRT_MODEL1D *mod1d);
void grt_static_topbound_RU(GRT_MODEL1D *mod1d);

/**
* 计算不同边界条件下底界面的静态反射系数RD，其中自由表面的公式见(6.3.12)
*
* @param[in,out] mod1d 模型结构体指针，结果保存在 M_bot
*
*/
void grt_static_botbound_RD(GRT_MODEL1D *mod1d);

/**
* 计算接收点位置的 P-SV 静态接收矩阵，将波场转为位移，公式(6.3.35)
Expand Down
19 changes: 19 additions & 0 deletions pygrt/C_extension/src/common/RT_matrix.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,31 @@
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>

#include "grt/common/const.h"
#include "grt/common/matrix.h"
#include "grt/common/RT_matrix.h"


void grt_reset_RT_matrix(RT_MATRIX *M)
{
RT_MATRIX *M0 = &(RT_MATRIX){
.RD = GRT_INIT_ZERO_2x2_MATRIX,
.RU = GRT_INIT_ZERO_2x2_MATRIX,
.TD = GRT_INIT_IDENTITY_2x2_MATRIX,
.TU = GRT_INIT_IDENTITY_2x2_MATRIX,
.RDL = 0.0,
.RUL = 0.0,
.TDL = 1.0,
.TUL = 1.0,
.invT = GRT_INIT_ZERO_2x2_MATRIX,
.invTL = 0.0,
.stats = GRT_INVERSE_SUCCESS
};
memcpy(M, M0, sizeof(*M));
}

void grt_recursion_RD(const RT_MATRIX *M1, const RT_MATRIX *M2, RT_MATRIX *M)
{
grt_recursion_RD_PSV(M1, M2, M);
Expand Down
20 changes: 20 additions & 0 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 @@ -421,10 +421,30 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea
// 先指定负频率，仅填充弹性模量
grt_attenuate_mod1d(mod1d, -1);

// 设置一个默认边界条件
mod1d->topbound = GRT_BOUND_FREE;
mod1d->botbound = GRT_BOUND_HALFSPACE;

return mod1d;
}


void grt_set_mod1d_boundary(GRT_MODEL1D *mod1d, GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound)
{
mod1d->topbound = topbound;
mod1d->botbound = botbound;

size_t n = mod1d->n;
bool atBottom = mod1d->isrc+1 >= n || mod1d->ircv+1 >= n;
if (atBottom && botbound != GRT_BOUND_HALFSPACE){
GRTRaiseError(
"When source/recevier is at the bottom, the bottom layer must be a halfspace. "
"Please check your boundary condition setting."
);
}
}


void grt_get_model_diglen_from_file(const char *modelpath, size_t diglen[6]){
FILE *fp = GRTCheckOpenFile(modelpath, "r");
size_t len;
Expand Down
44 changes: 42 additions & 2 deletions pygrt/C_extension/src/dynamic/grt_greenfn.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,12 @@ typedef struct {
char *s_depsrc;
char *s_deprcv;
} D;
/** 顶层和底层的边界条件 */
struct {
bool active;
GRT_BOUND_TYPE topbound;
GRT_BOUND_TYPE botbound;
} B;
/** 波形时窗 */
struct {
bool active;
Expand Down Expand Up @@ -212,6 +218,9 @@ static void print_Ctrl(const GRT_MODULE_CTRL *Ctrl){
printf(format, "applyDCM", (Ctrl->C.applyDCM)? "true" : "false");
printf(format, "applyPTAM", (Ctrl->C.applyPTAM)? "true" : "false");

printf(format_size, "topbound", (size_t)Ctrl->B.topbound);
printf(format_size, "botbound", (size_t)Ctrl->B.botbound);

printf(format_size, "nt", Ctrl->N.nt);
printf(format_real, "dt", Ctrl->N.dt);
printf(format_real, "winT", Ctrl->N.winT);
Expand Down Expand Up @@ -289,7 +298,7 @@ printf("\n"
" -R<r1>,<r2>[,...] -O<outdir> [-H<f1>/<f2>] \n"
" [-L<length>] [-C[d|p|n]] [-E<t0>[/<v0>]] \n"
" [-K[+k<k0>][+s<ampk>][+e<keps>][+v<vmin>]]\n"
" [-P<nthreads>] [-Ge|v|h|s] \n"
" [-P<nthreads>] [-Ge|v|h|s] [-Bf|F|r|R|h|H]\n"
" [-S[<i1>,<i2>,...]] [-e] [-s]\n"
"\n\n"
"Options:\n"
Expand Down Expand Up @@ -397,6 +406,13 @@ printf("\n"
" <b3>: Horizontal Force (HF)\n"
" <b4>: Shear (DC)\n"
"\n"
" -Bf|F|r|R|h|H\n"
" Boundary condition of top layer (lowercase) and\n"
" bottom layer (uppercase).\n"
" f|F: Free boundary.\n"
" r|R: Rigid boundary.\n"
" h|H: Halfspace.\n"
"\n"
" -S[<i1>,<i2>,...]\n"
" Frequency (index) of statsfile in wavenumber\n"
" integration to be output, require 0 <= i <= nf-1,\n"
Expand Down Expand Up @@ -429,6 +445,9 @@ printf("\n"
/** 从命令行中读取选项，处理后记录到全局变量中 */
static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
// 先为个别参数设置非0初始值
Ctrl->B.topbound = GRT_BOUND_FREE;
Ctrl->B.botbound = GRT_BOUND_HALFSPACE;

Ctrl->N.zeta = GRT_GREENFN_N_ZETA;
Ctrl->N.upsample_n = GRT_GREENFN_N_UPSAMPLE;
Ctrl->H.freq1 = GRT_GREENFN_H_FREQ1;
Expand All @@ -441,7 +460,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
Ctrl->G.doDC = GRT_GREENFN_G_DC;

int opt;
while ((opt = getopt(argc, argv, ":M:D:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) {
while ((opt = getopt(argc, argv, ":M:D:B:N:O:H:L:C:E:K:R:S::P:G:esh")) != -1) {
switch (opt) {
// 模型路径，其中每行分别为
// 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs
Expand Down Expand Up @@ -471,6 +490,24 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
}
break;

// 顶层和底层的边界条件 -Bf|F|r|R|h|H
case 'B':
Ctrl->B.active = true;
if(strlen(optarg) == 0 || strlen(optarg) > 2) GRTBadOptionError(B, "");
for(size_t i = 0; i < strlen(optarg); ++i) {
switch(optarg[i]) {
case 'f': Ctrl->B.topbound = GRT_BOUND_FREE; break;
case 'r': Ctrl->B.topbound = GRT_BOUND_RIGID; break;
case 'h': Ctrl->B.topbound = GRT_BOUND_HALFSPACE; break;
case 'F': Ctrl->B.botbound = GRT_BOUND_FREE; break;
case 'R': Ctrl->B.botbound = GRT_BOUND_RIGID; break;
case 'H': Ctrl->B.botbound = GRT_BOUND_HALFSPACE; break;
default:
GRTBadOptionError(B, "unsupported -B%s.", optarg);
}
}
break;

// 点数,采样间隔,虚频率 -Nnt/dt[+w<zeta>][+n<scale>][+a]
case 'N':
Ctrl->N.active = true;
Expand Down Expand Up @@ -870,6 +907,9 @@ int greenfn_main(int argc, char **argv) {
}
GRT_MODEL1D *mod1d = Ctrl->M.mod1d;

// 边界条件
grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound);

// 当震源位于液体层中时，仅允许计算爆炸源对应的格林函数
// 程序结束前会输出对应警告
if(mod1d->Vb[mod1d->isrc]==0.0){
Expand Down
Loading